text:
mutable struct MemcacheClient
    host::AbstractString
    port::Integer
    sock::IO
    debug::Bool

    function MemcacheClient(host::AbstractString="localhost", port::Integer=11211)
        sock = connect(host, port)
        new(host, port, sock, false)
    end
end

close(mc::MemcacheClient) = close(mc.sock)
quit(mc::MemcacheClient) = mc_send(mc, "quit")

for f in ("set", "add", "replace", "append", "prepend")
    @eval ($(Symbol(f)))(mc::MemcacheClient, key::AbstractString, val; exptime::Integer=0, cas::Integer=-1, noreply::Bool=false) = _setval(mc, $f, key, val, exptime=exptime, cas=cas, noreply=noreply)
end

cas(mc::MemcacheClient, key::AbstractString, val, cas::Integer; exptime::Integer=0, noreply::Bool=false) = _setval(mc, "cas", key, val, exptime=exptime, cas=cas, noreply=noreply)

function _setval(mc::MemcacheClient, cmd::AbstractString, key::AbstractString, val; exptime::Integer=0, cas::Integer=-1, noreply::Bool=false)
    flags, buff = mc_serialize(val)
    args = [key, flags, exptime, length(buff)]
    (cas != -1) && push!(args, cas)
    noreply && push!(args, "noreply")

    mc_send(mc, cmd, args...)
    mc_send(mc, buff)

    resp_line = validate(mc_recv(mc), "STORED", "NOT_STORED", "EXISTS", "NOT_FOUND")
    (resp_line[1] == "STORED") && return
    error(resp_line[1])
end

function get(mc::MemcacheClient, key::AbstractString)
    mc_send(mc, "get", key)
    cont = true

    local result::Any = nothing
    result_found = false
    while cont
        resp_line = validate(mc_recv(mc), "VALUE", "END")
        if resp_line[1] == "VALUE"
            result_found = true
            kv = split(resp_line[2])
            flags = parse(Int, kv[2])
            bytes = parse(Int, kv[3])
            result = mc_deserialize(mc_recv(mc, bytes), flags)
        else
            cont = false
        end
    end
    result_found ? result : error("NOT_FOUND")
end

function get(mc::MemcacheClient, key::AbstractString...; cas::Bool=false)
    c = (cas != nothing) ? "gets" : "get"
    mc_send(mc, c, key...)

    cont = true
    result = Dict{AbstractString,Any}()
    while cont
        resp_line = validate(mc_recv(mc), "VALUE", "END")
        if resp_line[1] == "VALUE"
            kv = split(resp_line[2])
            key = kv[1]
            flags = parse(Int, kv[2])
            bytes = parse(Int, kv[3])
            casval = cas ? parse(Int, kv[4]) : 0
            val = mc_deserialize(mc_recv(mc, bytes), flags)

            result[key] = cas ? (val, casval) : val
        else
            cont = false
        end
    end
    result
end

function delete(mc::MemcacheClient, key::AbstractString; noreply::Bool=false)
    mc_send(mc, "delete", key, noreply ? "noreply" : nothing)
    noreply && return
    resp_line = validate(mc_recv(mc), "DELETED", "NOT_FOUND")
    (resp_line[1] == "NOT_FOUND") && error(resp_line[1])
    nothing
end

incr(mc::MemcacheClient, key::AbstractString, val::Integer; noreply::Bool=false) = incr_decr(mc, "incr", key, val, noreply=noreply)
decr(mc::MemcacheClient, key::AbstractString, val::Integer; noreply::Bool=false) = incr_decr(mc, "decr", key, val, noreply=noreply)
function incr_decr(mc::MemcacheClient, c::AbstractString, key::AbstractString, val::Integer; noreply::Bool=false)
    mc_send(mc, c, key, val, noreply ? "noreply" : nothing)
    noreply && return
    resp_line = mc_recv(mc)
    (resp_line[1] == "NOT_FOUND") && error(resp_line[1])
    parse(Int, resp_line[1])
end

function touch(mc::MemcacheClient, key::AbstractString, exp::Integer; noreply::Bool=false)
    mc_send(mc, "touch", key, exp, noreply ? "noreply" : nothing)
    noreply && return
    resp_line = validate(mc_recv(mc), "TOUCHED", "NOT_FOUND")
    (resp_line[1] == "NOT_FOUND") && error(resp_line[1])
    nothing
end

function flush_all(mc::MemcacheClient, delay::Integer=0; noreply::Bool=false)
    mc_send(mc, "flush_all", (delay > 0) ? delay : nothing)
    !noreply && validate(mc_recv(mc), "OK")
    nothing
end

function stats(mc::MemcacheClient, args...)
    mc_send(mc, "stats", args...)

    result = Dict{AbstractString, AbstractString}()

    cont = true
    while cont
        resp_line = validate(mc_recv(mc), "STAT", "END")
        if resp_line[1] == "STAT"
            kv = split(resp_line[2])
            result[kv[1]] = kv[2]
        else
            cont = false
        end
    end
    result
end

function version(mc::MemcacheClient)
    mc_send(mc, "version")
    validate(mc_recv(mc), "VERSION")[2]
end

function slabs_reassign(mc::MemcacheClient, src::Integer, dst::Integer)
    mc_send(mc, "slabs", "reassign", src, dst)
    resp_line = validate(mc_recv(mc), "OK", "BUSY", "BADCLASS", "NOSPARE", "NOTFULL", "UNSAFE", "SAME")
    (resp_line[1] == "OK") && return
    error(join(resp_line, ':'))
end

function slabs_automove(mc::MemcacheClient, mode::Integer)
    !(mode in [0, 1, 2]) && error("invalid mode $mode")
    mc_send(mc, "slabs", "automove", mode)
    validate(mc_recv(mc), "OK")
    nothing
end

##
# serialize and deserialize
# plain ascii format for numbers and strings so that they are interoperable with other implementations.
# julia serialization for other types.
const ser_pipe = PipeBuffer()
mc_serialize(val) = (mc_serialize(ser_pipe, val), take!(ser_pipe))
mc_serialize(s, val) = (serialize(s, val); 0)
mc_serialize(s, val::T) where T <: Integer = (print(s, val); 1)
mc_serialize(s, val::T) where T <: AbstractFloat = (print(s, val); 2)
mc_serialize(s, val::T) where T <: AbstractString = (print(s, val); 3)

function mc_deserialize(buff::Array{UInt8,1}, typ::Int)
    if 0 == typ
        return deserialize(IOBuffer(buff))
    elseif 1 == typ
        return parse(Int, String(buff))
    elseif 2 == typ
        return parse(Float64, String(buff))
    elseif 3 == typ
        return String(buff)
    end
    error("unknown data type $typ")
end

##
# memcache send/recv protocol
const CMD_DLM = "\r\n"

function mc_send(mc::MemcacheClient, cmd::AbstractString, args...)
    iob = IOBuffer()
    write(iob, cmd)
    for arg in args
        (arg == nothing) && continue
        write(iob, ' ')
        print(iob, arg)
    end
    write(iob, CMD_DLM)
    cmdline = take!(iob)
    write(mc.sock, cmdline)
    mc.debug && println("send: ", String(cmdline))
    nothing
end

function mc_send(mc::MemcacheClient, data::Array{UInt8})
    write(mc.sock, data)
    write(mc.sock, CMD_DLM)
    mc.debug && println("send: $(length(data)) bytes data")
    nothing
end

function mc_recv(mc::MemcacheClient)
    s = readline(mc.sock)
    mc.debug && println("recv: $s")
    s = rstrip(s)
    split(s, ' ', limit=2)
end

mc_recv(mc::MemcacheClient, len::Integer) = mc_recv(mc, Vector{UInt8}(undef, len))

function mc_recv(mc::MemcacheClient, data::Array{UInt8,1})
    read!(mc.sock, data)
    readline(mc.sock)
    mc.debug && println("recv: $(length(data)) bytes data")
    data
end

validate(resp::Array, valids...) = (resp[1] in valids) ? resp : error("unrecognized data from server: $(resp[1])")


text:
# -----------------------------------------------------
#
# Monte Carlo Integration Techniques for Calculating Error
#
# -----------------------------------------------------

#
# More effective than sparse techniques for higher dimensions
#

function monte_carlo(f::Function, D::Int; count = 1000)
    total = 0.0
    for iter in 1:count
		x = ntuple(i -> rand(), D)
        total += f(x)
	end
    total /= count
end

function monte_carlo2(f::Function, D::Int; batch = 50, Z = 1.0)

    means = Real[]

    while true
        points = Real[]
        for i in 1:batch
            x = ntuple(i -> rand(), D)
            push!(points, f(x))
        end
        if length(means) >1  && (abs(mean(points)-mean(means))/std(means) < Z)
            break
        end
        push!(means, mean(points[end-batch+1:end]))
    end
    return mean(means)
end

function mcerr(f::Function, g::Function, D::Int; count = 1000)
    return sqrt(monte_carlo(x->(f(x)-g(x))^2, D; count=count))
end

function mcerr2(f::Function, g::Function, D::Int; batch = 50, Z = 1.0)
	return sqrt(monte_carlo2(x->(f(x)-g(x))^2, D; batch = batch, Z = Z))
end

ext: jl


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 208f935e-031c-11eb-172e-bf3a2463720c
using TestImages

# ╔═╡ 676fcc70-fe6a-11ea-228a-8f3af1750c03

# ╔═╡ ec7bfd20-fe57-11ea-11ff-ad7d1063c690
md"""
## Load these Packages First - Needed for Images Manipulation
Let's use the Images.jl package to load an image and see what we get

"""

# ╔═╡ 3e2b6a10-02b9-11eb-2026-5b14f80e3d02
	md"""
	Code Syntax  -->  import Pkg --> 	Pkg.activate(mktempdir())
	"""

# ╔═╡ ec64cba0-fe57-11ea-0a43-296a76f4e41f
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ ec4ef9b0-fe57-11ea-3218-19d3eb36b629
begin
	Pkg.add(["Images", "ImageIO", "ImageMagick", "TestImages"])
	using Images
end

# ╔═╡ 212b3722-031c-11eb-1002-c1d8a1e1cff2
begin
	Pkg.add("PlutoUI")
	using PlutoUI
end

# ╔═╡ 64efe220-02b9-11eb-33a9-05fc5c343994
	md"""
	Code Syntax –> 	Pkg.add(["Images", "ImageIO", "ImageMagick"]) --> using Images
	"""

# ╔═╡ e72eb300-fe69-11ea-0295-235feccceff5
	md"""
	## Open image, from URL
	"""

# ╔═╡ ec39eb10-fe57-11ea-1a17-f5177a3a2bc4
# defines a variable called `url_variable`
# whose value is a string (written inside `"`):

Url_Variable = "https://emahesha.files.wordpress.com/2008/08/9indian_flag.jpg"

# ╔═╡ ec2578b0-fe57-11ea-21f9-3d06f046e6c7
Indi_Flag_file = download(Url_Variable, "Indi_Flag.jpg")  # download to a local file

# ╔═╡ ec0dd1fe-fe57-11ea-2bff-8f9c00037c98
Tricolor = load(Indi_Flag_file)

# ╔═╡ df83edf0-02ba-11eb-3f62-b9af925b068d
	md"""
	Now can call it using the variable only like "Trocolor"
	"""

# ╔═╡ ebe03252-fe57-11ea-3150-d3945169079b
Tricolor

# ╔═╡ eb561abe-fe57-11ea-083b-7926bd86e4e9
	md"""
	How to see the Structure of a Image !!
	"""

# ╔═╡ 51309bb0-fe6a-11ea-04ca-e3b11cbf31f1
typeof(Tricolor)

# ╔═╡ 25be91a0-02be-11eb-0e62-65deaa359b84
	md"""
	The Above is Explained as follow:
- It's an `Array`

- The `2` means that it has **2 dimensions** (a **matrix**)

- `RGBX{Normed{UInt8,8}}` is the type of object stored in the array

- A Julia object representing a colour

- RGB = Red, Green, Blue
	"""

# ╔═╡ 51180aa0-fe6a-11ea-38ed-e157030c325c
	md"""
	To See the Size of an Image
	"""

# ╔═╡ d62a5da0-02bc-11eb-37f2-6d919583d571
size(Tricolor)

# ╔═╡ d60d86d0-02bc-11eb-2056-8b51aea35279
	md"""
	To see a portion of an Image - Provide the range (1st Range is Top to Bottom, 2nd 		Range is Left to right)
	"""

# ╔═╡ d5f3bd3e-02bc-11eb-05d7-d32bb9a027fc
Tricolor[30:120, 65:110]

# ╔═╡ d5d9cca0-02bc-11eb-048f-9320b39c74cd
md" Getting pieces of an image"

# ╔═╡ d5bd9212-02bc-11eb-222e-9bed9f96ac54
begin
	(h, w) = size(Tricolor)
	Slice = Tricolor[(h ÷ 4):h, (w ÷ 10): (9w ÷ 9)]
	# `÷` is typed as \div <TAB>  -- integer division
end

# ╔═╡ d5a268f0-02bc-11eb-115b-13534cc0c95e
	md"""Size of the Peice we just created above:	"""

# ╔═╡ d58284e0-02bc-11eb-2b6b-9fc3b786a3c9
size(Slice)

