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

using Markdown
using InteractiveUtils

# ╔═╡ 5bb6cf00-9b75-4e74-bc75-6a059bc1b32d
# Učita se paket SymPy
using SymPy

# ╔═╡ efc1cf20-8276-11eb-39ba-bd62a4342e68
md"
# Integrali racionalnih funkcija
"

# ╔═╡ d52c1d09-9bee-4359-8e2c-a969a2f32ea3
# Ova naredba daje popis svih djelova paketa, uglavnom su to funkcije.
# varinfo(SymPy)

# ╔═╡ 8153dabd-5be6-4308-b7c5-1f4ed54722ab
begin
	# Izračunajmo nekoliko integrala.
	x=Sym("x")
	integrate(x^2*sin(x))
end

# ╔═╡ 10b872ad-f070-404c-9ec8-19f0bc2032c5
integrate(exp(-x^2))

# ╔═╡ 056c397f-7303-464e-9ec6-8cf893bbc98c
# Primjer rekurzivne formule. Ovo traje malo duže
integrate(1/(1+x^2)^13)

# ╔═╡ b29b5f68-8978-4b4b-9b29-4840a4097978
# Kraća rekurzivna formula
integrate(1/(1+x^2)^3)

# ╔═╡ 232dcad5-1984-4421-9c1a-001406096217
methods(diff)

# ╔═╡ a2819cd4-d0b3-476a-ab8c-40838db72d54
# Derivacija
diff(x^2*exp(sin(1/x)))

# ╔═╡ 74af42dd-9768-4e91-aa3b-b26c69d6fb68
# Racionalna funkcija trigonometrijskih funkcija (univerzalna trigonometrijska supstitucija)
integrate(1/((2+cos(x))*sin(x)))

# ╔═╡ de024987-9c44-4199-9085-b259980e94a5
md"""
#### Zadatak 1.5.e)
"""

# ╔═╡ 62c5d76c-9eef-40a9-b8c9-446720733f80
integrate(1/(4*sin(x)+3*cos(x)+5))

# ╔═╡ 379bc6e9-7a38-4f3e-8e4a-2ff80917880b
integrate(1/(sin(x)*(2+cos(x)-2*sin(x))))

# ╔═╡ 28dd8649-ed7c-47e5-b7ab-69e99dc0ce86
# Ovo traje beskonačno?
# integrate(cos(x)^3/(sin(x)^2+sin(x)))

# ╔═╡ 12904f05-67c3-4d8a-aefd-e9645e6ef0eb
# Ovo je isti zadatak kao gore, ali traje kratko. Koristi se jednostavnija supstitucija.
integrate((1-sin(x)^2)*cos(x)/(sin(x)^2+sin(x)))

# ╔═╡ bf1829c9-b2fb-4e3b-b92d-38a93513678b
# Ovo traje jako dugo i ne uspije izračunati nego vrati polazni integral
# integrate( (2*tan(x)+3)/(sin(x)^2+2*cos(x)^2))

# ╔═╡ 3732f5f7-ed44-4254-8b16-d9d860d2c04f
md"""
Uz supstituciju $t=\tan(x)$ integral postaje jednostavan:
"""

# ╔═╡ 232927bc-ac81-4225-801f-9da5e5842d6d
begin
	t=Sym("t")
	I₁=integrate( (2*t+3)/(t^2+2))
end

# ╔═╡ 13e53f50-8c7d-403c-85c5-848bf2ccbe94
begin
	# Vratimo supstituciju natrag
	I₂=convert(Function,I₁)
	I₂(tan(x))
end

# ╔═╡ d9ada9db-2651-4895-876b-704e4ed05bcb
md"""
#### Zadatak 1.6
"""

# ╔═╡ f2eb749b-51d6-4875-8e14-1e0f688cdfc7
# Ovo traje beskonačno!
# I₄=integrate((1+sinh(x))/((2+cosh(x))*(3+sinh(x))),x)
# Program treba našu pomoć!!!

# ╔═╡ 244208a5-b153-424b-8305-5b6002fd8186
md"""
__Univerzalna hiperbolna supstitucija__

$$\begin{aligned}
t&=\tanh(\frac{x}{2}), \quad x=2\mathop{\mathrm{atanh}}(t),\quad  dx=\displaystyle\frac{2}{1-t^2}dt,\\
\sinh(x)&=\displaystyle\frac{2t}{1-t^2}, \quad \cosh(x)=\displaystyle\frac{1+t^2}{1-t^2}.\end{aligned}$$

"""

# ╔═╡ 78acc71b-b314-4799-b633-4ab5b7958496
I₄=(1+sinh(x))/((2+cosh(x))*(3+sinh(x)))

# ╔═╡ 463c1ac7-e503-4507-8aca-d671acf04d36
# supstitucija za sinh(x)
I₅=subs(I₄,sinh(x),(2*t)/(1-t^2))

# ╔═╡ 1b47139d-6c06-4227-91d2-1556d2f3a9d2
# Supstitucija za cosh(x)
I₆=subs(I₅,cosh(x),(1+t^2)/(1-t^2))

# ╔═╡ 646c21be-c589-43b6-b648-3d0dea8353db
# Pomnožimo s dx
I₇=I₆*2/(1-t^2)

# ╔═╡ 34320cf9-d978-46ca-ace2-ea638cc602d7
# Integriramo racionalnu funkciju
I₈=integrate(I₇,t)

# ╔═╡ 6e0e6146-fefa-4967-938e-1755a0a04607
# Vratimo supstituciju natrag
I₉=subs(I₈,t,tanh(x/2))

# ╔═╡ a817a9d1-f9d2-4fb9-ad16-2477a0a6dea6

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"

