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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ b415cf9a-fdb2-4070-a040-24ba4e7f5992
using Luxor, PlutoUI, Lindenmayer

# ╔═╡ a64aaa10-668c-11ec-3bc8-d9a3c11680f1
md"# First steps

with `Luxor.jl`, grids, some randomisation.

"

# ╔═╡ 16b2c558-8c14-4280-b609-06e8a390cc90
begin
	w = 400
	h = 400
	theOrigin = Point(0,0)
end

# ╔═╡ 0cc28bc8-fd18-403c-b1c5-e4c41d9f3c06
#@bind go Button("New colour")
@bind go Clock(0.5, true)

# ╔═╡ 8f08651a-2f4f-4ff7-adb2-bd6606b272f6
begin
	go # makereactive
	hue = rand(3);
end

# ╔═╡ 06e04f59-09e0-43de-b504-42b7cc6b752e
begin
	Drawing(w,h, "my-drawing.svg")
	origin()
	sethue(hue...)
	circle(theOrigin, 20, :fill)
	finish()
	preview()
end

# ╔═╡ 7c0f281a-cd4a-4b4d-9e78-5ecea65306c1
@bind shift Slider(-2:0.2:+3; show_value=true)

# ╔═╡ 2ffcabe8-c644-4383-8411-db204a8a3995
grid = Luxor.GridHex(theOrigin, 25, w, h)

# ╔═╡ 1a3fb2d1-d152-45be-b076-648ae90288b8
begin
	Drawing(w, h, "my-drawing.svg")
	#origin()
	for n in 1:200
		setcolor(hue...,0.5)

		#sethue(hue...)
		p = nextgridpoint(grid)
		s = rand([2,5,10,20], 1)
		jitter = shift.*rand(2)
		translate(jitter...)
		circle(p, s..., :fill)
	end
		finish()
	preview()
end

# ╔═╡ 7ae59e0f-0418-4ff6-9b41-aeba77a7dafc
md"## Drawing snowflakes

Simple version, getting 6-fold symmetry first, then elaborations
"

# ╔═╡ d0173b90-e89b-4484-8ff9-760354d02e19
@bind snow Button("Draw a snowflake")

# ╔═╡ a225706b-d889-4b9e-b6cf-2e248c74ee43
function drawArms(x0, x1, y; sz=10, color="white", lw=1.0)

	x0; # a point

	for iy in y
		isz = sz .* rand([0.2, 0.5, 1],1)[1]
		println(iy)
		sethue("white")
    	zb = between(x0, x1, iy)
		#circle(zb, 2, :fill)
		translate(zb)
		sethue("lightblue")
		rotate(π/3)
		line(x0, Point.(x0.x, isz), :stroke)
		rotate(-2π/3)
		line(x0, Point.(x0.x, isz), :stroke)
		rotate(π/3)
	end
	translate(x0)
end

# ╔═╡ 00588385-8862-46dd-87c4-fc2d2b4b540c
function drawSnowFlakelet(z0, z1; locs = [0.1, 0.3], color="white", lw=1.0)

	# draw one arm of snowflake:
	sethue(color)
	setline(lw)
	line(z0, z1, :stroke)
	sethue("red")
	rule(z0, 0) # a rule through z0
	drawArms(z0, z1, locs; sz=distance(z0, z1)./3,lw=0.5)
end

# ╔═╡ bc4f093a-9345-400f-93c5-33c38513e472
begin
	snow
	Drawing(w,h, "snow-flakes.svg")
	origin()
	background("black")
	sethue("white")
	sz = 120;
	n = 6;  #rand(2:10,1)[1]
	theLocs = [rand(n)...]
	for a in π/6:π/3:2π
		origin()
		rotate(a)
		drawSnowFlakelet(O, Point.(0, sz); locs = theLocs)
	end
	finish()
	preview()
end

# ╔═╡ abb44521-d176-4c80-a588-f90a4d75f8de
some_system  = LSystem(Dict(
    "F" => "+F+F+t",
    "G" => "F-G"),
    "F")

# ╔═╡ b02d7f97-b0ea-4430-8802-1bda35d550e8
@bind it Slider(1:10; show_value=true)

# ╔═╡ 56dad40c-0814-4d82-b9e6-b831755a7337
drawLSystem(some_system,
    forward     = 10,
    turn        = 60,
    iterations  = it,
    startingx   = 0,
    startingy   = 0,
	width=200,
	height=200,
    filename    = "some_system.svg",
    showpreview = true)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Lindenmayer = "d9db2020-d13e-4bf8-b13e-b12a3214b5ef"
Luxor = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Lindenmayer = "~0.2.0"
Luxor = "~2.18.0"
PlutoUI = "~0.7.27"
"""

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

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

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

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

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

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Cairo]]
deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"]
git-tree-sha1 = "d0b3f8b4ad16cb0a2988c6788646a5e6a17b6b1b"
uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
version = "1.0.5"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

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

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

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

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

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

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "2db648b6712831ecb333eae76dbfd1c156ca13bb"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.2"

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

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+2"

[[Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.1"

[[Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

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

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

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

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

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

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[Juno]]
deps = ["Base64", "Logging", "Media", "Profile"]
git-tree-sha1 = "07cb43290a840908a771552911a6274bc6c072c7"
uuid = "e5e0dc1b-0480-54bc-9374-aad01c23163d"
version = "0.8.4"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

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

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Librsvg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pango_jll", "Pkg", "gdk_pixbuf_jll"]
git-tree-sha1 = "25d5e6b4eb3558613ace1c67d6a871420bfca527"
uuid = "925c91fb-5dd6-59dd-8e8c-345e74382d89"
version = "2.52.4+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[Lindenmayer]]
deps = ["Colors", "Luxor"]
git-tree-sha1 = "5b64275624a79e78185fe0ab0ac6d3da646abba9"
uuid = "d9db2020-d13e-4bf8-b13e-b12a3214b5ef"
version = "0.2.0"

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

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

[[Luxor]]
deps = ["Base64", "Cairo", "Colors", "Dates", "FFMPEG", "FileIO", "Juno", "Random", "Rsvg"]
git-tree-sha1 = "6d663c668fd508e55889195f51bb13a442e951a7"
uuid = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc"
version = "2.18.0"

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

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

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

[[Media]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "75a54abd10709c01f1b86b84ec225d26e840ed58"
uuid = "e89f7d12-3494-54d1-8411-f7d8b9ae1f27"
version = "0.5.0"

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

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

[[NaNMath]]
git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.6"

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

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Pango_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9bc1871464b12ed19297fbc56c4fb4ba84988b0d"
uuid = "36c8627f-9965-5494-a995-c6b170f724f3"
version = "1.47.0+0"

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

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

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

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

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

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

[[Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.2.0"

[[Rsvg]]
deps = ["Cairo", "Glib_jll", "Librsvg_jll"]
git-tree-sha1 = "3d3dc66eb46568fb3a5259034bfc752a0eb0c686"
uuid = "c4c386cf-5103-5370-be45-f3a111cca3b8"
version = "1.0.0"

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

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

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

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

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

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

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

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

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[gdk_pixbuf_jll]]
deps = ["Artifacts", "Glib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg", "Xorg_libX11_jll", "libpng_jll"]
git-tree-sha1 = "c23323cd30d60941f8c68419a70905d9bdd92808"
uuid = "da03df04-f53b-5353-a52f-6a8b0620ced0"
version = "2.42.6+1"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

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

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

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"
"""

# ╔═╡ Cell order:
# ╟─a64aaa10-668c-11ec-3bc8-d9a3c11680f1
# ╟─b415cf9a-fdb2-4070-a040-24ba4e7f5992
# ╟─16b2c558-8c14-4280-b609-06e8a390cc90
# ╠═06e04f59-09e0-43de-b504-42b7cc6b752e
# ╟─0cc28bc8-fd18-403c-b1c5-e4c41d9f3c06
# ╟─8f08651a-2f4f-4ff7-adb2-bd6606b272f6
# ╟─7c0f281a-cd4a-4b4d-9e78-5ecea65306c1
# ╟─2ffcabe8-c644-4383-8411-db204a8a3995
# ╟─1a3fb2d1-d152-45be-b076-648ae90288b8
# ╟─7ae59e0f-0418-4ff6-9b41-aeba77a7dafc
# ╠═bc4f093a-9345-400f-93c5-33c38513e472
# ╟─d0173b90-e89b-4484-8ff9-760354d02e19
# ╟─00588385-8862-46dd-87c4-fc2d2b4b540c
# ╟─a225706b-d889-4b9e-b6cf-2e248c74ee43
# ╟─abb44521-d176-4c80-a588-f90a4d75f8de
# ╠═b02d7f97-b0ea-4430-8802-1bda35d550e8
# ╟─56dad40c-0814-4d82-b9e6-b831755a7337
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: schluppeck/generative-art

max_issues_repo_name: schluppeck/generative-art

max_forks_repo_name: schluppeck/generative-art


text:
"""
Represents a trajector matrix without using any more
memory than the source data.

https://en.wikipedia.org/wiki/Singular_spectrum_analysis#SSA_as_a_model-free_tool

Techinally this will work for for source data of arbitrary dimensions, however,
under most use cases it'll either be a 1d array or a 2d array
(ie: each time step as a set of features)
"""
type TrajectoryMatrix{T,N} <: AbstractArray{T,N}
    data::AbstractArray{T}
    window::Int
    nsteps::Int
    step_size::Int
    nrows::Int
    ncols::Int
    trans::Bool
end

function TrajectoryMatrix{T}(data::AbstractArray{T}, window; trans=false)
    src_size = size(data)
    nsteps = src_size[1]
    step_size = isempty(src_size[2:end]) ? 1 : prod(src_size[2:end])

    ncols = window * step_size
    nrows = nsteps - window + 1

    if trans
        # Reverse the nrows and ncols if this is to be a transposed TrajectoryMatrix
        return TrajectoryMatrix{T,2}(data, window, nsteps, step_size, ncols, nrows, trans)
    else
        return TrajectoryMatrix{T,2}(data, window, nsteps, step_size, nrows, ncols, trans)
    end
end

Base.size{T}(matrix::TrajectoryMatrix{T}) = (matrix.nrows, matrix.ncols)
Base.linearindexing(::Type{TrajectoryMatrix}) = Base.LinearFast()

function compute_index{T}(matrix::TrajectoryMatrix{T}, i::Int, j::Int)
    @fastmath begin
        # Base case
        if !matrix.trans
            # If we think of the source data as a matrix
            # the step is the row
            step = i + floor((j-1) / matrix.step_size)

            # the sub_step is the column
            sub_step = ((j-1) % matrix.step_size) + 1

            # info("Internal Index: ( $step, $sub_step )")
            idx = round(Int, (size(matrix.data, 1) * (sub_step-1)) + step)

            return idx
        else
            # If we think of the source data as a matrix
            # the step is the row
            step = j + floor((i-1) / matrix.step_size)

            # the sub_step is the column
            sub_step = ((i-1) % matrix.step_size) + 1

            # info("Internal Index: ( $step, $sub_step )")
            idx = round(Int, (size(matrix.data, 1) * (sub_step-1)) + step)

            return idx
        end
    end
end

function Base.getindex{T}(matrix::TrajectoryMatrix{T}, i::Int, j::Int)
    return matrix.data[compute_index(matrix, i, j)]
end