# ╔═╡ 108d3f70-0363-11eb-1fae-6713f3fef863
	md"""Original Size was	"""

# ╔═╡ 09aadb90-0363-11eb-165a-23513f34ae99
size(Tricolor)

# ╔═╡ 0b2d32e0-02c0-11eb-00d7-171f805c6e2d
md"## Manipulating matrices to Manipulate Images
An image is just a matrix, so we can manipulate *matrices* to manipulate the *image*
"

# ╔═╡ 0b139060-02c0-11eb-0a67-6920aee627ca
	md"""
	Lets Play with the Variable 'Slice' we created
	"""

# ╔═╡ 0af4bdc0-02c0-11eb-30aa-21865c74c6f3
[Slice Slice]

# ╔═╡ 0ad96d8e-02c0-11eb-1546-3534868d4894
	md"""	The Above keeps the image next to it - the below will put four of them creating an Kaliedoscope	"""

# ╔═╡ 0abf7cf0-02c0-11eb-0003-fb78d7a76dc0
[
 Slice                   reverse(Slice, dims=2)
 reverse(Slice, dims=1)  reverse(reverse(Slice, dims=1), dims=2)
]

# ╔═╡ 0aa58c50-02c0-11eb-1608-d1ec9abe94c6
	md"""	Create  Copy of Image (Good Habbit, Before Altering an Image)	"""

# ╔═╡ 0a8a3c20-02c0-11eb-15e2-6d41dab0deee
Dhwaj = copy(Tricolor)

# ╔═╡ 0a6ca200-02c0-11eb-0619-c1f7c9c51fb3
md"""## Syntax for three main colors - Red Green Blue"
red = RGB(1, 0, 0)   |
green = RGB(0, 1, 0)  |
blue = RGB(0, 0, 1)
"""

# ╔═╡ 0a519ff0-02c0-11eb-31e9-433a6762996c
red = RGB(1, 0, 0)

# ╔═╡ b273ab10-031a-11eb-2068-0796513e64c7
green = RGB(0, 1, 0)

# ╔═╡ 96c2e480-031a-11eb-29c3-b97cef4ae152
blue = RGB(0, 0, 1)

# ╔═╡ 0a3676d0-02c0-11eb-0510-e7536424f8b7
	md"""## Experiment - Changing Color
By Using Loop

Note that for loops do not return anything (or, rather, they return nothing)"""

# ╔═╡ 0a192ace-02c0-11eb-1e3f-b5ff3970f610
for i in 1:20
	for j in 1:20
		Dhwaj[i, j] = red
	end
end

# ╔═╡ 09c9ab40-02c0-11eb-2e14-294da8b643b4
for i in 70:100
	for j in 20:40
		Dhwaj[i, j] = blue
	end
end

# ╔═╡ 09ad70b0-02c0-11eb-1982-7f86a152b7b5
Dhwaj

# ╔═╡ 09922080-02c0-11eb-2d95-83e9fbe98d7a
	md"""## Element-wise operations: "Broadcasting"

- Julia provides powerful technology for operating element by element: **broadcasting**

- Adding "`.`" applies an operation element by element
"""

# ╔═╡ 09798f6e-02c0-11eb-09e0-3feb2b5a9048
begin
	Dhwaj2 = copy(Tricolor)
	Dhwaj2[70:90, 70:110] .= RGB(0, 1, 1)
	Dhwaj2
end

# ╔═╡ 095bce40-02c0-11eb-309a-7b4fda24cd8c
md"## Modifying the whole image at once
**TBD Left this Section for Now**
- We can use the same trick to modify the whole image at once

- Let's **redify** the image

- We define a **function** that turns a colour into just its red component
"

# ╔═╡ 0940cc30-02c0-11eb-252c-37c67ae8008e
md"**TBD Decimate not working - Need to add some package or import something**"

# ╔═╡ 21609f00-031c-11eb-2c4e-c7f936df961f
Low_Res_Dhwaj = decimate(Slice, 5)

# ╔═╡ 2146ae60-031c-11eb-10cf-6d0f7a5519ff
md"## Experiments come alive with interaction
- We start to get a feel for things when we can **experiment**!
- Next Sniped we Load --> Pkg.add(PlutoUI-quotes were used?) --> using PlutoUI (to make interaction work)"

# ╔═╡ 2110d150-031c-11eb-3ed4-4fc761239518
@bind repeat_count Slider(1:5, show_value=true)

# ╔═╡ 20f53300-031c-11eb-346c-555291725aa3
repeat(Dhwaj2, repeat_count, repeat_count)

# ╔═╡ 20d994b0-031c-11eb-2907-453907e22de5
md"## 10/5/2020
- From Here - https://github.com/JuliaImages/ImageInTerminal.jl
- TestImages was added all the way above"

# ╔═╡ 207888f2-031c-11eb-1800-5b0d98072ef8
testimage("cameraman")

# ╔═╡ 205c4e60-031c-11eb-1a04-4ba53bbe3cc4

# ╔═╡ 2044cec2-031c-11eb-38c8-d3f806df0607

# ╔═╡ 202ba170-031c-11eb-1cb2-5926b842e1f1

# ╔═╡ 20135e80-031c-11eb-1e42-459731477aa8

# ╔═╡ 1ff83560-031c-11eb-19f0-83458298f9ec

# ╔═╡ 1fde44c0-031c-11eb-33c8-87df48b53597

# ╔═╡ d527bab0-02bc-11eb-3ef2-61490bf35dab
	md"""
	Below, change values and will give you range of colors
	"""

# ╔═╡ 50fd2fa0-fe6a-11ea-22ea-93ec0773ebd3
RGBX(0.9, 0.9, 0.9)

# ╔═╡ 50e00ab0-fe6a-11ea-169e-89ce2d5d5e8b
	md"""
	Filler
	"""

# ╔═╡ 50c508a0-fe6a-11ea-0686-8391123b8455

# ╔═╡ 50a54ba0-fe6a-11ea-25e9-b148070c37f5
	md"""
	Filler
	"""

# ╔═╡ 505bbf82-fe6a-11ea-183e-65470af0bffe

# ╔═╡ Cell order:
# ╠═676fcc70-fe6a-11ea-228a-8f3af1750c03
# ╟─ec7bfd20-fe57-11ea-11ff-ad7d1063c690
# ╠═3e2b6a10-02b9-11eb-2026-5b14f80e3d02
# ╠═ec64cba0-fe57-11ea-0a43-296a76f4e41f
# ╟─64efe220-02b9-11eb-33a9-05fc5c343994
# ╠═ec4ef9b0-fe57-11ea-3218-19d3eb36b629
# ╟─e72eb300-fe69-11ea-0295-235feccceff5
# ╠═ec39eb10-fe57-11ea-1a17-f5177a3a2bc4
# ╠═ec2578b0-fe57-11ea-21f9-3d06f046e6c7
# ╠═ec0dd1fe-fe57-11ea-2bff-8f9c00037c98
# ╟─df83edf0-02ba-11eb-3f62-b9af925b068d
# ╟─ebe03252-fe57-11ea-3150-d3945169079b
# ╟─eb561abe-fe57-11ea-083b-7926bd86e4e9
# ╠═51309bb0-fe6a-11ea-04ca-e3b11cbf31f1
# ╟─25be91a0-02be-11eb-0e62-65deaa359b84
# ╟─51180aa0-fe6a-11ea-38ed-e157030c325c
# ╠═d62a5da0-02bc-11eb-37f2-6d919583d571
# ╟─d60d86d0-02bc-11eb-2056-8b51aea35279
# ╟─d5f3bd3e-02bc-11eb-05d7-d32bb9a027fc
# ╟─d5d9cca0-02bc-11eb-048f-9320b39c74cd
# ╠═d5bd9212-02bc-11eb-222e-9bed9f96ac54
# ╟─d5a268f0-02bc-11eb-115b-13534cc0c95e
# ╠═d58284e0-02bc-11eb-2b6b-9fc3b786a3c9
# ╟─108d3f70-0363-11eb-1fae-6713f3fef863
# ╟─09aadb90-0363-11eb-165a-23513f34ae99
# ╟─0b2d32e0-02c0-11eb-00d7-171f805c6e2d
# ╟─0b139060-02c0-11eb-0a67-6920aee627ca
# ╟─0af4bdc0-02c0-11eb-30aa-21865c74c6f3
# ╟─0ad96d8e-02c0-11eb-1546-3534868d4894
# ╟─0abf7cf0-02c0-11eb-0003-fb78d7a76dc0
# ╟─0aa58c50-02c0-11eb-1608-d1ec9abe94c6
# ╠═0a8a3c20-02c0-11eb-15e2-6d41dab0deee
# ╠═0a6ca200-02c0-11eb-0619-c1f7c9c51fb3
# ╠═0a519ff0-02c0-11eb-31e9-433a6762996c
# ╠═b273ab10-031a-11eb-2068-0796513e64c7
# ╠═96c2e480-031a-11eb-29c3-b97cef4ae152
# ╟─0a3676d0-02c0-11eb-0510-e7536424f8b7
# ╠═0a192ace-02c0-11eb-1e3f-b5ff3970f610
# ╠═09c9ab40-02c0-11eb-2e14-294da8b643b4
# ╠═09ad70b0-02c0-11eb-1982-7f86a152b7b5
# ╟─09922080-02c0-11eb-2d95-83e9fbe98d7a
# ╠═09798f6e-02c0-11eb-09e0-3feb2b5a9048
# ╟─095bce40-02c0-11eb-309a-7b4fda24cd8c
# ╟─0940cc30-02c0-11eb-252c-37c67ae8008e
# ╠═21609f00-031c-11eb-2c4e-c7f936df961f
# ╟─2146ae60-031c-11eb-10cf-6d0f7a5519ff
# ╠═212b3722-031c-11eb-1002-c1d8a1e1cff2
# ╠═2110d150-031c-11eb-3ed4-4fc761239518
# ╠═20f53300-031c-11eb-346c-555291725aa3
# ╟─20d994b0-031c-11eb-2907-453907e22de5
# ╠═208f935e-031c-11eb-172e-bf3a2463720c
# ╠═207888f2-031c-11eb-1800-5b0d98072ef8
# ╠═205c4e60-031c-11eb-1a04-4ba53bbe3cc4
# ╠═2044cec2-031c-11eb-38c8-d3f806df0607
# ╠═202ba170-031c-11eb-1cb2-5926b842e1f1
# ╠═20135e80-031c-11eb-1e42-459731477aa8
# ╠═1ff83560-031c-11eb-19f0-83458298f9ec
# ╠═1fde44c0-031c-11eb-33c8-87df48b53597
# ╟─d527bab0-02bc-11eb-3ef2-61490bf35dab
# ╠═50fd2fa0-fe6a-11ea-22ea-93ec0773ebd3
# ╠═50e00ab0-fe6a-11ea-169e-89ce2d5d5e8b
# ╠═50c508a0-fe6a-11ea-0686-8391123b8455
# ╟─50a54ba0-fe6a-11ea-25e9-b148070c37f5
# ╠═505bbf82-fe6a-11ea-183e-65470af0bffe


text:
"""
    get_neg_binomial_prn(p, r, shape=1; seed=nothing)

Generate a random variable from a Negative Binomial(`p`, `r`) distribution. Optionally you can set a specific seed.

# Examples

```julia-repl
julia> get_neg_binomial_prn(.2, 2)
 16

julia> get_neg_binomial_prn(.2, 2, seed=42)
 6

```

# References

Walk, C. Handbook on statistical distributions for experimentalists. 2007.
"""
function get_neg_binomial_prn(p::Real, r::Int; seed::Union{Int, Nothing}=nothing)
    check_p(p)
    n = convert(Int, r/p * 10)
    U = bernoulli_rng(p, n, seed=seed)
    X = sum(cumsum(U, dims=1) .< r) + 1   # TODO CHECK THIS AGAINST ANOTHER METHOD (e.g. convolution)
    return X
end