[compat]
SymPy = "~1.0.52"
"""

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

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

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

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

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "a325370b9dd0e6bf5656a6f1a7ae80755f8ccc46"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.7.2"

[[CommonEq]]
git-tree-sha1 = "d1beba82ceee6dc0fce8cb6b80bf600bbde66381"
uuid = "3709ef60-1bee-4518-9f2f-acd86f176c50"
version = "0.2.0"

[[CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

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

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

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

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

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

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

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

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

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

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

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

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

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

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

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

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

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

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

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

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

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

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[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"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

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

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

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

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

[[RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

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

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

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

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

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

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.7.0"

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

[[SymPy]]
deps = ["CommonEq", "CommonSolve", "LinearAlgebra", "Markdown", "PyCall", "RecipesBase", "SpecialFunctions"]
git-tree-sha1 = "1ef257ecbcab8058595a68ca36a6844b41babcbd"
uuid = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
version = "1.0.52"

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

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

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

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

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

[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

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

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

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

# ╔═╡ Cell order:
# ╟─efc1cf20-8276-11eb-39ba-bd62a4342e68
# ╠═5bb6cf00-9b75-4e74-bc75-6a059bc1b32d
# ╠═d52c1d09-9bee-4359-8e2c-a969a2f32ea3
# ╠═8153dabd-5be6-4308-b7c5-1f4ed54722ab
# ╠═10b872ad-f070-404c-9ec8-19f0bc2032c5
# ╠═056c397f-7303-464e-9ec6-8cf893bbc98c
# ╠═b29b5f68-8978-4b4b-9b29-4840a4097978
# ╠═232dcad5-1984-4421-9c1a-001406096217
# ╠═a2819cd4-d0b3-476a-ab8c-40838db72d54
# ╠═74af42dd-9768-4e91-aa3b-b26c69d6fb68
# ╟─de024987-9c44-4199-9085-b259980e94a5
# ╠═62c5d76c-9eef-40a9-b8c9-446720733f80
# ╠═379bc6e9-7a38-4f3e-8e4a-2ff80917880b
# ╠═28dd8649-ed7c-47e5-b7ab-69e99dc0ce86
# ╠═12904f05-67c3-4d8a-aefd-e9645e6ef0eb
# ╠═bf1829c9-b2fb-4e3b-b92d-38a93513678b
# ╟─3732f5f7-ed44-4254-8b16-d9d860d2c04f
# ╠═232927bc-ac81-4225-801f-9da5e5842d6d
# ╠═13e53f50-8c7d-403c-85c5-848bf2ccbe94
# ╟─d9ada9db-2651-4895-876b-704e4ed05bcb
# ╠═f2eb749b-51d6-4875-8e14-1e0f688cdfc7
# ╟─244208a5-b153-424b-8305-5b6002fd8186
# ╠═78acc71b-b314-4799-b633-4ab5b7958496
# ╠═463c1ac7-e503-4507-8aca-d671acf04d36
# ╠═1b47139d-6c06-4227-91d2-1556d2f3a9d2
# ╠═646c21be-c589-43b6-b648-3d0dea8353db
# ╠═34320cf9-d978-46ca-ace2-ea638cc602d7
# ╠═6e0e6146-fefa-4967-938e-1755a0a04607
# ╠═a817a9d1-f9d2-4fb9-ad16-2477a0a6dea6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: ivanslapnicar/Matematika

max_issues_repo_name: ivanslapnicar/Matematika

max_forks_repo_name: ivanslapnicar/Matematika


text:
using Test
using PlayingCards
using TexasHoldem

@testset "Game: play! (3 Bot5050's)" begin
    for n in 1:n_fuzz
        play!(Game(ntuple(i->Player(Bot5050(), i), 3)))
    end
end

@testset "Game: play! (10 Bot5050's)" begin
    for n in 1:n_fuzz_10_players
        play!(Game(ntuple(i->Player(Bot5050(), i), 10)))
    end
end

@testset "Game: tournament! (2 Bot5050's)" begin
    for n in 1:n_fuzz
        tournament!(Game(ntuple(i->Player(Bot5050(), i; bank_roll = 6), 2)))
    end
end

@testset "Game: tournament! (3 Bot5050's)" begin
    for n in 1:n_fuzz_3_players
        tournament!(Game(ntuple(i->Player(Bot5050(), i; bank_roll = 6), 3)))
    end
end

@testset "Game: tournament! (10 Bot5050's)" begin
    for n in 1:n_fuzz_10_players
        tournament!(Game(ntuple(i->Player(Bot5050(), i; bank_roll = 6), 10)))
    end
end

ext: jl


text:
function Base.show(io::IO, m::MIME"text/plain", t::Template)
    println(io, "Template:")
    foreach(fieldnames(Template)) do n
        n === :plugins || println(io, repeat(' ', 2), n, ": ", show_field(getfield(t, n)))
    end
    if isempty(t.plugins)
        print(io, "  plugins: None")
    else
        print(io, repeat(' ', 2), "plugins:")
        foreach(sort(t.plugins; by=string)) do p
            println(io)
            show(IOContext(io, :indent => 4), m, p)
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", p::T) where T <: Plugin
    indent = get(io, :indent, 0)
    print(io, repeat(' ', indent), nameof(T))
    ns = fieldnames(T)
    isempty(ns) || print(io, ":")
    foreach(ns) do n
        println(io)
        print(io, repeat(' ', indent + 2), n, ": ", show_field(getfield(p, n)))
    end
end

show_field(x) = repr(x)
if Sys.iswindows()
    show_field(x::AbstractString) = replace(repr(contractuser(x)), "\\\\" => "\\")
else
    show_field(x::AbstractString) = repr(contractuser(x))
end


text:
function randinit(
  n::Integer, # Number of elements to initialize
  m::Integer, # Number of codebooks
  h::Integer) # Number of elements per codebook
  B = convert( Matrix{Int16}, rand(1:h, m, n) );
  return B
end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
# Solution to Q4e (implement LogSoftMax) using expression derived in Q4a-c
import Base

type LogSoftMax{T}
    # Buffer of the outputs of fprop, needed during bprop.
    logσ::Matrix{T}

    LogSoftMax() = new(zeros(T, 0, 0))
end

function fprop!{T}(sm::LogSoftMax{T}, x::Matrix{T})
    # Subtract the largest input for stability (Solution to Q4a).
    x = x .- maximum(x, 2)

    # Compute the simplified log-softmax formula.
    return sm.logσ = x .- log.(sum(exp.(x), 2))
end

function bprop!{T}(sm::LogSoftMax{T}, ∂y::Matrix{T})
    # Return the gradient wrt. the input, see Q4b and Q4c.
    σ = exp.(sm.logσ)
    return ∂y - σ .* sum(∂y, 2)
end

params(::LogSoftMax) = ()
grads(::LogSoftMax) = ()

Base.show(io::IO, l::LogSoftMax) = print(io, "LogSoftMax")


text: abstract type AggregateState end


text:
module Day06

using DataStructures: counter

export day06

const TIME_TO_REPRODUCE = 7
const TIME_TO_MATURITY = 9

function readInput(inputfile="inputs/input06.txt")
    return parse.(UInt8, split(readline(inputfile), ","))
end

# Let L(x,t) be the population resulting starting with a SINGLE lanternfish with counter x
# after t days
# Use relations for L(x,t)
# - L(x, t) = L(0, t-x)
# - L(0, k) = 1 (for k ≤ 0)
# - L(0, t) = L(0, t-TIME_TO_REPRODUCE) + L(0, t-TIME_TO_MATURITY)

################################## more Julian solution ##################################
## a bit slower than the lower one
## but that's just due to the use of DataStructures.counter
## without that the 2 solution run in pretty much identical time

struct LaternfishPopulationPrediction
    precomputed::Vector{Int64}
    upto::Int64
    function LaternfishPopulationPrediction(upto)
        new(_precompute_lanternfish_popuplation(upto), upto)
    end
end

function _precompute_lanternfish_popuplation(upto)
    populations = Vector{Int64}(undef, upto+TIME_TO_MATURITY)
    populations[1:TIME_TO_MATURITY] .= 1
    for index in TIME_TO_MATURITY .+ (1:upto)
        populations[index] = populations[index-TIME_TO_REPRODUCE] + populations[index-TIME_TO_MATURITY]
    end
    return populations
end

function Base.getindex(pop::LaternfishPopulationPrediction, time)
    time <= pop.upto || error("time larger than simulated! ($time > $(pop.upto))")
    return pop.precomputed[time+9]
end

function population(counts, time, precomputed=nothing)
    precomputed === nothing && (precomputed = LaternfishPopulationPrediction(time))
    total = 0
    for (ini, count) in counts
        total += precomputed[time-ini]*count
    end
    return total
end

function day06(inputdata = readInput("inputs/input06.txt"))
    pop = LaternfishPopulationPrediction(256)
    counts = counter(inputdata)
    return population(counts, 80, pop), population(counts, 256, pop)
end

######################## using Memoize.jl and optimizing counter ########################

using Memoize

L(x,t) = _L(t-x)
@memoize function _L(t)
    if t <= 0
        return 1
    else
        return _L(t-TIME_TO_REPRODUCE) + _L(t-TIME_TO_MATURITY)
    end
end

_L(256) # warm cache

function population_alt(counts, time)
    total = 0
    for (ini, count) in counts
        total += L(ini-1, time)*count
    end
    return total
end

function day06_alt(inputdata = readInput("inputs/input06.txt"))
    counts = enumerate(_mycount(inputdata)) # ~4x faster than counter
    return population_alt(counts, 80), population_alt(counts, 256)
end

function _mycount(inputdata)
    ret = zeros(Int64, maximum(inputdata)+1)
    for i in inputdata
        ret[i+1] += 1
    end
    return ret
end

end #module


text:
function ConjugateGradients(d,operators,parameters;Niter=10,mu=0,tol=1.0e-15)
    # Conjugate Gradients following Algorithm 2 from Scales, 1987.
    # The user provides an array of linear operators. Ensure linear operator(s) pass the dot product.

    cost = Float64[]
    r = copy(d)
    g = LinearOperator(r,operators,parameters,adj=true)
    m = zeros(g)
    s = copy(g)
    gamma = InnerProduct(g,g)
    gamma00 = gamma
    cost0 = InnerProduct(r,r)
    push!(cost,1.0)
    for iter = 1 : Niter
	t = LinearOperator(s,operators,parameters,adj=false)
	delta = InnerProduct(t,t) + mu*InnerProduct(s,s)
	if delta <= tol
#	    println("delta reached tolerance, ending at iteration ",iter)
	    break;
	end
	alpha = gamma/delta
	m = m + alpha*s
	r = r - alpha*t
	g = LinearOperator(r,operators,parameters,adj=true)
	g = g - mu*m
	gamma0 = copy(gamma)
	gamma = InnerProduct(g,g)
        cost1 = InnerProduct(r,r) + mu*InnerProduct(m,m)
        push!(cost,cost1/cost0)
	beta = gamma/gamma0
	s = beta*s + g
	if (sqrt(gamma) <= sqrt(gamma00) * tol)
	    println("tolerance reached, ending at iteration ",iter)
	    break;
	end
    end

    return m, cost
end

function ConjugateGradients(m::String,d::String,operators,parameters,cost_file::String;Niter=10,mu=0,tol=1.0e-15)
    # Conjugate Gradients following Algorithm 2 from Scales, 1987.
    # The user provides an array of linear operators. Ensure linear operator(s) pass the dot product.

    cost = Float64[]
    rand_string = string(round(Int,rand()*100000))
    g = join(["tmp_CG_g_",rand_string])
    s = join(["tmp_CG_s_",rand_string])
    r = join(["tmp_CG_r_",rand_string])
    t = join(["tmp_CG_t_",rand_string])
    SeisCopy(d,r)
    fp = open(cost_file,"w")
    @compat write(fp,join(["began execution at: ",Libc.strftime(time()),"\n"]))
    close(fp)
    LinearOperator(g,r,operators,parameters,adj=true)
    SeisCopy(g,s)
    SeisCopy(g,m)
    CGStep(m,g,a=0.0,b=0.0)
    gamma = InnerProduct(g,g)
    gamma00 = gamma
    cost0 = InnerProduct(r,r)
    push!(cost,1.0)
    fp = open(cost_file,"a")
    write(fp,join([string(cost[1]),"\n"]))
    close(fp)
    for iter = 1 : Niter
	LinearOperator(s,t,operators,parameters,adj=false)
	delta = InnerProduct(t,t) + mu*InnerProduct(s,s)
	alpha = gamma/delta
	CGStep(m,s,a=1.0,b=alpha)
	CGStep(r,t,a=1.0,b=-alpha)
	LinearOperator(g,r,operators,parameters,adj=true)
	CGStep(g,m,a=1.0,b=-mu)
	gamma0 = copy(gamma)
	gamma = InnerProduct(g,g)
        cost1 = InnerProduct(r,r) + mu*InnerProduct(m,m)
        push!(cost,cost1/cost0)
	fp = open(cost_file,"a")
	write(fp,join([string(cost[iter+1]),"\n"]))
	close(fp)
	beta = gamma/gamma0
	CGStep(s,g,a=beta,b=1.0)
	if (sqrt(gamma) <= sqrt(gamma00) * tol)
	    println("tolerance reached, ending at iteration ",iter)
	    break;
	end
    end
    SeisRemove(g)
    SeisRemove(s)
    SeisRemove(r)
    SeisRemove(t)

end

function ConjugateGradients(m::Array{String,1},d::Array{String,1},operators,parameters,cost_file::String;Niter=10,mu=[0.0;0.0],tol=1.0e-15)
    # Conjugate Gradients following Algorithm 2 from Scales, 1987.
    # The user provides an array of linear operators. Ensure linear operator(s) pass the dot product.

    cost = Float64[]
    rand_string = string(round(Int,rand()*100000))
    g = [join(["tmp_CG_g1_",rand_string]);join(["tmp_CG_g2_",rand_string])]
    s = [join(["tmp_CG_s1_",rand_string]);join(["tmp_CG_s2_",rand_string])]
    r = [join(["tmp_CG_r1_",rand_string]);join(["tmp_CG_r2_",rand_string])]
    t = [join(["tmp_CG_t1_",rand_string]);join(["tmp_CG_t2_",rand_string])]
    SeisCopy(d,r)
    fp = open(cost_file,"w")
    @compat write(fp,join(["began execution at: ",Libc.strftime(time()),"\n"]))
    close(fp)
    LinearOperator(g,r,operators,parameters,adj=true)
    SeisCopy(g,s)
    SeisCopy(g,m)
    CGStep(m,g,a=[0.0;0.0],b=[0.0;0.0])
    gamma = InnerProduct(g,g)
    gamma00 = gamma
    cost0 = InnerProduct(r,r)
    push!(cost,1.0)
    fp = open(cost_file,"a")
    write(fp,join([string(cost[1]),"\n"]))
    close(fp)
    for iter = 1 : Niter
	LinearOperator(s,t,operators,parameters,adj=false)
	delta = InnerProduct(t,t) + mu[1]*InnerProduct(s[1],s[1]) + mu[2]*InnerProduct(s[2],s[2])
	alpha = gamma/delta
	CGStep(m,s,a=[1.0;1.0],b=[alpha;alpha])
	CGStep(r,t,a=[1.0;1.0],b=[-alpha;-alpha])
	LinearOperator(g,r,operators,parameters,adj=true)
	CGStep(g,m,a=[1.0;1.0],b=[-mu[1];-mu[2]])
	gamma0 = copy(gamma)
	gamma = InnerProduct(g,g)
        cost1 = InnerProduct(r,r) + mu[1]*InnerProduct(m[1],m[1]) + mu[2]*InnerProduct(m[2],m[2])
        push!(cost,cost1/cost0)
	fp = open(cost_file,"a")
	write(fp,join([string(cost[iter+1]),"\n"]))
	close(fp)
	beta = gamma/gamma0
	CGStep(s,g,a=[beta;beta],b=[1.0;1.0])
	if (sqrt(gamma) <= sqrt(gamma00) * tol)
	    println("tolerance reached, ending at iteration ",iter)
	    break;
	end
    end
    SeisRemove(g);
    SeisRemove(s);
    SeisRemove(r);
    SeisRemove(t);

end


text:
const invpio2 =  6.36619772367581382433e-01 # 0x3FE45F30, 0x6DC9C883 */
const pio2_1  =  1.57079632673412561417e+00 # 0x3FF921FB, 0x54400000 */
const pio2_1t =  6.07710050650619224932e-11 # 0x3DD0B461, 0x1A626331 */
const pio2_2  =  6.07710050630396597660e-11 # 0x3DD0B461, 0x1A600000 */
const pio2_2t =  2.02226624879595063154e-21 # 0x3BA3198A, 0x2E037073 */
const pio2_3  =  2.02226624871116645580e-21 # 0x3BA3198A, 0x2E000000 */
const pio2_3t =  8.47842766036889956997e-32 # 0x397B839A, 0x252049C1 */

function cody_waite_2c_pio2(x, signed_k)
    z = x - signed_k*pio2_1;
    y1 = z - signed_k*pio2_1t;
    y2 = (z - y1) - signed_k*pio2_1t;
    TwicePrecision(y1, y2)
end

function cody_waite_ext_pio2(x, xʰ⁺)
    fn = x*invpio2+0x1.8p52
    fn = fn-0x1.8p52
    n  = Int32(fn)
    r  = x-fn*pio2_1
    w  = fn*pio2_1t # 1st round good to 85 bit
    j  = xʰ⁺>>20
    y1 = r-w
    high = highword(y1)
    i = j-((high>>20)&0x7ff)
    if i>16  # 2nd iteration needed, good to 118
        t  = r
        w  = fn*pio2_2
        r  = t-w
        w  = fn*pio2_2t-((t-r)-w)
        y1 = r-w
        high = highword(y1)
        i = j-((high>>20)&0x7ff)
        if i>49 # 3rd iteration need, 151 bits acc
            t  = r # will cover all possible cases
            w  = fn*pio2_3
            r  = t-w
            w  = fn*pio2_3t-((t-r)-w)
            y1 = r-w
        end
    end
    y2 = (r-y1)-w
    return Int(n), TwicePrecision(y1, y2)
end


text:
import .ExpressionExplorer: external_package_names
import .PkgCompat
import .PkgCompat: select, is_stdlib

const tiers = [
	Pkg.PRESERVE_ALL,
	Pkg.PRESERVE_DIRECT,
	Pkg.PRESERVE_SEMVER,
	Pkg.PRESERVE_NONE,
]

const pkg_token = Token()

# This list appears multiple times in our codebase. Be sure to match edits everywhere.
function use_plutopkg(topology::NotebookTopology)
    !any(values(topology.nodes)) do node
        Symbol("Pkg.activate") ∈ node.references ||
        Symbol("Pkg.API.activate") ∈ node.references ||
        Symbol("Pkg.develop") ∈ node.references ||
        Symbol("Pkg.API.develop") ∈ node.references ||
        Symbol("Pkg.add") ∈ node.references ||
        Symbol("Pkg.API.add") ∈ node.references ||
        # https://juliadynamics.github.io/DrWatson.jl/dev/project/#DrWatson.quickactivate
        Symbol("quickactivate") ∈ node.references ||
        Symbol("@quickactivate") ∈ node.references ||
        Symbol("DrWatson.@quickactivate") ∈ node.references ||
        Symbol("DrWatson.quickactivate") ∈ node.references
    end
end

function external_package_names(topology::NotebookTopology)::Set{Symbol}
    union!(Set{Symbol}(), external_package_names.(c.module_usings_imports for c in values(topology.codes))...)
end

PkgCompat.project_file(notebook::Notebook) = PkgCompat.project_file(PkgCompat.env_dir(notebook.nbpkg_ctx))
PkgCompat.manifest_file(notebook::Notebook) = PkgCompat.manifest_file(PkgCompat.env_dir(notebook.nbpkg_ctx))

"""
```julia
sync_nbpkg_core(notebook::Notebook; on_terminal_output::Function=((args...) -> nothing))
```