function Base.getindex{T}(matrix::TrajectoryMatrix{T}, i::Int)
    idx = ind2sub(size(matrix), i)
    return matrix.data[compute_index(matrix, idx[1], idx[2])]
end

function Base.setindex!{T}(matrix::TrajectoryMatrix{T}, v::T, i::Int)
    idx = ind2sub(size(matrix), i)
    matrix.data[compute_index(matrix, idx[1], idx[2])] = v
end

function Base.setindex!{T}(matrix::TrajectoryMatrix{T}, v::T, i::Int, j::Int)
    return matrix.data[compute_index(matrix, i, j)] = v
end

function Base.transpose{T}(m::TrajectoryMatrix{T})
    if m.trans
        return TrajectoryMatrix(m.data, m.window; trans=false)
    else
        return TrajectoryMatrix(m.data, m.window; trans=true)
    end
end

function Base.ctranspose{T}(m::TrajectoryMatrix{T})
    data = map(i -> conj(i), m.data)

    if m.trans
        return TrajectoryMatrix(data, m.window; trans=false)
    else
        return TrajectoryMatrix(data, m.window; trans=true)
    end
end

function traj{T}(data::AbstractArray{T}, window; compressed=false)
    traj = TrajectoryMatrix(data, window)

    if compressed
        return traj
    else
        result = Array{T}(size(traj))
        for i in eachindex(traj)
            result[i] = traj[i]
        end

        return result
    end
end

ext: jl


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct ConceptMap
    data::Dict{String, AbstractConcept}
end

struct Explainable
    conjunction::String
    explainable_id::Int64
end

struct Explainables
    relations::Optional{Dict{String,Explainable}}
    attributes::Optional{Dict{String,Explainable}}
    ownerships::Optional{Dict{Tuple{String},Explainable}}
end

function Explainable(explainable::Proto.Explainable)
   return Explainable(explainable.conjunction, explainable.id)
end

Explainables() = Explainables(nothing, nothing, nothing)

function Explainables(explainables::Proto.Explainables)
    relations = Dict{String,Explainable}()
    for (var, explainable) in explainables.relations.items
        relations[var] = Explainable(explainable)
    end
    attributes = Dict{String,Explainable}()
    for (var, explainable) in explainables.attributes.items
        attributes[var] = Explainable(explainable)
    end
    ownerships = Dict{String,Explainable}()
    for (var, owned_map) in explainables.ownerships.items
        for (owned, explainable) in owned_map.owned.items
            ownerships[(var, owned)] = Explainable(explainable)
        end
    end
    return Explainables(relations, attributes, ownerships)
end

concepts(cm::ConceptMap) = values(cm.data)

Base.getindex(cm::ConceptMap, key::String) = cm.data[key]

function Base.show(io::IO, cm::ConceptMap)
    num_entries = length(cm.data)
    print(io, "ConceptMap:\n")
    for (i, k) in enumerate(sort(collect(keys(cm.data))))
        #TODO: decide which type of printing looks better. For me (Frank) it looks better
        #with a line break after each line.
        #i != num_entries && print(io, "\n")
        print(io, lpad("$i", 5), ". ", k, " => ", cm.data[k], "\n")
    end
end

function ConceptMap(res::Proto.ConceptMap)
    return ConceptMap(Dict(key => instantiate(concept) for (key, concept) in res.map))
end

function ConceptMap(array_concept_maps::Vector{Proto.ConceptMap})
    # return map(entry->ConceptMap(entry), array_concept_maps)
    return ConceptMap.(array_concept_maps)
end


text:
# Retrieve name of example and output directory
if length(ARGS) != 2
    error("please specify the name of the example and the output directory")
end
const EXAMPLE = ARGS[1]
const OUTDIR = ARGS[2]

# Activate environment
# Note that each example's Project.toml must include Literate as a dependency
using Pkg: Pkg
const EXAMPLEPATH = joinpath(@__DIR__, "..", "examples", EXAMPLE)
Pkg.activate(EXAMPLEPATH)
Pkg.instantiate()
using Literate: Literate

# Convert to markdown and notebook
const SCRIPTJL = joinpath(EXAMPLEPATH, "script.jl")
Literate.markdown(SCRIPTJL, OUTDIR; name=EXAMPLE, documenter=false, execute=true)


text:
module PredictMDExtra # Begin module PredictMDExtra

__precompile__(true)

include(joinpath("registry_url_list.jl"))
include(joinpath("package_directory.jl"))
include(joinpath("version.jl"))

include(joinpath("package_list.jl"))

include(joinpath("import_required_packages.jl"))

include(joinpath("import_all.jl"))

include(joinpath("welcome.jl"))
include(joinpath("init.jl"))

end # End module PredictMDExtra


text:
"""
Computes the cronbach's α
"""
function cronalpha(items)
    items = Matrix(items)
    vcmat = cov(items) # Variance-Covariance Matrix
    k = size(items, 2) # no of items
    σ2_Xi = tr(vcmat) # individual-item variances
    σ2_X = sum(vcmat) # total-variance

    α = (k/(k-1)) * (1-σ2_Xi/σ2_X)
end


text:
using DifferentialEquations
using Roots
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

function condition(u,t,integrator) # Event when event_f(u,t) == 0
 u[1]
end

function affect!(integrator)
end
#cb1 = SavingCallback((u,t,integrator)->(u[4],u[2]), saved_values)

cb2= ContinuousCallback(condition, affect!)
Energy=0.112012;
H=(2*Energy)^2
h=sqrt(H)
n_iter=100
t_end=1000

# H1(p,q,t)= ( (q[2]-q[1])^2+(p[2]-p[1])^2 )*( (q[2]+q[1])^2+(p[2]+p[1])^2 )
# H2(p,q,t)= (p[2]^4+2*p[2]^2*(q[1]^2-1)+(1+q[1]^2)^2 )*(q[2]^4+2*q[2]^2*(p[1]^2+1)+(p[1]^2-1)^2 )
# HAM(p,q,t)=H1(p,q,t)/H2(p,q,t)
Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))

ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    return
end
# refernce
# z=X+iP W=Q+iY . Plane Q=0, P=0
# z=u[1]+i u[2] ; w=u[3]+i u[4]

#Q=0, have P find Y
P_start=8e-3
P_end=1e-1
P=range(P_start,stop=P_end,length=n_iter)
Y=zeros(n_iter)
TEST1_0=zeros(n_iter)
TEST1_f=zeros(n_iter)
u1max=zeros(n_iter)
u2max=zeros(n_iter)
u3max=zeros(n_iter)
u4max=zeros(n_iter)
YofP_find(PP)=sqrt((h+PP^2-h*PP^2 )./(1+h-h*PP^2))
axis([ -.2,.2,-.02, .02]);

# figure(dpi=3600)
for k=1:n_iter
    #  u0=zeros(4)
    # Y[k]=YofP_find(P[k])
    Pval=P[k]
    Y_find1(y)=Hamil(0,y,0,Pval)-H
    u0=zeros(4)
    Y[k]=find_zero(Y_find1,.01, maxeval=100, tol=1e-15)
    u0[1]=0 #X
    u0[2]=P[k] #P
    u0[3]=0 #Q
    u0[4]=Y[k] #Y
    TEST1_0[k]=H_test(u0)/H
    prob = ODEProblem(Eq_of_M,u0,(0., t_end))
    t,A=solve(prob, Vern9(),reltol=1e-15,abstol=1e-12,maxiters=1e15, callback=cb2,save_end=false,save_everystep=false)
    plot(A[3,:],A[2,:],"k,")
    plot(A[3,:],-A[2,:],"k,")
    plot(-A[3,:],A[2,:],"k,")
    plot(-A[3,:],-A[2,:],"k,")

    #testing if energy is conserved at end
    uf=zeros(4)
    uf[1]=A[1,end]
    uf[2]=A[2,end]
    uf[3]=A[3,end]
    uf[4]=A[4,end]
    TEST1_f[k]=H_test(uf)/H

    #finding maxes for exit criteria
    u1max[k]=

 end
#P=0, have Q find Y

Q_start=5e-3
Q_end=1e-1
Q=range(Q_start,stop=Q_end,length=n_iter)
Y=zeros(n_iter)
TEST2=zeros(n_iter)
for k=1:n_iter
    Y_find(YY)=Hamil(0,YY,Q[k],0)-H
    u0=zeros(4)
    Y[k]=find_zero(Y_find,.01, maxeval=100,tol =1e-15)
    u0[1]=0 #X
    u0[2]=0 #P
    u0[3]=Q[k] #Q
    u0[4]=Y[k] #Y
    TEST2[k]=H_test(u0)
    prob = ODEProblem(Eq_of_M,u0,(0., t_end))
    t,A=solve(prob, Vern9(),reltol=1e-15,abstol=1e-12,maxiters=1e15, callback=cb2,save_end=false,save_everystep=false)
    plot(A[3,:],A[2,:],"k,")
    plot(A[3,:],-A[2,:],"k,")
    plot(-A[3,:],A[2,:],"k,")
    plot(-A[3,:],-A[2,:],"k,")
end
# savefig("PSS11.png")


text:
include("trackingcell.jl");
include("lineage.jl");

data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF";
ret_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200315_result";
#img_16_2 = load(File(format"TIFF","$data_dir/20200316/HE7-11-1-80uw-PWM_1_s2.ome.tiff"));
img_15_1 = load(File(format"TIFF","$data_dir/20200315/HE7-11-1-80uw-const_1_s1.ome.btf.tiff"));
img_edge, mask_edge, mask_clear = split_cell_LoG(img_15_1, size(img_15_1)[3]÷20);
save("$ret_dir/img_edge.tiff", img_edge);

nucleus_pos = [ component_centroids(label_components(mask_clear[:,:,i])) for i in 1:size(mask_clear)[3] ];

paths = generate_path(nucleus_pos);

for i in 1:length(paths)
	save("$ret_dir/s1_c$cell.tiff", N0f16.(export_cell(paths[i], img_15_l, mask_edge)));
end


text:
# Defaults for many of the interface functions are defined here as fallback.

## Coords
# Four options in SF, xy, xyz, xym, xyzm
const default_coord_names = (:X, :Y, :Z, :M)

coordnames(geom) = default_coord_names[1:ncoord(geom)]

x(geom) = getcoord(geom, findfirst(coordnames(geom), :X))
y(geom) = getcoord(geom, findfirst(coordnames(geom), :Y))
z(geom) = getcoord(geom, findfirst(coordnames(geom), :Z))
m(geom) = getcoord(geom, findfirst(coordnames(geom), :M))

## Points
npoint(c::AbstractCurve, geom) = ngeom(c, geom)
getpoint(c::AbstractCurve, geom, i) = getgeom(c, geom, i)

## Polygons
nring(p::AbstractPolygon, geom) = ngeom(p, geom)
getring(p::AbstractPolygon, geom, i) = getgeom(p, geom, i)
getexterior(p::AbstractPolygon, geom) = getring(p, geom, 1)
nhole(p::AbstractPolygon, geom) = nring(p, geom) - 1
gethole(p::AbstractPolygon, geom, i) = getring(p, geom, i + 1)

nlinestring(::AbstractMultiLineString, geom) = ngeom(p, geom)
nlinestring(::AbstractMultiLineString, geom, i) = getgeom(p, geom, i)

npolygon(::AbstractMultiPolygon, geom) = ngeom(p, geom)
getpolygon(::AbstractMultiPolygon, geom, i) = getgeom(p, geom, i)