"""
    neg_binomial_rng(p, r, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Negative Binomial(`p`, `r`) distribution. Optionally you can set a specific seed.

# Notes

The Negative Binomial distribution is given:

``f(x,p,r) = \\binom{x-1}{r-1} (1-p)^{x-r} p^r \\quad x = 0,1,\\dots, n``

# Examples

```julia-repl
julia> neg_binomial_rng(.5, 2)
1-element Vector{Float64}:
 3.0

julia> neg_binomial_rng(.5, 5, 5)
5-element Vector{Float64}:
  8.0
 10.0
  8.0
 13.0
 10.0

julia> neg_binomial_rng(.5, 2, (2,2))
2×2 Matrix{Float64}:
 3.0  4.0
 4.0  2.0

```

# References

Walk, C. Handbook on statistical distributions for experimentalists. 2007.
"""
function neg_binomial_rng(p::Real, r::Int, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    check_p(p)
    X = zeros(shape)
    X .= get_neg_binomial_prn.(p, r, seed=seed)
    return X
end

"""
    conv_neg_binomial_rng(p, r, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Negative Binomial(`p`, `r`) distribution. Optionally you can set a specific seed.

# Notes

The Negative Binomial distribution is given:

``f(x,p,r) = \\binom{x-1}{r-1} (1-p)^{x-r} p^r \\quad x = 0,1,\\dots, n``

Uses a convolution algorithm to generate random variables, which is slightly slower than [`neg_binomial_rng`](@ref).

# Examples

```julia-repl
julia> conv_neg_binomial_rng(.4, 5, 1)
1×1 Matrix{Int64}:
 8

julia> conv_neg_binomial_rng(.4, 5, 5)
5×1 Matrix{Int64}:
 11
 14
 8
 10
 13

julia> conv_neg_binomial_rng(.4, 5, (2,2))
2×2×1 Array{Int64, 3}:
[:, :, 1] =
 20  11
 7  10
```

# References

Law, A. Simulation modeling and analysis, 5th Ed. McGraw Hill Education, Tuscon, 2013.
"""
function conv_neg_binomial_rng(p::Real, r::Int, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    check_p(p)
    Y = geometric_rng(p, (shape..., r), seed=seed)
    X = sum(Y, dims=ndims(Y))
    return X
end


text:
using ICGPracticalSixSDW
using Test

@testset "ICGPracticalSixSDW.jl" begin
    # Write your tests here.
end


text:
using Test
using MianalyzerBackend
const Histogram = MianalyzerBackend.Histogram

@testset "Histogram" begin
  @testset ".fit" begin
    @test Histogram.fit([1, 2, 3, 4, 5]) == [1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0; 1 0 1 0 1 0 1 0 1]
  end
end


text:
"""
    probplot(fm::fittedEVA)

Probability plot
"""
function probplot(fm::fittedEVA)::Plot

    if getcovariatenumber(fm.model) > 0
        df = probplot_std_data(fm)
        plotTitle = "Residual Probability Plot"
    else
        df = probplot_data(fm)
        plotTitle = "Probability Plot"
    end

    return plot(df, x=:Model, y=:Empirical, Geom.point, Geom.abline(color="red", style=:dash),
        Guide.xlabel("Model"), Guide.ylabel("Empirical"), Guide.title(plotTitle),
        Theme(discrete_highlight_color=c->nothing))

end

"""
    qqplot(fm::fittedEVA)

Quantile-Quantile plot
"""
function qqplot(fm::fittedEVA)::Plot

    if getcovariatenumber(fm.model) > 0
        df = qqplot_std_data(fm)
        plotTitle = "Residual Quantile Plot"
    else
        df = qqplot_data(fm)
        plotTitle = "Quantile Plot"
    end

    return plot(df, x=:Model, y=:Empirical, Geom.point, Geom.abline(color="red", style=:dash),
        Guide.xlabel("Model"), Guide.ylabel("Empirical"), Guide.title(plotTitle),
        Theme(discrete_highlight_color=c->nothing))

end

"""
    returnlevelplot(fm::fittedEVA)

Return level plot
"""
function returnlevelplot(fm::fittedEVA)::Plot

    if getcovariatenumber(fm.model) > 0
        @warn "this graphic is not optimized for non-stationary model; plot ignored."
        return plot()
    else
        df = returnlevelplot_data(fm)
        l1 = layer(df, x=:Period, y=:Level,Geom.line, Theme(default_color="red", line_style=[:dash]))
        l2 = layer(df, x=:Period, y=:Data, Geom.point)
        return plot(l1,l2, Scale.x_log10, Guide.xlabel("Return Period"), Guide.ylabel("Return Level"),
            Guide.title("Return Level Plot"), Theme(discrete_highlight_color=c->nothing))
    end
end

"""
    histplot(fm::fittedEVA)

Histogram plot
"""
function histplot(fm::fittedEVA)::Plot

    if getcovariatenumber(fm.model) > 0
        df = histplot_std_data(fm)
        plotTitle = "Residual Density Plot"
    else
        df = histplot_data(fm)
        plotTitle = "Density Plot"
    end

    h = layer(df[:h], x = :Data, Geom.histogram(bincount=df[:nbin], density=true))
    d = layer(df[:d], x = :DataRange, y = :Density, Geom.line, Theme(default_color="red") )

    return plot(d,h, Coord.cartesian(xmin = df[:xmin], xmax = df[:xmax]), Guide.xlabel("Data"), Guide.ylabel("Density"), Guide.title(plotTitle))

end

"""
    diagnosticplots(fm::fittedEVA)

Diagnostic plots
"""
function diagnosticplots(fm::fittedEVA)::Gadfly.Compose.Context

    f1 = probplot(fm)
    f2 = qqplot(fm)
    f3 = histplot(fm)

    if getcovariatenumber(fm.model) > 0
        f4 = plot()
    else
        f4 = returnlevelplot(fm)
    end

    return gridstack([f1 f2; f3 f4])
end

"""
    mrlplot(y::Vector{<:Real}, steps::Int = 100)

Mean residual plot

See also [`mrlplot_data`](@ref).
"""
function mrlplot(y::Vector{<:Real}, steps::Int=100)::Plot

    df = mrlplot_data(y, steps)

    p = plot(df, x = :Threshold, y = :mrl, ymin = :lbound, ymax = :ubound,
        Geom.line, Geom.ribbon, Guide.ylabel("Mean Residual Life"))

    return p
end


text:
using AutoLandmarking
using Test
using JLD

@testset "AutoLM.jl" begin
    image_data, imnames = load_imgs(string(Base.@__DIR__, "/test_data"), (128,128,128), false)
    landmark_data = read_landmarks(string(Base.@__DIR__, "/test_data"), 10, "@1")
    solution_load_imgs = JLD.load(string(Base.@__DIR__, "/test_solution/solution_load_imgs.jld"))["solution"]
    solution_read_landmarks = JLD.load(string(Base.@__DIR__, "/test_solution/solution_read_landmarks.jld"))["solution"]
    @test solution_load_imgs == image_data
    @test solution_read_landmarks == landmark_data

    aligned, aligned_lm, reco = align_principal(image_data, landmark_data, 192)
    resized, lm_resized = resize_relevant(aligned, aligned_lm, 128)
    lm_om_volume = landmark_to_surface(resized, lm_resized, 10)
    solution_landmark_to_surface = JLD.load(string(Base.@__DIR__, "/test_solution/solution_landmark_to_surface.jld"))["solution"]
    @test solution_landmark_to_surface == lm_om_volume

    mirrored, lm_mirrored = mirror_vol(resized, lm_om_volume)
    solution_mirror_vol = JLD.load(string(Base.@__DIR__, "/test_solution/solution_mirror_vol.jld"))["solution"]
    @test solution_mirror_vol == mirrored

    flipped, lm_flipped = flip_volume_side(mirrored, lm_mirrored)
    solution_flip_volume_side = JLD.load(string(Base.@__DIR__, "/test_solution/solution_flip_volume_side.jld"))["solution"]
    @test solution_flip_volume_side == flipped

    tilted, lm_tilted = rotate_volumes(flipped, lm_flipped, 30)
    solution_rotate_volumes = JLD.load(string(Base.@__DIR__, "/test_solution/solution_rotate_volumes.jld"))["solution"]
    solution_lm_rotate_volumes = JLD.load(string(Base.@__DIR__, "/test_solution/solution_lms_rotate_volumes.jld"))["solution"]
    @test solution_lm_rotate_volumes == lm_tilted
    @test solution_rotate_volumes == tilted

    depthmaps = depth_map_all_sides(tilted)
    solution_depthmaps = JLD.load(string(Base.@__DIR__, "/test_solution/solution_depthmaps.jld"))["solution"]
    @test solution_depthmaps == depthmaps

end

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
file = "./src/Day2/input.txt"

using CSV, DataFrames

df = CSV.read(file, DataFrame, header = false)
rename!(df, [:Column1 => :Range, :Column2 => :Letter, :Column3 => :Password])
df.Min = parse.(Int, first.(split.(df.Range, "-")))
df.Max = parse.(Int, last.(split.(df.Range, "-")))
df.Letter = replace.(df.Letter, ":" =>"")

function checkpassword(password, letter, Min, Max)
    val = count( i -> (i == letter[1]), password)
    if (val >= Min) & (val <= Max)
        return 1
    else
        return 0
    end
end

vals = checkpassword.(df.Password, df.Letter, df.Min, df.Max)

sum(vals)

function checkpassword2(password, letter, id1, id2)
    l = letter[1] # string to char
    v1 = password[id1] == l ? 1 : 0
    v2 = password[id2] == l ? 1 : 0
    if v1 + v2 == 1
        return 1
    else
        return 0
    end
end

vals2 = checkpassword2.(df.Password, df.Letter, df.Min, df.Max)
sum(vals2)

## Sukera's solution

function parseInput(file)
    data = readlines(joinpath(dirname(@__FILE__), file))
    d = split.(data, ": ")
    map(d) do (x,pw)
        counts, char = split(x)
        low,high = parse.(Int, split(counts, '-'))
        ((low,high,char[1]),pw)
    end
end

function solve(input)
    count(input) do ((low,high,char),pw)
        low <= count(==(char),pw) <= high
    end
end

function solve2(input)
    count(input) do ((low,high,char),pw)
        (pw[low] == char) ⊻ (pw[high] == char)
    end
end

function main(file="test")
    data = parseInput(file)
    solve(data) |> println
    solve2(data) |> println
end

main("input.txt")


text:
export DummyConnection

using Sockets
using DataStructures

struct DummyConnection <: WebIO.AbstractConnection
    inbox::Queue{Any}

    DummyConnection() = new(Queue{Any}())
end

Base.take!(c::DummyConnection) = dequeue!(c.inbox)
Sockets.send(c::DummyConnection, data) = enqueue!(c.inbox, data)

"""
    open_window()

Open a window, optionally showing it when the BLINK_DEBUG environment variable
is set (to allow for seeing what happens in the Electron console).
"""
function open_window()
    if haskey(ENV, "BLINK_DEBUG")
        w = Window(Dict(:show => true))
        @js w localStorage.debug = "*"
        opentools(w)
        return w
    end
    return Window(Dict(:show => false))
end


text:
import MCMCChains: Chains

function create_a3d(noofsamples, noofvariables, noofchains)
   a3d = fill(0.0, noofsamples, noofvariables, noofchains)
   a3d
 end

 function insert_chain!(a3d, chain, posterior, trans)
   for i in 1:size(a3d, 1)
     a3d[i,:,chain] = inverse(trans, posterior[i])
   end
 end

 function insert_chain!(a3d, chain, posterior)
   for i in 1:size(a3d, 1)
     a3d[i,:,chain] = posterior[i, :]
   end
 end

 function create_mcmcchains(a3d, cnames;start=1)
   Chains(a3d, cnames; start=start)
 end

 function create_mcmcchains(a3d, cnames, sections::Dict{Symbol, Vector{String}};
   start=1)
   Chains(a3d, cnames, sections; start=start)
 end

"""
Convert DynamcHMC samples to a chain
* `posterior`: an array of NamedTuple consisting of mcmc samples
"""
function nptochain(posterior,tune)
    Np = length(vcat(posterior[1]...))+1 #include lf_eps
    Ns = length(posterior)
    a3d = Array{Float64,3}(undef,Ns,Np,1)
    ϵ=tune.ϵ
    for (i,post) in enumerate(posterior)
        temp = Float64[]
        for p in post
            push!(temp,values(p)...)
        end
        push!(temp,ϵ)
        a3d[i,:,1] = temp'
    end
    parameter_names = getnames(posterior)
    push!(parameter_names,"lf_eps")
    chns = MCMCChains.Chains(a3d,parameter_names,
        Dict(:internals => ["lf_eps"]))
    return chns
end

function getnames(post)
    nt = post[1]
    Np =length(vcat(nt...))
    parm_names = fill("",Np)
    cnt = 0
    for (k,v) in pairs(nt)
        N = length(v)
        if isa(v,Array)
            for i in 1:N
                cnt += 1
                parm_names[cnt] = string(k,"[",i,"]")
            end
        else
            cnt+=1
            parm_names[cnt] = string(k)
        end
    end
    return parm_names
end

export
  create_a3d,
  insert_chain!,
  nptochain,
  create_mcmcchains


text:
# ============================================================
# Output Methods
# ============================================================

function KitBase.plot_line(KS::AbstractSolverSet, uq::AbstractUQ, sol::AbstractSolution1D)

    x = deepcopy(KS.pSpace.x[1:KS.pSpace.nx])
    flowMean = zeros(KS.pSpace.nx, size(sol.w[1], 1))
    flowVar = similar(flowMean)

    for i = 1:KS.pSpace.nx
        flowMean[i, 1, 1] = mean(ctr[i].prim[1, :, 1], uq.op)
        flowMean[i, 1, 2] = mean(ctr[i].prim[1, :, 2], uq.op)
        flowVar[i, 1, 1] = var(ctr[i].prim[1, :, 1], uq.op)
        flowVar[i, 1, 2] = var(ctr[i].prim[1, :, 2], uq.op)
        for j = 2:4
            for k = 1:2
                flowMean[i, j, k] = mean(ctr[i].prim[j, :, k], uq.op)
                flowVar[i, j, k] = var(ctr[i].prim[j, :, k], uq.op)
            end
        end
        flowMean[i, 5, 1] = mean(lambda_tchaos(ctr[i].prim[5, :, 1], KS.gas.mi, uq), uq.op)
        flowVar[i, 5, 1] = var(lambda_tchaos(ctr[i].prim[5, :, 1], KS.gas.mi, uq), uq.op)
    end

    xlabel("X")
    ylabel("Expectation")
    legend("Ni", "Ne", "Ui", "Ue", "Ti", "Te")
    p1 = plot(x, flowMean[:, 1, 1])
    p1 = oplot(x, flowMean[:, 1, 2])
    p1 = oplot(x, flowMean[:, 2, 1])
    p1 = oplot(x, flowMean[:, 2, 2])
    p1 = oplot(x, flowMean[:, 5, 1])
    p1 = oplot(x, flowMean[:, 5, 2])
    display(p1)

    xlabel("X")
    ylabel("Expectation")
    legend("Ex", "Ey", "Ez", "Bx", "By", "Bz")
    p2 = plot(x, emMean[:, 1])
    p2 = oplot(x, emMean[:, 2])
    p2 = oplot(x, emMean[:, 3])
    p2 = oplot(x, emMean[:, 4])
    p2 = oplot(x, emMean[:, 5])
    p2 = oplot(x, emMean[:, 6])
    display(p2)

    xlabel("X")
    ylabel("Variance")
    legend("Ni", "Ne", "Ui", "Ue", "Ti", "Te")
    p3 = plot(x, flowVar[:, 1, 1])
    p3 = oplot(x, flowVar[:, 1, 2])
    p3 = oplot(x, flowVar[:, 2, 1])
    p3 = oplot(x, flowVar[:, 2, 2])
    p3 = oplot(x, flowVar[:, 5, 1])
    p3 = oplot(x, flowVar[:, 5, 2])
    display(p3)

    xlabel("X")
    ylabel("Variance")
    legend("Ex", "Ey", "Ez", "Bx", "By", "Bz")
    p4 = plot(x, emVar[:, 1])
    p4 = oplot(x, emVar[:, 2])
    p4 = oplot(x, emVar[:, 3])
    p4 = oplot(x, emVar[:, 4])
    p4 = oplot(x, emVar[:, 5])
    p4 = oplot(x, emVar[:, 6])
    display(p4)

end

function KitBase.write_jld(
    KS::AbstractSolverSet,
    ctr::AbstractArray{<:AbstractControlVolume,1},
    uq::AbstractUQ,
    t = 0::Real,
)

    strIter = string(t)
    fileOut = KS.outputFolder * "data/t=" * strIter * ".jld2"

    @save fileOut KS ctr uq t

end


text:
"""
    AbstractEvalMetric

The base class for all evaluation metrics. The sub-types should implement the following
interfaces:

* [`update!`](@ref)
* [`reset!`](@ref)
* [`get`](@ref)
"""
abstract AbstractEvalMetric

"""
    update!(metric, labels, preds)

Update and accumulate metrics.

# Arguments:
* `metric::AbstractEvalMetric`: the metric object.
* `labels::Vector{NDArray}`: the labels from the data provider.
* `preds::Vector{NDArray}`: the outputs (predictions) of the network.
"""
function update!{T <: AbstractEvalMetric}(metric :: T, labels :: Vector{NDArray}, preds :: Vector{NDArray})
  if length(labels) != length(preds)
    Base.warn_once(
      "The number of labels ($(length(labels))) does not correspond to the\
      number of outputs ($(length(preds))). The calculated metric might not be accuracte.")
  end
  for (label, pred) in zip(labels, preds)
    _update_single_output(metric, label, pred)
  end
end

"""
    reset!(metric)

Reset the accumulation counter.
"""
function reset!(metric :: AbstractEvalMetric)
  throw(MethodError(reset!, (typeof(metric),)))
end

import Base: get
"""
    get(metric)

Get the accumulated metrics.

Returns `Vector{Tuple{Base.Symbol, Real}}`, a list of name-value pairs.
For example, `[(:accuracy, 0.9)]`.
"""
function get(metric :: AbstractEvalMetric)
  throw(MethodError(get, (typeof(metric),)))
end

"""
    MultiMetric(metrics::Vector{AbstractEvalMetric})

Combine multiple metrics in one and get a result for all of them.

# Usage
To calculate both mean-squared error [`Accuracy`](@ref) and log-loss [`ACE`](@ref):
```julia
  mx.fit(..., eval_metric = mx.MultiMetric([mx.Accuracy(), mx.ACE()]))
```
"""
type MultiMetric <: mx.AbstractEvalMetric
    metrics :: Vector{mx.AbstractEvalMetric}
end

function update!(metric :: MultiMetric, labels :: Vector{NDArray}, preds :: Vector{NDArray})
    for m in metric.metrics
        update!(m, labels, preds)
    end
    return nothing
end

function reset!(metric :: MultiMetric)
    map(reset!, metric.metrics)
    return nothing
end

function get(metric :: MultiMetric)
    mapreduce(get, append!, metric.metrics)
end

"""
    Accuracy

Multiclass classification accuracy.

Calculates the mean accuracy per sample for softmax in one dimension.
For a multi-dimensional softmax the mean accuracy over all dimensions is calculated.
"""
type Accuracy <: AbstractEvalMetric
  acc_sum  :: Float64
  n_sample :: Int

  Accuracy() = new(0.0, 0)
end

function _update_single_output(metric :: Accuracy, label :: NDArray, pred :: NDArray)
  @nd_as_jl ro=(label,pred) begin
    # Samples are stored in the last dimension
    @assert size(label, ndims(label)) == size(pred, ndims(pred))

    if ndims(pred) == 4 # Multidimensional case
      # Reshape label to be of the same shape as pred.
      # Except for the third dimension where the predictions are stored.
      labels = reshape(label, size(pred, 1, 2)..., 1, size(pred, 4))

      for sample in 1:size(labels, 4)
        for j in 1:size(labels, 2)
          for i in 1:size(labels, 1)
            label = labels[i, j, 1, sample]
            klasses = view(pred, i, j, :, sample)
            klass = indmax(klasses) - 1 # Classes start at 0...k-1

            metric.acc_sum += klass == label
            metric.n_sample += 1
          end
        end
      end
    elseif ndims(pred) == 2 # 1-dimensional case
      for sample in 1:size(label, 1)
        klass = indmax(view(pred, :, sample)) - 1
        metric.acc_sum += klass == label[sample]
        metric.n_sample += 1
      end
    else
      error("Can't handle prediction with dimensions $(ndims(pred)).")
    end
  end
end

function get(metric :: Accuracy)
  return [(:accuracy, metric.acc_sum / metric.n_sample)]
end

function reset!(metric :: Accuracy)
  metric.acc_sum  = 0.0
  metric.n_sample = 0
end

"""
    MSE

Mean Squared Error. TODO: add support for multi-dimensional outputs.

Calculates the mean squared error regression loss in one dimension.
"""

type MSE <: AbstractEvalMetric
  mse_sum  :: Float64
  n_sample :: Int

  MSE() = new(0.0, 0)
end

function _update_single_output(metric :: MSE, label :: NDArray, pred :: NDArray)
  label = copy(label)
  pred  = copy(pred)

  n_sample = size(pred)[end]
  metric.n_sample += n_sample

  for i = 1:n_sample
    metric.mse_sum += (label[i] - pred[i])^2
  end
end

function get(metric :: MSE)
  return [(:MSE, metric.mse_sum / metric.n_sample)]
end

function reset!(metric :: MSE)
  metric.mse_sum  = 0.0
  metric.n_sample = 0
end

"""
    ACE

Calculates the averaged cross-entropy (logloss) for classification.
"""
type ACE <: AbstractEvalMetric
  ace_sum  :: Float64
  n_sample :: Int

  ACE() = new(0.0, 0)
end

function get(metric :: ACE)
  return [(:ACE, - metric.ace_sum / metric.n_sample)]
end

function reset!(metric :: ACE)
  metric.ace_sum = 0.0
  metric.n_sample = 0
end

function _update_single_output(metric :: ACE, label :: NDArray, pred :: NDArray)
  @nd_as_jl ro=(label,pred) begin
    # Samples are stored in the last dimension
    @assert size(label, ndims(label)) == size(pred, ndims(pred))
    if ndims(pred) == 4
      labels = reshape(label, size(pred, 1, 2)..., 1, size(pred, 4))
      for sample in 1:size(labels, 4)
        for j in 1:size(labels, 2)
          for i in 1:size(labels, 1)
            # Cross-entropy reduces to -(ln(p_1)*0 + ln(p_2)*1) for classification
            # Since we can only target labels right now this is the only thing we can do.
            target = Int(labels[i, j, 1, sample]) + 1 # klasses are 0...k-1 => julia indexing
            p_k = pred[i, j, target, sample]

            metric.ace_sum += log(p_k)
            metric.n_sample += 1
          end
        end
      end
    elseif ndims(pred) == 2 # 1-dimensional case
      for sample in 1:size(labels, 1)
        target = Int(labels[sample]) + 1
        p_k = pred[target, sample]
        metric.ace_sum += log(p_k)
        metric.n_sample += 1
      end
    else
      error("Can't handle prediction with dimensions $(ndims(pred)).")
    end
  end
end

"""
    MultiACE

Calculates the averaged cross-entropy per class and overall (see [`ACE`](@ref)).
This can be used to quantify the influence of different classes on the overall loss.
"""
type MultiACE <: AbstractEvalMetric
  aces  :: Vector{Float64}
  counts :: Vector{Int}

  MultiACE(nclasses) = new(Base.zeros(nclasses), Base.zeros(Int, nclasses))
end

function get(metric :: MultiACE)
  aces = [(Symbol("ACE_$(i-0)"), - metric.aces[i] / metric.counts[i]) for i in 1:length(metric.aces)]
  push!(aces, (:ACE, - Base.sum(metric.aces) / Base.sum(metric.counts)))
  return aces
end

function reset!(metric :: MultiACE)
  metric.aces = Base.zero(metric.aces)
  metric.counts = Base.zero(metric.counts)
end

function _update_single_output(metric :: MultiACE, label :: NDArray, pred :: NDArray)
  @nd_as_jl ro=(label,pred) begin
    # Samples are stored in the last dimension
    @assert size(label, ndims(label)) == size(pred, ndims(pred))

    if ndims(pred) == 4
      labels = reshape(label, size(pred, 1, 2)..., 1, size(pred, 4))
      for sample in 1:size(labels, 4)
        for j in 1:size(labels, 2)
          for i in 1:size(labels, 1)
            # Cross-entropy reduces to -(ln(p_1)*0 + ln(p_2)*1) for classification
            # Since we can only target labels right now this is the only thing we can do.
            target = Int(labels[i, j, 1, sample]) + 1 # klasses are 0...k-1 => julia indexing
            p_k = pred[i, j, target, sample]

            metric.aces[target] += log(p_k)
            metric.counts[target] += 1
          end
        end
      end
    elseif ndims(pred) == 2
      for sample in 1:size(label, 1)
        target = Int(label[sample]) + 1
        p_k = pred[target, sample]
        metric.aces[target] += log(p_k)
        metric.counts[target] += 1
      end
    else
      error("Can't handle prediction with dimensions $(ndims(pred)).")
    end
  end
end


text:
using MeshArrays, Plots, CSV, DataFrames, JLD

GCMGridSpec("LLC90")
GCMGridLoad()

"""
    CreateModelSamples(iFile,iVar,indx)

Example:
```
df, mp = CreateModelSamples(2,1,[])

df2, mp2 = CreateModelSamples(2,2,df.indx)
df = join(df,df2,on = [:indx,:lon,:lat])
mp = merge(mp,mp2)

save("mp.jld", mp)
mp2=load("mp.jld")

CSV.write("df.csv", df)
df2 = CSV.File("df.csv") |> DataFrame!
```
"""
function CreateModelSamples(iFile,iVar,indx)

    nSample=100

    #test cases
    if (iFile==1)
        InDir="occci-daily-llc90/"
        InFil="OC_CCI_L3S_Rrs_412_2010"
        InPrec=Float32
        varName="Rrs412"
    else #minimal test case
        InDir="GRID_LLC90/"
        InFil="Depth.data"
        InPrec=Float64
        varName="Depth"
    end
    #... later this should also treat model output
    #... and variable lists -> 1 csv file + n binaries

    #here should be a loops over
    #1) days
    #2) variables

    tmp=read_bin(InDir*InFil,InPrec)
    m=fill(1., MeshArrays.XC) + 0. * mask(MeshArrays.hFacC[:,:,1],NaN,0)
    tmp=m*tmp[:,:,1]
    tmp=mask(tmp,NaN,0)

    #extract random subset for each day
    lon=convert2gcmfaces(MeshArrays.XC)
    lat=convert2gcmfaces(MeshArrays.YC)
    df=convert2gcmfaces(tmp)
    isempty(indx) ? indx=rand(findall(df.>0.),nSample) : nothing
    #... store in a DataFrame(Row)
    #df = DataFrame(indx = indx, lon = lon[indx], lat = lat[indx], x = tmp[indx])
    df = DataFrame( [indx, lon[indx], lat[indx], df[indx] ],[:indx, :lon, :lat, Symbol(varName)])

    #get full fields for january 1
    mp=convert2array(tmp)
    mp=mp[:,62:239]
    mp=Dict( varName => circshift(mp,(142,0)))
    #... store to file
    #heatmap(tmp)

    #... add dataframe(row) to dataframe
    #... close loops

    #... write dataframe to file
    return df, mp

end

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
using AcuteML

# Example 3 - Template Rendering using Functions

# This method only uses functions that return string. You can build your desired string and call the function for rendering.

## create person function to store out html template
newTemplate("person", :function)

function person(;id, age, field, GPA, courses)

  # Build the taken courses section
  loopOut=""
  for course in courses
    loopOut = loopout * """ <taken-courses>$(course)</taken-courses>   """
  end

  # Append all the sections and variables together
  out = """
  <person id=$(id)>
    <age>$(age)</age>
    <study-field>$(field)</study-field>
    <GPA>$(GPA)</GPA>
    $loopout
  </person>
  """

  return out
end

# Call the function for rendering
out = person(
  id = "1",
  age = "24",
  field = "Mechanical Engineering",
  GPA = "4.5",
  courses = ["Artificial Intelligence", "Robotics"]
)

print(out)

# you can also write the output to a file:
file = open(filePath, "r"); print(file, out); close(file)

#################################################################
# Example 4 - Template Rendering using Files

# You can render variables into html/xml files. However, you can't have multiline control flow Julia code in this method.

# only to set path to current file
cd(@__DIR__)

# you can create a file and edit the file directly by using
newTemplate("person")

# Add the following html code to the generated html file
#=
<person id=$(id)>
  <age>$(age)</age>
  <study-field>$(field)</study-field>
  <GPA>$(GPA)</GPA>
  <taken-courses>$(courses[1])</taken-courses>
  <taken-courses>$(courses[2])</taken-courses>
</person>
=#

# Specify the template (or its path), and also the variables for rendering
out =render2file("person", false,
  id = 1,
  age = 24,
  field = "Mechanical Engineering",
  GPA = 4.5,
  courses = ["Artificial Intelligence", "Robotics"])

# you pass `true` as the 2nd argument to owerwrite person.html statically.


text:
#!/usr/bin/env julia

using DelimitedFiles: readdlm

import PyPlot
using LaTeXStrings
const plt = PyPlot
const mpl = PyPlot.matplotlib

struct Benchmark{Style}
    name     :: String
    filename :: String
    pyplot_style :: Style
    Benchmark(name, fname, style) = new{typeof(style)}(name, fname, style)
end

struct TimerData
    avg :: Vector{Float64}
    std :: Vector{Float64}
    min :: Vector{Float64}
    max :: Vector{Float64}
end

const MPI_TAG = Ref("IntelMPI.2019.9")

const STYLE_IDEAL = (color=:black, ls=:dotted, label="ideal")

function load_timings(bench::Benchmark, resolution)
    filename = joinpath("results", MPI_TAG[], "N$resolution", bench.filename)
    data = readdlm(filename, Float64, comments=true) :: Matrix{Float64}
    Nxyz = data[:, 1:3]
    @assert all(Nxyz .== resolution)
    procs = data[:, 4]
    proc_dims = data[:, 5:6]
    repetitions = data[:, 7]
    times = TimerData((data[:, j] for j = 8:11)...)
    (
        Nxyz = Nxyz,
        procs = procs,
        proc_dims = proc_dims,
        repetitions = repetitions,
        times = times,
    )
end

function plot_from_file!(ax, bench::Benchmark, resolution;
                         plot_ideal=false, error_bars=nothing)
    data = load_timings(bench, resolution)
    times = data.times
    t = times.avg

    ax.plot(data.procs, t; bench.pyplot_style...)

    colour = bench.pyplot_style.color
    if error_bars == :extrema
        ax.errorbar(data.procs, t; color=colour)
    elseif error_bars == :std
        δ = times.std ./ 2
        ax.fill_between(data.procs, t .- δ, t .+ δ; alpha=0.2, color=colour)
    end

    if plot_ideal
        plot_ideal_scaling!(ax, data, t)
        add_text_resolution!(ax, resolution, data.procs, t)
    end

    ax
end

function plot_ideal_scaling!(ax, data, t)
    p = data.procs
    t_ideal = similar(t)
    for n in eachindex(t)
        t_ideal[n] = t[1] * p[1] / p[n]
    end
    ax.plot(p, t_ideal; STYLE_IDEAL...)
end

function add_text_resolution!(ax, N, xs, ys)
    x = first(xs)
    y = first(ys)
    if N == 512
        kws = (ha = :left, va = :top)
        x *= 0.95
        y *= 0.65
    else
        kws = (ha = :right, va = :center)
        x *= 0.9
        y *= 1.02
    end
    ax.text(
        x, y, latexstring("$N^3");
        fontsize = "large",
        kws...
    )
end

function plot_lib_comparison!(ax, benchs, resolution)
    ax.set_xscale(:log, base=2)
    ax.set_yscale(:log, base=10)
    map(benchs) do bench
        plot_ideal = bench === first(benchs)
        plot_from_file!(
            ax, bench, resolution;
            plot_ideal,
            # error_bars = :std,
        )
    end
    ax
end

function legend_libs!(ax, benchs; with_ideal=false, outside=false)
    styles = Any[getfield.(benchs, :pyplot_style)...]
    labels = Any[getfield.(benchs, :name)...]
    if with_ideal
        push!(styles, STYLE_IDEAL)
        push!(labels, "Ideal")
    end
    kws = if outside
        (loc = "center left", bbox_to_anchor = (1.0, 0.5))
    else
        (loc = "lower left", )
    end
    draw_legend!(ax, styles, labels; frameon=false, kws...)
end

function draw_legend!(ax, styles, labels; kwargs...)
    leg = ax.get_legend()
    lines = [mpl.lines.Line2D(Float64[], Float64[]; style...)
             for style in styles]
    ax.legend(lines, labels; kwargs...)
    if leg !== nothing
        ax.add_artist(leg)
    end
    ax
end

# Wrap matplotlib's SVG writer.
struct SVGWriter{File<:IO} <: IO
    fh :: File
end

Base.isreadable(io::SVGWriter) = isreadable(io.fh)
Base.iswritable(io::SVGWriter) = iswritable(io.fh)
Base.isopen(io::SVGWriter) = isopen(io.fh)

function Base.write(io::SVGWriter, s::Union{SubString{String}, String})
    # We remove the image height and width from the header.
    # This way the SVG image takes all available space in browsers.
    p = "\"\\S+pt\""  # "316.8pt"
    pat = Regex("(<svg .*)height=$p (.*) width=$p")
    rep = replace(s, pat => s"\1\2")
    write(io.fh, rep)
end

function plot_timings()
    resolutions = (
        # 256,
        512,
        1024,
        2048,
    )
    style = (fillstyle = :none, ms = 7, markeredgewidth = 1.5, linewidth = 1.5)
    benchs = (
        Benchmark(
            "PencilFFTs (default)", "PencilFFTs_PI.dat",
            (style..., color = :C0, marker = :o, zorder = 5),
        ),
        Benchmark(
            "PencilFFTs (Alltoallv)", "PencilFFTs_PA.dat",
            (style..., color = :C2, marker = :s, zorder = 8,
             linestyle = :dashed, linewidth = 1.0),
        ),
        Benchmark(
            "P3DFFT", "P3DFFT2.dat",
            (style..., color = :C1, marker = :x, zorder = 4, markeredgewidth = 2),
        ),
    )

    fig = plt.figure(figsize = (6, 4.2) .* 1.1, dpi = 150)
    ax = fig.subplots()

    ax.set_xlabel("MPI processes")
    ax.set_ylabel("Time (milliseconds)")

    map(resolutions) do N
        plot_lib_comparison!(ax, benchs, N)
    end

    legend_libs!(ax, benchs, with_ideal=true, outside=false)

    # Show '1024' instead of '2^10'
    for axis in (ax.xaxis, ax.yaxis)
        axis.set_major_formatter(mpl.ticker.ScalarFormatter())
    end

    open("timing_comparison.svg", "w") do ff
        io = SVGWriter(ff)
        @time fig.savefig(io)
    end

    fig
end

plot_timings()


text:
function gen_iter(n, d)
    # based on https://twitter.com/evalparse/status/1107964924024635392
    iter = NTuple{n, Int}[]
    for x in partitions(d + n, n)
        x = x .- 1
        if all(x .<= d)
            ys = Set(multiset_permutations(x, n))
            for y in ys
                push!(iter, Tuple(y))
            end
        end
    end
    iter
end

function construct_iter_all(N::Int, order::Int)

    #Construct an ordered iterator going over all moments
    # sequentially in terms of order

    iters = NTuple{N, Int}[]
    for d in 0:order
        x = Base.sort(gen_iter(N, d), rev=true)
        push!(iters, x...)
    end

    iters

end

# Trim a string of form "(a, b, c, d, ...)" to "abcd..."
trim_key(expr) = filter(x -> !(isspace(x) || x == ')' || x== '(' || x==','), string(expr))

# Expand a symbolic expression (no binomial expansion)
expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) # apply distribution law
expand_mod = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) # distributes terms until no longer possible
flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) #
flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) #
expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod]))) # apply flatten and distribution until no longer possible

