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

using Markdown
using InteractiveUtils

# ╔═╡ c12ef471-40b9-4b3e-81b6-6d10fa4e5f04
begin
	bigbreak = HTML("<br>"^4);
	using PlutoUI
	using Test
end

# ╔═╡ 58863e50-8029-11ec-28da-35397c3e293f
md"""
# Loops
Topics:
1. `while` loops
2. `for` loops

## `while` loops
The syntax for a while is
```julia
while *condition*
    *loop body*
end
```
For example, we could use while to count or to iterate over an array.

"""

# ╔═╡ e75a73ad-c49a-4d0d-a548-e96ed41cc5cc
let
	n = 0
	while n < 10
	    n += 1
	    println(n)
	end
	n
end

# ╔═╡ bfe2dd76-54a3-479e-9f94-efaecde8266b
with_terminal() do
	myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

	i = 1
	while i <= length(myfriends)
	    friend = myfriends[i]
	    print("Hi $friend, it's great to see you!\n")
	    i += 1
	end
end

# ╔═╡ 35079b58-796c-4239-a201-cd6f55c16ca7
begin
end

# ╔═╡ ecd9ce51-6596-4a58-97d1-fbdee32ff1d8
bigbreak

# ╔═╡ 1e88489d-f123-4171-b406-47f6da914674
md"""
## `for` loops

The syntax for a for loop is
```julia
for *var* in *loop iterable*
    *loop body*
end
```
We could use a for loop to generate the same results as either of the examples above:
"""

# ╔═╡ 22ddbdaf-9483-417c-a702-317daa8c6537
with_terminal() do
	for n ∈ 1:10

		println(n)

	end
end

# ╔═╡ 3a38e1ef-7454-498d-af5d-876a9fa6e963
with_terminal() do
	myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

	for friend in myfriends
	    println("Hi $friend, it's great to see you!")
	end
end

# ╔═╡ c146bf2e-b890-482d-9b45-4fcfcf37eb93
Z = [10 30; 20 40]

# ╔═╡ 8e72a3b6-d276-46b9-b28e-3a9ab7530467
with_terminal() do

	for i in CartesianIndices(Z)
		println(i)
	end
end

# ╔═╡ f0390b57-a52e-4669-b992-2acd54f484da
bigbreak

# ╔═╡ f5a82c3e-47c6-4712-9c7e-7ff3fc1e5203
md"""
Now let's use for loops to create some addition tables, where the value of every entry is the sum of its row and column indices.

Note that we iterate over this array via column-major loops in order to get the best performance!

First, we initialize an array with zeros.

"""

# ╔═╡ 44c04d48-4e24-4034-b188-c92ba252ff22
begin
	m, n = 5, 5
	A = fill(0.0, (m, n))
end

# ╔═╡ 7a662dce-97fa-4c1c-ba60-51f6f147e808
begin
	for j in 1:n
	    for i in 1:m
	        A[i, j] = i + j
	    end
	end
	A
end

# ╔═╡ e2308376-31cd-4418-932c-8f166b4a316f
md"
Here's some syntactic sugar for the same nested for loop"

# ╔═╡ 0ae376b8-8451-46d7-9bb4-cc259a5b4557
B = fill(0, (m, n))

# ╔═╡ bf437047-2de0-4a5c-93d8-1ecab06dcda6
begin
	for j in 1:n, i in 1:m
	    B[i, j] = i + j
	end
	B
end

# ╔═╡ 4441001f-fd92-4f6b-9e5c-6daee5274d54
md"""The more "Julia" way to create this addition table would have been with an array comprehension."""

# ╔═╡ 2b9190e9-28c5-456f-97f4-b3b59b2d4957
C = [i + j for i in 1:m, j in 1:n]

# ╔═╡ 569ab1ca-623f-4569-b128-6a4fe4208d1c
bigbreak

# ╔═╡ d0ef189a-635e-4811-9c71-206e4a84fd81
md"""
### Generators

If you don't want to have the data in your memory, generators are your tool!
"""

# ╔═╡ e7afb4e9-f843-4113-ad67-f2eef26fa443
G = ( 10*i + j   for i in 1:5, j in 1:5    if i < j)

# ╔═╡ 5012a7a4-a37d-49df-9fcc-68a18295f5bd
G

# ╔═╡ bb5b8711-a987-4921-ac15-ffd3e71e903f
md"They can be iterated, which is enough for many utility functions, like
```julia
sum(), prod(), minimum(), maximum(), map(), broadcast(), ...
```
"

# ╔═╡ bcec978b-6da0-4720-a7fe-a65333233c2d
collect(G)

# ╔═╡ f19d75e2-a933-4534-b638-a724248671c4
sum(G)

# ╔═╡ 8d26881f-094c-42e6-95ea-93637646d210
minimum(G)

# ╔═╡ 50b79c47-09d1-41aa-8b60-5fc1dd114c29
map(x -> x^2, G)

# ╔═╡ d1782a41-8ea5-456c-a513-7fec9457c9ad
G .^ 2  # equivalent to the map command above

# ╔═╡ d0effe30-adf5-4cd2-9bd5-a99138a33e0e
with_terminal() do
	for i in 1:10
		print(i, " ")
	end
end

# ╔═╡ 99eefdf8-3ea5-41fb-92cc-5655b5e92fc2
bigbreak

# ╔═╡ f77c41cd-a36e-4ec0-968c-f1842b703032
md"""
## Exercises

4.1

Loop over integers between 1 and 100 and print their squares.
"""

# ╔═╡ 7021bc95-71bf-4861-98d2-a3f358e5aae6
with_terminal() do
	# write your code here
end

# ╔═╡ 33e4917e-ec98-44f5-85e4-b8f0c617f03c
bigbreak

# ╔═╡ e53d967c-d166-4208-a237-d4b692ab7431
md"""
4.2

Add to the code above a bit to create a dictionary, squares that holds integers and their squares as key, value pairs such that

```julia
squares[10] == 100
```
"""

# ╔═╡ 8a9b4ff0-1ebc-4546-9ac8-663c77d748c0
begin
	squares = missing
end

# ╔═╡ ebab307a-1a97-48da-89f2-940cb7936c6e
begin
	@assert squares[10] == 100
	@assert squares[11] == 121
end

# ╔═╡ d5ba15fa-255c-4bef-bb9d-d165232cd009
bigbreak

# ╔═╡ e7eeb88a-2e15-487a-b9f4-177b81d52b09
md"""
4.3

Use an array comprehension to create an an array `squares_arr` that stores the squares for all integers between 1 and 100.

"""

# ╔═╡ 5c13e9a8-491b-4486-ac35-c54ed91791c7
begin
	squares_arr = Int64[]
end

# ╔═╡ 01e7c0a1-9d4d-4f20-8625-6f8c08eeb418
begin
	@assert length(squares_arr) == 100
	@assert sum(squares_arr) == 338350
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[compat]
PlutoUI = "~0.7.32"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

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

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

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

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

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

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

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

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

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

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

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

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

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

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

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

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

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

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

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─c12ef471-40b9-4b3e-81b6-6d10fa4e5f04
# ╟─58863e50-8029-11ec-28da-35397c3e293f
# ╠═e75a73ad-c49a-4d0d-a548-e96ed41cc5cc
# ╠═bfe2dd76-54a3-479e-9f94-efaecde8266b
# ╠═35079b58-796c-4239-a201-cd6f55c16ca7
# ╟─ecd9ce51-6596-4a58-97d1-fbdee32ff1d8
# ╟─1e88489d-f123-4171-b406-47f6da914674
# ╠═22ddbdaf-9483-417c-a702-317daa8c6537
# ╠═3a38e1ef-7454-498d-af5d-876a9fa6e963
# ╠═c146bf2e-b890-482d-9b45-4fcfcf37eb93
# ╠═8e72a3b6-d276-46b9-b28e-3a9ab7530467
# ╟─f0390b57-a52e-4669-b992-2acd54f484da
# ╟─f5a82c3e-47c6-4712-9c7e-7ff3fc1e5203
# ╠═44c04d48-4e24-4034-b188-c92ba252ff22
# ╠═7a662dce-97fa-4c1c-ba60-51f6f147e808
# ╟─e2308376-31cd-4418-932c-8f166b4a316f
# ╠═0ae376b8-8451-46d7-9bb4-cc259a5b4557
# ╠═bf437047-2de0-4a5c-93d8-1ecab06dcda6
# ╟─4441001f-fd92-4f6b-9e5c-6daee5274d54
# ╠═2b9190e9-28c5-456f-97f4-b3b59b2d4957
# ╟─569ab1ca-623f-4569-b128-6a4fe4208d1c
# ╟─d0ef189a-635e-4811-9c71-206e4a84fd81
# ╠═e7afb4e9-f843-4113-ad67-f2eef26fa443
# ╠═5012a7a4-a37d-49df-9fcc-68a18295f5bd
# ╟─bb5b8711-a987-4921-ac15-ffd3e71e903f
# ╠═bcec978b-6da0-4720-a7fe-a65333233c2d
# ╠═f19d75e2-a933-4534-b638-a724248671c4
# ╠═8d26881f-094c-42e6-95ea-93637646d210
# ╠═50b79c47-09d1-41aa-8b60-5fc1dd114c29
# ╠═d1782a41-8ea5-456c-a513-7fec9457c9ad
# ╠═d0effe30-adf5-4cd2-9bd5-a99138a33e0e
# ╟─99eefdf8-3ea5-41fb-92cc-5655b5e92fc2
# ╟─f77c41cd-a36e-4ec0-968c-f1842b703032
# ╠═7021bc95-71bf-4861-98d2-a3f358e5aae6
# ╟─33e4917e-ec98-44f5-85e4-b8f0c617f03c
# ╟─e53d967c-d166-4208-a237-d4b692ab7431
# ╠═8a9b4ff0-1ebc-4546-9ac8-663c77d748c0
# ╠═ebab307a-1a97-48da-89f2-940cb7936c6e
# ╟─d5ba15fa-255c-4bef-bb9d-d165232cd009
# ╟─e7eeb88a-2e15-487a-b9f4-177b81d52b09
# ╠═5c13e9a8-491b-4486-ac35-c54ed91791c7
# ╠═01e7c0a1-9d4d-4f20-8625-6f8c08eeb418
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/3 - Loops.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/3 - Loops.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/3 - Loops.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
using LinearAlgebra

@testset "SimpleRectangle" begin
@testset "constructors" begin
    s = SimpleRectangle(1,2,3,4)
    @test typeof(s) == SimpleRectangle{Int}
    @test width(s) == 3
    @test height(s) == 4
    @test area(s) == 12
    @test xwidth(s) == 4
    @test yheight(s) == 6
    @test minimum(s) == Point(1,2)
    @test maximum(s) == Point(4,6)
    @test origin(s) == Point(1,2)
    h = HyperRectangle(s)
    @test typeof(h) == HyperRectangle{2,Int}
    @test h == HyperRectangle(1,2,3,4)
    h = HyperRectangle{3,Int}(s)
    @test h == HyperRectangle(1,2,0,3,4,0)

    @test SimpleRectangle(Vec(2,3)) == SimpleRectangle(0,0,2,3)
end

@testset "setops" begin
    r1 = SimpleRectangle(0,0,2,2)
    r2 = SimpleRectangle(2,2,2,2)
    r3 = SimpleRectangle(1,1,2,2)
    r4 = SimpleRectangle(0.5,0.5,3.,3.)
    @test intersect(r1,r2) == SimpleRectangle(2,2,0,0)
    @test intersect(r1,r3) == SimpleRectangle(1,1,1,1)
    @test intersect(r1,r4) == SimpleRectangle(0.5,0.5,1.5,1.5)
end

@testset "indexing" begin
    r = SimpleRectangle(0,0,2,2)
    a = Matrix(1.0I, 4, 4)
    @test a[r] == Matrix(1.0I, 2, 2)
    a[r] = Matrix(I, 2, 2)*2
    @test a == [2 0 0 0; 0 2 0 0; 0 0 1 0; 0 0 0 1]
end

end

ext: jl


text:
using MCMCBenchmarks, Test, Random