npoint(::Line, _) = 2
npoint(::Triangle, _) = 3
npoint(::Rectangle, _) = 4
npoint(::Quad, _) = 4
npoint(::Pentagon, _) = 5
npoint(::Hexagon, _) = 6

issimple(::AbstractCurve, geom) = allunique([getpoint(geom, i) for i in 1:npoint(geom) - 1]) && allunique([getpoint(geom, i) for i in 2:npoint(geom)])
isclosed(::AbstractCurve, geom) = getpoint(geom, 1) == getpoint(geom, npoint(geom))
isring(x::AbstractCurve, geom) = issimple(x, geom) && isclosed(x, geom)

# TODO Only simple if it's also not intersecting itself, except for its endpoints
issimple(::AbstractMultiCurve, geom) = all(i -> issimple(getgeom(geom, i)), 1:ngeom(geom))
isclosed(::AbstractMultiCurve, geom) = all(i -> isclosed(getgeom(geom, i)), 1:ngeom(geom))

issimple(::MultiPoint, geom) = allunique((getgeom(geom, i) for i in 1:ngeom(geom)))


text:
# Autogenerated wrapper script for SwiftShader_jll for x86_64-linux-gnu
export libvulkan

JLLWrappers.@generate_wrapper_header("SwiftShader")
JLLWrappers.@declare_library_product(libvulkan, "libvk_swiftshader.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libvulkan,
        "lib/libvulkan.so.1",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function project_recurse(trie::Trie{Any,ChoiceOrCallRecord},
                         selection::Selection)
    weight = 0.
    for (key, choice_or_call) in get_leaf_nodes(trie)
        if choice_or_call.is_choice
            if key in selection
                weight += choice_or_call.score
            end
        else
            subselection = selection[key]
            weight += project(choice_or_call.subtrace_or_retval, subselection)
        end
    end
    for (key, subtrie) in get_internal_nodes(trie)
        subselection = selection[key]
        weight += project_recurse(subtrie, subselection)
    end
    weight
end

function project(trace::DynamicDSLTrace, selection::Selection)
    project_recurse(trace.trie, selection)
end

project(trace::DynamicDSLTrace, ::EmptySelection) = trace.noise


text:
# CUBLAS.jl
#
# Julia interface to CUBLAS.  Inspired by CUDArt.jl and using Clang.jl.
#
# Author: Nick Henderson <nwh@stanford.edu>
# Created: 2014-08-27
# License: MIT
#
__precompile__(true)
module CUBLAS

importall Base.LinAlg.BLAS

using CUDArt
using CUDArt.rt.cudaStream_t

const BlasChar = Char #import Base.LinAlg.BlasChar
import Base.one
import Base.zero

export band, unband, bandex

include("util.jl")
include("libcublas_types.jl")

# get cublas status message
function statusmessage(status)
    if status == CUBLAS_STATUS_SUCCESS
        return "cublas success"
    end
    if status == CUBLAS_STATUS_NOT_INITIALIZED
        return "cublas not initialized"
    end
    if status == CUBLAS_STATUS_ALLOC_FAILED
        return "cublas alloc failed"
    end
    if status == CUBLAS_STATUS_INVALID_VALUE
        return "cublas invalid value"
    end
    if status == CUBLAS_STATUS_ARCH_MISMATCH
        return "cublas arch mismatch"
    end
    if status == CUBLAS_STATUS_MAPPING_ERROR
        return "cublas mapping error"
    end
    if status == CUBLAS_STATUS_EXECUTION_FAILED
        return "cublas execution failed"
    end
    if status == CUBLAS_STATUS_INTERNAL_ERROR
        return "cublas internal error"
    end
    if status == CUBLAS_STATUS_NOT_SUPPORTED
        return "cublas not supported"
    end
    if status == CUBLAS_STATUS_LICENSE_ERROR
        return "cublas license error"
    end
    return "cublas unknown status"
end

# error handling function
function statuscheck(status)
    if status == CUBLAS_STATUS_SUCCESS
        return nothing
    end
    # Because try/finally may disguise the source of the problem,
    # let's show a backtrace here
    warn("CUBLAS error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    println()
    throw(statusmessage(status))
end

# find the cublas library
const libcublas = Libdl.find_library(["libcublas"], ["/usr/local/cuda"])
if isempty(libcublas)
    error("CUBLAS library cannot be found. Please make sure that CUDA is installed")
end

include("libcublas.jl")

# setup cublas handle
const cublashandle = cublasHandle_t[0]

function __init__()
    cublasCreate_v2(cublashandle)
    # destroy cublas handle at julia exit
    atexit(()->cublasDestroy_v2(cublashandle[1]))
end

include("blas.jl")
include("highlevel.jl")

end # module


text:
# ---
# title: 705. Design HashSet
# id: problem705
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Design
# link: <https://leetcode.com/problems/design-hashset/description/>
# hidden: true
# ---
#
# Design a HashSet without using any built-in hash table libraries.
#
# To be specific, your design should include these functions:
#
#   * `add(value)`: Insert a value into the HashSet.
#   * `contains(value)` : Return whether the value exists in the HashSet or not.
#   * `remove(value)`: Remove a value in the HashSet. If the value does not exist in the HashSet, do nothing.
#
#
# **Example:**
#
#
#
#     MyHashSet hashSet = new MyHashSet();
#     hashSet.add(1);
#     hashSet.add(2);
#     hashSet.contains(1);    // returns true
#     hashSet.contains(3);    // returns false (not found)
#     hashSet.add(2);
#     hashSet.contains(2);    // returns true
#     hashSet.remove(2);
#     hashSet.contains(2);    // returns false (already removed)
#
#
#
# **Note:**
#
#   * All values will be in the range of `[0, 1000000]`.
#   * The number of operations will be in the range of `[1, 10000]`.
#   * Please do not use the built-in HashSet library.
#
#
## @lc code=start
using LeetCode

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


text:
abstract type AbstractSpectralElementSpace <: AbstractSpace end

Topologies.nlocalelems(space::AbstractSpectralElementSpace) =
    Topologies.nlocalelems(Spaces.topology(space))

local_geometry_data(space::AbstractSpectralElementSpace) = space.local_geometry

eachslabindex(space::AbstractSpectralElementSpace) =
    1:Topologies.nlocalelems(Spaces.topology(space))

function Base.show(io::IO, space::AbstractSpectralElementSpace)
    indent = get(io, :indent, 0)
    iio = IOContext(io, :indent => indent + 2)
    println(io, nameof(typeof(space)), ":")
    println(iio, " "^(indent + 2), space.topology)
    print(iio, " "^(indent + 2), space.quadrature_style)
end

topology(space::AbstractSpectralElementSpace) = space.topology
quadrature_style(space::AbstractSpectralElementSpace) = space.quadrature_style

"""
    SpectralElementSpace1D <: AbstractSpace

A one-dimensional space: within each element the space is represented as a polynomial.
"""
struct SpectralElementSpace1D{
    T,
    Q,
    GG <: Geometry.AbstractGlobalGeometry,
    LG,
    D,
} <: AbstractSpectralElementSpace
    topology::T
    quadrature_style::Q
    global_geometry::GG
    local_geometry::LG
    dss_weights::D
end

function SpectralElementSpace1D(
    topology::Topologies.IntervalTopology,
    quadrature_style,
)
    global_geometry = Geometry.CartesianGlobalGeometry()
    CoordType = Topologies.coordinate_type(topology)
    AIdx = Geometry.coordinate_axis(CoordType)
    FT = eltype(CoordType)
    nelements = Topologies.nlocalelems(topology)
    Nq = Quadratures.degrees_of_freedom(quadrature_style)

    LG = Geometry.LocalGeometry{AIdx, CoordType, FT, SMatrix{1, 1, FT, 1}}
    local_geometry = DataLayouts.IFH{LG, Nq}(Array{FT}, nelements)
    quad_points, quad_weights =
        Quadratures.quadrature_points(FT, quadrature_style)

    for elem in 1:nelements
        local_geometry_slab = slab(local_geometry, elem)
        for i in 1:Nq
            ξ = quad_points[i]
            # TODO: we need to massage the coordinate points because the grid is assumed 2D
            vcoords = Topologies.vertex_coordinates(topology, elem)
            x = Geometry.linear_interpolate(vcoords, ξ)
            ∂x∂ξ =
                (
                    Geometry.component(vcoords[2], 1) -
                    Geometry.component(vcoords[1], 1)
                ) / 2
            J = abs(∂x∂ξ)
            WJ = J * quad_weights[i]
            local_geometry_slab[i] = Geometry.LocalGeometry(
                x,
                J,
                WJ,
                Geometry.AxisTensor(
                    (
                        Geometry.LocalAxis{AIdx}(),
                        Geometry.CovariantAxis{AIdx}(),
                    ),
                    ∂x∂ξ,
                ),
            )
        end
    end
    dss_weights = copy(local_geometry.J)
    dss_weights .= one(FT)
    dss_1d!(dss_weights, dss_weights, local_geometry, topology, Nq)
    dss_weights = one(FT) ./ dss_weights

    return SpectralElementSpace1D(
        topology,
        quadrature_style,
        global_geometry,
        local_geometry,
        dss_weights,
    )
end

nlevels(space::SpectralElementSpace1D) = 1

const IntervalSpectralElementSpace1D =
    SpectralElementSpace1D{<:Topologies.IntervalTopology{<:Meshes.IntervalMesh}}

"""
    SpectralElementSpace2D <: AbstractSpace

A two-dimensional space: within each element the space is represented as a polynomial.
"""
struct SpectralElementSpace2D{
    T,
    Q,
    GG <: Geometry.AbstractGlobalGeometry,
    LG,
    D,
    IS,
    BS,
} <: AbstractSpectralElementSpace
    topology::T
    quadrature_style::Q
    global_geometry::GG
    local_geometry::LG
    dss_weights::D
    internal_surface_geometry::IS
    boundary_surface_geometries::BS
end

"""
    SpectralElementSpace2D(topology, quadrature_style)

Construct a `SpectralElementSpace2D` instance given a `topology` and `quadrature`.
"""
function SpectralElementSpace2D(topology, quadrature_style)
    domain = Topologies.domain(topology)
    if domain isa Domains.SphereDomain
        CoordType3D = Topologies.coordinate_type(topology)
        FT = Geometry.float_type(CoordType3D)
        CoordType2D = Geometry.LatLongPoint{FT} # Domains.coordinate_type(topology)
        global_geometry =
            Geometry.SphericalGlobalGeometry(topology.mesh.domain.radius)
    else
        CoordType2D = Topologies.coordinate_type(topology)
        FT = Geometry.float_type(CoordType2D)
        global_geometry = Geometry.CartesianGlobalGeometry()
    end
    AIdx = Geometry.coordinate_axis(CoordType2D)
    nelements = Topologies.nlocalelems(topology)
    Nq = Quadratures.degrees_of_freedom(quadrature_style)

    LG = Geometry.LocalGeometry{AIdx, CoordType2D, FT, SMatrix{2, 2, FT, 4}}

    local_geometry = DataLayouts.IJFH{LG, Nq}(Array{FT}, nelements)
    quad_points, quad_weights =
        Quadratures.quadrature_points(FT, quadrature_style)

    for elem in 1:nelements
        local_geometry_slab = slab(local_geometry, elem)
        for i in 1:Nq, j in 1:Nq
            if domain isa Domains.SphereDomain
                # compute the coordinate and partial derivative matrices for each quadrature point
                # Guba (2014))
                ξ = SVector(quad_points[i], quad_points[j])
                x = Geometry.spherical_bilinear_interpolate(
                    CoordType3D.(
                        Topologies.vertex_coordinates(topology, elem),
                    ),
                    ξ[1],
                    ξ[2],
                    global_geometry.radius,
                )
                u = Geometry.LatLongPoint(x, global_geometry)
                # [∂x1/∂ξ¹ ∂x1/∂ξ²
                #  ∂x2/∂ξ¹ ∂x2/∂ξ²
                #  ∂x3/∂ξ¹ ∂x3/∂ξ²]
                ∂x∂ξ = ForwardDiff.jacobian(ξ) do ξ
                    Geometry.components(
                        Geometry.spherical_bilinear_interpolate(
                            CoordType3D.(
                                Topologies.vertex_coordinates(topology, elem),
                            ),
                            ξ[1],
                            ξ[2],
                            global_geometry.radius,
                        ),
                    )
                end
                ϕ = u.lat
                λ = u.long
                # [∂u/∂x1 ∂u/∂x2 ∂u/∂x3
                #  ∂v/∂x1 ∂v/∂x2 ∂v/∂x3]
                # at the pole we orient u and v by taking the limit approaching
                # from the line λ == 0
                ∂u∂x = if ϕ == 90
                    # north pole => u axis is aligned with x2, v is aligned with -x1
                    @assert λ == 0
                    @SMatrix [
                        0 one(ϕ) 0
                        -one(ϕ) 0 0
                    ]
                elseif ϕ == -90
                    # south pole => u axis is aligned with x2, v is aligned with x1
                    @assert λ == 0
                    @SMatrix [
                        0 one(ϕ) 0
                        one(ϕ) 0 0
                    ]
                else
                    #=
                        # TODO: this might be more stable?
                        [
                            -sind(λ) cosd(λ) 0
                            -sind(ϕ)*cosd(λ) -sind(ϕ)*sind(λ) cosd(ϕ)
                        ]
                    =#
                    @SMatrix [
                        -sind(λ) cosd(λ) 0
                        0 0 1/cosd(ϕ)
                    ]
                end
                ∂u∂ξ = ∂u∂x * ∂x∂ξ
            else

                # this hard-codes a bunch of assumptions, and will unnecesarily duplicate data
                # e.g. where all metric terms are uniform over the space
                # alternatively: move local_geometry to a different object entirely, to support overintegration
                # (where the integration is of different order)
                ξ = SVector(quad_points[i], quad_points[j])
                u = Geometry.bilinear_interpolate(
                    CoordType2D.(
                        Topologies.vertex_coordinates(topology, elem),
                    ),
                    ξ[1],
                    ξ[2],
                )
                ∂u∂ξ = ForwardDiff.jacobian(ξ) do ξ
                    local x
                    x = Geometry.bilinear_interpolate(
                        CoordType2D.(
                            Topologies.vertex_coordinates(topology, elem),
                        ),
                        ξ[1],
                        ξ[2],
                    )
                    SVector(Geometry.component(x, 1), Geometry.component(x, 2))
                end
            end
            J = det(∂u∂ξ)
            WJ = J * quad_weights[i] * quad_weights[j]

            local_geometry_slab[i, j] = Geometry.LocalGeometry(
                u,
                J,
                WJ,
                Geometry.AxisTensor(
                    (
                        Geometry.LocalAxis{AIdx}(),
                        Geometry.CovariantAxis{AIdx}(),
                    ),
                    ∂u∂ξ,
                ),
            )
        end
    end

    # dss_weights = J ./ dss(J)
    dss_weights = copy(local_geometry.J)
    dss_2d!(dss_weights, local_geometry.J, local_geometry, topology, Nq)
    dss_weights .= local_geometry.J ./ dss_weights

    SG = Geometry.SurfaceGeometry{
        FT,
        Geometry.AxisVector{FT, Geometry.LocalAxis{AIdx}, SVector{2, FT}},
    }
    interior_faces = Topologies.interior_faces(topology)

    if quadrature_style isa Quadratures.GLL
        internal_surface_geometry =
            DataLayouts.IFH{SG, Nq}(Array{FT}, length(interior_faces))
        for (iface, (elem⁻, face⁻, elem⁺, face⁺, reversed)) in
            enumerate(interior_faces)
            internal_surface_geometry_slab =
                slab(internal_surface_geometry, iface)

            local_geometry_slab⁻ = slab(local_geometry, elem⁻)
            local_geometry_slab⁺ = slab(local_geometry, elem⁺)

            for q in 1:Nq
                sgeom⁻ = compute_surface_geometry(
                    local_geometry_slab⁻,
                    quad_weights,
                    face⁻,
                    q,
                    false,
                )
                sgeom⁺ = compute_surface_geometry(
                    local_geometry_slab⁺,
                    quad_weights,
                    face⁺,
                    q,
                    reversed,
                )

                @assert sgeom⁻.sWJ ≈ sgeom⁺.sWJ
                @assert sgeom⁻.normal ≈ -sgeom⁺.normal

                internal_surface_geometry_slab[q] = sgeom⁻
            end
        end

        boundary_surface_geometries =
            map(Topologies.boundary_tags(topology)) do boundarytag
                boundary_faces =
                    Topologies.boundary_faces(topology, boundarytag)
                boundary_surface_geometry =
                    DataLayouts.IFH{SG, Nq}(Array{FT}, length(boundary_faces))
                for (iface, (elem, face)) in enumerate(boundary_faces)
                    boundary_surface_geometry_slab =
                        slab(boundary_surface_geometry, iface)
                    local_geometry_slab = slab(local_geometry, elem)
                    for q in 1:Nq
                        boundary_surface_geometry_slab[q] =
                            compute_surface_geometry(
                                local_geometry_slab,
                                quad_weights,
                                face,
                                q,
                                false,
                            )
                    end
                end
                boundary_surface_geometry
            end
    else
        internal_surface_geometry = nothing
        boundary_surface_geometries = nothing
    end
    return SpectralElementSpace2D(
        topology,
        quadrature_style,
        global_geometry,
        local_geometry,
        dss_weights,
        internal_surface_geometry,
        boundary_surface_geometries,
    )
end

nlevels(space::SpectralElementSpace2D) = 1

const RectilinearSpectralElementSpace2D =
    SpectralElementSpace2D{<:Topologies.Topology2D{<:Meshes.RectilinearMesh}}

const CubedSphereSpectralElementSpace2D = SpectralElementSpace2D{
    <:Topologies.Topology2D{<:Meshes.AbstractCubedSphere},
}

function compute_surface_geometry(
    local_geometry_slab,
    quad_weights,
    face,
    q,
    reversed = false,
)
    Nq = length(quad_weights)
    @assert size(local_geometry_slab) == (Nq, Nq, 1, 1, 1)
    i, j = Topologies.face_node_index(face, Nq, q, reversed)

    local_geometry = local_geometry_slab[i, j]
    @unpack J, ∂ξ∂x = local_geometry

    # surface mass matrix
    n = if face == 4
        -J * ∂ξ∂x[1, :] * quad_weights[j]
    elseif face == 2
        J * ∂ξ∂x[1, :] * quad_weights[j]
    elseif face == 1
        -J * ∂ξ∂x[2, :] * quad_weights[i]
    elseif face == 3
        J * ∂ξ∂x[2, :] * quad_weights[i]
    end
    sWJ = norm(n)
    n = n / sWJ
    return Geometry.SurfaceGeometry(sWJ, n)
end

function variational_solve!(data, space::AbstractSpace)
    data .= RecursiveApply.rdiv.(data, space.local_geometry.WJ)
end

"""
    SpectralElementSpaceSlab <: AbstractSpace

A view into a `SpectralElementSpace2D` for a single slab.
"""
struct SpectralElementSpaceSlab{Q, G} <: AbstractSpectralElementSpace
    quadrature_style::Q
    local_geometry::G
end

const SpectralElementSpaceSlab1D =
    SpectralElementSpaceSlab{Q, DL} where {Q, DL <: DataLayouts.DataSlab1D}

const SpectralElementSpaceSlab2D =
    SpectralElementSpaceSlab{Q, DL} where {Q, DL <: DataLayouts.DataSlab2D}

nlevels(space::SpectralElementSpaceSlab1D) = 1
nlevels(space::SpectralElementSpaceSlab2D) = 1

function slab(space::AbstractSpectralElementSpace, v, h)
    SpectralElementSpaceSlab(
        space.quadrature_style,
        slab(space.local_geometry, v, h),
    )
end


text:
module FindComplexZeros

export countZeros, findZerosWithSubdivision, ZerosNearContourException

include("algorithm.jl")

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export radiobuttongroup

"""
    radiobuttongroup(;kwargs...)

A RadioButtonGroup component.
Radio Group Component
Keyword arguments:
- `legendText` (a list of or a singular dash component, string or number; optional): Provide a legend to the RadioButtonGroup input that you are
exposing to the user
- `radiosButtons` (optional): The radios buttons inside this group. radiosButtons has the following type: Array of lists containing elements 'label', 'value', 'id'.
Those elements have the following types:
  - `label` (String; optional)
  - `value` (String; optional)
  - `id` (String; optional)s
- `id` (String; optional): Specify a custom `id` for the radio group
- `value` (String; optional): The value selected
- `labelPosition` (a value equal to: 'left', 'right'; optional): Provide where label text should be placed
- `orientation` (a value equal to: 'horizontal', 'vertical'; optional): Provide where radio buttons should be placed
- `disabled` (Bool; optional): Specify whether the group is disabled
- `className` (String; optional): Provide an optional className to be applied to the container node
"""
function radiobuttongroup(; kwargs...)
        available_props = Symbol[:legendText, :radiosButtons, :id, :value, :labelPosition, :orientation, :disabled, :className]
        wild_props = Symbol[]
        return Component("radiobuttongroup", "RadioButtonGroup", "dash_carbon_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
"""
A simple Chatbot.

# Examples
```julia
julia> include("chatbot.jl")
Hi there! What is you name?: Julia
Hello, Julia! How are you?
Julia: I'm great!
Chatty: Cool!
Julia: What is your name?
Chatty: I'm Chatty
Julia: How are you Chatty?
Chatty: I'm happy as always
Julia: What is your hobby?
Chatty: I love to play the piano
Julia: exit
Good Bye!
```

# Contributors:
- [Sedat Aybars Nazlica](https://github.com/saybarsnazlica)
"""

# Chatbot type
struct Chatbot
    name::String
    mood::String
    hobby::String
end

# Parse the user input to the question for the chatbot
function parse_input(input::String)
    if isempty(input)
        return "nothing to say"
    end
    if occursin("name", input)
        input = "what is your name?"
    elseif occursin("hobby", input)
        input = "what is your hobby?"
    elseif occursin("how are", input)
        input = "how are you?"
    else
        input = "default"
    end
    return input
end

# Respond to the question
# Give a random response so the conversation feels more natural
function respond(bot::Chatbot, question::String)
    responses = Dict(
        "what is your name?" => ["I'm $(bot.name)", "Call me $(bot.name)", "My friends call me $(bot.name)"],
        "how are you?" => ["I feel $(bot.mood)", "I'm so $(bot.mood)", "I'm $(bot.mood) as always"],
        "what is your hobby?" => ["I $(bot.hobby)", "When I'm not busy, I $(bot.hobby)", "I love to $(bot.hobby)"],
        "nothing to say" => ["Hey, I can't hear you!", "Your silence means a lot", "Hmm..."],
        "default" => ["It's nice chatting with you", "Cool!", "Oh!"]
    )
    response = rand(responses[question])
    println("$(bot.name): ", response)
end

# Main function to chat
function chat()
    bot = Chatbot("Chatty", "happy", "play the piano")
    username = Base.prompt("Hi there! What is you name?")
    println("Hello, $(username)! How are you?")
    while true
        user_input = lowercase(Base.prompt("$(username)"))
        if user_input == "exit"
            println("Good Bye!")
            break
        else
            question = parse_input(user_input)
            respond(bot, question)
        end
    end
end

# Call the chat() function
chat()


text:
using Plots
using LinearAlgebra

x = 3 # True slope
a = -2:.25:2
b = a*x + randn(length(a)) # Add noise

p1 = plot(a, x*a, label = "True line") # True relationship
scatter!(a, b, label = "Noisy data") # Noisy measurements

U, S, V = svd(reshape(collect(a), length(a), 1))
xtilde = (V*inv(diagm(S))*U'*b)[1] # Least-square fit

plot!(a, xtilde*a, label = "Regression line") # Plot fit


text:
# n,m = map(x->parse(Int64,x),split(readline("rosalind_fibd.txt")," "))

function mortal_fibb(n::Int64,m::Int64)::Array{Int64,1}
	#initialse array
	rabbits = zeros(Int64,m)
	rabbits[end] = 1
	#println(rabbits)
	for _ in 2:n
		newborns = sum(rabbits[1:(end-1)])
		rabbits[1:(end-1)] = rabbits[2:end]
		rabbits[end] = newborns
		#println(rabbits)
	end
	return rabbits

end

sum(mortal_fibb(6,3))

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
module MACCOR

using DataFrames
using CSV
using Dates

include("clean_files.jl")

end # module


text:
struct SimpleOrbit{T,S,V,R} <: AbstractOrbit
    period::T
    t0::T
    b::S
    duration::T
    r_star::V

    b_norm::V
    speed::R
    half_period::T
    ref_time::T
end

"""
    SimpleOrbit(; period, duration, t0=0, b=0, r_star=1)

Circular orbit parameterized by the basic observables of a transiting system.

# Parameters
* `period` - The orbital period of the planets, nominally in days
* `duration` The duration of the transit, same units as `period`
* `t0` - The midpoint time of the reference transit, same units as `period`
* `b` - The impact parameter of the orbit
* `r_star` - The radius of the star, nominally in solar radii
"""
function SimpleOrbit(;period, duration, t0=zero(period), b=0.0, r_star=1.0)
    half_period = 0.5 * period
    duration > half_period && error("duration cannot be longer than half the period")
    b_norm = b * r_star
    speed = 2 * sqrt(r_star^2 - b_norm^2) / duration
    ref_time =  t0 - half_period
    # normalize time types
    period, t0, duration, half_period, ref_time = promote(period, t0, duration, half_period, ref_time)
    # normalize unitless types
    r_star, b_norm = promote(r_star, b_norm)
    SimpleOrbit(period, t0, b, duration, r_star, b_norm, speed, half_period, ref_time)
end

period(orbit::SimpleOrbit) = orbit.period
duration(orbit::SimpleOrbit) = orbit.duration

planet_position(orbit::SimpleOrbit, t) = relative_position(orbit, t)

function relative_position(orbit::SimpleOrbit, t)
    dt = (t - orbit.ref_time) % period(orbit) - orbit.half_period
    x = orbit.speed * dt
    y = orbit.b_norm
    z = abs(dt) < 0.5 * orbit.duration ? one(x) : -one(x)
    return SA[x, y, z]
end

function in_transit(orbit::SimpleOrbit, t, r; texp=0)
    dt = (t - orbit.ref_time) % period(orbit) - orbit.half_period
    tol = sqrt((r + orbit.r_star)^2 - orbit.b_norm^2) / orbit.speed
    tol += 0.5 * texp
    return abs(dt) < tol
end

function in_transit(orbit::SimpleOrbit, t; texp=0)
    dt = (t - orbit.ref_time) % period(orbit) - orbit.half_period
    tol = 0.5 * (orbit.duration + texp)
    return abs(dt) < tol
end

function flip(orbit::SimpleOrbit, r_planet)
    t0 = orbit.t0 + orbit.half_period
    b = orbit.b_norm / r_planet
    duration = 2 * sqrt(r_planet^2 + orbit.b_norm^2) / orbit.speed
    return SimpleOrbit(orbit.period,
                       t0,
                       b,
                       duration,
                       r_planet,
                       orbit.b_norm,
                       orbit.speed,
                       orbit.half_period,
                       orbit.t0)
end


text:
abstract type AbstractDimensionalityReduction <: RegressionModel end

"""
    projection(model::AbstractDimensionalityReduction)

Return the projection matrix of the model.
"""
projection(model::AbstractDimensionalityReduction) = error("'projection' is not defined for $(typeof(model)).")

"""
    reconstruct(model::AbstractDimensionalityReduction, y)

Return the model response (a.k.a. the dependent variable).
"""
reconstruct(model::AbstractDimensionalityReduction, y) = error("'reconstruct' is not defined for $(typeof(model)).")

abstract type LinearDimensionalityReduction <: AbstractDimensionalityReduction end

"""
    loadings(model::LinearDimensionalityReduction)

Return the model loadings (a.k.a. eigenvectors scaled up by the variances).
"""
loadings(model::LinearDimensionalityReduction) = error("'loadings' is not defined for $(typeof(model)).")

abstract type NonlinearDimensionalityReduction <: AbstractDimensionalityReduction end
abstract type LatentVariableDimensionalityReduction <: AbstractDimensionalityReduction end


text:
## Neuron name
# name = "spindle/whole_neurons/03a_spindle2aFI"
name = "spindle/whole_neurons/03b_spindle4aACC"
# name = "spindle/whole_neurons/03b_spindle7aACC"
# name = "spindle/whole_neurons/04b_spindle3aFI"
# name = "spindle/whole_neurons/06b_spindle8aACC"
# name = "spindle/whole_neurons/16o_spindle13aFI"
# name = "spindle/whole_neurons/16o_spindle13aFI.stl"
# name = "spindle/separated_neurons/03a_spindle2aFI_soma"
# name = "spindle/separated_neurons/03a_spindle2aFI_dendrites_1"
# name = "spindle/separated_neurons/03a_spindle2aFI_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle4aACC_dendrites_1"
# name = "spindle/separated_neurons/03b_spindle4aACC_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle4aACC_soma"
# name = "spindle/separated_neurons/03b_spindle6aACC_dendrites_1"
# name = "spindle/separated_neurons/03b_spindle6aACC_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle6aACC_soma"
# name = "spindle/separated_neurons/03b_spindle7aACC_dendrites_1"
# name = "spindle/separated_neurons/03b_spindle7aACC_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle7aACC_soma"
# name = "spindle/separated_neurons/04b_spindle3aFI_dendrites_1"
# name = "spindle/separated_neurons/04b_spindle3aFI_dendrites_2"
# name = "spindle/separated_neurons/04b_spindle3aFI_soma"
# name = "spindle/separated_neurons/19o_spindle14aFI_dendrites_2.stl"
# name = "pyramidal/separated_neurons/02b_pyramidal1aACC_dendrites"
# name = "pyramidal/separated_neurons/02b_pyramidal1aACC_soma"
# name = "pyramidal/whole_neurons/02b_pyramidal1aACC"
# name = "pyramidal/whole_neurons/02a_pyramidal2aFI"
# name = "pyramidal/whole_neurons/04a_pyramidal4aACC"
# name = "pyramidal/whole_neurons/04a_pyramidal5aACC"
# name = "pyramidal/whole_neurons/04b_pyramidal6aACC"
# name = "pyramidal/whole_neurons/04b_pyramidal6aFI"
# name = "pyramidal/whole_neurons/25o_pyramidal18aFI"

## Setup
setup = NeuronSetup(
    name = name,
    ecs_shape = "no_ecs",
    ecs_ratio = 0.3,
    # refinement = 0.5,
    D_in = 0.002 * I(3),
    D_out = 0.002 * I(3),
    D_ecs = 0.002 * I(3),
    T₂_in = Inf,
    T₂_out = Inf,
    T₂_ecs = Inf,
    ρ_in = 1.0,
    ρ_out = 1.0,
    ρ_ecs = 1.0,
    κ_in_out = 1e-3,
    κ_out_ecs = 1e-3,
    κ_in = 0.0,
    κ_out = 0.0,
    κ_ecs = 0.0,
)

experiment = Experiment(
    gradient = (
        directions = create_directions([1.0; 0.0; 0.0]),
        sequences = [PGSE(2500.0, 4000.0)],
        values = [1000.0],
        values_type = "b",
    ),
    btpde = (odesolver = QNDF(), reltol = 1e-4, abstol = 1e-6, nsave = 1),
    btpde_midpoint = (θ = 0.5, timestep = 5),
    mf = (length_scale = 3, neig_max = 400, ninterval = 500),
    hadc = (odesolver = QNDF(), reltol = 1e-4, abstol = 1e-6),
    karger = (odesolver = MagnusGL6(), timestep = 5),
)


text:
using JetPack, Jets, Test

@testset "JopDerivative, 2D, along first dimension" begin
    x = range(0, stop=4*pi, length=256)
    m = sin.(x) * ones(1,10)
    A = JopDerivative(JetSpace(Float64,256,10), delta=x[2]-x[1])
    d = A*m
    a = A'*d

    dim=1
    Rpre = CartesianIndices(size(m)[1:dim-1])
    Rpost = CartesianIndices(size(m)[dim+1:end])
    n = size(m, dim)

    d_expected = cos.(x) * ones(1,10)
    for j = 1:10, i = 5:length(x)-5
        @test isapprox(d[i,j], d_expected[i,j], rtol=1e-5)
    end

    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopDerivative, 3D, along first dimension" begin
    x = range(0, stop=4*pi, length=256)
    m = zeros(256,10,1)
    m[:,:,:] = sin.(x) * ones(1,10)
    A = JopDerivative(JetSpace(Float64,256,10,1), delta=x[2]-x[1])
    d = A*m
    a = A'*d

    d_expected = cos.(x) * ones(1,10)
    for j = 1:10, i = 5:length(x)-5
        @test isapprox(d[i,j,1], d_expected[i,j,1], rtol=1e-5)
    end

    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end


text:
"""
        dematel(comparisonMat; threshold = nothing)

Apply Dematel for a given comparison matrix.

# Arguments:
 - `comparisonMat::DataFrame`: n × n matrix of input values in DataFrame format (for convenience) where n is the number of criteria. Values are either 0, 1, 2, 3, or 4 which indicate the level of importance of the ith criterion relative to the jth criterion.
 - `threshold::Union{Nothing, Float64}`: Optional threshold used in calculating values of influence matrix. It is calculated when the argument is omitted.

# Description
dematel() applies the Dematel method to calculate criteria weights, possibly for use in another multi-criteria decision making tool.

# Output
- `::DematelResult`: DematelResult object that holds many results including weights for each single criterion.

# Examples
```julia-repl
julia> K = [
        0 3 0 2 0 0 0 0 3 0;
        3 0 0 0 0 0 0 0 0 2;
        4 1 0 2 1 3 1 2 3 2;
        4 1 4 0 1 2 0 1 0 0;
        3 2 3 1 0 3 0 2 0 0;
        4 1 4 4 0 0 0 1 1 3;
        3 0 0 0 0 2 0 0 0 0;
        3 0 4 3 2 3 1 0 0 0;
        4 3 2 0 0 1 0 0 0 2;
        2 1 0 0 0 0 0 0 3 0
    ];

julia> dmat = makeDecisionMatrix(K);
julia> result = dematel(dmat);

julia> result.weights
10-element Array{Float64,1}:
 0.1686568559124561
 0.07991375718719543
 0.14006200243438863
 0.10748052790517183
 0.08789022388276985
 0.12526272598854982
 0.03067915023486491
 0.10489168834828348
 0.092654758940811
 0.06250830916550884
```

# References
Celikbilek Yakup, Cok Kriterli Karar Verme Yontemleri, Aciklamali ve Karsilastirmali
Saglik Bilimleri Uygulamalari ile. Editor: Muhlis Ozdemir, Nobel Kitabevi, Ankara, 2018
"""
function dematel(comparisonMat::DataFrame; threshold::Union{Nothing,Float64}=nothing)::DematelResult

    n, _ = size(comparisonMat)

    csums =  colsums(comparisonMat)
    rsums =  rowsums(comparisonMat)

    largest = maximum(vcat(csums, rsums))

    # K = convert(Array{Float64,2}, comparisonMat)
    K = Matrix{Float64}(comparisonMat)
    ND = K ./ largest

    T = ND * inv(I(n) - ND)

    c = rowsums(T)
    r = colsums(T)

    c_plus_r  = c .+ r
    c_minus_r = c .- r

    if threshold === nothing
        threshold = sum(T) / length(T)
    end

    influence_matrix = zeros(Float64, n, n)
    for i in 1:n
        for j in 1:n
            if T[i, j] > threshold
                influence_matrix[i, j] = 1
    end
        end
    end

    w_row = sqrt.((c_plus_r.^2.0) .+ (c_minus_r.^2.0))
    w = w_row ./ sum(w_row)

    result = DematelResult(
        comparisonMat,
        ND,
        T,
        c,
        r,
        c_plus_r,
        c_minus_r,
        threshold,
        influence_matrix,
        w
    )

    return result
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
module SurfaceCouplingTests

using Test
using Gridap
using Gridap.Arrays
using Gridap.FESpaces
import Gridap: ∇
using LinearAlgebra: tr, ⋅

# Analytical functions

u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
∇u(x) = TensorValue( 2*x[1], 4*x[2], -2*x[1], zero(x[1]) )
Δu(x) = VectorValue( 6, -2 )

p(x) = x[1] + 3*x[2]
∇p(x) = VectorValue(1,3)

s(x) = -Δu(x)
f(x) = -Δu(x) + ∇p(x)
g(x) = tr(∇u(x))

∇(::typeof(u)) = ∇u
∇(::typeof(p)) = ∇p

# Geometry + Integration

n = 20
mesh = (n,n)
domain = 2 .* (0,1,0,1) .- 1
order = 1
model = CartesianDiscreteModel(domain, mesh)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,5])
add_tag_from_tags!(labels,"neumann",[6,7,8])

trian = Triangulation(model)

const R = 0.4

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

cell_to_coods = get_cell_coordinates(trian)
cell_to_is_solid = collect1d(apply(is_in,cell_to_coods))
cell_to_is_fluid = Vector{Bool}(.! cell_to_is_solid)

trian_solid = Triangulation(model, cell_to_is_solid)
trian_fluid = Triangulation(model, cell_to_is_fluid)

order = 2

degree = 2*order
quad = CellQuadrature(trian,degree)
quad_solid = CellQuadrature(trian_solid,degree)
quad_fluid = CellQuadrature(trian_fluid,degree)

btrian = BoundaryTriangulation(model,labels,"neumann")
bdegree = 2*order
bquad = CellQuadrature(btrian,bdegree)
n = get_normal_vector(btrian)

# This returns a SkeletonTriangulation whose normal vector
# goes outwards to the fluid domain.
trian_Γ = InterfaceTriangulation(model,cell_to_is_fluid)
n_Γ = get_normal_vector(trian_Γ)
quad_Γ = CellQuadrature(trian_Γ,bdegree)

# FESpaces

V = TestFESpace(
  model=model,
  valuetype=VectorValue{2,Float64},
  reffe=:QLagrangian,
  order=order,
  conformity =:H1,
  dirichlet_tags="dirichlet")

Q = TestFESpace(
  triangulation=trian_fluid,
  valuetype=Float64,
  order=order-1,
  reffe=:PLagrangian,
  conformity=:L2)

U = TrialFESpace(V,u)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

# FE problem

function a_solid(x,y)
  u,p = x
  v,q = y
  inner(∇(v),∇(u))
end

function l_solid(y)
  v,q = y
  v⋅s
end

function a_fluid(x,y)
  u,p = x
  v,q = y
  inner(∇(v),∇(u)) - (∇⋅v)*p + q*(∇⋅u)
end

function l_fluid(y)
  v,q = y
  v⋅f + q*g
end

function l_Γn_fluid(y)
  v,q = y
  v⋅(n⋅∇u) - (n⋅v)*p
end

# Pressure drop at the interface
function l_Γ(y)
  v,q = y
  - mean(n_Γ⋅v)*p
end

t_Ω_solid = AffineFETerm(a_solid,l_solid,trian_solid,quad_solid)
t_Ω_fluid = AffineFETerm(a_fluid,l_fluid,trian_fluid,quad_fluid)
t_Γn_fluid = FESource(l_Γn_fluid,btrian,bquad)
t_Γ = FESource(l_Γ,trian_Γ,quad_Γ)

op = AffineFEOperator(X,Y,t_Ω_solid,t_Ω_fluid,t_Γn_fluid,t_Γ)
uh, ph = solve(op)

# Visualization

ph_fluid = restrict(ph, trian_fluid)

eu = u - uh
ep = p - ph
ep_fluid = p - ph_fluid

#writevtk(trian_fluid,"trian_fluid",cellfields=["ph"=>ph_fluid, "ep"=>ep_fluid])
#
#writevtk(trian,"trian", cellfields=["uh" => uh, "ph"=> ph, "eu"=>eu, "ep"=>ep])

# Errors

l2(v) = v⋅v
h1(v) = v⋅v + inner(∇(v),∇(v))

eu_l2 = sqrt(sum(integrate(l2(eu),trian,quad)))
eu_h1 = sqrt(sum(integrate(h1(eu),trian,quad)))
ep_l2 = sqrt(sum(integrate(l2(ep_fluid),trian_fluid,quad_fluid)))

tol = 1.0e-9
@test eu_l2 < tol
@test eu_h1 < tol
@test ep_l2 < tol

end # module


text:
module DisposableResourceAllocation

using JuMP
using Gurobi

using LinearAlgebra
using MathOptInterface

export disposable_resource_allocation

function disposable_resource_allocation(
		initial_supply::Array{<:Real,1},
		supply::Array{<:Real,2},
		demand::Array{<:Real,2},
		adj_matrix::BitArray{2};
		objective::Symbol=:shortage,
		send_new_only::Bool=true,
		send_wait_period::Int=10,
		m::Real=1e-5,
		verbose::Bool=false,
)
	N, T = size(supply)
	@assert(size(initial_supply, 1) == N)
	@assert(size(demand, 1) == N)
	@assert(size(demand, 2) == T)
	@assert(size(adj_matrix, 1) == N)
	@assert(size(adj_matrix, 2) == N)
	@assert(objective in [:shortage, :overflow])

	model = Model(Gurobi.Optimizer)
	if !verbose set_silent(model) end

	@variable(model, sent[1:N,1:N,1:T] >= 0)
	@variable(model, used[1:N,1:T] >= 0)
	@variable(model, obj_dummy[1:N,1:T] >= 0)

	@objective(model, Min, sum(obj_dummy) + m*sum(sent))

	if send_new_only
		@constraint(model, [t=1:T],
			sum(sent[:,:,t], dims=2) .<= max.(0, supply[:,t])
		)
	end

	@constraint(model, [t=1:T],
		sum(sent[:,:,t], dims=2) .+ used[:,t] .<=
			initial_supply
			.+ sum(supply[:,1:t], dims=2)
			.- sum(used[:,1:t-1], dims=2)
			.- sum(sent[:,:,1:t-1], dims=[2,3])
			.+ sum(sent[:,:,1:t-1], dims=[1,3])
	)

	for i = 1:N
		for j = 1:N
			if ~adj_matrix[i,j]
				@constraint(model, sum(sent[i,j,:]) .== 0)
			end
		end
	end

	if send_wait_period > 0
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,i,t]), sum(sent[i,:,t:min(t+send_wait_period,T)])] in MOI.SOS1([1.0, 1.0])
		)
	end

	flip_sign = (objective == :shortage) ? 1 : -1
	@constraint(model, [i=1:N,t=1:T],
		obj_dummy[i,t] >= flip_sign * (demand[i,t] - used[i,t])
	)

	optimize!(model)
	return model
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
immutable XCubed <: Function end
@inline (::XCubed)(x) = x^3
derivative(::XCubed) = dXCubed()