function define_μ(iter::AbstractVector, iv::Union{Sym,Num})

    indices = map(trim_key, iter)

    μs = OrderedDict{eltype(iter), Any}()
    for (i, idx) in enumerate(iter)
        if sum(idx) == 0
            μs[idx] = 1
        else
            sym_name = Symbol('μ', join(map_subscripts(indices[i])))
            sym_raw = Sym{FnType{Tuple{Any}, Real}}(sym_name)
            term_raw = Term{Real}(sym_raw, [iv])
            μs[idx] = setmetadata(term_raw, Symbolics.VariableSource,
                                  (:momentclosure, sym_name))
        end
    end

    μs

end

define_μ(N::Int, order::Int, iv::Union{Sym,Num}) = define_μ(construct_iter_all(N, order), iv)

function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order))
    @parameters t
    return define_μ(iter, value(t))
end

function define_M(iter::AbstractVector, iv::Union{Sym,Num})

    indices = map(trim_key, iter)

    Ms = OrderedDict{eltype(iter), Any}()
    for (i, idx) in enumerate(iter)
        if sum(idx) == 0
            Ms[idx] = 1
        elseif sum(idx) == 1
            Ms[idx] = 0
        else
            sym_name = Symbol('M', join(map_subscripts(indices[i])))
            sym_raw = Sym{FnType{Tuple{Any}, Real}}(sym_name)
            term_raw = Term{Real}(sym_raw, [iv])
            Ms[idx] = setmetadata(term_raw, Symbolics.VariableSource,
                                  (:momentclosure, sym_name))
        end
    end

    Ms