@testset "Poisson Regression Tests " begin
    Random.seed!(641)
    path = pathof(MCMCBenchmarks)
    include(joinpath(path,
      "../../Models/Hierarchical_Poisson/Hierarhical_Poisson_Models.jl"))
    Nd = 1
    Nreps = 1
    Ns = 10
    a0 = 1.0
    a1 = .5
    a0_sig = .3
    ProjDir = @__DIR__
    cd(ProjDir)
    samplers = (
      CmdStanNUTS(CmdStanConfig, ProjDir),
      AHMCNUTS(AHMCpoisson, AHMCconfig),
      #DHMCNUTS(sampleDHMC)
      )
    options = (Nsamples=2000, Nadapt=1000, delta=.8, Nd=Nd, Ns=Ns)
    results = benchmark(samplers, simulatePoisson, Nreps;options...)
    #Models run slow and these are noisy
    @test results[!,:a0_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ a0 atol = .6
    @test results[!,:a0_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ a0 atol = .6
    #@test results[!,:a0_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ a0 atol = .6
    @test results[!,:a1_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ a1 atol = .6
    @test results[!,:a1_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ a1 atol = .6
    #@test results[!,:a1_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ a1 atol = .6
    @test results[!,:a0_sig_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ a1 atol = .6
    @test results[!,:a0_sig_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ a1 atol = .6
    #@test results[!,:a0_sig_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ a1 atol = .6
    @test results[!,:a0_sampler_rhat][1] ≈ 1 atol = .03
    @test results[!,:a1_sampler_rhat][1] ≈ 1 atol = .03
    @test results[!,:a0_sig_sampler_rhat][1] ≈ 1 atol = .03
    isdir("tmp") && rm("tmp", recursive=true)
end


text:
function VectorizationBase.vectorizable(ScA::ChunkedArray{E,M,T,N,Np2}) where {E,M,T,N,Np2}
    VectorizedChunkedArray{E,M,T,N,Np2}(pointer(ScA.data), size(ScA.data))
end

const ChunkedVector{E,M,T} = ChunkedArray{E,M,T,1,2}
const ChunkedMatrix{E,M,T} = ChunkedArray{E,M,T,2,3}

function chunked_ind_base_quote(T)
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    quote
        $(Expr(:meta, :inline))
        i1 = 1+((i[1]-1) & $(W-1))
        i2 = 1+((i[1]-1) >> $Wshift)
    end
end

function CcA_setindex_quote(N, L, T)
    q = chunked_ind_base_quote(T)
    inds = [:(i[$j]) for j ∈ 2:N]
    for j ∈ 1:L
        push!(q.args, Expr(:call, :setindex!, :(ScA.data), :(v[LinearStorage(), $j]), :i1, inds..., j,  :i2))
    end
    q
end

@generated function Base.setindex!(ScA::ChunkedArray{E,M,T,N,Np2}, v::T,
                                   # i::Vararg{<:Integer,N}) where {T,E,M,N,Np2}
                                   i::Vararg{<:Integer,N}) where {E,M,T,N,Np2}
    CcA_setindex_quote(N, type_length(T), E)
end

@generated function Base.setindex!(ScA::ChunkedArray{E,M,T,N,Np2}, v::T,
                                    # i::CartesianIndex{N}) where {T,E,M,N,Np2}
                                    i::MultiDimIndex{N}) where {E,M,T,N,Np2}
    CcA_setindex_quote(N, type_length(T), E)
end

@generated function Base.getindex(ScA::ChunkedArray{E,M,T,1,3}, i::Integer, j::Integer) where {E,M,T}
    q = chunked_ind_base_quote(E)
    # inds = [:(i[$n]) for n in 2:N]
    ind_expr = [Expr(:call, :getindex, :(ScA.data), :i1, k, :i2) for k ∈ 1:type_length(T)]
    push!(q.args, construct_expr(T, ind_expr))
    q
end
@generated function Base.getindex(ScA::ChunkedArray{E,M,T,N,Np2}, i::Vararg{<:Integer,N}) where {E,M,T,N,Np2}
    q = chunked_ind_base_quote(E)
    inds = [:(i[$n]) for n in 2:N]
    ind_expr = [Expr(:call, :getindex, :(ScA.data), :i1, inds..., j, :i2) for j ∈ 1:type_length(T)]
    push!(q.args, construct_expr(T, ind_expr))
    q
end

@generated function Base.getindex(ScA::ChunkedArray{E,M,T,N,Np2},
                        i::MultiDimIndex{N}) where {E,M,T,N,Np2}
    q = chunked_ind_base_quote(E)
    inds = [:(i[$n]) for n in 2:N]
    ind_expr = [Expr(:call, :getindex, :(ScA.data), :i1, inds..., j, :i2) for j ∈ 1:type_length(T)]
    push!(q.args, construct_expr(T, ind_expr))
    q
end

@noinline function calc_stride(sv::Core.SimpleVector)
    N = length(sv)
    X = Vector{Int}(undef,N)
    X[1] = 1
    for n in 2:N
        X[n] = X[n-1] * (sv[n])::Int
    end
    X
end

@noinline function chunked_array_view_inds(V, N, T)
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    inds = Vector{Expr}(undef, N)
    if V.parameters[1] == Colon
        ind1 = :(i[1])
    else
        ind1 = :(ScA.view[1])
    end
    inds[1] = :( 1 + (($ind1-1) & $(W-1)) )
    i = 1
    for n ∈ 2:N
        if V.parameters[n] == Colon
            i += 1
            inds[n] = :(i[$i])
        else
            inds[n] = :(ScA.view[$n])
        end
    end
    inds, :(1 + (($ind1-1) >> $Wshift))
end

@generated function Base.getindex(ScA::ChunkedArrayView{E,M,T,new_N,N,Np2,V},
                        i::Vararg{<:Any,new_N}) where {E,M,T,new_N,N,Np2,V}
    inds, indslast = chunked_array_view_inds(V, N, E)
    # The strategy is to calculate the first index, as well as the stride between indices.
    # After the first index is calculated, the remaining indices can be calcualted much more
    # efficiently via the stide.
    q = quote
        $(Expr(:meta, :inline))
        stride = $(Expr(:call, :*, [:(ScA.full_size[$n]) for n ∈ 1:N]...))
        ind_1 = sub2ind(ScA.full_size, $(Expr(:tuple, inds..., 1, indslast)))
    end
    for j ∈ 2:type_length(T)
        push!(q.args, :($(Symbol(:ind_,j)) = ind_1 + stride * $(j-1)))
    end
    push!(q.args, construct_expr(T, [:(unsafe_load(ScA.ptr, $(Symbol(:ind_,j)))) for j ∈ 1:type_length(T)]))
    q
end
@generated function Base.setindex!(ScA::ChunkedArrayView{E,M,T,new_N,N,Np2,V}, v,
                        i::Vararg{<:Any,new_N}) where {E,M,T,new_N,N,Np2,V}
    inds, indslast = chunked_array_view_inds(V, N, E)
    q = quote
        $(Expr(:meta, :inline))
        stride = $(Expr(:call, :*, [:(ScA.full_size[$n]) for n ∈ 1:N]...))
        ind_1 = sub2ind(ScA.full_size, $(Expr(:tuple, inds..., 1, indslast)))
        unsafe_store!(ScA.ptr, v[LinearStorage(), 1], ind_1)
    end
    for j ∈ 2:type_length(T)
        push!(q.args, quote
            unsafe_store!(ScA.ptr, v[LinearStorage(), $j], ind_1 + stride * $(j-1))
        end)
    end
    push!(q.args, :(v))
    q
end
@generated function Base.setindex!(
                    ScA::ChunkedArrayView{E,M,T,new_N,N,Np2,V}, v,
                    i::Vararg{<:Any,new_N}
                ) where {E,M,S,T<:PaddedMatrices.AbstractConstantFixedSizeArray{S,E,M},new_N,N,Np2,V}
    inds, indslast = chunked_array_view_inds(V, N, E)
    q = quote
        $(Expr(:meta, :inline))
        stride = $(Expr(:call, :*, [:(ScA.full_size[$n]) for n ∈ 1:N]...))
        ind_1 = sub2ind(ScA.full_size, $(Expr(:tuple, inds..., 1, indslast)))
        # unsafe_store!(ScA.ptr, v[1], ind_1)
    end
    sv = S.parameters
    for (j,ci) ∈ enumerate(ntuple(i -> 1:sv[i], Val(M)))
        push!(q.args, quote
            unsafe_store!(ScA.ptr, v[$ci], ind_1 + stride * $(j-1))
        end)
    end
    push!(q.args, :(v))
    q
end
@inline function Base.view(ScA::ChunkedArray{E,M,T,N,Np2}, i::Vararg{<:Any,N}) where {E,M,T,N,Np2}
    ChunkedArrayView(ScA, i)
end
@generated function ChunkedArrayView(ScA::ChunkedArray{E,M,T,N,Np2,U}, i::V) where {E,M,T,N,Np2,V,U}
    new_N = 0
    new_size = Expr(:tuple)
    for n ∈ 1:N
        if V.parameters[n] == Colon
            new_N += 1
            push!(new_size.args, :(full_size[$n]))
        end
    end
    quote
        full_size = size(ScA.data)
        ChunkedArrayView{$E,$M,$T,$new_N,$N,$Np2,$V,$U}(pointer(ScA.data), $new_size, full_size, i, ScA.mask)
    end
end

# @inline function SArray{Tuple{S},T,N,L}(vs::Vararg{SIMDPirates.SVec{W,T},L}) where {S,T,N,L,W}
#     SArray{S}(vs...)
# end
@inline function Base.:+(i::Integer, v::VectorizedChunkedArray{E,M,T,N,Np2}) where {E,M,T,N,Np2}
    VectorizedChunkedArray{E,M,T,N,Np2}(v.ptr + sizeof(T) * i, v.fullsize)
end
@inline function Base.:+(v::VectorizedChunkedArray{E,M,T,N,Np2}, i::Integer) where {E,M,T,N,Np2}
    VectorizedChunkedArray{E,M,T,N,Np2}(v.ptr + sizeof(T) * i, v.fullsize)
end
@generated function SIMDPirates.vload(::Type{V}, vScA::VectorizedChunkedArray{E,M,T,N,Np2}) where {E,M,T,N,Np2,W, V <: Union{SIMDPirates.SVec{W,E},SIMDPirates.Vec{W,E}}}
    W_full, Wshift_full = VectorizationBase.pick_vector_width_shift(E)
    quote
        $(Expr(:meta,:inline))
        $(construct_expr(T, [Expr(:call, :vload, V, :(vScA.ptr + $(W_full*sizeof(E)*j))) for j ∈ 0:type_length(T)-1]))
    end
end
@generated function SIMDPirates.vload(
                    ::Type{V}, vScA::VectorizedChunkedArray{E,M,T,N,Np2}
                ) where {E,M,S,T<:PaddedMatrices.AbstractFixedSizeArray{S},N,Np2,W, V <: Union{SIMDPirates.SVec{W,E},SIMDPirates.Vec{W,E}}}
    W_full, Wshift_full = VectorizationBase.pick_vector_width_shift(E)
    @assert W == W_full
    sv = S.parameters
    R = calc_stride(sv)
    :(PaddedMatrices.PtrArray{$S,Vec{$W_full,$E},$(length(sv)),$R,$(PaddedMatrices.simple_vec_prod(sv)),true}(vScA.ptr))
end
@generated function vload!(
                A::PaddedMatrices.AbstractMutableFixedSizeArray{S,NTuple{W,Core.VecElement{E}},N3,P2,L2},
                vScA::VectorizedChunkedArray{E,M,ConstantFixedSizeArray{S,E,N3,P1,L1},N,Np2}
            ) where {E,M,N,Np2,W,S,P1,L1,P2,L2,N3}
    W_full, Wshift_full = VectorizationBase.pick_vector_width_shift(E)
    quote
        $(Expr(:meta,:inline))
        # $(construct_expr(T, [Expr(:call, :vload, V, :(vScA.ptr + $(W_full*sizeof(E)*j))) for j ∈ 0:type_length(T)-1]))
        ptrvScA = vScA.ptr
        ptrA = pointer(A)
        @inbounds for l ∈ 0:$(L2-1)
            SIMDPirates.vstore!(
                ptrA + $(W_full*sizeof(E))*l,
                SIMDPirates.vload(NTuple{$W,Core.VecElement{$E}}, ptrvScA + $(W_full*sizeof(E))*l)
            )
        end
        nothing
    end
end

function chunked_vload_quote(N, W, E, T, V)
    W_full, Wshift_full = VectorizationBase.pick_vector_width_shift(E)
    inds = [:(i[$n]) for n in 2:N]
    quote
        $(Expr(:meta,:inline))
        $(construct_expr(T, [Expr(:call, :vload, V, :(vScA.ptr),
            Expr(:call, :sub2ind, :(vScA.fullsize), Expr(:tuple, 1, inds..., j, 1 + :( (i[1]-1) >> $Wshift_full ) ))) for j in 1:type_length(T)]))
    end
end

@generated function SIMDPirates.vload(::Type{V}, vScA::VectorizedChunkedArray{E,M,T,N,Np2},
                                                # i::CartesianIndex{N}) where {E,M,T,N,Np2,W}
                                                i::Vararg{<:Integer,N}) where {E,M,T,N,Np2,W, V <: Union{SIMDPirates.SVec{W,E},SIMDPirates.Vec{W,E}}}
    chunked_vload_quote(N, W, E, T, V)
end
@generated function SIMDPirates.vload(::Type{V}, vScA::VectorizedChunkedArray{E,M,T,N,Np2},
                                                # i::Vararg{<:Integer,N}) where {E,M,T,N,Np2,W}
                                                i::MultiDimIndex{N}) where {E,M,T,N,Np2,W, V <: Union{SIMDPirates.SVec{W,E},SIMDPirates.Vec{W,E}}}
    chunked_vload_quote(N, W, E, T, V)
end
"""
For vload(::Type{<:SVec}, ::VectorizedChunkedArray), we assume the loads are alligned with chunk boundaries.
"""
@generated function SIMDPirates.vload(::Type{V}, vScA::VectorizedChunkedArray{E,M,T,N,Np2},
                                                # i::Vararg{<:Integer,N}) where {E,M,T,N,Np2,W}
                                                i::Integer) where {E,M,T,N,Np2,W, V <: Union{SIMDPirates.SVec{W,E},SIMDPirates.Vec{W,E}}}
    W_full, Wshift_full = VectorizationBase.pick_vector_width_shift(E)
    WE = W_full * sizeof(E)
    TL = type_length(T)
    quote
        $(Expr(:meta,:inline))
        ind = $(TL*sizeof(E)) * (i-1)
        $(construct_expr(T, [Expr(:call, :vload, V, :(vScA.ptr + ind + $(WE*j))) for j in 0:TL-1]))
    end
end

function Base.copyto!(ScA::ChunkedArray{E,M,T,N,Np2}, A::AbstractArray{T}) where {E,M,T,N,Np2}
    @boundscheck size(ScA) == size(A) || throwdimensionmismatcherror(size(ScA), size(A))
    @inbounds for i in eachindex(A)
        ScA[i] = A[i]
    end
    ScA
end

@generated function ChunkedArray(A::AbstractArray{T,N}) where {E,M,T <: AbstractArray{E,M},N}
    isbitstype(T) || thrownotisbitserror(T)
    TL = type_length(T)
    W, Wshift = VectorizationBase.pick_vector_width_shift(E)
    U = VectorizationBase.mask_type(T)
    quote
        L = size(A,1)
        mask = VectorizationBase.mask_from_remainder($T, L & $(W-1))
        ScA = ChunkedArray{$E,$M,$T,$N,$(N+2),$U}(
            Array{$E}(undef, $(Expr(:tuple, W, [:(size(A,$n)) for n ∈ 2:N]..., TL, :((L+$(W-1)) >> $Wshift) )) ),
            size(A), mask
        )
        @inbounds copyto!(ScA, A)
    end
end
@generated function ChunkedArray(A::AbstractArray{T,N}) where {E,M,S,T <: PaddedMatrices.AbstractFixedSizeArray{S,E,M},N}
    isbitstype(T) || thrownotisbitserror(T)
    TL = type_length(T)
    W, Wshift = VectorizationBase.pick_vector_width_shift(E)
    U = VectorizationBase.mask_type(T)
    sv = S.parameters
    X = calc_stride(sv)
    quote
        L = size(A,1)
        mask = VectorizationBase.mask_from_remainder($T, L & $(W-1))
        ScA = ChunkedArray{$E,$M,PaddedMatrices.ConstantFixedSizeArray{$S,$E,$M,$(Tuple{X...}),$(PaddedMatrices.simple_vec_prod(sv))},$N,$(N+2),$U}(
            Array{$E}(undef, $(Expr(:tuple, W, [:(size(A,$n)) for n ∈ 2:N]..., TL, :((L+$(W-1)) >> $Wshift) )) ),
            size(A), mask
        )
        @inbounds copyto!(ScA, A)
    end
end
@generated function ChunkedArray{E,M,T,N,Np2}(::UndefInitializer, I::Vararg{<:Integer,N}) where {E,M,T <: AbstractArray{E,M},N,Np2}
    isbitstype(T) || thrownotisbitserror(T)
    N + 2 == Np2 || throw("N + 1 = $(N + 1) != Np2 = $Np2")
    # @show T
    TL = type_length(T)
    W, Wshift = VectorizationBase.pick_vector_width_shift(E)
    U = VectorizationBase.mask_type(T)
    quote
        L = I[1]
        mask = VectorizationBase.mask_from_remainder($T, L & $(W-1))
        ChunkedArray{$E,$M,$T,$N,$Np2,$U}(
            Array{$E}(undef, $(Expr(:tuple, W, [:(I[$n]) for n ∈ 2:N]..., TL, :(I[1] >> $Wshift))) ),
            I, mask
        )
    end
end

# @generated function Base.size(ScA::ChunkedArray{E,M,T,N,Np2}) where {E,M,T,N,Np2}
#     quote
#         $(Expr(:meta, :inline))
#         s = size(ScA.data)
#         @inbounds $(Expr(:tuple, :(s[1]*s[end]), [:(s[$n]) for n ∈ 2:N]...))
#     end
# end
@inline Base.size(ScA::ChunkedArray) = ScA.size
@inline function Base.size(ScA::ChunkedArray{E,M,T,N}, n::Integer) where {E,M,T,N}
    @boundscheck n > N && ThrowBoundsError("$n > N, the number of dimensions.")
    @inbounds ScA.size[n]
end
# @inline Base.size(ScA::ChunkedArray{E,M,T,N,Np2}, n::Integer) where {E,M,T,N,Np2} = n == 1 ? size(ScA.data, 1)*size(ScA.data, Np2) : size(ScA.data, n)
@inline function Base.length(ScA::ChunkedArray{E,M,T,N,Np2}) where {E,M,T,N,Np2}
    s = ScA.size
    @inbounds l = s[1]
    @inbounds for n ∈ 2:N
        l *= s[n]
    end
    l
end


text:
function _svd_(A::AbstractSymMatrix)
    T = eltype(A)
    S = vtype(A)

    u_blocks = SortedDict{Sector{S, 2}, Matrix{T}}()
    s_blocks = SortedDict{Sector{S, 2}, Diagonal{Float64, Vector{Float64}}}()
    v_blocks = SortedDict{Sector{S, 2}, Matrix{T}}()

    midl = SortedDict{S, Int}()
    midr = SortedDict{S, Int}()
    fact = svd(zeros(T,1,1), full=false)
    for (sect,blk) in A.blocks
        c1, c2 = sect
        try
            fact = svd(blk, full=false)
        catch e
            @warn "error $e"
            fact = svd(blk, full=false, alg=LinearAlgebra.QRIteration())
        end
        u_blocks[Sector(c1, c1)] = fact.U
        s_blocks[Sector(c1, c2)] = Diagonal(fact.S)
        v_blocks[Sector(c2, c2)] = fact.Vt
        midl[c1] = length(fact.S)
        midr[c2] = length(fact.S)
    end

    Vl = VectorSpace{S}(midl, A.space[1].isdual)
    Vr = VectorSpace{S}(midr, A.space[2].isdual)

    u  = SymMatrix{S,T}(zero(S), (A.space[1], dual(Vl)), u_blocks)
    s  = SymDiagonal{S,Float64}(A.charge, (Vl, Vr), s_blocks)
    v = SymMatrix{S,T}(zero(S), (dual(Vr), A.space[2]), v_blocks)

    u, s, v
end

function _svd_(A::AbstractArray)
    u,s,v = svd(A, full=false)
    # note that iterating the component produces U,S,V instead of U,S,Vt!
    u, Diagonal(s), v'
end

"""
    fermionswapgate(l1, l2)

Fot the two given legs, `l1` and `l2` make the fermionic swap gate
which is `-1` if the two legs has charge odd and `+1` otherwise.

The output legs are l2_afterX, l1_afterX, l2, l1
"""
function fermionswapgate(l1::VectorSpace{S}, l2::VectorSpace{S}) where{S}

    space = (l2, l1, dual(l2), dual(l1))

    sects, sizes = _allsectorsandsizes(0, space)
    blocks = SortedDict{Sector{S, 4}, Array{Float64, 4}}()

    for index in eachindex(sects)
        c2, c1, c2X, c1X = sects[index]
        if c1 == c1X && c2 == c2X
            d2, d1 = sizes[index][1:2]
            s = isodd(c1) && isodd(c2) ? -1 : +1
            blocks[sects[index]] =  s * reshape(I(d1*d2), d2, d1, d2, d1)
        else
            blocks[sects[index]] = zeros(sizes[index])
        end
    end
    SymTensor(zero(S), space, blocks)
end

# function isrightisometry(A::SymTensor)
#     N = rank(A)

# end

function isrightisometry(A::SymMatrix)
    if eye(eltype(A), A.space[1]) ≈ contract(A, (1, -1), dual(A), (2, -1))
        return true
    else
        println(contract(A, (1, -1), dual(A), (2, -1)))
    end
    false
end

function isleftisometry(A::SymMatrix)
    eye(eltype(A), dual(A.space[2])) ≈ contract(dual(A), (-1, 1), A, (-1, 2))
end

isunitary(A::SymMatrix) = isrightisometry(A) && isleftisometry(A)

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
module Glutils

using ModernGL
using StaticArrays

include("utils.jl")
include("uniforms.jl")
include("shaders.jl")
include("buffers.jl")
include("textures.jl")

export @fragment_srt
export @geometry_srt
export @map
export @program_str
export @uniforms
export @uniforms!
export @vertex_srt
export @g_str

export arraybuffer
export bindbuffer
export bufferdata
export compile_shader
export compiled
export genbuffer
export getuniform
export indexbuffer
export linked
export program_infolog
export shader_infolog
export shaderiv
export sprogram
export SymString
export texture2D
export uniform!
export vertexarray
export vertexattrib

include("mat.jl"); using .Mat

# workaround (for vs-code) to get julia intellisense
# on files that are not included in the package's
# entry point
if (false)
    include("../examples/runexamples.jl")
end

end # Glutils


text:
# To do: Replace Length with any quantity, scale with scale_to_pt
# Write a macro to extend all curve-like functions which take input not a Point.
"""
    circle(p::Point, r::Length, action)
    circle(p::Point; r = missing, d = missing, action = :stroke)

    where

action can be :nothing, :fill, :stroke, :fillstroke, :fillpreserve, :strokepreserve, :clip.
    The default without keyword arguments is :nothing.
"""
circle(p::Point, r::Length, action) = circle(p, scale_to_pt(r), action)
function circle(p::Point; r = missing, d = missing, action = :stroke)
    @assert ismissing(r) + ismissing(d) == 1
    r = ismissing(d) ? r : d / 2
    circle(p, r, action)
end
"ellipse(p::Point, w::Q, h::Q; action = :stroke)  where {Q<:Quantity}"
function ellipse(p::Point, w::Q, h::Q; action = :stroke)  where {Q<:Quantity}
    ellipse(p, scale_to_pt(w), scale_to_pt(h), action)
end
"squircle(center::Point, hradius::Q, vradius::Q; action=:stroke, kwargs...)  where {Q<:Quantity}"
function squircle(center::Point, hradius::Q, vradius::Q; action = :stroke, kwargs...) where {Q<:Quantity}
    squircle(center, scale_to_pt(hradius), scale_to_pt(vradius), action;  kwargs...)
end

"rect(p::Point, w::Q, h::Q; action=:stroke, kwargs...)  where {Q<:Quantity}"
function rect(p::Point, w::Q, h::Q; action=:stroke, kwargs...)  where {Q<:Quantity}
    rect(p, scale_to_pt(w), scale_to_pt(h), action; kwargs...)
end
#= Full list of curve functions:
 circle, center3pts, ellipse, squircle,
arc, carc, arc2r, carc2r, sector, pie, curve, circlepath,
hypotrochoid, epitrochoid, spiral, intersection2circles,
intersectioncirclecircle, circlepointtangent
=#

"""
    trace_diminishing(center, vx, vy)

center is local origo, a Point, typically equal to global O.
vx, vy are vectors of corresponding coordinate quantities.

Plots a polyline until NaN or end of vectors.
The first segment is drawn with opacity 1, decreasing to one step before translucent(zero)
at the last segment
"""
function trace_diminishing(center, vx, vy)
    gsave()
    n = length(vx)
    opacity = 1.0
    xo, yo  = vx[1], vy[1]
    for (x, y) in zip(vx[2:end], vy[2:end])
        isnan(x) && break
        isnan(y) && break
        setopacity(opacity)
        move(center + (xo, yo))
        line(center + (x, y))
        do_action(:stroke)
        xo, yo = x, y
        opacity -= 1.0 / (n - 1)
    end
    grestore()
end

"""
    trace_rotate_hue(center, vx, vy; rotatehue_degrees_total = 270°)

center is local origo, a Point, typically equal to global O.
vx, vy are vectors of corresponding coordinate quantities.
rotatehue_degrees_total = 270° is the amount of hue rotation from the current
hue. It varies gradually. The hue is reset afterwards. If 360°, the starting
color and the ending color will be the same.

Plots a polyline until NaN or end of vectors.
The first segment is drawn with opacity 1, decreasing to one step before translucent(zero)
at the last segment
"""
function trace_rotate_hue(center, vx, vy; rotatehue_degrees_total = 270°)
    gsave()
    n = length(vx)
    startcolo = get_current_RGB()
    xo, yo  = vx[1], vy[1]
    count = 0
    for (x, y) in zip(vx[2:end], vy[2:end])
        isnan(x) && break
        isnan(y) && break
        move(center + (xo, yo))
        line(center + (x, y))
        do_action(:stroke)
        count += 1
        rotatedeg = count * rotatehue_degrees_total / n
        sethue(rotate_hue(startcolo, rotatedeg))
        xo, yo = x, y
    end
    grestore()
end

setline(x::Length) = setline(scale_to_pt(x))


text:
"""
Returns a bitarray for sites to be selected from the first (false) and second (true) homologous chromosome.
"""
function crossing_overs(nsites::Int, ncrossing_overs::Int)
  output = falses(nsites)
  if ncrossing_overs == 0
    return output
  elseif ncrossing_overs ≥ nsites
    output[1:2:end] .= true
    return output
  end
  breaking_points = sample(1:(nsites-1), ncrossing_overs, replace=false, ordered=true)
  last = true
  counter = 0
  for site in 1:nsites
    if counter < ncrossing_overs && site == breaking_points[counter+1]
      counter += 1
      output[site] = last
      last = !last
    else
      output[site] = last
    end
  end
  return output
end

"""
Returns gametes for epistasisMat, pleiotropyMat, and q.

A ametes includes `cross_overs` sites from one homologous chr and the rest from another.the corresponding column of the `epistasisMat` and `pleiotropyMat` matrices.
Each gamete is half of `epistasisMat` and `pleiotropyMat` (column-wise).
"""
function create_gamete(agent, cross_overs, nsites, first::Bool)
  indices1 = 1:nsites
  indices2 = indices1 .+ nsites
  if !first
    indices1, indices2 = indices2, indices1
  end
  epistasisMat_gamete = agent.epistasisMat[:, indices1]
  epistasisMat_gamete[:, cross_overs] .= agent.epistasisMat[:, indices2][:, cross_overs]

  pleiotropyMat_gamete = agent.pleiotropyMat[:, indices1]
  pleiotropyMat_gamete[:, cross_overs] .= agent.pleiotropyMat[:, indices2][:, cross_overs]

  q_gamete = agent.q[indices1]
  q_gamete[cross_overs] .= agent.q[indices2][cross_overs]

  return epistasisMat_gamete, pleiotropyMat_gamete, q_gamete
end

"""
Adds new individual(s) to the model as offsprings of `ag1` and `ag2`.
"""
function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
  species = ag1.species
  nsites = model.ngenes[species]

  if model.recombination[species] == 0
    nco1, nco2 = 0, 0
  else
    nco1, nco2 = rand(model.recombination[species], 2)
  end

  cross_overs1 = crossing_overs(nsites, nco1)
  cross_overs2 = crossing_overs(nsites, nco2)
  gametes1 = create_gamete(ag1, cross_overs1, nsites, rand((true, false)))
  gametes2 = create_gamete(ag2, cross_overs2, nsites, rand((true, false)))

  sex = rand((true, false))
  interaction_history = MVector{model.nspecies, Int}(fill(0, model.nspecies))

  # Merge gametes
  nsites2 = nsites*2
  episMat = MMatrix{nsites2, nsites2}(hcat(gametes1[1], gametes2[1]))
  pleioMat = MMatrix{model.nphenotypes[species],nsites2}(hcat(gametes1[2], gametes2[2]))
  q = MVector{nsites2}(vcat(gametes1[3], gametes2[3]))
  abph = get_abiotic_phenotype(species, episMat, pleioMat, q, model)
  bph = get_biotic_phenotype(species, episMat, pleioMat, q, model)
  W = abiotic_fitness(abph, species, ag1.pos, model)
  W = adjust_fitness(W, ag1.species, model)
  initial_energy = model.initial_energy[species]

  offspring = add_agent!(ag1.pos, model, ag1.species, bph, abph, episMat, pleioMat, q, 0, sex, interaction_history, initial_energy, W, true)
  return offspring
end

"""
Sexual reproduction for diploid individuals
"""
function reproduce!(ag1::Ind, ag2::Ind, model::ABM)
  reproduction_success = 1.0 - phenotypic_distance(ag1, ag2, model)
  growth_rate = model.growthrates[ag1.species]
  nchildren = rand(Poisson(reproduction_success * growth_rate))
  for c in 1:nchildren
    offspring = create_one_offspring(ag1, ag2, model)
    mutate!(offspring, model)
  end
end

"""
Asexual reproduction for the haploid
"""
function reproduce!(agent::Ind, model::ABM)
  if model.ploidy[agent.species] == 1 && in_reproduction_age(agent, model)
    growth_rate = model.growthrates[agent.species]
    W = agent.W >= 0.0 ? agent.W : 0.0
    nchildren = rand(Poisson(growth_rate * W))
    nchildren == 0 && return
    interaction_history = deepcopy(agent.interaction_history)
    interaction_history[1:end] .= 0
    for c in 1:nchildren
      offspring = add_agent!(agent.pos, model, agent.species, agent.biotic_phenotype, agent.abiotic_phenotype, agent.epistasisMat, agent.pleiotropyMat, agent.q, 0, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true)
      mutate!(offspring, model)
    end
  end
end

function in_reproduction_age(agent, model)
  agent.age >= model.repro_start[agent.species] && agent.age <= model.repro_end[agent.species]
end


text:
using Ballistics

function Δrange(v0::Float64, θ::Float64, target::Types.Target2D, trajectory::Function)

vx = v0*cosd(θ)
vy = v0*sind(θ)
p = scenario.createProjectile(0.03, 0.04,[vx,vy], [0.0,0.0])
trajectory(p,target)
Δrange = target.position[2]-p.position[2]

return Δrange
end

function bissection(target::Types.Target2D; θmin = 0.0, θmax = 45.0, ϵ = 0.001 )

#t2D = scenario.createTarget(0.2,[200.0,50.0])

#global θmin = 0.0
#global θmax = 45.0
v0 = 900.0
#ϵ = 0.001
#global fmin = Δrange(inputs.v0, θmin,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
#global fmax = Δrange(inputs.v0, θmax,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
#global error = min(abs(fmin),abs(fmax))
fmin = Δrange(inputs.v0, θmin,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
fmax = Δrange(inputs.v0, θmax,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
error = min(abs(fmin),abs(fmax))

if sign(fmin)==sign(fmax)
  println("Bisection method fails.")
else

  while error>ϵ

    #global θmin
    #global θmax
    #global fmin
    #global fmax
    #global error
    #global θ = 1/2*(θmin+θmax)
    θ = 1/2*(θmin+θmax)
    #println(error)

    f = Δrange(inputs.v0, θ,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
    error = abs(f)
    if f*fmin < 0
            θmax = θ
        else f*fmax < 0
            θmin = θ
      end
    end

  end
  return θ

end


text:
module NMF

include("NMF_Exact.jl")
include("NMF_EM.jl")
include("NMF_VB.jl")
include("NMF_MC.jl")
include("BackwardKernel.jl")

using .NMF_Exact, .NMF_EM, .NMF_VB, .NMF_MC, .BackwardKernel

export generate, log_marginal
export standard_EM, dual_EM
export standard_VB, online_VB
export Particle, smc_weight, particle_mcmc
export EventQueue, sum, length, eltype, iterate

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
import NLopt
import TransformVariables
const TV = TransformVariables
import ReverseDiff

export NLoptim

Base.@kwdef struct NLoptim{N,G,T} <: ComradeSoss.AbstractOptimizer
    nlopt::N
    grado::G
    transform::T
    dimension::Int
end

function NLoptimRev(lj::Soss.ConditionalModel; alg=:LD_LBFGS, xtol=1e-8, maxeval::Int=50_000)
    #Constuct our objective function and gradients
    tc = xform(lj)
    dim = tc.dimension

    function f(x)
        p, logjac = TV.transform_and_logjac(tc, x)
        return logdensity(lj, p) + logjac
    end

    grado = ReverseDiff.GradientTape(f, rand(dim))
    ctape = ReverseDiff.compile(grado)
    res = DiffResults.GradientResult(zeros(dim))
    function gradf!(x, grad)
        ReverseDiff.gradient!(res, ctape, x)
        grad .= DiffResults.gradient(res)
        return DiffResults.value(res)
    end
    g = zeros(dim)
    gradf!(zeros(dim), g)
    println("A derivative took ")
    @time gradf!(zeros(dim), g)
    println("Adjust your expectations accordingly")

    srange = [(-7.0,7.0) for _ in 1:dim]
    nlopt = NLopt.Opt(alg, length(srange))
    NLopt.lower_bounds!(nlopt, first.(srange))
    NLopt.upper_bounds!(nlopt, last.(srange))
    NLopt.max_objective!(nlopt, gradf!)
    NLopt.xtol_rel!(nlopt, xtol)
    NLopt.maxeval!(nlopt, maxeval)
    return NLoptim(nlopt, grado, tc, dim)
end

function NLoptim(lj::Soss.ConditionalModel; alg=:LD_LBFGS, xtol=1e-8, maxeval::Int=50_000, chunksize::Int=5)
    #Constuct our objective function and gradients
    tc = xform(lj)
    dim = tc.dimension

    function f(x)
        p, logjac = TV.transform_and_logjac(tc, x)
        return logdensity(lj, p) + logjac
    end

    grado = GradObj(f, rand(dim); chunksize=chunksize)
    res = DiffResults.GradientResult(zeros(dim))
    function gradf!(x, grad)
        ForwardDiff.gradient!(res, f, x, grado.cfg)
        grad .= DiffResults.gradient(res)
        return DiffResults.value(res)
    end
    g = zeros(dim)
    gradf!(zeros(dim), g)
    println("A derivative took ")
    @time gradf!(zeros(dim), g)
    println("Adjust your expectations accordingly")

    srange = [(-7.0,7.0) for _ in 1:dim]
    nlopt = NLopt.Opt(alg, length(srange))
    NLopt.lower_bounds!(nlopt, first.(srange))
    NLopt.upper_bounds!(nlopt, last.(srange))
    NLopt.max_objective!(nlopt, gradf!)
    NLopt.xtol_rel!(nlopt, xtol)
    NLopt.maxeval!(nlopt, maxeval)
    return NLoptim(nlopt, grado, tc, dim)
end

function optimize(opt::NLoptim, start=nothing)
    dim = opt.dimension
    nlopt = opt.nlopt
    x0 = zeros(dim)
    lower = nlopt.lower_bounds
    upper = nlopt.upper_bounds
    if isnothing(start)
        x0 = lower .+ rand(dim).*(upper .- lower)
    else
        x0 .= TV.inverse(opt.transform, start)
    end

    (minf,minx,ret) = NLopt.optimize(nlopt, x0)
    println("NLOpt stopped because: $ret")
    return TV.transform(opt.transform, minx), minf
end


text:
const CorrelatedOUDouble = CorrelatedOU+CorrelatedOU
n = 10
Δ = 0.5
@testset "memoryallocation" begin
    @testset "Acv2EIStorage" begin
        x = randn(ComplexF64,3,n); y = copy(x)
        acvstore = Acv2EIStorage(y,reinterpret(SHermitianCompact{2,eltype(y),3}, vec(y)),ones(length(x)),FFTW.ESTIMATE)
        @test fft(x,2) == acvstore.planned_fft*acvstore.allocatedarray
    end
    @testset "Sdf2AcvStorage" begin
        x = randn(ComplexF64,3,n); y = copy(x)
        sdfstore = Sdf2AcvStorage(y,FFTW.ESTIMATE)
        @test ifft(x,2) == sdfstore.planned_ifft*sdfstore.allocatedarray
    end
    @testset "Acv2EIStorageUni" begin
        x = randn(ComplexF64,n); y = copy(x)
        acvstore = Acv2EIStorageUni(y,ones(length(x)),FFTW.ESTIMATE)
        @test fft(x,1) == acvstore.planned_fft*acvstore.allocatedarray
    end
    @testset "Sdf2AcvStorageUni" begin
        x = randn(ComplexF64,n); y = copy(x)
        sdfstore = Sdf2AcvStorageUni(y,FFTW.ESTIMATE)
        @test ifft(x,1) == sdfstore.planned_ifft*sdfstore.allocatedarray
    end
    @testset "EIstorage_function" begin
        @test EIstorage_function(CorrelatedOU,n,ones(2n)) isa Acv2EIStorage
        @test EIstorage_function(OU,n,ones(2n))           isa Acv2EIStorageUni
        @test EIstorage_function(Matern2D,n,ones(2n))     isa Sdf2EIStorage
        @test EIstorage_function(OUUnknown{1},n,ones(2n))  isa Sdf2EIStorageUni
    end
    @testset "EIstorage_gradient" begin
        @test EIstorage_gradient(CorrelatedOU,n,ones(2n)) isa Acv2EIStorage
        @test EIstorage_gradient(OU,n,ones(2n))           isa Acv2EIStorageUni
        @test EIstorage_gradient(Matern2D,n,ones(2n))     isa Sdf2EIStorage
        @test EIstorage_gradient(OUUnknown{1},n,ones(2n)) isa Sdf2EIStorageUni
    end
    @testset "EIstorage_hessian" begin
        @test EIstorage_hessian(CorrelatedOU,n,ones(2n)) isa Acv2EIStorage
        @test EIstorage_hessian(OU,n,ones(2n))           isa Acv2EIStorageUni
        @test EIstorage_hessian(Matern2D,n,ones(2n))     isa Sdf2EIStorage
        @test EIstorage_hessian(OUUnknown{1},n,ones(2n)) isa Sdf2EIStorageUni
    end
    @testset "FreqAcvEst" begin
        freq = FreqAcvEst(3,Δ)
        @test freq.Ωₘ ≈ [0.0,4.1887902047863905,-4.1887902047863905]
        @test freq.n == 3
        @test freq.Δ == Δ
        @test length(FreqAcvEst(Matern2D,n,Δ).Ωₘ) > n
    end
    @testset "LagsEI" begin
        @test LagsEI(3,Δ).lags ≈ [0.0,Δ,2Δ,-3Δ,-2Δ,-Δ]
    end
    @testset "encodetimescale" begin
        @test encodetimescale(OU,n,Δ) isa LagsEI
        @test encodetimescale(Matern2D,n,Δ) isa FreqAcvEst
    end
    @testset "extract_S" begin
        @test WhittleLikelihoodInference.extract_S(allocate_memory_EI_F(CorrelatedOUDouble,n,Δ,taper=nothing)) isa Acv2EIStorage
        @test WhittleLikelihoodInference.extract_S(allocate_memory_EI_F(CorrelatedOUUnknown{1},n,Δ,taper=nothing)) isa Sdf2EIStorage
    end
    @testset "allocate_memory_EI_F" begin
        @test allocate_memory_EI_F(CorrelatedOU,n,Δ,taper=nothing) isa PreallocatedEI
        add_alloc_func = allocate_memory_EI_F(CorrelatedOUDouble,n,Δ,taper=nothing)
        @test add_alloc_func isa AdditiveStorage
        @test add_alloc_func.store1 isa PreallocatedEI
        @test add_alloc_func.store2 isa PreallocatedEI
        @test add_alloc_func.npar1 == 3
    end
    @testset "allocate_memory_EI_FG" begin
        @test allocate_memory_EI_FG(CorrelatedOU,n,Δ,taper=nothing) isa PreallocatedEIGradient
        add_alloc_grad = allocate_memory_EI_FG(CorrelatedOUDouble,n,Δ,taper=nothing)
        @test add_alloc_grad isa AdditiveStorage
        @test add_alloc_grad.store1 isa PreallocatedEIGradient
        @test add_alloc_grad.store2 isa PreallocatedEIGradient
        @test add_alloc_grad.npar1 == 3
    end
    @testset "allocate_memory_EI_FGH" begin
        @test allocate_memory_EI_FGH(CorrelatedOU,n,Δ,taper=nothing) isa PreallocatedEIHessian
        add_alloc_hess = allocate_memory_EI_FGH(CorrelatedOUDouble,n,Δ,taper=nothing)
        @test add_alloc_hess isa AdditiveStorage
        @test add_alloc_hess.store1 isa PreallocatedEIHessian
        @test add_alloc_hess.store2 isa PreallocatedEIHessian
        @test add_alloc_hess.npar1 == 3
    end
    @testset "sdfstorage_function" begin
        @test sdfstorage_function(OU,n) isa SdfStorageUni
        @test sdfstorage_function(CorrelatedOU,n) isa SdfStorage
    end
    @testset "sdfstorage_gradient" begin
        @test sdfstorage_gradient(OU,n) isa SdfStorageUni
        @test sdfstorage_gradient(CorrelatedOU,n) isa SdfStorage
    end
    @testset "sdfstorage_hessian" begin
        @test sdfstorage_hessian(OU,n) isa SdfStorageUni
        @test sdfstorage_hessian(CorrelatedOU,n) isa SdfStorage
    end
    @testset "allocate_memory_sdf_F" begin
        @test allocate_memory_sdf_F(CorrelatedOU,n,3Δ) isa PreallocatedSdf
        add_alloc_func = allocate_memory_sdf_F(CorrelatedOUDouble,n,3Δ)
        @test add_alloc_func isa AdditiveStorage
        @test add_alloc_func.store1 isa PreallocatedSdf
        @test add_alloc_func.store2 isa PreallocatedSdf
        @test add_alloc_func.npar1 == 3
    end
    @testset "allocate_memory_sdf_FG" begin
        @test allocate_memory_sdf_FG(CorrelatedOU,n,3Δ) isa PreallocatedSdfGradient
        add_alloc_grad = allocate_memory_sdf_FG(CorrelatedOUDouble,n,3Δ)
        @test add_alloc_grad isa AdditiveStorage
        @test add_alloc_grad.store1 isa PreallocatedSdfGradient
        @test add_alloc_grad.store2 isa PreallocatedSdfGradient
        @test add_alloc_grad.npar1 == 3
    end
    @testset "allocate_memory_sdf_FGH" begin
        @test allocate_memory_sdf_FGH(CorrelatedOU,n,Δ) isa PreallocatedSdfHessian
        add_alloc_hess = allocate_memory_sdf_FGH(CorrelatedOUDouble,n,Δ)
        @test add_alloc_hess isa AdditiveStorage
        @test add_alloc_hess.store1 isa PreallocatedSdfHessian
        @test add_alloc_hess.store2 isa PreallocatedSdfHessian
        @test add_alloc_hess.npar1 == 3
    end
    @testset "unpack" begin
        x = randn(ComplexF64,3,n)
        y1,y2,y3,y4,y5,y6 = fill(copy(x),6)
        h1 = reinterpret(SHermitianCompact{2,eltype(y1),3}, vec(y1))
        h2 = reinterpret(SHermitianCompact{2,eltype(y2),3}, vec(y2))
        h3 = reinterpret(SHermitianCompact{2,eltype(y3),3}, vec(y3))
        acvstore    = Acv2EIStorage(y1,h1,ones(length(x)),FFTW.ESTIMATE)
        sdfstore    = Sdf2EIStorage(y2,y2,h2,ones(length(x)),FFTW.ESTIMATE)
        sdfonly     = SdfStorage(y3,h3)
        acvstoreuni = Acv2EIStorageUni(y4,ones(length(x)),FFTW.ESTIMATE)
        sdfstoreuni = Sdf2EIStorageUni(y5,y5,ones(length(x)),FFTW.ESTIMATE)
        sdfonlyuni  = SdfStorageUni(y6)
        @test WhittleLikelihoodInference.unpack(acvstore)    === h1
        @test WhittleLikelihoodInference.unpack(sdfstore)    === h2
        @test WhittleLikelihoodInference.unpack(sdfonly )    === h3
        @test WhittleLikelihoodInference.unpack(acvstoreuni) === y4
        @test WhittleLikelihoodInference.unpack(sdfstoreuni) === y5
        @test WhittleLikelihoodInference.unpack(sdfonlyuni)  === y6
    end
    @testset "taper2timedomainkernel" begin
        @test WhittleLikelihoodInference.taper2timedomainkernel(ones(4)./2,4) == WhittleLikelihoodInference.taper2timedomainkernel(nothing,4)
    end
end


text:
export get_coin_list, get_coin_latest, get_historical_daily, get_mining_calculator_data

function get_coin_list(; extraParams::String="", sign::Bool=false, api_key::String="")
    url = build_blockchain_data_url("available_coin_list", extraParams=extraParams, sign=sign, api_key=api_key)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)

        datadict = jobj["Data"]
        df = DataFrame(id=Int[], data_available_from_timestamp= Int[], data_available_from=DateTime[], key=String[], symbol=String[], partner_symbol=String[], )

        for k in keys(datadict)
            value = datadict[k]
            push!(df.id, value["id"])
            push!(df.data_available_from_timestamp, value["data_available_from"])
            push!(df.data_available_from, unix2datetime(value["data_available_from"]))
            push!(df.key, k)
            push!(df.symbol, value["symbol"])
            push!(df.partner_symbol, value["partner_symbol"])
        end

        return df

    catch e
        return e
    end

end

function get_coin_latest(fsym::String; extraParams::String="", sign::Bool=false, api_key::String="")
    url = build_blockchain_data_url("latest", fsym=fsym, extraParams=extraParams, sign=sign, api_key=api_key)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]
        return datadict

    catch e
        return e
    end
end

"""
    get_historical_daily(fsym::String; limit::Int=30, toTs::String="", extraParams::String="", sign::String="", api_key::String=api_key)

1. fsym [Required]: The cryptocurrency symbol of interest [ Min length - 1] [ Max length - 10]
2. limit: The number of data points to return. [ Min - 1] [ Max - 2000] [ Default - 30]
3. toTs: Returns historical data before that timestamp. If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}
4. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
5. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_historical_daily(fsym::String; limit::Int=30, toTs::String="", extraParams::String="", sign::String="", api_key::String="")
    url = build_blockchain_data_url("historical_daily", fsym=fsym, limit=limit, toTs=toTs, extraParams=extraParams, sign=sign, api_key=api_key)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]["Data"]

        df = DataFrame()
        nrows = size(datadict)[1]
        for i in 1:nrows
            push!(df, datadict[i], cols=:union)
        end
        df.datetime = unix2datetime.(df.time)
        return df

    catch e
        return e
    end
end

"""
    get_mining_calculator_data(fsyms::Array{String,1}, tsyms::Array{String,1}, extraParams::String, sign::Bool, api_key::String)

    Get the price, the total supply, the hash rate, the block reward and the block time. Based on these values you can calculate the profitability of a mining rig.
1. fsyms [Required]: Comma separated cryptocurrency symbols list [ Min length - 1] [ Max length - 300]
2. tsyms [Required]: Comma separated cryptocurrency symbols list to convert into [ Min length - 1] [ Max length - 100]
3. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
4. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_mining_calculator_data(fsyms::Array{String,1}, tsyms::Array{String,1}; extraParams::String="", sign::Bool=false, api_key::String)
    url = build_blockchain_data_url("mining_calculator_data", fsyms=fsyms, tsyms=tsyms, extraParams=extraParams, sign=sign, api_key=api_key)
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]
        return datadict

    catch e
        return e
    end
end


text:
# ---
# title: 1529. Bulb Switcher IV
# id: problem1529
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/bulb-switcher-iv/description/>
# hidden: true
# ---
#
# There is a room with `n` bulbs, numbered from `0` to `n-1`, arranged in a row
# from left to right. Initially all the bulbs are **turned off**.
#
# Your task is to obtain the configuration represented by `target` where
# `target[i]` is '1' if the i-th bulb is turned on and is '0' if it is turned
# off.
#
# You have a switch to flip the state of the bulb, a flip operation is defined
# as follows:
#
#   * Choose **any** bulb (index `i`) of your current configuration.
#   * Flip each bulb from index `i` to `n-1`.
#
# When any bulb is flipped it means that if it is 0 it changes to 1 and if it is
# 1 it changes to 0.
#
# Return the **minimum** number of flips required to form `target`.
#
#
#
# **Example 1:**
#
#
#
#     Input: target = "10111"
#     Output: 3
#     Explanation: Initial configuration "00000".
#     flip from the third bulb:  "00000" -> "00111"
#     flip from the first bulb:  "00111" -> "11000"
#     flip from the second bulb:  "11000" -> "10111"
#     We need at least 3 flip operations to form target.
#
# **Example 2:**
#
#
#
#     Input: target = "101"
#     Output: 3
#     Explanation: "000" -> "111" -> "100" -> "101".
#
#
# **Example 3:**
#
#
#
#     Input: target = "00000"
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: target = "001011101"
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * `1 <= target.length <= 10^5`
#   * `target[i] == '0'` or `target[i] == '1'`
#
#
## @lc code=start
using LeetCode

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


text:
module FosterParameters

    using GlobalDecls
    using CareData: CPIINDEX
    using Parameters

    export get_default_params, Params, Minimum_Allowances, AllowancesDict, makeAllowances, getAllowances
    export DEFAULT_UPRATE_FACTOR

    export ContributionType, actual_contribution, no_contribution, benefits_only, flat_rate, all_people, hb_only
    export PaymentType, actual_payment, min_payment, age_16_17
    export UpratingType, no_uprating, by_benefits, by_cpi

    @with_kw  mutable struct Minimum_Allowances
        babies :: Real
        pre_primary :: Real
        primary :: Real
        age_11_15 :: Real
        age_16_17 :: Real
    end

    function Base.:*( a::Minimum_Allowances, x::Real )
        a.babies *= x
        a.pre_primary *= x
        a.primary *= x
        a.age_11_15 *= x
        a.age_16_17 *= x
        a
    end

    AllowancesDict = Dict{Region,Minimum_Allowances}

    function makeAllowances( year :: Integer ) :: AllowancesDict
        d = AllowancesDict()
        if year <= 2016 # don't really need to calculate further back
            d[london] = Minimum_Allowances( 142.0, 145.0, 163.0, 184.0, 216.0 )
            d[se] = Minimum_Allowances( 136.0, 140.0, 156.0, 177.0, 208.0 )
            d[rest_of_england] = Minimum_Allowances( 123.0, 126.0, 139.0, 159.0, 185.0 )
        elseif year == 2017
                d[london] = Minimum_Allowances( 144.0, 147.0, 165.0, 187.0, 216.0 )
                d[se] = Minimum_Allowances( 138.0, 142.0, 158.0, 179.0, 211.0 )
                d[rest_of_england] = Minimum_Allowances( 125.0, 128.0, 141.0, 161.0, 188.0 )
        elseif year == 2018
            d[london] = Minimum_Allowances( 146.0, 149.0, 168.0, 190.0, 222.0 )
            d[se] = Minimum_Allowances( 140.0, 144.0, 160.0, 182.0, 214.0 )
            d[rest_of_england] = Minimum_Allowances( 127.0, 130.0, 143.0, 164.0, 191.0 )
        elseif year == 2019
            d[london] = Minimum_Allowances( 149.0, 152.0, 171.0, 193.0, 226.0 )
            d[se] = Minimum_Allowances( 143.0, 147.0, 163.0, 185.0, 218.0 )
            d[rest_of_england] = Minimum_Allowances( 129.0, 132.0, 146.0, 164.0, 194.0 )
        elseif year == 2020
            d[london] = Minimum_Allowances( 152.0, 155.0, 174.0, 197.0, 231.0 )
            d[se] = Minimum_Allowances( 146.0, 150.0, 166.0, 189.0, 222.0 )
            d[rest_of_england] = Minimum_Allowances( 132.0, 135.0, 149.0, 170.0, 198.0 )
        elseif year > 2020
            d[london] = Minimum_Allowances( 152.0, 155.0, 174.0, 197.0, 231.0 )* CPIINDEX[year]
            d[se] = Minimum_Allowances( 146.0, 150.0, 166.0, 189.0, 222.0 )* CPIINDEX[year]
            d[rest_of_england] = Minimum_Allowances( 132.0, 135.0, 149.0, 170.0, 198.0 )*CPIINDEX[year]
        end
        d
    end

    DEFAULT_UPRATE_FACTOR =
        makeAllowances( 2020 )[london].age_16_17/
        makeAllowances( 2019 )[london].age_16_17

    @enum ContributionType actual_contribution no_contribution  benefits_only flat_rate all_people hb_only
    @enum PaymentType actual_payment min_payment age_16_17
    @enum UpratingType no_uprating by_benefits by_cpi

    mutable struct Params
        name            :: AbstractString
        yp_contrib_type :: ContributionType
        contrib_hb      :: ContributionType
        yp_contrib      :: Real
        payment         :: PaymentType
        uprating        :: UpratingType
        taper           :: Array{Real}
        fee             :: Array{Real}
        extra_grant     :: Real
        prop_fees_deleted :: Real
    end

    function get_default_params()
        Params(
            "default parameters",
            actual_contribution,
            actual_contribution,
            0.0,
            actual_payment,
            no_uprating,
            [],
            [],
            0.0,
            0.0 )
    end

    function getAllowances( reg :: Region, year :: Integer ) :: Minimum_Allowances
        # println( "getAllowances; reg=$reg year=$year")
        return makeAllowances( year )[reg]
    end

end


text:
using LightGraphs

"""Get a set of all remaining keys in plans"""
remaining_keys(plans) = Set(vcat([[c for c in plan if c ∈ 'a':'z'] for plan in plans]...))

"""Get a dictionary of keychar => doorcoord for doors in plan"""
doorcoords(plan) = Dict(lowercase(plan[coord]) => coord for coord in CartesianIndices(plan) if plan[coord] ∈ 'A':'Z')

"""Build a graph with locked doors not included. Return graph, nodes, nodemap, keynodes, doornodes"""
function build_graph(plan)
    nodes = [coord for coord in CartesianIndices(plan) if plan[coord] != '#' && plan[coord] ∉ 'A':'Z']
    nodemap = Dict(coord => i for (i, coord) in enumerate(nodes))
    graph = SimpleGraph(length(nodes))
    keynodes = []
    doornodes = Dict(lowercase(plan[coord]) => coord for coord in CartesianIndices(plan) if plan[coord] ∈ 'A':'Z')
    for (i, coord) in enumerate(nodes)
        # Skip walls
        if plan[coord] == '#' continue end
        # Everything else is walkable, although we'll need keys to access the door nodes
        if plan[coord] != '@' && plan[coord] != '.'
            push!(keynodes, i)
        end
        # Add east and south edges (other nodes will add north and west for us)
        for offset in [CartesianIndex(0, 1), CartesianIndex(1, 0)]
            if (coord + offset) ∈ keys(nodemap) add_edge!(graph, i, nodemap[coord + offset]) end
        end
    end
    return graph, nodes, nodemap, keynodes, doornodes
end

"""Collect all keys in the plans recursively, returning the shortest distance to collect them"""
function collect_keys(plans, current_coords, cached_solutions=Dict())
    # Do we have an easy cached answer?
    search_state = (current_coords, remaining_keys(plans))
    if search_state in keys(cached_solutions) return cached_solutions[search_state] end
    # Have we actually collected all the keys?
    if isempty(last(search_state)) return 0 end

    path_lengths = []
    for (p, (plan, current_coord)) in enumerate(zip(plans, current_coords))
        # We need to exhaustively search all key pickups that are currently possible
        graph, nodes, nodemap, keynodes, doornodes = build_graph(plan)
        # Have we finished this plan?
        if length(keynodes) == 0 continue end
        possible_paths = []
        collected = []
        ds = dijkstra_shortest_paths(graph, nodemap[current_coord])
        for ki in keynodes
            if ds.pathcounts[ki] > 0
                path = enumerate_paths(ds, ki)
                # Skip any paths with more than one key on - the path to the first key will do just fine
                if length(keynodes ∩ path) > 1 continue end
                push!(possible_paths, path)
                push!(collected, plan[nodes[ki]])
            end
        end

        # For each potential path with a single key, unlock the corresponding door, then recursively search for the next move
        for path in possible_paths
            collected_key_coord = nodes[last(path)]
            # Clone list of all plans, then replace this one with a copy after moving along this path
            new_plans = copy(plans)
            new_plan = copy(plan)
            new_plans[p] = new_plan
            new_plan[collected_key_coord] = '.'
            # Figure out which plan the door was in and remove it
            for (op, otherplan) in enumerate(plans)
                dc = doorcoords(otherplan)
                if plan[collected_key_coord] ∈ keys(dc)
                    if op != p
                        # We have not yet copied this plan - clone it
                        new_plans[op] = copy(otherplan)
                    end
                    new_plans[op][dc[plan[collected_key_coord]]] = '.'
                    break
                end
            end
            new_coords = copy(current_coords)
            new_coords[p] = collected_key_coord
            push!(path_lengths, length(path) - 1 + collect_keys(new_plans, new_coords, cached_solutions))
        end
    end

    # Cache this solution in case we have another path ending at the same place in future
    cached_solutions[search_state] = minimum(path_lengths)
    if length(cached_solutions) % 100 == 0 println("Cached solutions: $(length(cached_solutions))") end
    return minimum(path_lengths)
end

function main()
    # Part 1: build graph of corridors, keys and doors, then solve for shortest path to get all keys
    # Sounds a bit travelling salesperson-y.
    plan = hcat(open(f -> collect.(readlines(f)), "day18_input.txt")...)
    shortest_dist = collect_keys([plan], [indexin('@', plan)[1]])
    println("Part 1: Shortest distance: $shortest_dist")

    # Part 2: this time there are 4 vaults with robots, one moving at a time. Patch and split the map.
    centre = (size(plan) .+ 1) .÷ 2
    replacement = ['@' '#' '@'; '#' '#' '#'; '@' '#' '@']
    plan[(centre[1]-1):(centre[1]+1), (centre[2]-1):(centre[2]+1)] = replacement
    subplans = [
        plan[1:(centre[1]-1), 1:(centre[2]-1)],
        plan[1:(centre[1]-1), (centre[2]+1):end],
        plan[(centre[1]+1):end, 1:(centre[2]-1)],
        plan[(centre[1]+1):end, (centre[2]+1):end]
    ]

    shortest_dist = collect_keys(subplans, [indexin('@', plan)[1] for plan in subplans])
    println("Part 2: Shortest distance: $shortest_dist")
end
main()

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
module MusicTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Music" => Ahorn.EntityPlacement(
        Maple.MusicTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.MusicTrigger)
    return Dict{String, Any}(
        "track" => sort(collect(keys(Maple.Songs.songs)))
    )
end

end


text:
# Fannie Mae

# download("http://rapidsai-data.s3-website.us-east-2.amazonaws.com/notebook-mortgage-data/mortgage_2000.tgz", "c:/data/mortgage_2000.tgz")

# un-tar and uncompress the file
# ;tar zxvg c:/data/mortgage_2000.tgz

# uncomment for debugging
dirpath = "C:/data/Performance_All/"
largest_file = "Performance_2000Q4.txt"
outpath = "c:/data/jdf-bench/Performance_2000Q4.txt"
data_label =  "Fannie Mae Performance 2000Q4"
# delim = ','
# header = true

include("C:/Users/RTX2080/git/JDF/benchmarks/benchmarks.jl")
#@timegen_benchmark("c:/data/AirOnTimeCSV/", "airOT199302.csv", "c:/data/jdf-bench/airOT199302.csv", "Air On Time 199302")
@time gen_benchmark(dirpath, largest_file, outpath, data_label, delim = '|', header = false);

sizedf = sum_file_size(outpath)

using StatsPlots

p = plot(
    sizedf.pkg,
    sizedf.fs/1024^3,
    linetype = :bar,
    ylab = "Size (GB)",
    legend = false,
    title = "On-disk file Size for various formats\n $data_label data")
savefig(p, joinpath(outpath, largest_file*"_filesize.png"))


text:
using Test
using Flux
using LinearAlgebra
using NeuralArithmetic
using NeuralArithmetic: weights, gate
using Random

Random.seed!(0)

include("nau.jl")
include("nmu.jl")
include("npu.jl")

include("nac.jl")
include("nalu.jl")
include("inalu.jl")

using DiffEqFlux
include("fast_nau.jl")
include("fast_npu.jl")

@testset "show functions" begin
    layers = [
        "NPU","RealNPU","NaiveNPU","NAU","NMU","NAC","NALU",
        "FastNPU","FastRealNPU","FastNAU"
    ]
    for l in layers
        m = eval(Symbol(l))(2,2)
        @test repr(m) == "$l(in=2, out=2)"
    end
end


text:
include("help.jl")
function mark_kross_max(r::Robot)
    for side in (Nord, West, Sud, Ost)
        num_steps = putmarkers_1!(r,side)
        movements!(r,inverse(side), num_steps)
    end
    putmarker!(r)
end

function putmarkers_1!(r::Robot,side::HorizonSide)
    num_steps=0
    while move_if_possible!(r, side) == true
        putmarker!(r)
        num_steps += 1
    end
    return num_steps
end
movements!(r::Robot, side::HorizonSide, num_steps::Int) = for _ in 1:num_steps move!(r,side) end

    movements!(r::Robot, side::HorizonSide) = while isborder(r,side)==false move!(r,side) end

movements!(r::Robot, side::HorizonSide, num_steps::Int) =
for _ in 1:num_steps
    move_if_possible!(r,side)
end
go_left(side::HorizonSide) =  HorizonSide(mod(Int(side)+1, 4))

function move_if_possible!(r::Robot, direct_side::HorizonSide)::Bool
    orthogonal_side = go_left(direct_side)
    reverse_side = inverse(orthogonal_side)
    num_steps=0
    while isborder(r,direct_side) == true
        if isborder(r, orthogonal_side) == false
            move!(r, orthogonal_side)
            num_steps += 1
        else
            break
        end
    end

    if isborder(r,direct_side) == false
        move!(r,direct_side)
        while isborder(r,reverse_side) == true
            move!(r,direct_side)
        end
        result = true
    else
        result = false
    end
    for _ in 1:num_steps
        move!(r,reverse_side)
    end
    return result
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
"""
    isweekend(dt) :: Bool

Returns `true` for Saturdays or Sundays.
Returns `false` otherwise.
"""
@inline isweekend(dt::Dates.Date) :: Bool = signbit(5 - Dates.dayofweek(dt))

"""
    isweekday(dt) :: Bool

Returns `true` for Monday to Friday.
Returns `false` otherwise.
"""
@inline isweekday(dt::Dates.Date) :: Bool = signbit(Dates.dayofweek(dt)  - 6)

"""
    isbday(calendar, dt) :: Bool

Returns `false` for weekends or holidays.
Returns `true` otherwise.
"""
function isbday(hc::HolidayCalendar, dt::Dates.Date) :: Bool
    if _getcachestate(hc)
        return isbday(_getholidaycalendarcache(hc), dt)
    else
        return !(isweekend(dt) || isholiday(hc, dt))
    end
end

@inline isbday(calendar, dt) :: Bool = isbday(convert(HolidayCalendar, calendar), dt)

"""
    tobday(calendar, dt; [forward=true]) :: Dates.Date

Adjusts `dt` to next Business Day if it's not a Business Day.
If `isbday(dt)`, returns `dt`.
"""
function tobday(hc::HolidayCalendar, dt::Dates.Date; forward::Bool = true) :: Dates.Date
    if isbday(hc, dt)
        return dt
    else
        increment = forward ? 1 : -1
        next_date = dt + Dates.Day(increment)

        while !isbday(hc, next_date)
            next_date += Dates.Day(increment)
        end
    end

    return next_date
end

tobday(calendar, dt; forward::Bool = true) = tobday(convert(HolidayCalendar, calendar), dt; forward=forward)

"""
    advancebdays(calendar, dt, bdays_count) :: Dates.Date

Increments given date `dt` by `bdays_count`.
Decrements it if `bdays_count` is negative.
`bdays_count` can be a `Int`, `Dates.Day`, `Vector{Int}`, `Vector{Dates.Day}` or a `UnitRange`.

Computation starts by next Business Day if `dt` is not a Business Day.
"""
function advancebdays(hc::HolidayCalendar, dt::Dates.Date, bdays_count::Int) :: Dates.Date
    result = tobday(hc, dt)

    # does nothing
    if bdays_count == 0
        return result
    end

    # if bdays_count is positive, goes forward. Otherwise, goes backwards.
    increment = bdays_count > 0 ? +1 : -1

    num_iterations = abs(bdays_count)

    while num_iterations > 0
        result += Dates.Day(increment)

        # Looks for previous / next Business Day
        while !isbday(hc, result)
            result += Dates.Day(increment)
        end

        num_iterations += -1
    end

    return result
end

advancebdays(hc::HolidayCalendar, dt::Dates.Date, bdays_count::Dates.Day) = advancebdays(hc, dt, Dates.value(bdays_count))

const BDaysCountType = Union{Int, Dates.Day}

function advancebdays(calendar, dt, bdays_count::Union{T, Vector{T}, AbstractRange}) where {T<:BDaysCountType}
    advancebdays(convert(HolidayCalendar, calendar), convert(Dates.Date, dt), bdays_count)
end

"""
    bdayscount(calendar, dt0, dt1) :: Int

Counts the number of Business Days between `dt0` and `dt1`.
Returns `Int`.

Computation is always based on next Business Day if given dates are not Business Days.
"""
function bdayscount(hc::HolidayCalendar, dt0::Dates.Date, dt1::Dates.Date) :: Int
    if _getcachestate(hc)
        return bdayscount(_getholidaycalendarcache(hc), dt0, dt1)
    else
        dt0 = tobday(hc, dt0)
        dt1 = tobday(hc, dt1)
        inc = dt0 <= dt1 ? +1 : -1

        result = 0
        while dt0 != dt1
            dt0 = advancebdays(hc, dt0, inc)
            result += inc
        end

        return result
    end
end

bdayscount(calendar, dt0::Dates.Date, dt1::T) where {T<:Union{Dates.Date, Vector{Dates.Date}}} = bdayscount(convert(HolidayCalendar, calendar), dt0, dt1)
bdayscount(calendar, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date}) = bdayscount(convert(HolidayCalendar, calendar), dt0, dt1)

"""
    bdays(calendar, dt0, dt1) :: Dates.Day

Counts the number of Business Days between `dt0` and `dt1`.
Returns instances of `Dates.Day`.

Computation is always based on next Business Day if given dates are not Business Days.
"""
bdays(hc::HolidayCalendar, dt0::Dates.Date, dt1::Dates.Date) = Dates.Day(bdayscount(hc, dt0, dt1))
bdays(calendar, dt0::Dates.Date, dt1::T) where {T<:Union{Dates.Date, Vector{Dates.Date}}} = bdays(convert(HolidayCalendar, calendar), dt0, dt1)
bdays(calendar, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date}) = bdays(convert(HolidayCalendar, calendar), dt0, dt1)

"""
    firstbdayofmonth(calendar, dt) :: Dates.Date
    firstbdayofmonth(calendar, yy, mm) :: Dates.Date

Returns the first business day of month.
"""
firstbdayofmonth(calendar, dt::Dates.Date) = tobday(calendar, Dates.firstdayofmonth(dt))

"""
    lastbdayofmonth(calendar, dt) :: Dates.Date
    lastbdayofmonth(calendar, yy, mm) :: Dates.Date

Returns the last business day of month.
"""
lastbdayofmonth(calendar, dt::Dates.Date) = tobday(calendar, Dates.lastdayofmonth(dt), forward=false)
firstbdayofmonth(calendar, yy::T, mm::T) where {T<:Integer} = firstbdayofmonth(calendar, Dates.Date(yy, mm, 1))
firstbdayofmonth(calendar, yy::Dates.Year, mm::Dates.Month) = firstbdayofmonth(calendar, Dates.Date(yy, mm, Dates.Day(1)))
lastbdayofmonth(calendar, yy::T, mm::T) where {T<:Integer} = lastbdayofmonth(calendar, Dates.Date(yy, mm, 1))
lastbdayofmonth(calendar, yy::Dates.Year, mm::Dates.Month) = lastbdayofmonth(calendar, Dates.Date(yy, mm, Dates.Day(1)))


text:
#!/usr/bin/env julia

using SynPUF, Test

function tests()
  @testset "Download of SynPUF" begin
    using DataDeps
    using SynPUF
    SynPUF.register()
    path = datadep"synpuf_01"
    @test path == endswith(path, "datadeps/synpuf_01")
  end
end

tests()


text:
module FolsomData

    using ProjAssistant
    @gen_sub_proj

    include("data.jl")

    function __init__()
        _populate_bundle()
        @create_proj_dirs
    end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    AbstractMaxima

A collection of maxima (e.g. block maxima)
"""
abstract type AbstractMaxima end

"""
    collect(maxima)

Collect `maxima` values using appropriate model.
"""
Base.collect(m::AbstractMaxima) = error("not implemented")

#------------------
# IMPLEMENTATIONS
#------------------
"""
    BlockMaxima(x, n)

Maxima obtained by splitting the data `x` into blocks of size `n`.
"""
struct BlockMaxima{V<:AbstractVector} <: AbstractMaxima
  x::V
  n::Int
end

function Base.collect(m::BlockMaxima)
  result = Vector{Float64}()
  for i in 1:m.n:length(m.x)-m.n+1
    xs = skipmissing(view(m.x, i:i+m.n-1))

    xmax = -Inf
    for x in xs
      x > xmax && (xmax = x)
    end

    isinf(xmax) || push!(result, xmax)
  end

  result
end

"""
    PeakOverThreshold(x, u)

Maxima obtained by thresholding the data `x` with threshold `u`.
"""
struct PeakOverThreshold{V<:AbstractVector} <: AbstractMaxima
  x::V
  u::Float64
end

Base.collect(m::PeakOverThreshold) = filter(x -> !ismissing(x) && x > m.u, m.x)


text:
# bootstrapping traits
const avg_rate = 0.05
const max_rate = 1.0

const AVG_HAZARD_RATE = 0.01
const MAX_HAZARD_RATE = 1.0

# Discount bootstrap trait
type Discount <: BootstrapTrait end
type HazardRate <: BootstrapTrait end

# Type aliases
# typealias InterpolatedHazardRateCurve{P, T} InterpolatedDefaultProbabilityCurve{P, HazardRate}

initial_value(::Discount) = 1.0
max_iterations(::Discount) = 100

initial_value(::HazardRate) = AVG_HAZARD_RATE
max_iterations(::HazardRate) = 30

function guess{T <: TermStructure}(::Discount, i::Int, ts::T, valid::Bool)
  if valid
    # return previous iteration value
    return ts.data[i]
  end

  if i == 2
    # first pillar
    return 1.0 / (1.0 + avg_rate * ts.times[2])
  end

  r = -log(ts.data[i - 1]) / ts.times[i - 1]
  return exp(-r * ts.times[i])
end

function min_value_after{T <: TermStructure}(::Discount, i::Int, ts::T, valid::Bool)
  if valid
    return ts.data[end] / 2.0
  end

  dt = ts.times[i] - ts.times[i - 1]
  return ts.data[i - 1] * exp( -max_rate * dt)
end

function max_value_after{T <: TermStructure}(::Discount, i::Int, ts::T, ::Bool)
  dt = ts.times[i] - ts.times[i - 1]
  return ts.data[i - 1] * exp(max_rate * dt)
  #return ts.data[i - 1]
end

function update_guess!{T <: TermStructure}(::Discount, i::Int, ts::T, discount::Float64)
  ts.data[i] = discount
  return ts
end

function guess(::HazardRate, i::Int, ts::TermStructure, valid::Bool)
  if valid
    return ts.data[i]
  end

  if i == 2 # first pillar
    return AVG_HAZARD_RATE
  end

  # extrapolate
  d = ts.dates[i]
  return hazard_rate(ts, d)
end

function min_value_after(::HazardRate, i::Int, ts::TermStructure, valid::Bool)
  if valid
    r = min(ts.data[1], ts.data[end])
    return r / 2.0
  end

  return eps()
end

function max_value_after(::HazardRate, i::Int, ts::TermStructure, valid::Bool)
  if valid
    r = max(ts.data[1], ts.data[end])
    return r * 2.0
  end

  return MAX_HAZARD_RATE
end

function update_guess!(::HazardRate, i::Int, ts::TermStructure, rate::Float64)
  ts.data[i] = rate
  if i == 2
    ts.data[1] = rate # first point handled as well
  end

  return ts
end


text:
using StatisticalMachineLearning
using Documenter

DocMeta.setdocmeta!(StatisticalMachineLearning, :DocTestSetup, :(using StatisticalMachineLearning); recursive=true)

makedocs(;
    modules=[StatisticalMachineLearning],
    authors="KronosTheLate",
    repo="https://github.com/KronosTheLate/StatisticalMachineLearning.jl/blob/{commit}{path}#{line}",
    sitename="StatisticalMachineLearning.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://KronosTheLate.github.io/StatisticalMachineLearning.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/KronosTheLate/StatisticalMachineLearning.jl",
    devbranch="master",
)


text:
"""
    Functions in this file read the subfind output.

"""

using ProgressMeter
using Base.Threads

"""
    struct SubfindHeader

Contains the data of the `HEAD` block in the subfind output

# Fields
| Name                                 | Meaning                                                                                |
| :----------------------------------  | :------------------------------------------------------------------------------------- |
| `nhalos::Int32`                      | number of halos in the output file                                                     |
| `nsubhalos::Int32`                   | number of subhalos in the output file                                                  |
| `nfof::Int32`                        | number of particles in the FoF                                                         |
| `ngroups::Int32`                     | number of large groups in the output file                                              |
| `time::Float64`                      | time / scale factor of the simulation                                                  |
| `z::Float64`                         | redshift of the simulation                                                             |
| `tothalos::UInt32`                   | total number of halos over all output files                                            |
| `totsubhalos::UInt32`                | total number of subhalos over all output files                                         |
| `totfof::UInt32`                     | total number of particles in the FoF                                                   |
| `totgroups::UInt32`                  | 1 if simulation was run with cooling, else 0                                           |
| `num_files::Int32`                   | number of files over which subfind data is distributed                                 |
| `boxsize::Float64`                   | total size of the simulation box                                                       |
| `omega_0::Float64`                   | Omega matter                                                                           |
| `omega_l::Float64`                   | Omega dark enery                                                                       |
| `h0::Float64`                        | little h                                                                               |
| `flag_doubleprecision::Int32`        | 1 if snapshot is in double precision, else 0                                           |
| `flag_ic_info::Int32`                | 1 if initial snapshot file contains an info block, else 0                              |

"""
struct SubfindHeader <: AbstractGadgetHeader
    nhalos::Int32                       # number of halos in the output file
    nsubhalos::Int32                    # number of subhalos in the output file
    nfof::Int32                         # number of particles in the FoF
    ngroups::Int32                      # number of large groups in the output file
    time::Float64                       # time / scale factor of the simulation
    z::Float64                          # redshift of the simulation
    tothalos::UInt32                    # total number of halos over all output files
    totsubhalos::UInt32                 # total number of subhalos over all output files
    totfof::UInt32                      # total number of particles in the FoF
    totgroups::UInt32                   # total number of large groups over all output files
    num_files::Int32                    # number of files
    boxsize::Float64                    # total size of the simulation box
    omega_0::Float64                    # Omega matter
    omega_l::Float64                    # Omega dark enery
    h0::Float64                         # little h
    flag_doubleprecision::Int32         # 1 if snapshot is in double precision, else 0
    flag_ic_info::Int32
end

"""
    struct HaloID
        file::Int64
        id::Int64
    end

Stores the subfile that contains the halo in `file` and the position in the block in `id`.
"""
struct HaloID
    file::Int64
    id::Int64
end

"""
    read_subfind_header(filename::String)

Reads the header of a subfind file or file base (without .0, .1, etc.) into a SubfindHeader struct.
"""
function read_subfind_header(filename::String)

    filename = select_file(filename, 0)
    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize == 8
        swap = 0
        snap_format = 2
    elseif blocksize == 256
        swap = 0
        snap_format = 1
    else
        blocksize = bswap(blocksize)
        if blocksize == 8
            swap = 1
            snap_format = 2
        elseif blocksize == 256
            swap = 1
            snap_format = 1
        else
            println("incorrect file format encountered when reading header of", filename)
        end
    end

    #println("Reading snapshot format: ", snap_format)

    if snap_format == 2
        seek(f, 16)
        skip_line = read(f, Int32)
    end

    nhalos = read(f, Int32)
    nsubhalos = read(f, Int32)
    nfof = read(f, Int32)
    ngroups = read(f, Int32)

    # rest of numbers and massarr are obsolete
    DUMMY = read!(f, Array{Int32,1}(undef,2))
    DUMMY = read!(f, Array{Float64,1}(undef,6))

    time = read(f, Float64)
    z = read(f, Float64)

    # flags not needed
    DUMMY = read(f, Int32)
    DUMMY = read(f, Int32)

    tothalos = read(f, UInt32)
    totsubhalos = read(f, UInt32)
    totfof = read(f, UInt32)
    totgroups = read(f, UInt32)

    DUMMY = read!(f, Array{Int32,1}(undef,2))
    DUMMY = read(f, Int32)

    num_files = read(f, Int32)
    boxsize = read(f, Float64)
    omega_0 = read(f, Float64)
    omega_l = read(f, Float64)
    h0 = read(f, Float64)

    DUMMY = read(f, Int32)
    DUMMY = read(f, Int32)
    DUMMY = read!(f, Array{UInt32,1}(undef,6))
    DUMMY = read(f, Int32)

    flag_doubleprecision = read(f, Int32)
    flag_ic_info = read(f, Int32)

    close(f)

    return SubfindHeader(nhalos, nsubhalos, nfof, ngroups,
                         time, z,
                         tothalos, totsubhalos, totfof, totgroups,
                         num_files,
                         boxsize, omega_0, omega_l, h0,
                         flag_doubleprecision,
                         flag_ic_info)
end

"""
    read_subfind(filename::String, blockname::String)

Reads a block of a subfind file.
"""
function read_subfind(filename::String, blockname::String)

    # read the info block
    info = read_info(filename)

    blocknames = getfield.(info, :block_name)
    ind = findfirst(==(blockname), blocknames)

    # the the block is not contained in the file throw and error
    if isnothing(ind)
        error("Block $blockname not present!")
    end

    # get the relevant entry
    info_selected = info[ind]

    # blocks are type specific so we can use this to make our life easier
    parttype = findfirst(==(1), info_selected.is_present) - 1

    return read_block(filename, blockname, info = info_selected, parttype = parttype)
end


text:
#!/usr/bin/env julia

const CPU_FLOPS = peakflops()
#Start Test Script
using FiniteElementDiffEq, DiffEqDevTools, DiffEqProblemLibrary, DiffEqPDEBase
using Base.Test

tic()

# Internals
println("Assembly Tests")
@time @test include("internals/assembly_tests.jl")

#Heat
println("Finite Element Heat Dt Tests")
@time @test include("heat/femheat_dtconvergence_tests.jl")
println("Finite Element Heat Dx Tests")
@time @test include("heat/femheat_dxconvergence_tests.jl")
println("Finite Element Heat Method Tests")
@time @test include("heat/femheat_methods_tests.jl")
println("Finite Element Nonlinear Heat Methods Tests")
@time @test include("heat/femheat_nonlinearmethods_tests.jl")
println("Finite Element Nonlinear System Heat Tests")
@time @test include("heat/femheat_system_tests.jl")
println("Heat Animation Test")
@time @test include("heat/femheat_animation_tests.jl")
println("Stochastic Heat Animation Test")
@time @test include("heat/femheat_stochasticanimation_tests.jl")

#Poisson
println("Finite Element Poisson Convergence Test")
@time @test include("poisson/fempoisson_convergence_tests.jl")
println("Finite Element Nonlinear Poisson Tests")
@time @test include("poisson/fempoisson_nonlinear_tests.jl")
println("Finite Element Nonlinear System Poisson Tests")
@time @test include("poisson/fempoisson_nonlinearsystem_tests.jl")
println("Finite Element Stochastic Poisson")
@time @test include("poisson/fempoisson_stochastic_tests.jl")
println("Finite Element Poisson")
@time @test include("poisson/fempoisson_linear_tests.jl")

println("Units Tests")
@time @test include("internals/units_tests.jl")

toc()


text:
export Companion

immutable Companion{T} <: AbstractArray{T, 2}
    c :: Vector{T}
end
#From polynomial

using Polynomials
#Generate companion matrix from a polynomial

function Companion(P::Poly)
   n = length(P)
   c = Array(Number,n-1)
   d=P.a[n]
   for i=1:n-1
       c[i]=P.a[i]/d
   end
   Companion(c)
end

#Basic property computations
size(C::Companion, r::Int) = (r==1 || r==2) ? length(C.c) :
    throw(ArgumentError("Companion is of rank 2"))

function size(C::Companion)
    n = length(C.c)
    n, n
end

#XXX Inefficient but works
getindex(C::Companion, i, j) = getindex(full(C), i, j)
isassigned(C::Companion, i, j) = isassigned(full(C), i, j)

function full{T}(C::Companion{T})
    M = zeros(T, size(C)...)
    M[:,end]=-C.c
    for i=1:size(C,1)-1
    	M[i+1, i] = one(T)
	end
    M
end

#Linear algebra stuff
function A_mul_B!{T}(C::Companion{T}, b::Vector{T})
	x = b[end]
	y = -C.c[1]*x
	b[2:end] = b[1:end-1]-C.c[2:end]*x
	b[1] = y
    b
end
*{T}(C::Companion{T}, b::Vector{T}) = A_mul_B!(C, copy(b))

function A_mul_B!{T}(A::Matrix{T}, C::Companion{T})
	v = Array(T, size(A,1))
	for i=1:size(A,1)
		v[i] =(A[i,:]*-C.c)[1]
	end
	for i=1:size(A,1), j=1:size(A,2)-1
		A[i,j] = A[i,j+1]
	end
	A[:,end] = v
	A
end
*{T}(A::Matrix{T}, C::Companion{T}) = A_mul_B!(copy(A), C)

function inv{T}(C::Companion{T})
	M = zeros(T, size(C)...)
    for i=1:size(C,1)-1
    	M[i, i+1] = one(T)
	end
	d = M[end, 1] = -one(T)/C.c[1]
	M[1:end-1, 1] = d*C.c[2:end]
    M
end


text:
"""
Register citations with correct bibliography
"""
abstract type RegisterCitations <: Builder.DocumentPipeline end
Selectors.order(::Type{RegisterCitations}) = 4.0
Selectors.matcher(::Type{RegisterCitations}) = true
Selectors.strict(::Type{RegisterCitations}) = false
function Selectors.runner(::Type{RegisterCitations}, doc::Documents.Document)
    Utilities.log(doc, "registering citations.")
    register_citations(doc)
end

function register_citations(doc::Documents.Document)
    for (src, page) in doc.internal.pages
        empty!(page.globals.meta)
        for element in page.elements
            register_citations(page.mapping[element], page, doc)
        end
    end
end

function register_citations(element, page, doc)
    walk_and_replace(page.globals.meta, element) do link
        register_citation(link, page.globals.meta, page, doc)
    end
end

register_citation(other, meta, page, doc) = nothing

function register_citation(link::Markdown.Link, meta, page, doc)
    if occursin(r"^@citation", link.url)
        if link.url == "@citation" && length(link.text) == 1
            citationkey = link.text[1]
            alttext = nothing
        elseif occursin(r"^@citation .+", link.url)
            citationkey = split(link.url)[2]
            alttext = link.text[1]
        else
            Utilities.warn("could not find citation key. Giving up.")
            return false
        end

        # Find the correct bibnode if one exists
        settings = Documents.getplugin(doc, CitationPlugin)
        for bibnode in settings.bibnodes
            if length(bibnode.pages) == 0 || page.source in bibnode.pages
                push!(bibnode.keys, citationkey)
                return CitationNode(citationkey, alttext, bibnode)
            end
        end

        Utilities.warn("citation $citationkey has no bibliography.")
    end

    return nothing
end


text:
struct BDMBasis{T,M,P} <: Space{T}
    geo::M
    fns::Vector{Vector{Shape{T}}}
    pos::Vector{P}
end

refspace(s::BDMBasis{T}) where {T} = BDMRefSpace{T}()
subset(s::BDMBasis,I) = BDMBasis(s.geo, s.fns[I], s.pos[I])

function brezzidouglasmarini(mesh)
    edges = skeleton(mesh, 1)
    cps = cellpairs(mesh, edges, dropjunctionpair=true)
    ids = findall(x -> x>0, cps[2,:])
    brezzidouglasmarini(mesh, cps[:,ids])
end

function brezzidouglasmarini(mesh, cellpairs::Array{Int,2})

    @warn "brezzidouglasmarini(mesh, cellpairs) assumes mesh is oriented"

    @assert size(cellpairs,1) == 2

    T = coordtype(mesh)
    P = vertextype(mesh)

    S = Shape{T}
    F = Vector{Shape{T}}

    nf = 2*size(cellpairs,2)
    fns = Vector{F}(undef, nf)
    pos = Vector{P}(undef, nf)

    for i in axes(cellpairs)[2]
        c1, c2 = cellpairs[:,i]
        cell1 = cells(mesh)[c1]
        cell2 = cells(mesh)[c2]
        e1, e2 = getcommonedge(cell1, cell2)
        @assert e1*e2 < 0
        e1, e2 = abs(e1), abs(e2)
        fns[2*(i-1)+1] = [ S(c1, 2*(e1-1)+1 ,+1.0), S(c2, 2*(e2-1)+2,-1.0)]
        fns[2*(i-1)+2] = [ S(c1, 2*(e1-1)+2 ,+1.0), S(c2, 2*(e2-1)+1,-1.0)]

        v1 = cell1[mod1(e1+1,3)]
        v2 = cell1[mod1(e1+2,3)]
        edge = simplex(mesh.vertices[[v1,v2]]...)
        cntr = cartesian(center(edge))
        pos[2*(i-1)+1] = cntr
        pos[2*(i-1)+2] = cntr
    end

    BDMBasis(mesh, fns, pos)
end