immutable dXCubed <: Function end
@inline (::dXCubed)(x) = 3x^2
derivative(::dXCubed) = ddXCubed()

immutable ddXCubed <: Function end
@inline (::ddXCubed)(x) = 6x

immutable XSquaredLogX <: Function end
@inline (::XSquaredLogX)(x) = x^2 * log(x)
derivative(::XSquaredLogX) = dXSquaredLogX()

immutable dXSquaredLogX <: Function end
@inline (::dXSquaredLogX)(x) = 2x*log(x) + x
derivative(::dXSquaredLogX) = ddXSquaredLogX()

immutable ddXSquaredLogX <: Function end
@inline (::ddXSquaredLogX)(x) = 2*log(x) + 3


text:
@testset "MapReduce-style computation" begin
    run_with_job("Filling") do job
        println(typeof(Base.fill(1.0, 2048)))
        x = BanyanArrays.fill(10.0, 2048)
        println(typeof(x))
        x = map(e -> e / 10, x)
        println(typeof(x))
        res = sum(x)

        res = collect(res)
        @test typeof(res) == Float64
        @test res == 2048
    end

    run_with_job("Multiple evaluations apart") do job
        x = BanyanArrays.fill(10.0, 2048)
        x = map(e -> e / 10, x)
        res1 = collect(sum(x)) # Note: failed here with "key :val_6HTGdt08_idx_0 not found"
        res2 = collect(minimum(x))

        @test typeof(res1) == Float64
        @test res1 == 2048
        @test typeof(res2) == Float64
        @test res2 == 1.0
    end

    run_with_job("Multiple evaluations together") do job
        x = BanyanArrays.fill(10.0, 2048)
        x = map(e -> e / 10, x)
        res1 = sum(x)
        res2 = minimum(x)

        res1 = collect(res1)
        res2 = collect(res2)
        @test typeof(res1) == Float64
        @test res1 == 2048
        @test typeof(res2) == Float64
        @test res2 == 1.0
    end

    run_with_job("Simple computing") do job
        for _ in 1:8
            # NOTE: This also tests simple writing to and reading from local disk
            x = BanyanArrays.fill(10.0, 2048)
            # x = map(e -> e / 10, x)
            @show typeof(x)
            write_to_disk(x)
            # write_to_disk(x)
            @show typeof(x)
            sleep(15)
            @show typeof(x)
            # NOTE: The only reason why we're not putting `collect(x)` inside the
            # the `@test` is because `@test` will catch exceptions and prevent the
            # job from getting destroyed when an exception occurs and we can't keep
            # running this test if the job ends
            x_collect = collect(x)
            @test x_collect == Base.fill(10.0, 2048)
        end
    end

    run_with_job("Computing") do job
        # NOTE: This also tests simple writing to and reading from local disk
        x = BanyanArrays.fill(10.0, 2048)
        x = map(e -> e / 10, x)
        @show typeof(x)
        write_to_disk(x)
        write_to_disk(x)
        @show typeof(x)
        # NOTE: The only reason why we're not putting `collect(x)` inside the
        # the `@test` is because `@test` will catch exceptions and prevent the
        # job from getting destroyed when an exception occurs and we can't keep
        # running this test if the job ends
        x_collect = collect(x)
        @test x_collect == Base.fill(1.0, 2048)
        @show typeof(x)
        write_to_disk(x)
        x_collect = collect(x)
        @test x_collect == Base.fill(1.0, 2048)
        x_collect = collect(x)
        @test x_collect == Base.fill(1.0, 2048)
    end

    run_with_job("Re-computing") do job
        x = BanyanArrays.fill(10.0, 2048)
        x_sum = reduce(+, x)
        x = map(e -> e / 10, x)
        write_to_disk(x)
        write_to_disk(x_sum)
        x_sum_collect = collect(x_sum)
        @test x_sum_collect == 10.0 * 2048
        write_to_disk(x_sum)
        x_collect = collect(x)
        @show length(x_collect)
        @test x_collect == Base.fill(1.0, 2048)
        collect(x_sum)
        x_sum_collect = collect(x_sum)
        @test x_sum_collect == 10.0 * 2048
    end

    run_with_job("Map with multiple values") do job
        a = BanyanArrays.fill(10.0, 2048)
        b = BanyanArrays.fill(10.0, 2048)
        c = a + b
        c_sum_collect = collect(sum(c))
        @test c_sum_collect == 2048 * 10.0 * 2
    end

    run_with_job("Complex dependency graphs") do job
        # Here we test more complex dependency graphs where some values are destroyed

        x = BanyanArrays.fill(10.0, 2048)
        y = BanyanArrays.fill(10.0, 2048)
        a = BanyanArrays.fill(10.0, 2048)
        x += y
        x += a
        y_sum_collect = collect(sum(y))
        @test y_sum_collect == 2048 * 10.0
        a = nothing
        x_sum_collect = collect(sum(x))
        @test x_sum_collect == 2048 * 10.0 * 3
        y = nothing
        z = x + x
        z_sum_collect = collect(sum(z))
        @test z_sum_collect == 2048 * 10.0 * 6
        x_sum = sum(x)
        x=nothing
        x_sum_collect = collect(x_sum)
        @test x_sum_collect == 2048 * 10.0 * 3
    end

    run_with_job("Multiple arrays") do job
        x1 = BanyanArrays.fill(10.0, 2048)
        x2 = BanyanArrays.fill(10.0, 2048)
        res = map((a, b) ->  a * b, x1, x2)

        res_sum_collect = collect(sum(res))
        @test res_sum_collect == 204_800.0
        res_minimum_collect = collect(minimum(res))
        @test res_minimum_collect == 100.0
    end

    run_with_job("2D arrays") do job
        x1 = BanyanArrays.fill(1.0, (2048, 2048))
        x2 = BanyanArrays.fill(2.0, (2048, 2048))
        res = map((a, b) ->  a * b, x1, x2)
        res += BanyanArrays.ones((2048, 2048))

        res_sum_collect = collect(sum(res))
        @test res_sum_collect == 3.0 * 2048 * 2048
        res_maximum_collect = collect(maximum(res))
        @test res_maximum_collect == 3.0
    end

    # TODO: Re-enable this test once we ensure that we can write out small
    # enough datasets without unnecessary batching
    # run_with_job("String arrays") do job
    #     x1 = BanyanArrays.fill("hello\n", 2048)
    #     x2 = deepcopy(x1)
    #     x3 = BanyanArrays.fill("world\n", 2048)
    #     res = map(*, x1, x2, x3)
    #     res_lengths = map(s -> length(s), res)

    #     res_lengths_minimum_collect = collect(minimum(res_lengths))
    #     @test res_lengths_minimum_collect == 18

    #     # TODO: Support writing string arrays for this to work
    #     # This is unnecessary but will cache `res` on disk
    #     # write_to_disk(res)
    #     # res_collect = collect(res)
    #     # @test res_collect == BanyanArrays.fill("hello\nhello\nworld\n", 2048)

    #     # TODO: Test this once we implement a merging function for
    #     # variable-sized reductions
    #     # res_minimum_collect = collect(minimum(res))
    #     # @test res_minimum_collect == "hello\nhello\nworld\n"

    #     # TODO: Test this once we implement a merging function for
    #     # variable-sized reductions
    #     # x = BanyanArrays.fill("hi\n", 8)
    #     # res = reduce(*, x)

    #     # res_collect = collect(res)
    #     # @test res_collect == "hi\nhi\nhi\nhi\nhi\nhi\nhi\nhi\n"
    # end

    # TODO: Test HDF5 from URL and from S3
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
"""
MCMC iteration for pairwise alignment
"""
function simuiter(iter, p, q1, q2, L, tau, times, kappa, alpha, beta,
                  powera, dist, dist_min, best_match, match, thin, cut)

    match -= 1;
    best_match -= 1;

    Ixout = Array(Int64, 2*times);
    Ioriginal = Array(Int64, L+1);
    increment = 0.0; n_dist = 0.0; o_dist = 0.0; adjustcon = 0.0;
    ration = 0.0; prob = 0.0; u = 0.0; logpost = 0.0;
    increment_int = 0; newj = 0; tempnew = 0; tempold = 0; tempx = 0;
    newmatchvec = Array(Float64, 3);
    oldmatchvec = Array(Float64, 3);
    idenmatchvec = Array(Float64, 3);
    pnewvec = Array(Float64, 2);
    poldvec = Array(Float64, 2);
    internew = Array(Float64, 2*times);
    interold = Array(Float64, 2*times);
    qt1_5_interx = Array(Float64, 2*times);
    qt2_5_internew = Array(Float64, 2*times);
    qt2_5_interold = Array(Float64, 2*times);
    diff_ynew = Array(Float64, 2*times);
    diff_yold = Array(Float64, 2*times);
    original = Array(Float64, L+1);
    idy = Array(Float64, p);
    scalevec = Array(Float64, p);
    qt_5_fitted = Array(Float64, p);
    kappa_collect = Array(Float64, iter);
    log_collect = Array(Float64, iter);
    dist_collect = Array(Float64, iter);
    Irow = collect(1:p)-1;
    row = float(Irow);
    scale1 = Array(Float64, L);
    match_collect = Array(Float64, round(Integer,iter/thin), L+1);

    q2LLlen = (p-1)*times+1;
    q2LL = Array(Float64, q2LLlen);
    tempspan1 = collect(1:q2LLlen)-1;
    temp_span2 = float(tempspan1);
    timesf = float(times);
    q2LL_time = (temp_span2*(1/timesf));
    q2L_time1 = collect(1:p)-1;
    q2L_time2 = float(q2L_time1);
    q2LL = approx(q2L_time2, q2, q2LL_time);

    bar = Progress(iter, 1)   # minumum update interval: 1 second
    for j = 1:iter
        for i = 2:L
            if ((match[i+1]-match[i-1])>2)
                increment = rand(Normal(0, tau));
                increment_int = round(Integer,round(increment));
                if (increment_int == 0)
                    increment_int = (increment>0) ? 1 : -1;
                end
                newj = match[i] + increment_int;
                if (newj < match[i+1] && newj > match[i-1])
                    newmatchvec[1] = match[i-1];
                    newmatchvec[2] = newj;
                    newmatchvec[3] = match[i+1];
                    oldmatchvec[1] = match[i-1];
                    oldmatchvec[2] = match[i];
                    oldmatchvec[3] = match[i+1];
                    idenmatchvec[1] = times*(i-2);
                    idenmatchvec[2] = times*(i-1);
                    idenmatchvec[3] = times*i;
                    Ixout = collect(1:2*times)+times*(i-2)-1;
                    xout = float(Ixout);
                    internew = approx(idenmatchvec, newmatchvec, xout);
                    interold = approx(idenmatchvec, oldmatchvec, xout);
                    interx = copy(xout);
                    interynew = copy(internew);
                    push!(interynew, match[i+1]);
                    interyold = copy(interold);
                    push!(interyold, match[i+1]);
                    diff_ynew = interynew[2:2*times+1]-interynew[1:2*times];
                    diff_yold = interyold[2:2*times+1]-interyold[1:2*times];
                    for ll = 1:2*times
                        tempx = round(Integer, interx[ll]);
                        qt1_5_interx[ll] = q1[tempx+1];
                        internew[ll] = (internew[ll]>(p-1)) ? (p-1) : internew[ll];
                        tempnew = round(Integer, times*internew[ll]+1);
                        qt2_5_internew[ll] = q2LL[tempnew]*sqrt(diff_ynew[ll]);
                        interold[ll] = (interold[ll]>(p-1)) ? (p-1) : interold[ll];
                        tempold = round(Integer, times*interold[ll]+1);
                        qt2_5_interold[ll] = q2LL[tempold]*sqrt(diff_yold[ll]);
                    end
                    n_dist = norm(qt1_5_interx - qt2_5_internew)^2/p;
                    o_dist = norm(qt1_5_interx - qt2_5_interold)^2/p;
                    pnewvec = (newmatchvec[2:3] - newmatchvec[1:2])/p;
                    poldvec = (oldmatchvec[2:3] - oldmatchvec[1:2])/p;
                    adjustcon = exp((powera-1)*(sum(log(pnewvec))-
                                    sum(log(poldvec))));
                    ratio = adjustcon*exp(kappa*o_dist-kappa*n_dist);
                    prob = (ratio < 1) ? ratio : 1;
                    u = rand(Uniform(0,1));
                    match[i] = (u<prob) ? newj : match[i];
                end
            end
        end

        Ioriginal = (collect(1:L+1)-1)*times;
        original = float(Ioriginal);
        idy = round(approx(original, match, row));
        for ii = 1:L
            scale1[ii] = sqrt((match[ii+1]-match[ii])/times);
        end
        for kk = 1:p
            idy[kk] = (idy[kk]<p) ? idy[kk] : (p-1);
            scalevec[kk] = scale1[ceil(Integer,kk/times)];
            qt_5_fitted[kk] = scalevec[kk]*q2[round(Integer,idy[kk])+1];
        end
        dist = norm(q1-qt_5_fitted)^2/p;
        dist_collect[j] = dist;
        if (dist < dist_min)
            best_match = copy(match);
            dist_min = dist;
        end
        if (mod(j,thin)==0)
            match_collect[ceil(Integer, j/thin),:] = match;
        end
        kappa = rand(Gamma(p/2+alpha, 1/(dist+beta)));
        kappa_collect[j] = kappa;
        logpost = (p/2+alpha)*log(kappa)-kappa*(beta+dist);
        log_collect[j] = logpost;
        next!(bar);
    end
    match_collect += 1;
    best_match += 1;

    out = Dict("best_match" => best_match, "match_collect" => match_collect,
               "dist_collect" => dist_collect, "kappafamily" => kappa_collect,
               "log_posterior" => log_collect, "dist_min" => dist_min);
    return out