end

define_M(N::Int, order::Int, iv::Union{Sym,Num}) = define_M(construct_iter_all(N, order), iv)

function define_M(N::Int, order::Int, iter = construct_iter_all(N, order))
    @parameters t
    return define_M(iter, value(t))
end

function extract_variables(eqs::Array{Equation, 1}, μ, M=[])

    vars = vcat(values(μ)..., values(M)...)
    # extract variables from rhs of each equation
    eq_vars = unique(vcat(get_variables.(eqs)...))
    # need this as get_variables does not extract var from `Differential(t)(var(t))`
    diff_vars = [var_from_nested_derivative(eq.lhs)[1] for eq in eqs]
    # filter out the unique ones
    eq_vars = unique(vcat(eq_vars..., diff_vars...))
    # this should preserve the correct ordering

    intersect!(vars, eq_vars)

end

## Set of functions to deconstruct polynomial propensities ##

#=
    Consider a polynomial propensity function: a₁ = x^2 + c₁*y²*x,
    where x(t) and y(t) are the molecule numbers variables and c₁ is reaction parameter
    First we need to split a₁ into separate terms, i.e. x^2 and c₁*y²*x. Then we
    determine the independent multiplication factor in each term (1 and c₁). Finally,
    we obtain the power each variable is raised to in each term (x to power 2 in term 1;
    x to power 1 and y to power 2 in term 2). Having this information we can proceed
    in constructing raw moment equations. IF any propensity function is non-polynomial
    then the function `polynomial_propensities` will throw an error.
=#

function isconstant(expr::Union{Symbolic, Number}, vars::Vector, iv::Sym) :: Bool

    # Check that the given expression does NOT depend on the given variables `vars` (expr is constant wrt. vars)
    # A variable here is defined as a function of the independent variable `iv`, e.g,  X(t) is variable, where t
    # is the independent variable

    istree(expr) || return true

    args = arguments(expr)
    if length(args) == 1
        if isequal(args[1], iv)
            !isvar(expr, vars)
        else
            isconstant(args[1], vars, iv)
        end
    else
        all(arg -> isconstant(arg, vars, iv), args)
    end

end

isvar(x::Symbolic, vars::Vector) = any(var -> isequal(x, var), vars)

function extract_pwr!(powers::Vector, expr::Symbolic, smap::Dict, vars::Vector)
    args = arguments(expr)
    if isvar(args[1], vars) && args[2] isa Int && args[2] >= 0
        idx = smap[args[1]]
        if powers[idx] != 0
            error(args[1], " occurring multiple times in a monomial is unexpected. In ", expr)
        else
            powers[idx] = args[2]
        end
    else
        error("Unexpected term: ", expr)
    end
end

function extract_mul(expr::Symbolic, smap::Dict, vars::Vector, iv::Sym)

    powers = zeros(Int, length(smap))
    factors = []

    for arg in arguments(expr)
        if isconstant(arg, vars, iv)
            push!(factors, arg)
        elseif operation(arg) == ^
            extract_pwr!(powers, arg, smap, vars)
        elseif isvar(arg, vars)
            powers[smap[arg]] = 1
        else
            error("Unexpected operation: ", operation(arg), " in ", expr)
        end
    end

    factor = isempty(factors) ? 1 : prod(factors)
    powers, factor