Update the notebook package environment to match the notebook's code. This will:
- Add packages that should be added (because they are imported in a cell).
- Remove packages that are no longer needed.
- Make sure that the environment is instantiated.
- Detect the use of `Pkg.activate` and enable/disabled nbpkg accordingly.
"""
function sync_nbpkg_core(notebook::Notebook; on_terminal_output::Function=((args...) -> nothing))

    👺 = false

    use_plutopkg_old = notebook.nbpkg_ctx !== nothing
    use_plutopkg_new = use_plutopkg(notebook.topology)

    if !use_plutopkg_old && use_plutopkg_new
        @debug "Started using PlutoPkg!! HELLO reproducibility!"

        👺 = true
        notebook.nbpkg_ctx = PkgCompat.create_empty_ctx()
    end
    if use_plutopkg_old && !use_plutopkg_new
        @debug "Stopped using PlutoPkg 💔😟😢"

        no_packages_loaded_yet = (
            notebook.nbpkg_restart_required_msg === nothing &&
            notebook.nbpkg_restart_recommended_msg === nothing &&
            all(PkgCompat.is_stdlib, keys(PkgCompat.project(notebook.nbpkg_ctx).dependencies))
        )
        👺 = !no_packages_loaded_yet
        notebook.nbpkg_ctx = nothing
    end

    if notebook.nbpkg_ctx !== nothing
        PkgCompat.mark_original!(notebook.nbpkg_ctx)

        old_packages = String.(keys(PkgCompat.project(notebook.nbpkg_ctx).dependencies))
        new_packages = String.(external_package_names(notebook.topology)) # search all cells for imports and usings

        removed = setdiff(old_packages, new_packages)
        added = setdiff(new_packages, old_packages)

        iolistener = let
            busy_packages = notebook.nbpkg_ctx_instantiated ? added : new_packages
            IOListener(callback=(s -> on_terminal_output(busy_packages, s)))
        end

        # We remember which Pkg.Types.PreserveLevel was used. If it's too low, we will recommend/require a notebook restart later.
        local used_tier = Pkg.PRESERVE_ALL

        if !isready(pkg_token)
            println(iolistener.buffer, "Waiting for other notebooks to finish Pkg operations...")
            trigger(iolistener)
        end

        can_skip = isempty(removed) && isempty(added) && notebook.nbpkg_ctx_instantiated

        if !can_skip
            return withtoken(pkg_token) do
                PkgCompat.refresh_registry_cache()

                if !notebook.nbpkg_ctx_instantiated
                    notebook.nbpkg_ctx = PkgCompat.clear_stdlib_compat_entries(notebook.nbpkg_ctx)
                    PkgCompat.withio(notebook.nbpkg_ctx, IOContext(iolistener.buffer, :color => true)) do
                        withinteractive(false) do
                            try
                                Pkg.resolve(notebook.nbpkg_ctx)
                            catch e
                                @warn "Failed to resolve Pkg environment. Removing Manifest and trying again..." exception=e
                                reset_nbpkg(notebook; keep_project=true, save=false, backup=false)
                                Pkg.resolve(notebook.nbpkg_ctx)
                            end
                        end
                    end
                end

                to_remove = filter(removed) do p
                    haskey(PkgCompat.project(notebook.nbpkg_ctx).dependencies, p)
                end
                if !isempty(to_remove)
                    @debug to_remove
                    # See later comment
                    mkeys() = Set(filter(!is_stdlib, [m.name for m in values(PkgCompat.dependencies(notebook.nbpkg_ctx))]))
                    old_manifest_keys = mkeys()

                    Pkg.rm(notebook.nbpkg_ctx, [
                        Pkg.PackageSpec(name=p)
                        for p in to_remove
                    ])

                    # We record the manifest before and after, to prevent recommending a reboot when nothing got removed from the manifest (e.g. when removing GR, but leaving Plots), or when only stdlibs got removed.
                    new_manifest_keys = mkeys()

                    # TODO: we might want to upgrade other packages now that constraints have loosened? Does this happen automatically?
                end

                # TODO: instead of Pkg.PRESERVE_ALL, we actually want:
                # "Pkg.PRESERVE_DIRECT, but preserve exact verisons of Base.loaded_modules"

                to_add = filter(PkgCompat.package_exists, added)

                if !isempty(to_add)
                    @debug to_add
                    startlistening(iolistener)

                    PkgCompat.withio(notebook.nbpkg_ctx, IOContext(iolistener.buffer, :color => true)) do
                        withinteractive(false) do
                            # We temporarily clear the "semver-compatible" [deps] entries, because Pkg already respects semver, unless it doesn't, in which case we don't want to force it.
                            notebook.nbpkg_ctx = PkgCompat.clear_auto_compat_entries(notebook.nbpkg_ctx)

                            try
                                for tier in [
                                    Pkg.PRESERVE_ALL,
                                    Pkg.PRESERVE_DIRECT,
                                    Pkg.PRESERVE_SEMVER,
                                    Pkg.PRESERVE_NONE,
                                ]
                                    used_tier = tier

                                    try
                                        Pkg.add(notebook.nbpkg_ctx, [
                                            Pkg.PackageSpec(name=p)
                                            for p in to_add
                                        ]; preserve=used_tier)

                                        break
                                    catch e
                                        if used_tier == Pkg.PRESERVE_NONE
                                            # give up
                                            rethrow(e)
                                        end
                                    end
                                end
                            finally
                                notebook.nbpkg_ctx = PkgCompat.write_auto_compat_entries(notebook.nbpkg_ctx)
                            end

                            # Now that Pkg is set up, the notebook process will call `using Package`, which can take some time. We write this message to the io, to notify the user.
                            println(iolistener.buffer, "\e[32m\e[1mLoading\e[22m\e[39m packages...")
                        end
                    end

                    @debug "PlutoPkg done"
                end

                should_instantiate = !notebook.nbpkg_ctx_instantiated || !isempty(to_add) || !isempty(to_remove)
                if should_instantiate
                    startlistening(iolistener)
                    PkgCompat.withio(notebook.nbpkg_ctx, IOContext(iolistener.buffer, :color => true)) do
                        @debug "Instantiating"

                        # Pkg.instantiate assumes that the environment to be instantiated is active, so we will have to modify the LOAD_PATH of this Pluto server
                        # We could also run the Pkg calls on the notebook process, but somehow I think that doing it on the server is more charming, though it requires this workaround.
                        env_dir = PkgCompat.env_dir(notebook.nbpkg_ctx)
                        pushfirst!(LOAD_PATH, env_dir)

                        # update registries if this is the first time
                        PkgCompat.update_registries(notebook.nbpkg_ctx)
                        # instantiate without forcing registry update
                        PkgCompat.instantiate(notebook.nbpkg_ctx; update_registry=false)

                        @assert LOAD_PATH[1] == env_dir
                        popfirst!(LOAD_PATH)
                    end
                    notebook.nbpkg_ctx_instantiated = true
                end

                stoplistening(iolistener)

                return (
                    did_something=👺 || (
                        should_instantiate || (use_plutopkg_old != use_plutopkg_new)
                    ),
                    used_tier=used_tier,
                    # changed_versions=Dict{String,Pair}(),
                    restart_recommended=👺 || (
                        (!isempty(to_remove) && old_manifest_keys != new_manifest_keys) ||
                        used_tier != Pkg.PRESERVE_ALL
                    ),
                    restart_required=👺 || (
                        used_tier ∈ [Pkg.PRESERVE_SEMVER, Pkg.PRESERVE_NONE]
                    ),
                )
            end
        end
    end
    return (
        did_something=👺 || (use_plutopkg_old != use_plutopkg_new),
        used_tier=Pkg.PRESERVE_ALL,
        # changed_versions=Dict{String,Pair}(),
        restart_recommended=👺 || false,
        restart_required=👺 || false,
    )
end

"""
```julia
sync_nbpkg(session::ServerSession, notebook::Notebook; save::Bool=true)
```