end

"""
MCMC iteration for group alignment
"""
function itermatch(iter, n, m, mu_5, match_matrix, qt_matrix,
                   qt_fitted_matrix, L, tau, times, kappa, alpha, beta,
                   powera, best_vec, dist_vec, best_match_matrix, mu_prior,
                   var_const, sumdist, thin, mu_q, mu_q_standard, logmax,
                   burnin, AVG)

    burnin -= 1;
    match_matrix -= 1;
    best_match_matrix -= 1;

    scale1 = Array(Float64, L);
    increment = 0.0; adjustcon = 0.0; ratio = 0.0; prob = 0.0; u = 0.0;
    n_dist = 0.0; o_dist = 0.0; dist = 0.0; logpost = 0.0; SigmaVar = 0.0;
    rescale = 0.0;
    increment_int = 0; newj = 0; tempnew = 0; tempold = 0; tempx = 0;
    newmatchvec = Array(Float64, 3);
    oldmatchvec = Array(Float64, 3);
    idenmatchvec = Array(Float64, 3);
    qt2_5 = Array(Float64, m);
    match = Array(Float64, L+1);
    interynew = Array(Float64, 2*times);
    interyold = Array(Float64, 2*times);
    interx = Array(Float64, 2*times);
    scalevec = Array(Float64, m);
    qt_5_fitted = Array(Float64, m);
    Ixout = Array(Int64, 2*times);
    Ioriginal = Array(Int64, L+1);
    xout = Array(Float64, 2*times);
    internew = Array(Float64, 2*times);
    interold = Array(Float64, 2*times);
    mu_5_interx = Array(Float64, 2*times);
    qt2_5_internew = Array(Float64, 2*times);
    qt2_5_interold = Array(Float64, 2*times);
    diff_ynew = Array(Float64, 2*times);
    diff_yold = Array(Float64, 2*times);
    pnewvec = Array(Float64, 2);
    poldvec = Array(Float64, 2);
    original = Array(Float64, L+1);
    idy = Array(Float64, m);
    Irow = collect(1:m)-1;
    row = float(Irow);
    karcher_res = Array(Float64, m+1);
    tempkarcher_res = Array(Float64, m+1);
    revscalevec = Array(Float64, L*times);
    kappa_collect = Array(Float64, iter);
    log_collect = Array(Float64, iter);
    MAP = Array(Float64, m);
    Mean = Array(Float64, m);
    qt2_5_warp = Array(Float64, m);
    mu_5_warp = Array(Float64, m);
    tempmatch_matrix = Array(Float64, L+1, n);
    bayes_warps = zeros(L+1, n);

    q2LLlen = (m-1)*times+1;
    q2LL_mat = Array(Float64, q2LLlen, n);
    tempspan1 = collect(1:q2LLlen)-1;
    temp_span2 = float(tempspan1);
    timesf = float(times);
    q2LL_time = (temp_span2*(1/timesf));
    q2L_time1 = collect(1:m)-1;
    q2L_time2 = float(q2L_time1);

    for LLL in 1:n
        q2LL_mat[:, LLL] = approx(q2L_time2, qt_matrix[:, LLL], q2LL_time);
    end

    bar = Progress(iter, 1);
    for j = 1:iter
        for t = 1:n
            qt2_5 = qt_matrix[:, t];
            match = match_matrix[:, t];
            for i = 2:L
                if ((match[i+1]-match[i-1])>2)
                    increment = rand(Normal(0, tau));
                    increment_int = round(Integer, round(increment));
                    if (increment_int == 0)
                        break
                        # increment_int = (increment>0)?(1):(-1);
                    end
                    newj = match[i] + increment_int;
                    if (newj < match[i+1] && newj > match[i-1])
                        newmatchvec[1] = match[i-1];
                        newmatchvec[2] = newj;
                        newmatchvec[3] = match[i+1];
                        oldmatchvec[1] = match[i-1];
                        oldmatchvec[2] = match[i];
                        oldmatchvec[3] = match[i+1];
                        idenmatchvec[1] = times*(i-2);
                        idenmatchvec[2] = times*(i-1);
                        idenmatchvec[3] = times*i;
                        Ixout = collect(1:2*times)+times*(i-2)-1;
                        xout = float(Ixout);
                        internew = approx(idenmatchvec, newmatchvec, xout);
                        interold = approx(idenmatchvec, oldmatchvec, xout);
                        interx = copy(xout);
                        interynew = copy(internew);
                        push!(interynew, match[i+1]);
                        interyold = copy(interold);
                        push!(interyold, match[i+1]);
                        diff_ynew = interynew[2:(2*times+1)]-interynew[1:2*times];
                        diff_yold = interyold[2:(2*times+1)]-interyold[1:2*times];
                        for ll = 1:2*times
                            tempx = round(Integer, interx[ll]);
                            mu_5_interx[ll] = mu_5[tempx+1];
                            internew[ll] = (internew[ll]>(m-1)) ? (m-1) : internew[ll];
                            tempnew = round(Integer, times*internew[ll]+1);
                            tempnew = (tempnew > q2LLlen) ? q2LLlen : tempnew;
                            qt2_5_internew[ll] = q2LL_mat[tempnew, t] *
                                sqrt(diff_ynew[ll]);
                            interold[ll] = (interold[ll]>(m-1)) ? (m-1) : interold[ll];
                            tempold = round(Integer, times*interold[ll]+1);
                            tempold = (tempold > q2LLlen) ? q2LLlen : tempold;
                            qt2_5_interold[ll] = q2LL_mat[tempold, t] *
                                sqrt(diff_yold[ll]);
                        end
                        n_dist = norm(mu_5_interx - qt2_5_internew)^2/m;
                        o_dist = norm(mu_5_interx - qt2_5_interold)^2/m;
                        pnewvec = (newmatchvec[2:3] - newmatchvec[1:2])/m;
                        poldvec = (oldmatchvec[2:3] - oldmatchvec[1:2])/m;
                        adjustcon = exp((powera-1)*(sum(log(pnewvec))-
                                        sum(log(poldvec))));
                        ratio = adjustcon*exp(kappa*o_dist-kappa*n_dist);
                        prob = (ratio < 1) ? ratio : 1;
                        u = rand(Uniform(0,1));
                        match[i] = (u<prob) ? newj : match[i];
                    end
                end
            end
            Ioriginal = (collect(1:L+1)-1)*times;
            original = float(Ioriginal);
            idy = approx(original, match, row);
            qt2_5_warp = approx(q2L_time2, qt2_5, idy);

            for ii = 1:L
                scale1[ii] = sqrt((match[ii+1]-match[ii])/times);
            end

            for kk = 1:m
                scalevec[kk] = scale1[ceil(Integer, kk/times)];
                qt_5_fitted[kk] = scalevec[kk]*qt2_5_warp[kk];
            end

            dist = norm(mu_5-qt_5_fitted)^2/m;
            match_matrix[:, t] = match;
            dist_vec[t] = dist;
            qt_fitted_matrix[:, t] = qt_5_fitted;
            if ((j/thin)>=burnin)
                bayes_warps[:, t] = bayes_warps[:, t] + match./AVG;
            end
        end
        sumdist[j] = sum(dist_vec);
        if (sumdist[j] < sum(best_vec))
            best_vec = copy(dist_vec);
            best_match_matrix = copy(match_matrix);
        end
        if (mod(j,thin)==0)
            mu_q[ceil(Integer, j/thin),:] = mu_5;
            karcher_res, revscalevec = findkarcherinv(match_matrix+1, times);
            mu_5_warp = approx(q2L_time2, vec(mu_5), karcher_res-1);
            for qq in 1:m
                mu_q_standard[ceil(Integer, j/thin), qq] = revscalevec[qq]*mu_5_warp[qq];
            end
        end
        kappa = rand(Gamma(n*m/2+alpha, 1/(sum(dist_vec)+beta)));
        kappa_collect[j] = kappa;
        logpost = (n*m/2+alpha)*log(kappa)-kappa*(beta+sum(dist_vec));
        log_collect[j] = logpost;
        if (logpost > logmax)
            logmax = logpost;
            MAP = copy(mu_5);
        end
        Mean = (var_const*m)/(m+2*kappa*n*var_const)*(2*kappa*n*mean(qt_fitted_matrix,2)/m + mu_prior/var_const);
        SigmaVar = (var_const*m)/(m+2*kappa*n*var_const);
        mu_5 = Mean + rand(Normal(0, sqrt(SigmaVar)), m);
        rescale = sqrt(m/sum(mu_5.^2));
        mu_5 *= rescale;
        next!(bar);
    end

    best_match_matrix += 1;
    bayes_warps += 1;

    out = Dict("mu_q_standard" => mu_q_standard, "mu_q" => mu_q,
               "best_match_matrix" => best_match_matrix,
               "sumdist" => sumdist, "kappafamily" => kappa_collect,
               "log_posterior" => log_collect, "dist" => sum(best_vec),
               "MAP" => MAP, "bayes_warps" => bayes_warps);
    return out