end

function polynomial_propensities(a::AbstractArray, iv::Sym, smap::Dict)

    vars = [x for (x,_) in Base.sort(collect(smap), by=x->x[2])]
    N = length(vars)

    all_factors = copy.(fill([], size(a)))
    all_powers = copy.(fill(Vector{Int}[], size(a)))

    for (rind, expr) in enumerate(a)

        expr = expand(expr)

        if isconstant(expr, vars, iv)

            push!(all_factors[rind], expr)
            push!(all_powers[rind], zeros(Int, N))

        elseif isvar(expr, vars)

            push!(all_factors[rind], 1)
            push!(all_powers[rind], map(v -> isequal(expr, v), vars))

        elseif operation(expr) == ^ #Symbolics.Pow

            try
                powers = zeros(Int, N)
                extract_pwr!(powers, expr, smap, vars)
                push!(all_factors[rind], 1)
                push!(all_powers[rind], powers)
            catch e
                error("Propensity function ", expr, " is non-polynomial? \n" * string(e))
            end

        elseif operation(expr) == * #Symbolics.Mul

            try
                powers, factor = extract_mul(expr, smap, vars, iv)
                push!(all_powers[rind], powers)
                push!(all_factors[rind], factor)
            catch e
                error("Propensity function ", expr, " is non-polynomial?\n" * string(e))
            end

        elseif operation(expr) == + #Symbolics.Add

            for term in arguments(expr)
                if isconstant(term, vars, iv)
                    push!(all_factors[rind], term)
                    push!(all_powers[rind], zeros(Int, N))
                else
                    try
                        op = operation(term)
                        if op == ^
                            powers = zeros(Int, N)
                            extract_pwr!(powers, term, smap, vars)
                            push!(all_factors[rind], 1)
                            push!(all_powers[rind], powers)
                        elseif op == *
                            powers, factor = extract_mul(term, smap, vars, iv)
                            push!(all_factors[rind], factor)
                            push!(all_powers[rind], powers)
                        elseif isvar(term, vars)
                            powers = zeros(Int, N)
                            powers[smap[term]] = 1
                            push!(all_factors[rind], 1)
                            push!(all_powers[rind], powers)
                        else
                            error("Unexpected operation: ", op, " in ", term)
                        end
                    catch e
                        error("Propensity function ", expr, " is non-polynomial?\n" * string(e))
                    end
                end
            end

        else

            error("Expression $expr could not be parsed correctly!")

        end

    end

    max_power = maximum(sum.(vcat(all_powers...)))
    all_factors, all_powers, max_power

end

#=
function degree(p::Num, sym::AbstractVector)
    p = value(p)
    sym = Set(value.(sym))
    if p isa Number
        return 0
    elseif p ∈ sym
        return 1
    elseif p isa Symbolic
        return degree(p, sym)
    end
end
degree(p::Sym, sym::Set) = Int(p ∈ sym)
degree(p::Sym, sym::AbstractVector) = degree(p, Set(value.(sym)))
degree(p::Term, sym::Set) = Int(p ∈ sym)
degree(p::Term, sym::AbstractVector) = degree(p, Set(value.(sym)))
degree(p::Symbolics.Add, sym::AbstractVector) = degree(p, Set(value.(sym)))
degree(p::Symbolics.Mul, sym::AbstractVector) = degree(p, Set(value.(sym)))
degree(p::Symbolics.Pow, sym::AbstractVector) = degree(p, Set(value.(sym)))

poly_subs(ex::Union{Num, Number, Symbolics.Mul, Symbolics.Add}, subs::AbstractDict, ps::AbstractArray, flag::Bool = false) = poly_subs(ex, subs, Set(ps), flag)
poly_subs(ex::Num, subs::AbstractDict, ps::Set = Set(), flag::Bool = false) = poly_subs(value(expand(ex)), subs, ps, flag)
poly_subs(ex::Number, ::AbstractDict, ps::Set = Set(), ::Bool = false) = ex
poly_subs(ex::Sym, subs::AbstractDict, ::Set, ::Bool = false) = haskey(subs, ex) ? subs[ex] : ex
poly_subs(ex::Term, subs::AbstractDict, ::Set, ::Bool = false) = haskey(subs, ex) ? subs[ex] : ex
poly_subs(ex::Symbolics.Pow, subs::AbstractDict, ::Set, flag::Bool = false) = substitute(ex, subs)
poly_subs(ex::Symbolics.Add, subs::AbstractDict, ps::Set = Set(), flag::Bool = false) = sum(poly_subs(term, subs, ps, flag) for term in arguments(ex))

function poly_subs(ex::Symbolics.Mul, subs::AbstractDict, ps = Set(), flag::Bool = false)
    mono = 1
    coeff = ex.coeff
    if flag # flag controls whether there are variables not being substituted -> those are shifted to the coefficient
        for (base, exponent) in pairs(ex.dict)
            if base ∉ ps && haskey(subs, base) # very hacky but does the job ...
                mono *= base^exponent
            else
                coeff *= base^exponent
            end
        end
    else
        for (base, exponent) in pairs(ex.dict)
            if base ∉ ps
                mono *= base^exponent
            else
                coeff *= base^exponent
            end
        end
    end
    return haskey(subs, mono) ? coeff*subs[mono] : ex
end

function taylor_expand(f::Union{Num, Term, Symbolics.Mul, Symbolics.Add, Symbolics.Pow}, vars::AbstractVector, ref_point::AbstractVector, order::Int, aux_flag::Bool = true)
    @assert order >= 0 "Order of Taylor expansion must be non-negative"
    @assert length(ref_point) == length(vars) "Variables and reference point must have the same dimension"
    N = length(vars)
    if aux_flag
        Δ = vars .- ref_point
    else
        Δ = vars
    end
    sub_ref_point = Dict(vars[i] => ref_point[i] for i in 1:N)
    f_taylor = substitute(f, sub_ref_point)
    iter_all = filter(x -> sum(x) > 0, MomentClosure.construct_iter_all(N, order))
    for iter in iter_all
        f_taylor += substitute(derivative(f, vars, iter), sub_ref_point)*prod(Δ .^ iter)/MomentClosure.fact(iter)
    end
    return f_taylor
end

function derivative(f::Union{Num, Term, Symbolics.Mul, Symbolics.Add, Symbolics.Pow}, vars::AbstractVector, order::Tuple)
    for i in 1:length(vars)
        if order[i] != 0
            f = nth_differential(f, vars[i], order[i])
        end
    end
    return f
end

function nth_differential(f::Union{Num, Term, Symbolics.Mul, Symbolics.Add, Symbolics.Pow}, var::Union{Num, Sym, Term}, n::Int)
    for i in 1:n
        f = expand_derivatives(Differential(var)(f), true)
    end
    return f
end
=#


text:
include("Model.jl")
include("Routine.jl")
include("SystemState.jl")
include("DataDist.jl")

mutable struct System
    model::Model
    routine::CEM
    data::Vector{Array}
    times::Vector{Float64}
    state::SystemState

    function System(model::Model, folder::String; routine::CEM = CEM(), state::SystemState = SystemState())
        ds = loadData(folder)
        da = arData(ds)
        times = Float64.(da[1][:,1])
        new(model, routine, da, times, state)
    end
end

function Base.show(io::IO, s::System)
    print("\nSystem:", "\n\t", s.model, "\n\tRoutine: ", s.routine, "\n\tState: ", s.state)
end

function resetSystem!(sys::System)
    sys.state = SystemState()
end

include("Path.jl")
include("Ensemble.jl")

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
module InfAlleles
include("../src/aliases.jl")
include("../src/trial_result.jl")
include("../src/inf_alleles.jl")
include("../src/propsel.jl")
include("../src/watterson.jl")
include("../src/entropy.jl")
end

#using InfAlleles

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
@doc raw"""
    AbstractQuasiNewtonDirectionUpdate

An abstract representation of an Quasi Newton Update rule to determine the next direction
given current [`QuasiNewtonOptions`](@ref).

All subtypes should be functors, i.e. one should be able to call them as `H(M,x,d)` to compute a new direction update.
"""
abstract type AbstractQuasiNewtonDirectionUpdate end

@doc raw"""
    AbstractQuasiNewtonUpdateRule

Specify a type for the different [`AbstractQuasiNewtonDirectionUpdate`](@ref)s.
"""
abstract type AbstractQuasiNewtonUpdateRule end