In addition to the steps performed by [`sync_nbpkg_core`](@ref):
- Capture terminal outputs and store them in the `notebook`
- Update the clients connected to `notebook`
- `try` `catch` and reset the package environment on failure.
"""
function sync_nbpkg(session, notebook; save::Bool=true)
	try
		pkg_result = withtoken(notebook.executetoken) do
			function iocallback(pkgs, s)
				notebook.nbpkg_busy_packages = pkgs
				for p in pkgs
					notebook.nbpkg_terminal_outputs[p] = s
				end
                update_nbpkg_cache!(notebook)
				send_notebook_changes!(ClientRequest(session=session, notebook=notebook))
			end
			sync_nbpkg_core(notebook; on_terminal_output=iocallback)
		end

		if pkg_result.did_something
			@debug "PlutoPkg: success!" pkg_result

			if pkg_result.restart_recommended
				@debug "PlutoPkg: Notebook restart recommended"
				notebook.nbpkg_restart_recommended_msg = "Yes, something changed during regular sync."
			end
			if pkg_result.restart_required
				@debug "PlutoPkg: Notebook restart REQUIRED"
				notebook.nbpkg_restart_required_msg = "Yes, something changed during regular sync."
			end

			notebook.nbpkg_busy_packages = String[]
            update_nbpkg_cache!(notebook)
			send_notebook_changes!(ClientRequest(session=session, notebook=notebook))
			save && save_notebook(session, notebook)
		end
	catch e
		bt = catch_backtrace()
		old_packages = try String.(keys(PkgCompat.project(notebook.nbpkg_ctx).dependencies)); catch; ["unknown"] end
		new_packages = try String.(external_package_names(notebook.topology)); catch; ["unknown"] end
		@warn """
		PlutoPkg: Failed to add/remove packages! Resetting package environment...
		""" PLUTO_VERSION VERSION old_packages new_packages exception=(e, bt)
		# TODO: send to user

		error_text = sprint(showerror, e, bt)
		for p in notebook.nbpkg_busy_packages
            old = get(notebook.nbpkg_terminal_outputs, p, "")
			notebook.nbpkg_terminal_outputs[p] = old * "\n\n\nPkg error!\n\n" * error_text
		end
		notebook.nbpkg_busy_packages = String[]
        update_nbpkg_cache!(notebook)
		send_notebook_changes!(ClientRequest(session=session, notebook=notebook))

		# Clear the embedded Project and Manifest and require a restart from the user.
		reset_nbpkg(notebook; keep_project=false, save=save)
		notebook.nbpkg_restart_required_msg = "Yes, because sync_nbpkg_core failed. \n\n$(error_text)"
        notebook.nbpkg_ctx_instantiated = false
        update_nbpkg_cache!(notebook)
		send_notebook_changes!(ClientRequest(session=session, notebook=notebook))

		save && save_notebook(session, notebook)
	end
end

function writebackup(notebook::Notebook)
    backup_path = backup_filename(notebook.path)
    Pluto.readwrite(notebook.path, backup_path)

    @info "Backup saved to" backup_path

    backup_path
end

function reset_nbpkg(notebook::Notebook; keep_project::Bool=false, backup::Bool=true, save::Bool=true)
    backup && save && writebackup(notebook)

    if notebook.nbpkg_ctx !== nothing
        p = PkgCompat.project_file(notebook)
        m = PkgCompat.manifest_file(notebook)
        keep_project || (isfile(p) && rm(p))
        isfile(m) && rm(m)

        notebook.nbpkg_ctx = PkgCompat.load_ctx(PkgCompat.env_dir(notebook.nbpkg_ctx))
    else
        notebook.nbpkg_ctx = use_plutopkg(notebook.topology) ? PkgCompat.create_empty_ctx() : nothing
    end

    save && save_notebook(notebook)
end

function update_nbpkg_core(notebook::Notebook; level::Pkg.UpgradeLevel=Pkg.UPLEVEL_MAJOR, on_terminal_output::Function=((args...) -> nothing))
    if notebook.nbpkg_ctx !== nothing
        PkgCompat.mark_original!(notebook.nbpkg_ctx)

        old_packages = String.(keys(PkgCompat.project(notebook.nbpkg_ctx).dependencies))

        iolistener = let
            # we don't know which packages will be updated, so we send terminal output to all installed packages
            IOListener(callback=(s -> on_terminal_output(old_packages, s)))
        end

        # We remember which Pkg.Types.PreserveLevel was used. If it's too low, we will recommend/require a notebook restart later.
        local used_tier = Pkg.PRESERVE_ALL

        if !isready(pkg_token)
            println(iolistener.buffer, "Waiting for other notebooks to finish Pkg operations...")
            trigger(iolistener)
        end

        return withtoken(pkg_token) do
            PkgCompat.refresh_registry_cache()

            if !notebook.nbpkg_ctx_instantiated
                notebook.nbpkg_ctx = PkgCompat.clear_stdlib_compat_entries(notebook.nbpkg_ctx)
                PkgCompat.withio(notebook.nbpkg_ctx, IOContext(iolistener.buffer, :color => true)) do
                    withinteractive(false) do
                        try
                            Pkg.resolve(notebook.nbpkg_ctx)
                        catch e
                            @warn "Failed to resolve Pkg environment. Removing Manifest and trying again..." exception=e
                            reset_nbpkg(notebook; keep_project=true, save=false, backup=false)
                            Pkg.resolve(notebook.nbpkg_ctx)
                        end
                    end
                end
            end

            startlistening(iolistener)

            PkgCompat.withio(notebook.nbpkg_ctx, IOContext(iolistener.buffer, :color => true)) do
                # We temporarily clear the "semver-compatible" [deps] entries, because it is difficult to update them after the update 🙈. TODO
                notebook.nbpkg_ctx = PkgCompat.clear_auto_compat_entries(notebook.nbpkg_ctx)

                try
                    ###
                    Pkg.update(notebook.nbpkg_ctx; level=level)
                    ###
                finally
                    notebook.nbpkg_ctx = PkgCompat.write_auto_compat_entries(notebook.nbpkg_ctx)
                end
            end

            stoplistening(iolistener)

            🐧 = !PkgCompat.is_original(notebook.nbpkg_ctx)
            (
                did_something=🐧,
                restart_recommended=🐧,
                restart_required=🐧,
            )
        end
    end
    (
        did_something=false,
        restart_recommended=false,
        restart_required=false,
    )
end

function update_nbpkg(session, notebook::Notebook; level::Pkg.UpgradeLevel=Pkg.UPLEVEL_MAJOR, backup::Bool=true, save::Bool=true)
    if backup && save
        bp = writebackup(notebook)
    end

    try
		pkg_result = withtoken(notebook.executetoken) do
            original_outputs = deepcopy(notebook.nbpkg_terminal_outputs)
			function iocallback(pkgs, s)
				notebook.nbpkg_busy_packages = pkgs
				for p in pkgs
                    original = get(original_outputs, p, "")
					notebook.nbpkg_terminal_outputs[p] = original * "\n\n" * s
				end
                update_nbpkg_cache!(notebook)
				send_notebook_changes!(ClientRequest(session=session, notebook=notebook))
			end
			update_nbpkg_core(notebook; level=level, on_terminal_output=iocallback)
		end

		if pkg_result.did_something
			if pkg_result.restart_recommended
				@debug "PlutoPkg: Notebook restart recommended"
				notebook.nbpkg_restart_recommended_msg = "Yes, something changed during regular update_nbpkg."
			end
			if pkg_result.restart_required
				@debug "PlutoPkg: Notebook restart REQUIRED"
				notebook.nbpkg_restart_required_msg = "Yes, something changed during regular update_nbpkg."
			end
		else
            isfile(bp) && rm(bp)
        end
	finally
		notebook.nbpkg_busy_packages = String[]
        update_nbpkg_cache!(notebook)
		send_notebook_changes!(ClientRequest(session=session, notebook=notebook))
		save && save_notebook(session, notebook)
	end
end

nbpkg_cache(ctx::Union{Nothing,PkgContext}) = ctx === nothing ? Dict{String,String}() : Dict{String,String}(
    x => string(PkgCompat.get_manifest_version(ctx, x)) for x in keys(PkgCompat.project(ctx).dependencies)
)

function update_nbpkg_cache!(notebook::Notebook)
    notebook.nbpkg_installed_versions_cache = nbpkg_cache(notebook.nbpkg_ctx)
    notebook
end

function is_nbpkg_equal(a::Union{Nothing,PkgContext}, b::Union{Nothing,PkgContext})::Bool
    if (a isa Nothing) != (b isa Nothing)
        the_other = something(a, b)

        ptoml_contents = PkgCompat.read_project_file(the_other)
        the_other_is_empty = isempty(strip(ptoml_contents))

        if the_other_is_empty
            # then both are essentially 'empty' environments, i.e. equal
            true
        else
            # they are different
            false
        end
    elseif a isa Nothing
        true
    else
        ptoml_contents_a = strip(PkgCompat.read_project_file(a))
        ptoml_contents_b = strip(PkgCompat.read_project_file(b))

        if ptoml_contents_a == ptoml_contents_b == ""
            true
        else
            mtoml_contents_a = strip(PkgCompat.read_project_file(a))
            mtoml_contents_b = strip(PkgCompat.read_project_file(b))

            (ptoml_contents_a == ptoml_contents_b) && (mtoml_contents_a == mtoml_contents_b)
        end
    end
end

const is_interactive_defined = isdefined(Base, :is_interactive) && !Base.isconst(Base, :is_interactive)
function withinteractive(f::Function, value::Bool)
    old_value = isinteractive()
    @static if is_interactive_defined
        Core.eval(Base, :(is_interactive = $value))
    end
    try
        f()
    finally
        @static if is_interactive_defined
            Core.eval(Base, :(is_interactive = $old_value))
        end
    end
end

"A polling system to watch for writes to an IOBuffer. Up-to-date content will be passed as string to the `callback` function."
Base.@kwdef struct IOListener
    callback::Function
    buffer::IOBuffer=IOBuffer()
    interval::Real=1.0/60
    running::Ref{Bool}=Ref(false)
    last_size::Ref{Int}=Ref(-1)
end
function trigger(listener::IOListener)
    new_size = listener.buffer.size
    if new_size > listener.last_size[]
        listener.last_size[] = new_size
        new_contents = String(listener.buffer.data[1:new_size])
        listener.callback(new_contents)
    end
end
function startlistening(listener::IOListener)
    if !listener.running[]
        listener.running[] = true
        @async while listener.running[]
            trigger(listener)
            sleep(listener.interval)
        end
    end
end
function stoplistening(listener::IOListener)
    if listener.running[]
        listener.running[] = false
        trigger(listener)
    end
end


text:
function upgradeToBot(token = token)
    isempty(token) && error(
        "Need to set the token via `setDefaultToken` or pass it as an argument",
    )
    headers = Dict("Authorization" => "Bearer $token")
    r = HTTP.request(
        "POST",
        "https://lichess.org/api/bot/account/upgrade",
        headers,
        "";
        status_exception = false,
    )
    if r.status != 200
        @warn(r.status)
        @warn(String(r.body))
        @warn "Couldn't convert account"
    else
        @info("Successfully upgraded account to Bot account")
    end
end

function setDefaultToken!(tok)
    global token, defaultHeader
    token = tok
    defaultHeader = Dict("Authorization" => "Bearer $token")
    myId =
        isempty(token) ? "" :
        JSON.parse(String(HTTP.get(baseUrl * "/account", defaultHeader).body))["id"]
    nothing
end

function countbits(n)
    c = 0
    while n != 0
        c += 1
        n &= n - 1
    end
    c
end

function isempty(p::Piece)
    p == EMPTY
end

function randomenginegame(d = 3; evalFunc = simpleEval)
    game = Game()
    while !isterminal(game)
        move = LichessBot.negamax(
            game,
            d,
            true,
            -typemax(Int),
            typemax(Int);
            evalFunc = evalFunc,
        )
        domove!(game, move)
        print("\r", ply(game))
    end
    println()
    game
end


text:
export parse_sequence

function parse_sequence(sequence::String)
    if length(sequence) == 1
        return parse_single_char_sequence(Char(sequence[1]))
    end

    (c, state) = iterate(sequence)
    if c == '\e'
        return parse_esc_leaded_sequence(sequence, state)
    end

    return PasteEvent(sequence) # fallback
end

function parse_single_char_sequence(sequence::Char, is_alt=false)
    if is_alt
        adjoint_key = ALT
        adjoint_ctrl = CTRL_ALT
    else
        adjoint_key = NO_ADJOINT
        adjoint_ctrl = CTRL
    end

    if sequence == '\n' || sequence == '\r' # ENTER
        return KeyPressedEvent(ENTER, adjoint_key)
    elseif sequence == '\e' # ESC
        return KeyPressedEvent(ESC, adjoint_key)
    elseif sequence == '\t' # TAB
        return KeyPressedEvent(TAB, adjoint_key)
    elseif sequence == '\x7F' #BACKSPACE
        return KeyPressedEvent(BACKSPACE, adjoint_key)
    elseif sequence == '\0' # NULL
        return KeyPressedEvent(NULL, adjoint_key)
    elseif any(isequal(sequence), Char(0x01):Char(0x1A)) # CTRL + key
        return KeyPressedEvent(Char(UInt8(sequence)-0x01+UInt8('a')), adjoint_ctrl)
    elseif any(isequal(sequence), Char(0x0C):Char(0x1F)) # CTRL + key
        return KeyPressedEvent(Char(UInt8(sequence)-0x1C+UInt8('4')), adjoint_ctrl)
    else # single char
        return KeyPressedEvent(sequence, adjoint_key)
    end
end

function parse_esc_leaded_sequence(sequence::String, state::Int)
    (c, state) = iterate(sequence, state)
    if c == 'O' # F1 - F4
        return parse_xterm_f1_2_f4(sequence, state)
    elseif c == '[' # Some CSI sequence
        return parse_csi(sequence, state)
    elseif iterate(sequence, state) === nothing # ALT + key
        return parse_single_char_sequence(c, true)
    end

    return PasteEvent(sequence) # fallback
end

function parse_xterm_f1_2_f4(sequence::String, state::Int)
    next = iterate(sequence, state)
    (next === nothing) && (return PasteEvent(sequence)) # fallback

    (c, state) = next
    if any(isequal(c), 'P':'S')
        return KeyPressedEvent(SpecialKey(UInt8(c)-UInt8('P')))
    end

    return PasteEvent(sequence) # fallback
end

function parse_csi(sequence::String, state::Int)
    next = iterate(sequence, state)
    (next === nothing) && (return PasteEvent(sequence)) # fallback

    if sequence[end] == '~' # vt sequence
        code_sequence = sequence[state:(end-1)]
        if ';' in code_sequence # with adjoint keys
            # +------------------------------------------------------------------+
            # | the form of the sequence: "\e[<code>;<adjoint_code>~" and code=1:35 |
            # +------------------------------------------------------------------+
            code, adjoint_code = tryparse.(Int, split(code_sequence, ';'))
            (code === nothing || adjoint_code === nothing) && (return PasteEvent(sequence))  # fallback
            adjoint = parse_adjoint_key_code(adjoint_code)
            (adjoint == -1) && (return PasteEvent(sequence))  # fallback
        else # without adjoint keys
            # +------------------------------------------------------+
            # | the form of the sequence: "\e[<code>~" and code=1:35 |
            # +------------------------------------------------------+
            code = tryparse(Int, code_sequence)
            (code === nothing) && (return PasteEvent(sequence))  # fallback
            adjoint = NO_ADJOINT
        end

        return parse_vt_code(sequence, code, adjoint)
    else # xterm sequence
        code_sequence = sequence[state:end]
        if ';' in code_sequence # with adjoint keys
            # +-----------------------------------------------------------------+
            # | the form of the sequence: "\e[1;<adjoint_code><code>" and code=A:Z |
            # +-----------------------------------------------------------------+
            id, code_sequence = split(code_sequence, ';')
            (id != "1" || length(code_sequence) != 2) && (return PasteEvent(sequence))  # fallback

            # determing key code
            c = code_sequence[2]
            if any(isequal(c), 'P':'S') # F1 - F4 with adjoint keys have the same form
                code = Int(c)-Int('P') # code=1:3, matches the enum F1 - F4
            else
                code = Int(c)
            end

            # construct AdjointKey array
            adjoint_code = tryparse(Int, string(code_sequence[1]))
            (adjoint_code === nothing) && (return PasteEvent(sequence)) # fallback
            adjoint = parse_adjoint_key_code(adjoint_code)
            (adjoint == -1) && (return PasteEvent(sequence))  # fallback
        else # without adjoint keys (not including F1 - F4)
            # +----------------------------------------------------+
            # | the form of the sequence: "\e[<code>" and code=A:Z |
            # +----------------------------------------------------+
            (length(code_sequence) != 1) && (return PasteEvent(sequence))  # fallback

            # determing key code
            code = Int(code_sequence[1])

            # construct AdjointKey
            adjoint = NO_ADJOINT
        end

        return parse_xterm_code(sequence, code, adjoint)
    end
end

function parse_vt_code(sequence::String, code::Int, adjoint=NO_ADJOINT)
    if code in 1:6 # HOME INSERT DELETE END PAGEUP PAGEDOWN
        enum_bias = 11
        return KeyPressedEvent(SpecialKey(code+enum_bias), adjoint)
    elseif code == 7 # HOME
        return KeyPressedEvent(HOME, adjoint)
    elseif code == 8 # END
        return KeyPressedEvent(END, adjoint)
    elseif code in 11:15 # F1 - F5
        return KeyPressedEvent(SpecialKey(code-11), adjoint)
    elseif code in 17:21 # F6 - f10
        return KeyPressedEvent(SpecialKey(code-12), adjoint)
    elseif code in 23:24 # F11 - F12
        return KeyPressedEvent(SpecialKey(code-13), adjoint)
    end

    return PasteEvent(sequence) # fallback
end

function parse_xterm_code(sequence::String, code::Int, adjoint=NO_ADJOINT)
    if code in Int('A'):Int('D') || code == Int('Z') # UP DOWN RIGHT LEFT BACKTAB
        return KeyPressedEvent(SpecialKey(code), adjoint)
    elseif code == Int('F') # END
        return KeyPressedEvent(END, adjoint)
    elseif code == Int('G') # Keypad 5
        return KeyPressedEvent('5', adjoint)
    elseif code == Int('H') # HOME
        return KeyPressedEvent(HOME, adjoint)
    elseif code in 0:3 # F1 - F4 with adjoint keys have the same form
        return KeyPressedEvent(SpecialKey(code), adjoint)
    end

    return PasteEvent(sequence) # fallback
end

function parse_adjoint_key_code(code::Int)
    if  code in 2:8
        return AdjointKey(code)
    end

    return -1 # fallback
end


text:
#=
DataReader.jl
This script implements basic IO functionalities reading data from the files described at (https://www.nature.com/articles/sdata201889) into an in-memory data structure. Possible improvements include rethinking data structure use-cases, reading directly from GTFS files, investigating possible gains of using the MetaGraphs Julia package, and supporting multiple edges between nodes.
=#

export ReadNatureFeed, ReadGTFS

#importing all necessary packages
using DataFrames
using LightGraphs
using CSV
using MetaGraphs

#function gets as parameter a directory path containing the necessary files described in Nature's article
"""
	ReadNatureFeed(directory_path)::ProblemInstance