end


text:
using KnapsackLib
using GLPK
using Printf

# runs = 10000
# sizes = [ 5 ]
# max_vals = [ 10 ]
runs = 10
sizes = [ 10, 100, 500, 1000, 2000, 4000 ]
max_vals = [ 10, 100, 500, 1000, 2000 ]

times = []

for v in max_vals
    sizes_times = []
    for n in sizes
        algos_times = zeros(Float64, 4)
        for r = 1:runs
            data = genKnap(n, v; seed = r)

            ret1 = @timed solveKnapNaive(data)
            ret2 = @timed solveKnapModel(data, GLPK.Optimizer)
            ret3 = @timed solveKnapExpCore(data)
            ret4 = @timed solveKnapHeur(data)

            val1, sol1 = ret1.value
            val2, sol2 = ret2.value
            val3, sol3 = ret3.value

            algos_times[1] += ret1.time / runs
            algos_times[2] += ret2.time / runs
            algos_times[3] += ret3.time / runs
            algos_times[4] += ret4.time / runs

            w1 = sum(data.items[j].weight for j in sol1)
            w2 = sum(data.items[j].weight for j in sol2)
            w3 = sum(data.items[j].weight for j in sol3)

            if val1 != val2 || val2 != val3
                println(val1, " (", w1, ") x ", val2, " (", w2, ") x ", val3, " (", w3, ")")

                println(items)
                println(capacity)

                println(sol1)
                println(sol2)
                println(sol3)

                @assert 1 != 1
            end
        end
        push!(sizes_times, algos_times)
    end
    push!(times, sizes_times)
end

@printf "%10s " ""
for n in 1:length(sizes)
    @printf "%10d " sizes[n]
end
@printf "\n"

for v in 1:length(max_vals)
    for i in 1:4
        if i == 2
            @printf "%10d " max_vals[v]
        else
            @printf "%10s " ""
        end

        for n in 1:length(sizes)
            @printf "%10.7lf " times[v][n][i]
        end
        @printf "\n"
    end
end