@doc raw"""
    BFGS <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian BFGS update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{H}_k^\mathrm{BFGS}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{BFGS}_{k+1} = \widetilde{H}^\mathrm{BFGS}_k  + \frac{y_k y^{\mathrm{T}}_k }{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{H}^\mathrm{BFGS}_k s_k s^{\mathrm{T}}_k \widetilde{H}^\mathrm{BFGS}_k }{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{BFGS}_k s_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively.
"""
struct BFGS <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    InverseBFGS <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the inverse Riemanian BFGS update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{B}_k^\mathrm{BFGS}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{BFGS}_{k+1}  = \Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{s_k y^{\mathrm{T}}_k }{s^{\mathrm{T}}_k y_k}
\Bigr)
\widetilde{B}^\mathrm{BFGS}_k
\Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{y_k s^{\mathrm{T}}_k }{s^{\mathrm{T}}_k y_k}
\Bigr) + \frac{s_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively.

"""
struct InverseBFGS <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    DFP <: AbstractQuasiNewtonUpdateRule

indicates in an [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian DFP update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{H}_k^\mathrm{DFP}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{DFP}_{k+1} = \Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{y_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
\Bigr)
\widetilde{H}^\mathrm{DFP}_k
\Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{s_k y^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
\Bigr) + \frac{y_k y^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively.

"""
struct DFP <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    InverseDFP <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the inverse Riemanian DFP update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{B}_k^\mathrm{DFP}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{DFP}_{k+1} = \widetilde{B}^\mathrm{DFP}_k
+ \frac{s_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
- \frac{\widetilde{B}^\mathrm{DFP}_k y_k y^{\mathrm{T}}_k \widetilde{B}^\mathrm{DFP}_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{DFP}_k y_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively.

"""
struct InverseDFP <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    SR1 <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian SR1 update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{H}_k^\mathrm{SR1}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{SR1}_{k+1} = \widetilde{H}^\mathrm{SR1}_k
+ \frac{
  (y_k - \widetilde{H}^\mathrm{SR1}_k s_k) (y_k - \widetilde{H}^\mathrm{SR1}_k s_k)^{\mathrm{T}}
}{
(y_k - \widetilde{H}^\mathrm{SR1}_k s_k)^{\mathrm{T}} s_k
}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively.

This method can be stabilized by only performing the update if denominator is larger than
``r\lVert s_k\rVert_{x_{k+1}}\lVert y_k - \widetilde{H}^\mathrm{SR1}_k s_k \rVert_{x_{k+1}}``
for some ``r>0``. For more details, see Section 6.2 in [^NocedalWright2006]

[^NocedalWright2006]:
    > Nocedal, J., Wright, S.: Numerical Optimization, Second Edition, Springer, 2006.
    > doi: [10.1007/978-0-387-40065-5](https://doi.org/10.1007/978-0-387-40065-5)

# Constructor
    SR1(r::Float64=-1.0)

Generate the `SR1` update, which by default does not include the check (since the default sets ``t<0```)

"""
struct SR1 <: AbstractQuasiNewtonUpdateRule
    r::Float64
    SR1(r::Float64=-1.0) = new(r)
end

@doc raw"""
    InverseSR1 <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the inverse Riemanian SR1 update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{B}_k^\mathrm{SR1}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{SR1}_{k+1} = \widetilde{B}^\mathrm{SR1}_k
+ \frac{
  (s_k - \widetilde{B}^\mathrm{SR1}_k y_k) (s_k - \widetilde{B}^\mathrm{SR1}_k y_k)^{\mathrm{T}}
}{
  (s_k - \widetilde{B}^\mathrm{SR1}_k y_k)^{\mathrm{T}} y_k
}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively.

This method can be stabilized by only performing the update if denominator is larger than
``r\lVert y_k\rVert_{x_{k+1}}\lVert s_k - \widetilde{H}^\mathrm{SR1}_k y_k \rVert_{x_{k+1}}``
for some ``r>0``. For more details, see Section 6.2 in [^NocedalWright2006].

# Constructor
    InverseSR1(r::Float64=-1.0)

Generate the `InverseSR1` update, which by default does not include the check,
since the default sets ``t<0```.
"""
struct InverseSR1 <: AbstractQuasiNewtonUpdateRule
    r::Float64
    InverseSR1(r::Float64=-1.0) = new(r)
end

@doc raw"""
    Broyden <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian Broyden update is used in the Riemannian quasi-Newton method, which is as a convex combination of [`BFGS`](@ref) and [`DFP`](@ref).

We denote by ``\widetilde{H}_k^\mathrm{Br}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{Br}_{k+1} = \widetilde{H}^\mathrm{Br}_k
  - \frac{\widetilde{H}^\mathrm{Br}_k s_k s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k}{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k} + \frac{y_k y^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
  + φ_k s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k
  \Bigl(
        \frac{y_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{H}^\mathrm{Br}_k s_k}{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k}
  \Bigr)
  \Bigl(
        \frac{y_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{H}^\mathrm{Br}_k s_k}{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k}
  \Bigr)^{\mathrm{T}}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively, and ``φ_k`` is the Broyden factor which is `:constant` by default but can also be set to `:Davidon`.

# Constructor
    Broyden(φ, update_rule::Symbol = :constant)

"""
mutable struct Broyden <: AbstractQuasiNewtonUpdateRule
    φ::Float64
    update_rule::Symbol
end
Broyden(φ::Float64) = Broyden(φ, :constant)

@doc raw"""
    InverseBroyden <: AbstractQuasiNewtonUpdateRule

Indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian Broyden update
is used in the Riemannian quasi-Newton method, which is as a convex combination
of [`InverseBFGS`](@ref) and [`InverseDFP`](@ref).

We denote by ``\widetilde{H}_k^\mathrm{Br}`` the operator concatenated with a vector transport
and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{Br}_{k+1} = \widetilde{B}^\mathrm{Br}_k
 - \frac{\widetilde{B}^\mathrm{Br}_k y_k y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k}
   + \frac{s_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
 + φ_k y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k
 \Bigl(
     \frac{s_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{B}^\mathrm{Br}_k y_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k}
    \Bigr) \Bigl(
        \frac{s_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{B}^\mathrm{Br}_k y_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k}
 \Bigr)^{\mathrm{T}}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
\operatorname{grad}f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(\operatorname{grad}f(x_k)) ∈ T_{x_{k+1}} \mathcal{M},
```

respectively, and ``φ_k`` is the Broyden factor which is `:constant` by default but can also be set to `:Davidon`.

# Constructor
    InverseBroyden(φ, update_rule::Symbol = :constant)

"""
mutable struct InverseBroyden <: AbstractQuasiNewtonUpdateRule
    φ::Float64
    update_rule::Symbol
end
InverseBroyden(φ::Float64) = InverseBroyden(φ, :constant)

@doc raw"""
    QuasiNewtonOptions <: Options

These Quasi Newton [`Options`](@ref) represent any quasi-Newton based method and can be
used with any update rule for the direction.

# Fields
* `x` – the current iterate, a point on a manifold
* `gradient` – the current gradient
* `sk` – the current step
* `yk` the current gradient difference
* `direction_update` - an [`AbstractQuasiNewtonDirectionUpdate`](@ref) rule.
* `retraction_method` – an `AbstractRetractionMethod`
* `stop` – a [`StoppingCriterion`](@ref)

# See also
[`GradientProblem`](@ref)
"""
mutable struct QuasiNewtonOptions{
    P,
    T,
    U<:AbstractQuasiNewtonDirectionUpdate,
    SC<:StoppingCriterion,
    S<:Stepsize,
    RTR<:AbstractRetractionMethod,
    VT<:AbstractVectorTransportMethod,
} <: AbstractGradientOptions
    x::P
    gradient::T
    sk::T
    yk::T
    direction_update::U
    retraction_method::RTR
    stepsize::S
    stop::SC
    vector_transport_method::VT
end
function QuasiNewtonOptions(
    x::P,
    gradient::T,
    direction_update::U,
    stop::SC,
    stepsize::S;
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
) where {P,T,U<:AbstractQuasiNewtonDirectionUpdate,SC<:StoppingCriterion,S<:Stepsize}
    return QuasiNewtonOptions{
        P,T,U,SC,S,typeof(retraction_method),typeof(vector_transport_method)
    }(
        x,
        gradient,
        deepcopy(gradient),
        deepcopy(gradient),
        direction_update,
        retraction_method,
        stepsize,
        stop,
        vector_transport_method,
    )
end

@doc raw"""
    QuasiNewtonMatrixDirectionUpdate <: AbstractQuasiNewtonDirectionUpdate

These [`AbstractQuasiNewtonDirectionUpdate`](@ref)s represent any quasi-Newton update rule, where the operator is stored as a matrix. A distinction is made between the update of the approximation of the Hessian, ``H_k \mapsto H_{k+1}``, and the update of the approximation of the Hessian inverse, ``B_k \mapsto B_{k+1}``. For the first case, the coordinates of the search direction ``η_k`` with respect to a basis ``\{b_i\}^{n}_{i=1}`` are determined by solving a linear system of equations, i.e.

```math
\text{Solve} \quad \hat{η_k} = - H_k \widehat{\operatorname{grad}f(x_k)}
```

where ``H_k`` is the matrix representing the operator with respect to the basis ``\{b_i\}^{n}_{i=1}`` and ``\widehat{\operatorname{grad}f(x_k)}`` represents the coordinates of the gradient of the objective function ``f`` in ``x_k`` with respect to the basis ``\{b_i\}^{n}_{i=1}``.
If a method is chosen where Hessian inverse is approximated, the coordinates of the search direction ``η_k`` with respect to a basis ``\{b_i\}^{n}_{i=1}`` are obtained simply by matrix-vector multiplication, i.e.

```math
\hat{η_k} = - B_k \widehat{\operatorname{grad}f(x_k)}
```

where ``B_k`` is the matrix representing the operator with respect to the basis ``\{b_i\}^{n}_{i=1}`` and ``\widehat{\operatorname{grad}f(x_k)}`` as above. In the end, the search direction ``η_k`` is generated from the coordinates ``\hat{eta_k}`` and the vectors of the basis ``\{b_i\}^{n}_{i=1}`` in both variants.
The [`AbstractQuasiNewtonUpdateRule`](@ref) indicates which quasi-Newton update rule is used. In all of them, the Euclidean update formula is used to generate the matrix ``H_{k+1}`` and ``B_{k+1}``, and the basis ``\{b_i\}^{n}_{i=1}`` is transported into the upcoming tangent space ``T_{x_{k+1}} \mathcal{M}``, preferably with an isometric vector transport, or generated there.

# Fields
* `basis` – the basis.
* `matrix` – the matrix which represents the approximating operator.
* `scale` – indicates whether the initial matrix (= identity matrix) should be scaled before the first update.
* `update` – a [`AbstractQuasiNewtonUpdateRule`](@ref).
* `vector_transport_method` – an `AbstractVectorTransportMethod`

# See also
[`QuasiNewtonLimitedMemoryDirectionUpdate`](@ref)
[`QuasiNewtonCautiousDirectionUpdate`](@ref)
[`AbstractQuasiNewtonDirectionUpdate`](@ref)
"""
mutable struct QuasiNewtonMatrixDirectionUpdate{
    NT<:AbstractQuasiNewtonUpdateRule,
    B<:AbstractBasis,
    VT<:AbstractVectorTransportMethod,
    M<:AbstractMatrix,
} <: AbstractQuasiNewtonDirectionUpdate
    basis::B
    matrix::M
    scale::Bool
    update::NT
    vector_transport_method::VT
end
function QuasiNewtonMatrixDirectionUpdate(
    update::U,
    basis::B,
    m::M,
    ;
    scale::Bool=true,
    vector_transport_method::V=ParallelTransport(),
) where {
    U<:AbstractQuasiNewtonUpdateRule,
    M<:AbstractMatrix,
    B<:AbstractBasis,
    V<:AbstractVectorTransportMethod,
}
    return QuasiNewtonMatrixDirectionUpdate{U,B,V,M}(
        basis, m, scale, update, vector_transport_method
    )
end
function (d::QuasiNewtonMatrixDirectionUpdate{T})(
    p, o
) where {T<:Union{InverseBFGS,InverseDFP,InverseSR1,InverseBroyden}}
    return get_vector(
        p.M, o.x, -d.matrix * get_coordinates(p.M, o.x, o.gradient, d.basis), d.basis
    )
end
function (d::QuasiNewtonMatrixDirectionUpdate{T})(
    p, o
) where {T<:Union{BFGS,DFP,SR1,Broyden}}
    return get_vector(
        p.M, o.x, -d.matrix \ get_coordinates(p.M, o.x, o.gradient, d.basis), d.basis
    )
end

@doc raw"""
    QuasiNewtonLimitedMemoryDirectionUpdate <: AbstractQuasiNewtonDirectionUpdate

This [`AbstractQuasiNewtonDirectionUpdate`](@ref) represents the limited-memory Riemanian BFGS update, where the approximating  operator is represented by ``m`` stored pairs of tangent vectors ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}`` in the ``k``-th iteration.
For the calculation of the search direction ``η_k``, the generalisation of the two-loop recursion is used (see [^HuangGallivanAbsil2015]), since it only requires inner products and linear combinations of tangent vectors in ``T_{x_k} \mathcal{M}``. For that the stored pairs of tangent vectors ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}``, the gradient ``\operatorname{grad}f(x_k)`` of the objective function ``f`` in ``x_k`` and the positive definite self-adjoint operator

```math
\mathcal{B}^{(0)}_k[⋅] = \frac{g_{x_k}(s_{k-1}, y_{k-1})}{g_{x_k}(y_{k-1}, y_{k-1})} \; \mathrm{id}_{T_{x_k} \mathcal{M}}[⋅]
```

are used. The two-loop recursion can be understood as that the [`InverseBFGS`](@ref) update is executed ``m`` times in a row on ``\mathcal{B}^{(0)}_k[⋅]`` using the tangent vectors ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}``, and in the same time the resulting operator ``\mathcal{B}^{LRBFGS}_k [⋅]`` is directly applied on ``\operatorname{grad}f(x_k)``.
When updating there are two cases: if there is still free memory, i.e. ``k < m``, the previously stored vector pairs ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}`` have to be transported into the upcoming tangent space ``T_{x_{k+1}} \mathcal{M}``; if there is no free memory, the oldest pair ``\{ \widetilde{s}_{k−m}, \widetilde{y}_{k−m}\}`` has to be discarded and then all the remaining vector pairs ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m+1}^{k-1}`` are transported into the tangent space ``T_{x_{k+1}} \mathcal{M}``. After that we calculate and store ``s_k = \widetilde{s}_k = T^{S}_{x_k, α_k η_k}(α_k η_k)`` and ``y_k = \widetilde{y}_k``. This process ensures that new information about the objective function is always included and the old, probably no longer relevant, information is discarded.

# Fields
* `method` – the maximum number of vector pairs stored.
* `memory_s` – the set of the stored (and transported) search directions times step size ``\{ \widetilde{s}_i\}_{i=k-m}^{k-1}``.
* `memory_y` – set of the stored gradient differences ``\{ \widetilde{y}_i\}_{i=k-m}^{k-1}``.
* `ξ` – a variable used in the two-loop recursion.
* `ρ` – a variable used in the two-loop recursion.
* `scale` –
* `vector_transport_method` – a `AbstractVectorTransportMethod`

# See also
[`InverseBFGS`](@ref)
[`QuasiNewtonCautiousDirectionUpdate`](@ref)
[`AbstractQuasiNewtonDirectionUpdate`](@ref)

[^HuangGallivanAbsil2015]:
    > Huang, Wen and Gallivan, K. A. and Absil, P.-A., A Broyden Class of Quasi-Newton Methods for Riemannian Optimization,
    > SIAM J. Optim., 25 (2015), pp. 1660-1685.
    > doi: [10.1137/140955483](https://doi.org/10.1137/140955483)
"""
mutable struct QuasiNewtonLimitedMemoryDirectionUpdate{
    NT<:AbstractQuasiNewtonUpdateRule,T,VT<:AbstractVectorTransportMethod
} <: AbstractQuasiNewtonDirectionUpdate
    method::NT
    memory_s::CircularBuffer{T}
    memory_y::CircularBuffer{T}
    ξ::Vector{Float64}
    ρ::Vector{Float64}
    scale::Float64
    project::Bool
    vector_transport_method::VT
end
function QuasiNewtonLimitedMemoryDirectionUpdate(
    method::NT,
    ::T,
    memory_size::Int;
    scale::Bool=true,
    project=true,
    vector_transport_method::V=ParallelTransport(),
) where {NT<:AbstractQuasiNewtonUpdateRule,T,V<:AbstractVectorTransportMethod}
    return QuasiNewtonLimitedMemoryDirectionUpdate{NT,T,V}(
        method,
        CircularBuffer{T}(memory_size),
        CircularBuffer{T}(memory_size),
        zeros(memory_size),
        zeros(memory_size),
        scale,
        project,
        vector_transport_method,
    )
end
function (d::QuasiNewtonLimitedMemoryDirectionUpdate{InverseBFGS})(p, o)
    r = deepcopy(o.gradient)
    m = length(d.memory_s)
    m == 0 && return -r
    for i in m:-1:1
        d.ρ[i] = 1 / inner(p.M, o.x, d.memory_s[i], d.memory_y[i]) # 1 sk 2 yk
        d.ξ[i] = inner(p.M, o.x, d.memory_s[i], r) * d.ρ[i]
        r .= r .- d.ξ[i] .* d.memory_y[i]
    end
    r .= 1 / (d.ρ[m] * norm(p.M, o.x, last(d.memory_y))^2) .* r
    for i in 1:m
        r .= r .+ (d.ξ[i] - d.ρ[i] * inner(p.M, o.x, d.memory_y[i], r)) .* d.memory_s[i]
    end
    d.project && project!(p.M, r, o.x, r)
    return -r
end

@doc raw"""
    QuasiNewtonCautiousDirectionUpdate <: AbstractQuasiNewtonDirectionUpdate

These [`AbstractQuasiNewtonDirectionUpdate`](@ref)s represent any quasi-Newton update rule, which are based on the idea of a so-called cautious update. The search direction is calculated as given in [`QuasiNewtonMatrixDirectionUpdate`](@ref) or [`LimitedMemoryQuasiNewctionDirectionUpdate`]. But the update given in [`QuasiNewtonMatrixDirectionUpdate`](@ref) or [`LimitedMemoryQuasiNewctionDirectionUpdate`] is only executed if

```math
\frac{g_{x_{k+1}}(y_k,s_k)}{\lVert s_k \rVert^{2}_{x_{k+1}}} \geq \theta(\lVert \operatorname{grad}f(x_k) \rVert_{x_k}),
```

is satisfied, where ``\theta`` is a monotone increasing function satisfying ``\theta(0) = 0`` and ``\theta`` is strictly increasing at ``0``. If this is not the case, the corresponding update will be skipped, which means that for [`QuasiNewtonMatrixDirectionUpdate`](@ref) the matrix ``H_k`` or ``B_k`` is not updated, but the basis ``\{b_i\}^{n}_{i=1}`` is nevertheless transported into the upcoming tangent space ``T_{x_{k+1}} \mathcal{M}``, and for [`LimitedMemoryQuasiNewctionDirectionUpdate`] neither the oldest vector pair ``\{ \widetilde{s}_{k−m}, \widetilde{y}_{k−m}\}`` is discarded nor the newest vector pair ``\{ \widetilde{s}_{k}, \widetilde{y}_{k}\}`` is added into storage, but all stored vector pairs ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}`` are transported into the tangent space ``T_{x_{k+1}} \mathcal{M}``.
If [`InverseBFGS`](@ref) or [`InverseBFGS`](@ref) is chosen as update, then the resulting method follows the method of [^HuangAbsilGallivan2018], taking into account that the corresponding step size is chosen.

# Fields
* `update` – an [`AbstractQuasiNewtonDirectionUpdate`](@ref)
* `θ` – a monotone increasing function satisfying ``θ(0) = 0`` and ``θ`` is strictly increasing at ``0``.

# See also
[`QuasiNewtonMatrixDirectionUpdate`](@ref)
[`QuasiNewtonLimitedMemoryDirectionUpdate`](@ref)

[^HuangAbsilGallivan2018]:
    > Huang, Wen and Absil, P.-A and Gallivan, Kyle, A Riemannian BFGS Method Without Differentiated Retraction for Nonconvex Optimization Problems,
    > SIAM J. Optim., 28 (2018), pp. 470-495.
    > doi: [10.1137/17M1127582](https://doi.org/10.1137/17M1127582)
"""
mutable struct QuasiNewtonCautiousDirectionUpdate{U} <:
               AbstractQuasiNewtonDirectionUpdate where {
    U<:Union{QuasiNewtonMatrixDirectionUpdate,QuasiNewtonLimitedMemoryDirectionUpdate{T}}
} where {T<:AbstractQuasiNewtonUpdateRule}
    update::U
    θ::Function
end
function QuasiNewtonCautiousDirectionUpdate(
    update::U; θ::Function=x -> x
) where {
    U<:Union{QuasiNewtonMatrixDirectionUpdate,QuasiNewtonLimitedMemoryDirectionUpdate{T}}
} where {T<:AbstractQuasiNewtonUpdateRule}
    return QuasiNewtonCautiousDirectionUpdate{U}(update, θ)
end
(d::QuasiNewtonCautiousDirectionUpdate)(p, o) = d.update(p, o)

# access the inner vector transport method
function get_update_vector_transport(u::AbstractQuasiNewtonDirectionUpdate)
    return u.vector_transport_method
end
function get_update_vector_transport(u::QuasiNewtonCautiousDirectionUpdate)
    return get_update_vector_transport(u.update)
end


text:
"""
    onehot(lvl, x::Vector)

Create an indicator variable of `x` equal to `lvl`.
`lvl` and `x` can have missing values,
`ismissing` is used where necessary to ensure correct result.
"""
function onehot(lvl::Missing, x::Vector{T}) where T >: Missing
    return ismissing.(x)
end

function onehot(lvl::Missing, x::Vector)
    return zeros(Bool, size(x))
end

function onehot(lvl::T, x::Vector{Union{Missing,T}}) where T
    out = zeros(Bool, size(x))
    for (i,v) in enumerate(x)
        if !ismissing(v)
            out[i] = (v == lvl)
        end
    end

    return out
end

onehot(lvl::T, x::Vector{T}) where T = (x .== lvl)

"""
    onehot!(df::AbstractDataFrame, var::Symbol)

One hot encode unique values of `var` in `df`.
New variable name is constructed as `var_lvl`.
If `df.var` is string, it is standardized first via `strstd` before one hot encoding.
"""
function onehot!(df::AbstractDataFrame, var::Symbol)
    if eltype(df[!, var]) <: AbstractString
        vec = strstd(df[!, var])
    else
        vec = df[!, var]
    end

    for lvl in sort!(unique(vec))
        var_lvl = Symbol(var, '_', lvl)
        println("... creating $var_lvl")
        df[!, var_lvl] = onelvl(lvl, vec)
    end
end

"""
    onehot!(df::AbstractDataFrame, vars::Vector{Symbol})

One hot encode variables in `vars` in `df`.
"""
function onehot!(df::AbstractDataFrame, vars::Vector{Symbol})
    for var in unique(vars)
        onehot!(df, var)
    end
end


text:
module Wrap

using MacroTools: isexpr, prewalk
import Core: SSAValue, GotoNode
import Core.Compiler: CodeInfo, IRCode, CFG, ReturnNode,
  just_construct_ssa, compact!, OptimizationState, GotoIfNot, PhiNode, StmtRange

PhiNode(x, y) = PhiNode(Any[x...], Any[y...])

function Base.getindex(phi::PhiNode, b)
  j = findfirst(c -> c == b, phi.edges)
  return phi.values[j]
end

# SSA contruction (forked from Base for untyped code)

sparams(opt::OptimizationState) = VERSION > v"1.2-" ? Any[t.val for t in opt.sptypes] : Any[opt.sp...]

using ..Inner, ..IRTools
import ..Inner: IR, Variable, Statement, Branch, BasicBlock, TypedMeta, Meta, block!,
  unreachable, varmap, argument!, branch!, return!

vars(ex) = prewalk(x -> x isa SSAValue ? Variable(x.id) : x, ex)

Branch(x::GotoNode) = Branch(nothing, x.label, [])
Branch(x::GotoIfNot) = Branch(vars(x.cond), x.dest, [])
Branch(x::ReturnNode) = isdefined(x, :val) ? Branch(nothing, 0, [vars(x.val)]) : unreachable

function IRCode(meta::TypedMeta)
  opt = OptimizationState(meta.frame)
  Base.Meta.partially_inline!(meta.code.code, [], meta.method.sig, sparams(opt), 0, 0, :propagate)
  ir = just_construct_ssa(meta.code, copy(meta.code.code),
                          Int(meta.method.nargs)-1, opt)
  resize!(ir.argtypes, meta.method.nargs)
  return compact!(ir)
end

function branches_for!(ir, (from, to))
  brs = []
  for br in ir.blocks[from].branches
    br.block == to && push!(brs, br)
  end
  if isempty(brs)
    br = Branch(nothing, to, [])
    push!(ir.blocks[from].branches, br)
    push!(brs, br)
  end
  return brs
end

function rewrite_phis!(ir::IR, offset)
  for (v, st) in ir
    ex = st.expr
    ex isa PhiNode || continue
    to, = IRTools.Inner.blockidx(ir, v)
    bb = BasicBlock(to)
    push!(bb.args, v)
    push!(bb.argtypes, st.type)
    for (from, arg) in zip(ex.edges, ex.values), br in branches_for!(ir, from+offset=>to.id)
      push!(br.args, ex[from])
    end
    delete!(ir, v)
    ir.defs[v.id] = (to.id, -length(bb.args))
  end
  return ir
end

function IR(ir::IRCode)
  isempty(ir.new_nodes) || error("IRCode must be compacted")
  ir2 = IR(ir.linetable)
  defs = Dict()
  args = map(T -> argument!(ir2, T), ir.argtypes)
  offset = !isempty(ir.cfg.blocks[1].preds)
  offset && block!(ir2)
  for i = 1:length(ir.stmts)
    findfirst(==(i), ir.cfg.index) == nothing || block!(ir2)
    if ir.stmts[i] isa Union{GotoIfNot,GotoNode,ReturnNode}
      push!(ir2.blocks[end].branches, Branch(ir.stmts[i]))
    else
      x = push!(ir2, Statement(vars(ir.stmts[i]), ir.types[i], ir.lines[i]))
      defs[Variable(i)] = x
    end
  end
  ir2 = prewalk(x -> x isa Variable ? defs[x] :
                     x isa Core.Compiler.Argument ? Variable(x.n) :
                     x, ir2)
  rewrite_phis!(ir2, offset)
  return ir2 |> IRTools.renumber
end

unvars(ex) = prewalk(x -> x isa Variable ? SSAValue(x.id) : x, ex)

function IRCode(ir::IR)
  defs = Dict()
  stmts, types, lines = [], [], Int32[]
  index = Int[]
  for b in IRTools.blocks(ir)
    if b.id == 1
      for (i, arg) in enumerate(IRTools.arguments(b))
        defs[arg] = Core.Compiler.Argument(i)
      end
    else
      @assert isempty(BasicBlock(b).args)
    end
    for (v, st) in b
      defs[v] = Variable(length(stmts)+1)
      ex = varmap(x -> get(defs, x, x), st.expr) |> unvars
      push!(stmts, ex)
      push!(types, st.type)
      push!(lines, st.line)
    end
    for br in BasicBlock(b).branches
      if IRTools.isreturn(br)
        x = get(defs, br.args[1], br.args[1]) |> unvars
        push!(stmts, ReturnNode(x))
      elseif br == unreachable
        push!(stmts, ReturnNode())
      elseif br.condition == nothing
        push!(stmts, GotoNode(br.block))
      else
        cond = get(defs, br.condition, br.condition) |> unvars
        push!(stmts, GotoIfNot(cond, br.block))
      end
      push!(types, Any); push!(lines, 0)
    end
    push!(index, length(stmts)+1)
  end
  ranges = StmtRange.([1, index[1:end-1]...], index.-1)
  succs = map.(x -> x.id, IRTools.successors.(IRTools.blocks(ir)))
  preds = map.(x -> x.id, IRTools.predecessors.(IRTools.blocks(ir)))
  bs = Core.Compiler.BasicBlock.(ranges, preds, succs)
  cfg = CFG(bs, index)
  flags = [0x00 for _ in stmts]
  sps = VERSION > v"1.2-" ? [] : Core.svec()
  IRCode(stmts, types, lines, flags, cfg, ir.lines, ir.blocks[1].argtypes, [], sps)
end

function blockstarts(ci::CodeInfo)
  bs = Int[]
  terminator = false
  for i = 1:length(ci.code)
    ex = ci.code[i]
    if isexpr(ex, :gotoifnot)
      push!(bs, ex.args[2])
      terminator = true
    elseif isexpr(ex, GotoNode, :return)
      ex isa GotoNode && push!(bs, ex.label)
      i < length(ci.code) && push!(bs, i+1)
      terminator = false
    elseif terminator
      push!(bs, i)
      terminator = false
    end
  end
  return sort(unique(bs))
end

# TODO more informative names, while avoiding clashes.
slotname(ci, s) = Symbol(:_, s.id)

function IR(ci::CodeInfo, nargs::Integer; meta = nothing)
  bs = blockstarts(ci)
  ir = IR([ci.linetable...], meta = meta)
  _rename = Dict()
  rename(ex) = prewalk(ex) do x
    haskey(_rename, x) && return _rename[x]
    x isa Core.SlotNumber && return Inner.Slot(slotname(ci, x))
    return x
  end
  for i = 1:nargs
    _rename[Core.SlotNumber(i)] = argument!(ir)
  end
  for i = 1:length(ci.code)
    i in bs && block!(ir)
    ex = ci.code[i]
    if ex isa Core.NewvarNode
      continue
    elseif isexpr(ex, GotoNode)
      branch!(ir, findfirst(==(ex.label), bs)+1)
    elseif isexpr(ex, :gotoifnot)
      branch!(ir, findfirst(==(ex.args[2]), bs)+1,
              unless = rename(ex.args[1]))
    elseif isexpr(ex, :return)
      return!(ir, rename(ex.args[1]))
    else
      _rename[Core.SSAValue(i)] = push!(ir, IRTools.stmt(rename(ex), line = ci.codelocs[i]))
    end
  end
  return ir
end

function IR(meta::Union{Meta,TypedMeta}; slots = false, prune = true)
  if slots
    return IR(meta.code, meta.method.nargs, meta = meta)
  elseif meta isa Meta # TODO check this works for meta
    ir = IR(meta.code, meta.nargs, meta = meta) |> IRTools.ssa!
    if prune
      ir = ir |> IRTools.prune! |> IRTools.renumber
    end
    return ir
  else
    return IR(IRCode(meta))
  end
end

function IR(Ts::Type...; slots = false, prune = true)
  m = IRTools.meta(Tuple{Ts...})
  m == nothing && return
  IR(m, slots = slots, prune = prune)
end

end