Read data from a feed as specified in [this article](https://www.nature.com/articles/sdata201889)
"""
function ReadNatureFeed(directory_path)::ProblemInstance

	#checks if all needed files exist and loads them into dataframes:

	nodes_path = string(directory_path, "network_nodes.csv")
	edges_path = string(directory_path, "network_bus.csv") #loading the bus dataset for now
	temp_day_path = string(directory_path, "network_temporal_day.csv")
	stats_path = string(directory_path, "stats.csv")

	#why load all files into memory simultaneously?
	#reads network nodes into dataframe, if file exists. Else, raises error
	if isfile(nodes_path)
		nodes = CSV.read(nodes_path)
	else
		error(string("File ", nodes_path, " not found"))
	end

	#reads network connections into dataframe, if file exists. Else, raises error
	if isfile(edges_path)
		edges_df = CSV.read(edges_path)
	else
		error(string("File ", edges_path, " not found"))
	end

	# "", network_temporal_day
	#we extract, from the network_temporal_day file, infos about each route
	if isfile(temp_day_path)
		temp_day = CSV.read(temp_day_path)
	else
		error(string("File ", temp_day_path, " not found"))
	end

	if isfile(stats_path)
		stats_df = CSV.read(stats_path)
	else
		error(string("File ", stats_path, " not found"))
	end

	#existem arestas paralelas no grafo?
	edges_set = Set{Tuple{Int64, Int64}}()

	for row in eachrow(edges_df)
		edge = (row.from_stop_I, row.to_stop_I)
		if in(edge, edges_set)
			error("Error: parallel edges were detected in the network graph, but are not currently supported")
		else
			push!(edges_set, edge)
		end
	end

	#a identificação dos vértices é feita por inteiros, mas estes vem do arquivo e não necessariamente tem o índice < n_nodes

	#por isso, além de criar o grafo, crio um dicionário auxiliar que mapeia cada índice original ao seu novo índice
	#posso, posteriormente, criar outro dicionário para o processo reverso, caso seja necessário

	n_nodes = size(nodes, 1)
	graph = SimpleDiGraph(n_nodes)
	metagraph = MetaDiGraph(graph)
	new_indices = Dict{Int64, Int64}()

	i = 1
	for row in eachrow(nodes)
		new_indices[row.stop_I] = i

		#set properties for storing a stop's latitude and longitude
		set_props!(metagraph, i, Dict(:latitude => row.lat, :longitude => row.lon, :name => row.name) )

		i = i + 1
	end

	#important note: lightgraphs does not allow for multiple edges between two nodes

	#creates a problem instance given the above information

	n_edges = size(edges_df, 1)

	for row in eachrow(edges_df)
		#metagraph
		s,d = new_indices[row.from_stop_I], new_indices[row.to_stop_I] #source and destination nodes
		edge = Edge(s,d)
		add_edge!(metagraph, edge)
		set_prop!(metagraph, Edge(s,d), :weight, row.d * 1.0 + row.duration_avg * 0.0) #the weight of the edge is given as a function of the distance and the avg time

		#route_I_counts field: list (string)
		#as in https://www.nature.com/articles/sdata201889/tables/7
		#=
			A list of route_I's and the number of times each route has operated between two stops.
		For the network extracts, this data is formatted as a string where each element is written as ``route_I:count'' and different routes are separated by a comma.
		An example value for this field is thus ``1:3,2:131,10:93''. (...)Please note that the definition of a route varies across the cities provided,
		and that routes can have deviations from their main paths for instance when traveling to and from a depot.

			I use this field to estimate the minimum frequency requirements for the given edge. I heuristically assume the minimum frequency keep the same number of vehicles described by the network.
				This should guarantee that any optimization done at least mantains a similar level of service, for each edge
		=#
		str = row.route_I_counts
		sum = 0 # the sum of all vehicles covering this edge
		#route_coverage = Array{Int64,1}() #the routes id's covering this specific edge. ###### NOTE THAT, FOR NOW, THERE IS NO GUARANTEE THAT THE ID IS UNIQUE PER LINE (TO DO)########
		for m in eachmatch(r"(\d+)\:(\d+)", str)
			#push!(route_coverage, parse(Int64, m.captures[1]))
			sum += parse(Int64, m.captures[2])
		end

		set_prop!(metagraph, Edge(s,d), :base_freq, sum)

		#I'm removing these properties so I can place these kind of approximations in the ConstructModel function itself
		set_prop!(metagraph, Edge(s,d), :f_min, Int64(ceil(sum*0.5))) #setting it to floor may make requirements drop to zero on some edges
		set_prop!(metagraph, Edge(s,d), :f_max, Int64(sum*10))

		#line_coverage is a dict containing each line covering this line, mapping to is frequency
		#however, it isn't populated now because we haven't constructed the lines array yet
		set_prop!(metagraph, Edge(s,d), :line_coverage, Dict{Int64, Int64}())

	end

	#constructing the routes array

	#temporary, i just want to check if tripIds are unique
	all_trips = Set{Int64}()

	#this temporary dictionary records walks and their respective (id, frequency, cost, trip_id)
	#trip_id is later discarded, but it is useful for debugging
	walks_read = Dict{Array{Tuple{Int64, Int64}}, Tuple{Int64, Int64, Float64, Int64}}()

	routes_read = Set{Int64}()
	count = 0
	invalid_route = false
	i, row_count = 1, size(temp_day, 1)
	while i <= row_count # -1
		col = temp_day[i, :]

		#read the route in this column:
		routeId = col.route_I
		tripId = col.trip_I

		if in(tripId, all_trips)
			error("trip ids are not unique!" * string(tripId))
		else
			push!(all_trips, tripId)
		end

		#read the entire route,
		walk = Array{Tuple{Int64, Int64}, 1}()
		walk_cost = 0.0
		invalid_route = false
		while i <= row_count && tripId == temp_day[i, :trip_I] #we've not reached the end of the data-frame and we're still reading the same trip
			col = temp_day[i, :]

			#the edge is a tuple indicating its origin and destination nodes,
			#and we have to take the new node numeration as defined above
			s,d = new_indices[col.from_stop_I], new_indices[col.to_stop_I]

			#if a route isn't of the bus mode,
			#or if an edge in the route has not been added to the network, it means this route is not in the same mode of transport we've loaded
			#thus we should ignore it for now
			if col.route_type != 3 || ~has_edge(metagraph, s, d) #3 is the identifier for the 'bus' mode in a GTFS feed(https://developers.google.com/transit/gtfs/reference/)
				invalid_route = true
				i += 1
				continue # continue is used instead of break to make the iteration run to the end of this unwanted trip
			end

			edge = (s,d)
			push!(walk, edge)

			#increment walk's total cost. The cost per edge is given by a funcion of the edge's distance and time
			walk_cost += get_prop(metagraph, s, d, :weight)
			i += 1
		end

		if ~invalid_route
			#counts how many unique route identifiers there are
			count += 1
			push!(routes_read, routeId)

			#if the route I just read had not been read yet,
			if ~haskey(walks_read, walk)
				walks_read[walk] = (routeId, 1, walk_cost, tripId)
			else
				#increment the walk's frequency
				walks_read[walk] = (walks_read[walk][1], walks_read[walk][2] + 1, walks_read[walk][3], walks_read[walk][4])
			end
		end
	end

	#FILTERING WALKS:
	#it is common that some routes feature little variations during the day, making them register as an unique walk. We discard these instances
	#freq_threshold = 1 #possibly parametrize this later
	#walks_read = filter((pair) -> pair.second[2] > freq_threshold, walks_read) #filter to remove frequency values < freq_threshold

	#=
	#it is also common for different walks to be part of the same line(as in: one way and back). We merge these into one bigger walk
	walks_dict_perID = Dict{Int64, Tuple{Array{Tuple{Int64, Int64}, 1}, Int64, Float64}}() #dictionary of (walks, freq, cost), indexed per id
	for (k,v) in walks_read
		id = v[1]
		if haskey(walks_dict_perID, id) #if we've already seen this ID, but its a different walk,
			#try to merge the two walks into one
			println("Merging two walks. Id: $(id)")
			freq = v[2]
			if walks_dict_perID[id][2] != v[2]
				@warn "Merging two walks with different frequencies. Lowest frequency is kept"
				freq = min(walks_dict_perID[id][2], v[2])
			end
			walk = vcat(walks_dict_perID[id][1], k)
			cost = v[3] + walks_dict_perID[id][3]
			walks_dict_perID[id] = (walk, freq, cost)
		else
			walks_dict_perID[id] = (k, v[2], v[3])
		end
	end
	=#

	#possibly divide the frequency for a given time period (24h) and/or filter lines with too low a frequency

	#take every read walk, its frequency and create the lines array
	lines = Array{Line, 1}()
	for (k,v) in walks_read
		line = Line(v[1], k, v[2], v[3])
		push!(lines, line)
	end

	#println("unique route ids:" ,  length(routes_read))
	#println("Unique lines:", length(lines))
	#println("Count:", count)

	#iterate through the lines array, adding some redundant info on the graphs edges that refer back to each line
	# O(e*v), isn't great, but will contribute to avoiding constant O(e*v) checks in the future
	for i in 1:length(lines)
		line = lines[i]
		for (s, d) in line.walk
			dict = get_prop(metagraph, Edge(s, d), :line_coverage)
			dict[i] = line.frequency
			set_prop!(metagraph, Edge(s,d), :line_coverage, dict)
		end
	end

	#add some useful (for plotting) info to the metagraph
	set_prop!(metagraph, :buffer_center_lat, stats_df[1, :buffer_center_lat])
	set_prop!(metagraph, :buffer_center_lon, stats_df[1, :buffer_center_lon])
	set_prop!(metagraph, :buffer_radius_km, stats_df[1, :buffer_radius_km])

	#finally, construct new problem instance

	problemInstance = ProblemInstance(metagraph, lines)
	return problemInstance

#ends ReadProblemInstance function
end

"""
	ReadGTFS(directory_path, transportModeCode=3)::ProblemInstance

Read the data in an uncompressed GTFS feed into a ProblemInstance struct.

The transportModeCode parameter defines which mode of transportation should
be considered, with the default value = 3 indicating buses. Another common value would be 1 for trains.

"""
function ReadGTFS(directory_path::String; transportModeCode::Int64 = 3)::ProblemInstance
	#checks if all needed files exist and loads them into dataframes:

	nodes_path = string(directory_path, "stops.txt")
	routes_path = string(directory_path, "routes.txt")
	trips_path = string(directory_path, "trips.txt")
	stop_times_path = string(directory_path, "stop_times.txt")

	#why load all files into memory simultaneously?

	#reads network nodes into dataframe, if file exists. Else, raises error
	if !isfile(nodes_path)
		error("File " * nodes_path * " not found")
	end
	nodes = CSV.read(nodes_path)

	if !isfile(routes_path)
		error("File " * routes_path * " not found")
	end
	routes = CSV.read(routes_path)
	#a mapping of route ids to route_type
    route_types = Dict{Any, Int}()
    for row in eachrow(routes)
        if(row.route_id !== DataFrames.missing && row.route_type !== DataFrames.missing)
            route_types[row.route_id] = row.route_type
        end
    end

	if !isfile(trips_path)
		error("File " * trips_path * " not found")
	end
	trips_df = CSV.read(trips_path)
	#a mapping of trip_ids to route_ids
    trips_to_routes = Dict{Any, Any}()
    for row in eachrow(trips_df)
        if(row.trip_id !== DataFrames.missing && row.route_id !== DataFrames.missing)
            trips_to_routes[row.trip_id] = row.route_id
        end
    end

	n_nodes = size(nodes, 1)
	graph = SimpleDiGraph(n_nodes)
	metagraph = MetaDiGraph(graph)

    #a identificação dos vértices é feita por inteiros, mas estes vem do arquivo e não necessariamente tem o índice < n_nodes
	#por isso, além de criar o grafo, crio um dicionário auxiliar que mapeia cada índice original ao seu novo índice
	new_indices = Dict{Any, Int64}()
	i = 1
    #at the same time, we compute median latitude and longitude, it may later be used for plotting purposes
    lat_sum, lon_sum = 0.0, 0.0
	for row in eachrow(nodes)
		new_indices[row.stop_id] = i
		#set properties in metagraph for storing a stop's latitude and longitude
		set_props!(metagraph, i, Dict(:latitude => row.stop_lat, :longitude => row.stop_lon, :name => row.stop_name) )
		lat_sum += row.stop_lat
        lon_sum += row.stop_lon
        i = i + 1
	end
    set_prop!(metagraph, :buffer_center_lat, lat_sum / (i - 1))
    set_prop!(metagraph, :buffer_center_lon, lon_sum / (i - 1))

	#reads network connections into dataframe, if file exists. Else, raises error
	if !isfile(stop_times_path)
		error(string("File ", stop_times_path, " not found"))
	end
	stop_times_df = CSV.read(stop_times_path)

	#constructing the routes array AND the graph's edges

	#temporary, i just want to check if tripIds are unique
	all_trips = Set{Any}()

	#this temporary dictionary records walks and their respective (id, frequency, cost, trip_id)
	#trip_id is later discarded, but it is useful for debugging
	walks_read = Dict{Array{Tuple{Int64, Int64}, 1}, Tuple{String, Int64, Float64}}()

	routes_read = Set{Int64}()
	count = 0
	invalid_route = false
	i, row_count = 1, size(stop_times_df, 1)
	while i <= row_count
        col = stop_times_df[i, :]

        #read the route in this column:
        tripId = col.trip_id
        routeId = trips_to_routes[tripId]

        if in(tripId, all_trips)
            error("trip ids are not unique!" * string(tripId))
        else
            push!(all_trips, tripId)
        end

        #read the entire route,
        invalid_route = route_types[trips_to_routes[stop_times_df[i, :trip_id]]] != transportModeCode #compares with the given gtfs standard transport mode code

        stops = Array{Int64,1}() #the sequence of all stops in this trip
        while i <= row_count && tripId == stop_times_df[i, :trip_id] #we've not reached the end of the data-frame and we're still reading the same trip
            col = stop_times_df[i, :]
            #if the route type isnt of type bus, just skip it for now. But the while continues, as we just need to find its end
            if invalid_route
                i += 1
                continue
            end
            stop_id = new_indices[col.stop_id]
            push!(stops, stop_id)
            i += 1
        end

        #iterate through all of the trip's stops, constructing the edges,
        #only if there's more than one stop
        walk = Array{Tuple{Int64, Int64}, 1}()
        walk_cost = 0.0
        for i in 1:(length(stops) - 1)
            from = stops[i]
            to = stops[i+1]
            #creates edge
            edge = (from,to)
            push!(walk, edge)

            #increments walk's cost with calculated point to point distance
            lat1, lon1 = get_prop(metagraph, from, :latitude), get_prop(metagraph, from, :longitude)
            lat2, lon2 = get_prop(metagraph, to, :latitude), get_prop(metagraph, to, :longitude)
            d = GreatCircleDistance(lat1, lon1, lat2, lon2)
            walk_cost += d

            if !has_edge(metagraph, from,to)
                #adds edge to graph, also setting its weigth:
                edge = Edge(from, to)
                add_edge!(metagraph, edge)
                set_prop!(metagraph, edge, :weight, d)
            end

        end

        #if the route I just read had not been read yet,
        if ~haskey(walks_read, walk)
			if length(walk) > 0
				walks_read[walk] = (string(trips_to_routes[tripId]), 1, walk_cost)
			else
				@warn "Read walk with zero length?"
			end
        else
            #increment the walk's frequency
            walks_read[walk] = (walks_read[walk][1], walks_read[walk][2] + 1, walks_read[walk][3])
        end

    end
    #to-do - filtering walks: merging forw and back walks into one, frequency threshold, etc.

	#FILTERING WALKS:
	#it is common that some routes feature little variations during the day, making them register as an unique walk. We discard these instances
	#freq_threshold = 1 #possibly parametrize this later
	#walks_read = filter((pair) -> pair.second[2] > freq_threshold, walks_read) #filter to remove frequency values < freq_threshold

	#=
	#it is also common for different walks to be part of the same line(as in: one way and back). We merge these into one bigger walk
	walks_dict_perID = Dict{Int64, Tuple{Array{Tuple{Int64, Int64}, 1}, Int64, Float64}}() #dictionary of (walks, freq, cost), indexed per id
	for (k,v) in walks_read
		id = v[1]
		if haskey(walks_dict_perID, id) #if we've already seen this ID, but its a different walk,
			#try to merge the two walks into one
			println("Merging two walks. Id: $(id)")
			freq = v[2]
			if walks_dict_perID[id][2] != v[2]
				@warn "Merging two walks with different frequencies. Lowest frequency is kept"
				freq = min(walks_dict_perID[id][2], v[2])
			end
			walk = vcat(walks_dict_perID[id][1], k)
			cost = v[3] + walks_dict_perID[id][3]
			walks_dict_perID[id] = (walk, freq, cost)
		else
			walks_dict_perID[id] = (k, v[2], v[3])
		end
	end
	=#

    #take every read walk, its frequency and create the lines array
    lines = Array{Line, 1}()
	edges_freq = Dict{Tuple{Int64, Int64}, Int64}() #set the edge's min_freq and max_freq based on the lines array

	#line_coverage is a dict containing each line covering this line, mapping to is frequency,
	#and we should have one of them per edge
	edges_lineCoverage = Dict{Tuple{Int64, Int64}, Dict{Int64, Int64}}()

	i = 1
	for (k,v) in walks_read
        line = Line(v[1], k, v[2], v[3])
        push!(lines, line)

        for (s,d) in k
			if haskey(edges_freq, (s,d))
				edges_freq[(s,d)] = edges_freq[(s,d)] + v[2]

				edges_lineCoverage[(s,d)][i] = v[2]
            else
				edges_freq[(s,d)] = v[2]
				#println(typeof(edges_lineCoverage))
				edges_lineCoverage[(s,d)] = Dict{Int64,Int64}(i => v[2])
            end
		end
		i += 1

	end

    #set min_freq and max_freq for each edge
    for e in edges(metagraph)
        s = src(e)
        d = dst(e)
        if ~haskey(edges_freq, (s,d))
            error("An edge was included in the graph but its accumulated frequency is non existent.")
        else
            base_freq = edges_freq[(s,d)]
            set_prop!(metagraph, e, :base_freq, base_freq)
            set_prop!(metagraph, e, :f_min, Int64(max(1, floor(base_freq*0.5))))
			set_prop!(metagraph, e, :f_max, Int64(base_freq*10))

			set_prop!(metagraph, e, :line_coverage, edges_lineCoverage[(s,d)])

        end
    end

    #finally, construct new problem instance
    problemInstance = ProblemInstance(metagraph, lines)
    return problemInstance

#end read gtfs function
end


text:
module TraitSimulation
using GLM # this already defines some useful distribution and link types
using DataFrames # so we can test it
using StatsModels # useful distributions #lots more useful distributions
using LinearAlgebra
using Random

include("calculate_mean_vector.jl")

include("apply_inverse_link_new.jl")
export LogLink, IdentityLink, SqrtLink, ProbitLink, LogitLink, InverseLink, CauchitLink, CloglogLink

include("simulate_glm_trait_new.jl")
export PoissonResponse, NormalResponse, BinomialResponse, BernoulliResponse, GammaResponse, InverseGaussianResponse, TResponse, WeibullResponse #Exporting these from the Distributions package

#this is the main functionality of this package, to run the actual simulation now for the split up responsedist type and linkfunction type
function GLM_trait_simulation(mu, dist::ResponseDistribution, link::InverseLinkFunction)
  transmu = apply_inverse_link(mu, link)
  Simulated_Trait = simulate_glm_trait(transmu, dist)
  return(Simulated_Trait)
end

########

include("Multiple_traits.jl")

include("Model_Framework.jl")

"""
```
simulate(trait::GLMTrait)
```
this for GLM trait
"""
function simulate(trait::GLMTrait)
    simulated_trait = GLM_trait_simulation(trait.mu, trait.dist, trait.link)
    rep_simulation = DataFrame(trait1 = simulated_trait)
    return(rep_simulation)
end

function simulate(trait::GLMTrait, n_reps::Int64)
  rep_simulation = Vector{DataFrame}(undef, n_reps)
  for i in 1:n_reps
    rep_simulation[i] = simulate(trait) # store each data frame in the vector of dataframes rep_simulation
  end
    return(rep_simulation)
end

"""
```
simulate(trait::LMMTrait)
```
this for LMMtrait
"""
function simulate(trait::LMMTrait)
  rep_simulation = LMM_trait_simulation(trait.mu, trait.vc)
  return(rep_simulation)
end

function simulate(trait::LMMTrait, n_reps::Int64)
  rep_simulation = Vector{DataFrame}(undef, n_reps)
  for i in 1:n_reps
    rep_simulation[i] = simulate(trait)
  end
  return(rep_simulation)
end

export ResponseType, GLM_trait_simulation, mean_formula, VarianceComponent, append_terms!, LMM_trait_simulation
export GLMTrait, Multiple_GLMTraits, LMMTrait, simulate, @vc, vcobjtuple, SimulateMVN, SimulateMVN!, Aggregate_VarianceComponents!
export TResponse, WeibullResponse, PoissonResponse, NormalResponse, BernoulliResponse, BinomialResponse
export GammaResponse, InverseGaussianResponse, ExponentialResponse
export CauchitLink, CloglogLink, IdentityLink, InverseLink, LogitLink, LogLink, ProbitLink, SqrtLink
end #module

max_stars_repo_name: sarah-ji/Tutorials

max_issues_repo_name: sarah-ji/Tutorials

max_forks_repo_name: sarah-ji/Tutorials


text:
module FixedEffects

##############################################################################
##
## Dependencies
##
##############################################################################
using Base
using LinearAlgebra
using Distributed
using CategoricalArrays
using FillArrays
using Reexport
using CUDAapi
if has_cuda()
	try
		using CuArrays
		@eval has_cuarrays() = true
	catch
		@info "CUDA is installed, but CuArrays.jl fails to load"
		@eval has_cuarrays() = false
	end
else
	has_cuarrays() = false
end
@reexport using StatsBase # used for Weights
##############################################################################
##
## Exported methods and types
##
##############################################################################

export
group,
FixedEffect,
AbstractFixedEffectSolver,
solve_residuals!,
solve_coefficients!

##############################################################################
##
## Load files
##
##############################################################################
include("utils/lsmr.jl")
include("FixedEffect.jl")
include("FixedEffectSolvers/FixedEffectLinearMap.jl")
include("FixedEffectSolvers/AbstractFixedEffectSolver.jl")
include("FixedEffectSolvers/FixedEffectSolverLSMR.jl")
include("FixedEffectSolvers/FixedEffectSolverLSMRParallel.jl")
if has_cuarrays()
	include("FixedEffectSolvers/FixedEffectSolverLSMRGPU.jl")
end

end  # module FixedEffects


text:
struct Ranq1 {
Recommended generator for everyday use. The period is 	 1:8  10 19 . Calling conventions
same as Ran, above.
Ullong v;
Ranq1(Ullong j) : v(4101842887655102017LL) {
v ^= j;
v = int64();
}
inline Ullong int64() {
v ^= v >> 21; v ^= v << 35; v ^= v >> 4;
return v * 2685821657736338717LL;
}
inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
inline Uint int32() { return (Uint)int64(); }
};

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
"Return a vector of the indices defining the parameter or variable."
function getdims(model::Model, component::Symbol, name::Symbol)
    dims = getdimnames(model, component, name)
    Vector{Int64}(map(dim -> dim_count(model, dim), dims))
end

"""
Return the symbols representing each of the dimensions for this variable or parameter.
"""
function getdimnames(model::Model, component::Symbol, name::Symbol)
    if name in parameter_names(model, component)
        parameter_dimensions(model, component, name)
    else
        variable_dimensions(model, component, name)
    end
end

"""
Avoid the distinction between set_param! and update_param!
"""
function set_or_update_param!(model::Model, component::Symbol, name::Symbol, value)
    try
        update_param!(model, name, value)
    catch
        set_param!(model, component, name, value)
    end
end


text:
if !("." in LOAD_PATH) # for ease of local testing
    push!(LOAD_PATH, ".")
end

import Random
Random.seed!(1234)
using Test
import Thermodynamics
using SurfaceFluxes
const SF = SurfaceFluxes
const UF = SF.UniversalFunctions
using StaticArrays
import KernelAbstractions: CPU

import CLIMAParameters
const CP = CLIMAParameters
const CPP = CP.Planet
const APS = CP.AbstractEarthParameterSet
const CPSGS = CP.SubgridScale
const CPSGS = CP.SubgridScale
const TD = Thermodynamics

# TODO: add GPU tests back in, similar to Thermodynamics
# if get(ARGS, 1, "Array") == "CuArray"
#     using CUDA
#     import CUDAKernels: CUDADevice
#     ArrayType = CUDA.CuArray
#     CUDA.allowscalar(false)
#     device(::T) where {T <: CuArray} = CUDADevice()
# else
ArrayType = Array
device(::T) where {T <: Array} = CPU()
# end

@show ArrayType

using CLIMAParameters: AbstractEarthParameterSet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

@testset "SurfaceFluxes - Recovery Profiles" begin
    FT = Float32
    ρ_sfc = FT(1.15)
    ρ_in = FT(1.13)
    qt_sfc = FT(0.01)
    qt_in = FT(0.009)
    ## Discretisation altitude z
    z = ArrayType(FT[29.432779269303, 30.0497139076724, 31.6880000418153, 34.1873479240475])
    ## Virtual Potential Temperature at height z
    θ = ArrayType(FT[268.559120403867, 269.799228886728, 277.443023238556, 295.79192777341])
    ## Surface Pottemp
    θ_sfc = ArrayType(FT[273.42369841804, 272.551410044203, 278.638168565727, 298.133068766049])
    ## Roughness lengths
    z0 = ArrayType(FT[5.86144925739178e-05, 0.0001, 0.000641655193293549, 3.23383768877187e-05])
    ## Speed
    speed = ArrayType(FT[2.9693638452068, 2.43308757772094, 5.69418282305367, 9.5608693754561])
    ## Scale velocity and moisture
    u_star = ArrayType(FT[0.109462510724615, 0.0932942802513508, 0.223232887323184, 0.290918439028557])
    # No explicit buoyancy terms in ClimateMachine
    b_star = ArrayType([0.00690834676781433, 0.00428178089592372, 0.00121229800895103, 0.00262353784027441])

    κ = CPSGS.von_karman_const(param_set)
    for ii in 1:length(b_star)
        # Compute L_MO given u_star and b_star
        L_MO = u_star[ii]^2 / κ / b_star[ii]

        ts_sfc = TD.PhaseEquil_ρθq(param_set, ρ_sfc, θ_sfc[ii], qt_sfc)
        ts_in = TD.PhaseEquil_ρθq(param_set, ρ_in, θ[ii], qt_in)

        state_sfc = SF.SurfaceValues(FT(0), SVector{2, FT}(0, 0), ts_sfc)
        state_in = SF.InteriorValues(z[ii], SVector{2, FT}(speed[ii], 0), ts_in)

        # State containers
        kwargs = (; state_in, state_sfc, z0m = z0[ii], z0b = FT(0.001))
        sc = (
            SF.Fluxes{FT}(; kwargs..., shf = FT(0), lhf = FT(0)),
            SF.FluxesAndFrictionVelocity{FT}(; kwargs..., shf = FT(0), lhf = FT(0), ustar = u_star[ii]),
            SF.Coefficients{FT}(; kwargs..., Cd = FT(0.001), Ch = FT(0.001)),
            SF.ValuesOnly{FT}(; kwargs...),
        )
        uf = UF.Businger()
        for jj in 1:length(sc)
            u_scale_fd =
                SF.compute_physical_scale_coeff(param_set, sc[jj], L_MO, UF.MomentumTransport(), uf, SF.FDScheme())
            Δu_fd = u_star[ii] / u_scale_fd
            u_scale_fv =
                SF.compute_physical_scale_coeff(param_set, sc[jj], L_MO, UF.MomentumTransport(), uf, SF.FVScheme())
            Δu_fv = u_star[ii] / u_scale_fv
            @test (Δu_fd - Δu_fv) ./ Δu_fd * 100 <= FT(50)
        end
    end
end

include("test_profiles.jl")
include("test_universal_functions.jl")


text:
# Do monkey patches (via JuliaVenv) only in the "root" julia-venv
# process. `Base.compilecache` uses `Base.load_path_setup_code` to
# propagate so there is no need to setup `Base.DEPOT_PATH` manually
# inside child processes launched by `Base.compilecache`.
if get(ENV, "_JULIA_VENV_PROCESS", "") != "true"
    if !any(isequal(@__DIR__), Base.LOAD_PATH)
        pushfirst!(Base.LOAD_PATH, @__DIR__)
    end
    @eval using JuliaVenv
end
ENV["_JULIA_VENV_PROCESS"] = "true"

max_stars_repo_name: tkf/julia-venv

max_issues_repo_name: tkf/julia-venv

max_forks_repo_name: tkf/julia-venv


text:
# cat
A1data, A2data = [1 3; 2 4], [5 7; 6 8]

A1 = AxisArray(A1data, Axis{:Row}([:First, :Second]), Axis{:Col}([:A, :B]))
A2 = AxisArray(A2data, Axis{:Row}([:Third, :Fourth]), Axis{:Col}([:A, :B]))
@test isa(cat(1, A1, A2), AxisArray)
@test cat(1, A1, A2) == AxisArray(vcat(A1data, A2data),
                                  Axis{:Row}([:First, :Second, :Third, :Fourth]), Axis{:Col}([:A, :B]))

A2 = AxisArray(A2data, Axis{:Row}([:First, :Second]), Axis{:Col}([:C, :D]))
@test isa(cat(2, A1, A2), AxisArray)
@test cat(2, A1, A2) == AxisArray(hcat(A1data, A2data),
                                  Axis{:Row}([:First, :Second]), Axis{:Col}([:A, :B, :C, :D]))

A2 = AxisArray(A2data, Axis{:Row}([:First, :Second]), Axis{:Col}([:A, :B]))
@test isa(cat(3, A1, A2), AxisArray)
@test cat(3, A1, A2) == AxisArray(cat(3, A1data, A2data),
                                       Axis{:Row}([:First, :Second]), Axis{:Col}([:A, :B]),
                                       Axis{:page}(1:2))

A1 = AxisArray(A1data, :Row, :Col)
A2 = AxisArray(A2data, :Row, :Col)
@test_throws ArgumentError cat(2, A1, A2)
@test cat(3, A1, A2) == AxisArray(cat(3, A1data, A2data), :Row, :Col)

# merge
Adata, Bdata, Cdata = randn(4,4,2), randn(4,4,2), randn(4,4,2)
A = AxisArray(Adata, Axis{:X}([1,2,3,4]), Axis{:Y}([10.,20,30,40]), Axis{:Z}([:First, :Second]))
B = AxisArray(Bdata, Axis{:X}([3,4,5,6]), Axis{:Y}([30.,40,50,60]), Axis{:Z}([:First, :Second]))

ABdata = zeros(6,6,2)
ABdata[1:4,1:4,:] = Adata
ABdata[3:6,3:6,:] = Bdata
@test merge(A,B) == AxisArray(ABdata, Axis{:X}([1,2,3,4,5,6]), Axis{:Y}([10.,20,30,40,50,60]), Axis{:Z}([:First, :Second]))

AC = AxisArray(cat(3, Adata, Cdata), :X, :Y, :Z)
B2 = AxisArray(Bdata, :X, :Y, :Z)
@test merge(AC,B2) == AxisArray(cat(3, Bdata, Cdata), :X, :Y, :Z)

# join
ABdata = zeros(6,6,2,2)
ABdata[1:4,1:4,:,1] = Adata
ABdata[3:6,3:6,:,2] = Bdata
@test join(A,B) == AxisArray(ABdata, Axis{:X}([1,2,3,4,5,6]), Axis{:Y}([10.,20,30,40,50,60]), Axis{:Z}([:First, :Second]))
@test join(A,B, newaxis=Axis{:JoinAxis}([:A, :B])) == AxisArray(ABdata, Axis{:X}([1,2,3,4,5,6]), Axis{:Y}([10.,20,30,40,50,60]), Axis{:Z}([:First, :Second]), Axis{:JoinAxis}([:A, :B]))
@test join(A,B,method=:inner) == AxisArray(ABdata[3:4, 3:4, :, :], Axis{:X}([3,4]), Axis{:Y}([30.,40]), Axis{:Z}([:First, :Second]))
@test join(A,B,method=:left) == AxisArray(ABdata[1:4, 1:4, :, :], A.axes...)
@test join(A,B,method=:right) == AxisArray(ABdata[3:6, 3:6, :, :], B.axes...)
@test join(A,B,method=:outer) == join(A,B)