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

using Markdown
using InteractiveUtils

# ╔═╡ f42eb779-f2f2-43d7-996f-97ba181399f4
begin
	using Downloads, DataFrames, CSV, TimeZones, Dates, Statistics
	using HypertextLiteral, JSONTables
end

# ╔═╡ 8b3681b9-8a63-4168-a1f1-04444b24a8bf
md"# 🧗♀ O’BLOC
This notebook is used to analyse scrape data of the [O'BLOC Website](https://obloc.ch).

O'BLOC is a excellent climbing gym! If you are in Bern consider visiting them 😊"

# ╔═╡ 05a8266a-e30e-4e6f-bbd5-3d5c8114fca6
md"## DataFrame"

# ╔═╡ 5fa2251f-9aec-424b-b0fe-48ad15d74b7a
url = "https://raw.githubusercontent.com/ioboi/obloc-data/main/visitors.csv"

# ╔═╡ cc4c8310-5004-42fe-9cbb-0f33e7b034dd
Markdown.parse("This is the [source]($(url)) of the data.")

# ╔═╡ eb448f9c-add8-4c8f-b850-0bac4881c6f2
begin
	df = DataFrame(CSV.File(Downloads.download(url)))

	# Convert column datetime to "ZonedDateTime" and change timezone
	df.datetime = ZonedDateTime.(df.datetime, dateformat"yyyy-mm-ddTHH:MM:SS.sssz")
	df.datetime = astimezone.(df.datetime, tz"UTC+1")

	# Add columns
	df.date = Date.(df.datetime)
	df.dayofweek = Dates.dayofweek.(df.datetime)
	df.hour = Dates.hour.(df.datetime)

	# Filter out NaN entries
	df = df[.!isnan.(df.value), :]

	hourBetween(startHour::Int, endHour::Int)::BitVector =
		df.hour .>= startHour .&& df.hour .<= endHour

	# Day filters
	mo = Dates.dayofweek.(df.datetime) .== 1
	tue = Dates.dayofweek.(df.datetime) .== 2
	wed = Dates.dayofweek.(df.datetime) .== 3
	thu = Dates.dayofweek.(df.datetime) .== 4
	fri = Dates.dayofweek.(df.datetime) .== 5
	sat = Dates.dayofweek.(df.datetime) .== 6
	sun = Dates.dayofweek.(df.datetime) .== 7

	# Actual schedules
	monTueThu = (mo .|| tue .|| thu) .&& hourBetween(10, 23)
	wedFri = (wed .|| fri) .&& hourBetween(8, 23)
	satSun = (sat .|| sun) .&& hourBetween(9, 19)

	# TODO: Respect holidays in schedule filter

	# Filter schedules
	df = df[monTueThu .|| wedFri .|| satSun, :]
end

# ╔═╡ 2f412aa6-ab67-4ff8-8bf1-5282545ea012
md"## 📊 Overview"

# ╔═╡ 24221234-bc42-4da0-a5bd-146aaa8791b8
md"### Mean Occupancy throughout the Week"

# ╔═╡ fddf3b84-98e3-4571-bf47-611db0003679
begin
	dfMeanWeeks = combine(groupby(df, :dayofweek), :value => mean)
	dfMeanWeeks.dayofweek = Dates.dayname.(dfMeanWeeks.dayofweek)
	dfMeanWeeks.value_mean = dfMeanWeeks.value_mean * 100
	rename!(dfMeanWeeks, [:dayofweek, :value_mean] .=>  ["Day of Week", "ø %Occupancy"])
end

# ╔═╡ d3631006-c877-4241-8fbb-8b6c3ba0df9d
md"### What is the best time to go climbing?"

# ╔═╡ a6666bdf-8096-4895-a46d-2b68986482fe
begin
	groupedByHour = combine(groupby(df, :hour), :value => mean, renamecols=false)
	function generateBestTimeBarChart()
		data = arraytable(groupedByHour)
		@htl("""
		<script>
			const d3 = await import("https://cdn.skypack.dev/d3@7");
			const width = 600;
			const height = 400;
			const m = {left: 40, top: 40, right: 40, bottom: 40};

			const graphWidth = width - m.left - m.right;
			const graphHeight = height - m.top - m.bottom;

			const data = JSON.parse($(data));

			const svg = DOM.svg(width, height);

			const xScale = d3.scaleBand()
				.range([0, graphWidth])
				.domain(data.map(d => d.hour))
				.padding(0.05);

			const yScale = d3.scaleLinear()
				.range([graphHeight, 0])
				.domain([0, 1])

			d3.select(svg)
				.append("g")
				.attr("transform", `translate(\${m.left},\${m.top})`)
  				.call(d3.axisLeft(yScale));

			d3.select(svg)
				.append("g")
				.attr("transform", `translate(\${m.left},\${graphHeight+m.top})`)
  				.call(d3.axisBottom(xScale));

			d3.select(svg)
				.append("g")
				.attr("transform", `translate(\${m.left},\${m.top})`)
				.selectAll("rect")
				.data(data, (d) => d)
				.join("rect")
	      		.attr("x", (d) => xScale(d.hour))
	      		.attr("y", (d) => yScale(d.value))
	      		.attr("width", xScale.bandwidth())
	      		.attr("height", (d) => graphHeight - yScale(d.value))
				.attr("fill", "#006446");

			return svg;
		</script>
		""")
	end
	generateBestTimeBarChart()
end

# ╔═╡ 52813bc0-29b9-4494-a9a2-6437da451897
begin
	function topHour()
		x = sort(groupedByHour, [:value])
		x.value = x.value * 100
		x
	end
	topHour()
end

# ╔═╡ 7b498d91-63ed-4a70-b4d4-bdcc21ceacd4
md"### Heatmap Year
This heatmap shows the mean occupancy per calendar week and day.

Sadly I started scraping in calendar week 5😢
"

# ╔═╡ 06df7a66-81d9-4673-9902-b873e3a4fe2b
begin
	# Here be Dragons!
	function generateHeatmapYear()
		data = combine(groupby(df, :date), :value=>mean, renamecols=false)
		data.week = Dates.week.(data.date)
		data.day = Dates.dayofweek.(data.date)
		data = arraytable(data)
		@htl("""
		<script>
			const d3 = await import("https://cdn.skypack.dev/d3@7");
			const margins = {left: 40, top: 0, right: 20, bottom: 140};
			const data = JSON.parse($(data));

			const weeks = d3.range(d3.min(data.map(d=>d.week)), d3.max(data.map(d=>d.week))+1);
			const days = d3.range(1, 8);

			const DAYNAMES = ["Mo", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]

			const width = 600;
			const height = 400;
			const svg = DOM.svg(width, height);

			const xScale = d3.scaleBand()
	  			.range([ 0, width - margins.left - margins.right ])
				.domain(weeks)
				.paddingInner(0.05);

			const yScale = d3.scaleBand()
				.range([ 0, height - margins.top - margins.bottom ])
	  			.domain(days)
				.paddingInner(0.04);

			const dayScale = d3.scaleBand()
				.range([ 0, height - margins.top - margins.bottom ])
	  			.domain(DAYNAMES)
				.paddingInner(0.04);

			const color = d3.scaleLinear()
	    	.domain([0, d3.max(data.map(d => d.value))])
	    	.range(["#e5f5e0", "#31a354"]);

			d3.select(svg)
	  			.append("g")
				.style("font-size", 14)
				.attr("transform", `translate(\${margins.left-0.04}, 0)`)
	  			.call(d3.axisLeft(dayScale).tickSize(0))
				.select(".domain").remove();

			d3.select(svg)
	  			.append("g")
				.style("font-size", 14)
				.attr("transform", `translate(\${margins.left},\${height-margins.bottom})`)
	  			.call(d3.axisBottom(xScale).tickSize(0))
				.select(".domain").remove();

			d3.select(svg)
				.append("g")
				.attr("transform", `translate(\${margins.left},\${margins.top})`)
				.selectAll("rect")
				.data(data, (d) => d)
				.join("rect")
	      		.attr("x", (d) => xScale(d.week))
	      		.attr("y", (d) => yScale(d.day))
				.attr("rx", 4)
	      		.attr("ry", 4)
	      		.attr("width", xScale.bandwidth())
	      		.attr("height", yScale.bandwidth())
				.attr("fill", (d) => color(d.value));

			const gradient = DOM.uid();

			const grad = d3.select(svg)
				.append("linearGradient")
      			.attr("id", gradient.id)
      			.attr("x1", "0%")
    			.attr("y1", "0%")
    			.attr("x2", "100%")
    			.attr("y2", "0%");

			grad.append("stop")
					.attr("offset", "0%")
					.attr("stop-color", color(0));

			grad.append("stop")
					.attr("offset", "100%")
					.attr("stop-color", color(d3.max(data.map(d => d.value))));

			const legend = d3.select(svg)
				.append("g")
				.attr("transform", `translate(\${margins.left},\${height-margins.bottom + 1.5*yScale.bandwidth()})`)
				.style("font-size", 14);

			legend.append("text")
				.attr("x", "0")
				.attr("alignment-baseline", "middle")
				.text("Less");

			legend.append("rect")
				.attr("x", "30")
				.attr("y", "-10")
  				.attr("width", 100)
  				.attr("height", 20)
				.attr("fill", gradient);

			legend.append("text")
				.attr("x", "137")
				.attr("alignment-baseline", "middle")
				.text("More");

			return svg;
		</script>
		""")
	end
	generateHeatmapYear()
end

# ╔═╡ ae48bc41-89fc-4f71-bd91-2e93c9c673eb
begin
	# This cell is mainly used to display when the static html was generated.
	# It has no meaning when this notebook is used interactively. 😂
	generatedtimestap = now(tz"UTC+1")
	md"*Generated at $(generatedtimestap)*"
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
JSONTables = "b9914132-a727-11e9-1322-f18e41205b0b"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"

[compat]
CSV = "~0.10.3"
DataFrames = "~1.3.2"
HypertextLiteral = "~0.9.3"
JSONTables = "~1.0.3"
TimeZones = "~1.7.2"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

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

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

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

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "9310d9495c1eb2e4fa1955dd478660e2ecab1fbb"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.3"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

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

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

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

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

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

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

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

[[deps.ExprTools]]
git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.8"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "04d13bfa8ef11720c24e4d840c0033d145537df7"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.17"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

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

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

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

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "8c1f668b24d999fb47baf80436194fdccec65ad2"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.4"

[[deps.JSONTables]]
deps = ["JSON3", "StructTypes", "Tables"]
git-tree-sha1 = "13f7485bb0b4438bb5e83e62fcadc65c5de1d1bb"
uuid = "b9914132-a727-11e9-1322-f18e41205b0b"
version = "1.0.3"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

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

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

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

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

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

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

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

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

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

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

[[deps.Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

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

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

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

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

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

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

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

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

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

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

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

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

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

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

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

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

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

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

[[deps.StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

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

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

[[deps.TimeZones]]
deps = ["Dates", "Downloads", "InlineStrings", "LazyArtifacts", "Mocking", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "2d4b6de8676b34525ac518de36006dc2e89c7e2e"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.7.2"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

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

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

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

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

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

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

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

# ╔═╡ Cell order:
# ╟─f42eb779-f2f2-43d7-996f-97ba181399f4
# ╟─8b3681b9-8a63-4168-a1f1-04444b24a8bf
# ╟─05a8266a-e30e-4e6f-bbd5-3d5c8114fca6
# ╟─5fa2251f-9aec-424b-b0fe-48ad15d74b7a
# ╟─cc4c8310-5004-42fe-9cbb-0f33e7b034dd
# ╟─eb448f9c-add8-4c8f-b850-0bac4881c6f2
# ╟─2f412aa6-ab67-4ff8-8bf1-5282545ea012
# ╟─24221234-bc42-4da0-a5bd-146aaa8791b8
# ╟─fddf3b84-98e3-4571-bf47-611db0003679
# ╟─d3631006-c877-4241-8fbb-8b6c3ba0df9d
# ╟─a6666bdf-8096-4895-a46d-2b68986482fe
# ╟─52813bc0-29b9-4494-a9a2-6437da451897
# ╟─7b498d91-63ed-4a70-b4d4-bdcc21ceacd4
# ╟─06df7a66-81d9-4673-9902-b873e3a4fe2b
# ╟─ae48bc41-89fc-4f71-bd91-2e93c9c673eb
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: ioboi/obloc

max_issues_repo_name: ioboi/obloc

max_forks_repo_name: ioboi/obloc


text:
using SafeTestsets

@safetestset "SDDEProblem, solve" begin
    include("test_prob_sol.jl")
end

ext: jl


text:
"""
    AbstractAxis

An `AbstractVector` subtype optimized for indexing.
"""
abstract type AbstractAxis{I<:Integer,Inds} <: AbstractUnitRange{I} end

Base.keytype(::Type{T}) where {I,T<:AbstractAxis{I}} = I
Base.keys(axis::AbstractAxis) = eachindex(axis)

"""
    AbstractOffsetAxis{I,Inds}

Supertype for axes that begin indexing offset from one. All subtypes of `AbstractOffsetAxis`
use the keys for indexing and only convert to the underlying indices when
`to_index(::OffsetAxis, ::Integer)` is called (i.e. when indexing the an array
with an `AbstractOffsetAxis`. See [`OffsetAxis`](@ref), [`CenteredAxis`](@ref),
and [`IdentityAxis`](@ref) for more details and examples.
"""
abstract type AbstractOffsetAxis{I,Inds,F} <: AbstractAxis{I,Inds} end

"""
    IndexAxis

Index style for mapping keys to an array's parent indices.
"""
struct IndexAxis <: IndexStyle end

Base.valtype(::Type{T}) where {I,T<:AbstractAxis{I}} = I

Base.IndexStyle(::Type{T}) where {T<:AbstractAxis} = IndexAxis()

Base.parent(axis::AbstractAxis) = getfield(axis, :parent)

ArrayInterface.parent_type(::Type{T}) where {I,Inds,T<:AbstractAxis{I,Inds}} = Inds

Base.eachindex(axis::AbstractAxis) = static_first(axis):static_last(axis)

for f in (:(==), :isequal)
    @eval begin
        Base.$(f)(x::AbstractAxis, y::AbstractAxis) = $f(eachindex(x), eachindex(y))
        Base.$(f)(x::AbstractArray, y::AbstractAxis) = $f(x, eachindex(y))
        Base.$(f)(x::AbstractAxis, y::AbstractArray) = $f(eachindex(x), y)
        Base.$(f)(x::AbstractRange, y::AbstractAxis) = $f(x, eachindex(y))
        Base.$(f)(x::AbstractAxis, y::AbstractRange) = $f(eachindex(x), y)
        Base.$(f)(x::StaticRanges.GapRange, y::AbstractAxis) = $f(x, eachindex(y))
        Base.$(f)(x::AbstractAxis, y::StaticRanges.GapRange) = $f(eachindex(x), y)
        Base.$(f)(x::OrdinalRange, y::AbstractAxis) = $f(x, eachindex(y))
        Base.$(f)(x::AbstractAxis, y::OrdinalRange) = $f(eachindex(x), y)
    end
end

Base.allunique(a::AbstractAxis) = true

Base.empty!(axis::AbstractAxis) = set_length!(axis, 0)

@inline Base.in(x::Integer, axis::AbstractAxis) = !(x < first(axis) || x > last(axis))
@inline Base.length(axis::AbstractAxis) = length(parent(axis))

Base.pairs(axis::AbstractAxis) = Base.Iterators.Pairs(a, keys(axis))

# This is required for performing `similar` on arrays
Base.to_shape(axis::AbstractAxis) = length(axis)

Base.haskey(axis::AbstractAxis, key) = key in keys(axis)

@inline Base.axes(axis::AbstractAxis) = (Base.axes1(axis),)

@inline Base.axes1(axis::AbstractAxis) = copy(axis)

@inline Base.unsafe_indices(axis::AbstractAxis) = (axis,)

Base.isempty(axis::AbstractAxis) = isempty(parent(axis))
Base.empty(axis::AbstractAxis) = unsafe_reconstruct(axis, _empty(parent(axis)))
_empty(axis::AbstractAxis) = empty(axis)
_empty(axis::AbstractUnitRange) = One():Zero()

Base.sum(axis::AbstractAxis) = sum(eachindex(axis))

function ArrayInterface.can_change_size(::Type{T}) where {T<:AbstractAxis}
    return can_change_size(parent_type(T))
end

Base.collect(a::AbstractAxis) = collect(eachindex(a))

Base.step(axis::AbstractAxis) = oneunit(eltype(axis))

Base.step_hp(axis::AbstractAxis) = 1

Base.size(axis::AbstractAxis) = (length(axis),)

Base.:-(axis::AbstractAxis) = maybe_unsafe_reconstruct(axis, -eachindex(axis))

function Base.:+(r::AbstractAxis, s::AbstractAxis)
    indsr = axes(r, 1)
    if indsr == axes(s, 1)
        data = eachindex(r) + eachindex(s)
        axs = (unsafe_reconstruct(r, eachindex(data)),)
        return  AxisArray{eltype(data),ndims(data),typeof(data),typeof(axs)}(data, axs)
    else
        throw(DimensionMismatch("axes $indsr and $(axes(s, 1)) do not match"))
    end
end
function Base.:-(r::AbstractAxis, s::AbstractAxis)
    indsr = axes(r, 1)
    if indsr == axes(s, 1)
        data = eachindex(r) - eachindex(s)
        axs = (unsafe_reconstruct(r, eachindex(data); keys=keys(data)),)
        return  AxisArray{eltype(data),ndims(data),typeof(data),typeof(axs)}(data, axs)
    else
        throw(DimensionMismatch("axes $indsr and $(axes(s, 1)) do not match"))
    end
end

Base.UnitRange(axis::AbstractAxis) = UnitRange(eachindex(axis))
function Base.AbstractUnitRange{T}(axis::AbstractAxis) where {T<:Integer}
    if eltype(axis) <: T && !can_change_size(axis)
        return axis
    else
        return unsafe_reconstruct(axis, AbstractUnitRange{T}(parent(axis)); keys=keys(axis))
    end
end

Base.pop!(axis::AbstractAxis) = pop!(parent(axis))
Base.popfirst!(axis::AbstractAxis) = popfirst!(parent(axis))

# TODO check for existing key first
function push_key!(axis::AbstractAxis, key)
    grow_last!(parent(axis), 1)
    return nothing
end

function pushfirst_axis!(axis::AbstractAxis, key)
    grow_last!(parent(axis), 1)
    return nothing
end

function popfirst_axis!(axis::AbstractAxis)
    shrink_last!(parent(axis), 1)
    return nothing
end

Base.lastindex(a::AbstractAxis) = last(a)
Base.last(axis::AbstractAxis) = last(parent(axis))

ArrayInterface.known_last(::Type{T}) where {T<:AbstractAxis} = known_last(parent_type(T))

###
### first
###
Base.firstindex(axis::AbstractAxis) = first(axis)
Base.first(axis::AbstractAxis) = first(parent(axis))

ArrayInterface.known_first(::Type{T}) where {T<:AbstractAxis} = known_first(parent_type(T))

# This is different than how most of Julia does a summary, but it also makes errors
# infinitely easier to read when wrapping things at multiple levels or using Unitful keys
Base.summary(io::IO, axis::AbstractAxis) = show(io, axis)

function reverse_keys(axis::AbstractAxis, newinds::AbstractUnitRange)
    return Axis(reverse(keys(axis)), newinds; checks=NoChecks)
end

###
### maybe_unsafe_reconstruct
###
maybe_unsafe_reconstruct(axis, inds::Integer; kwargs...) = eltype(axis)(inds)
function maybe_unsafe_reconstruct(axis, inds::AbstractArray; keys=nothing)
    if known_step(inds) === 1 && eltype(inds) <: Integer
        return unsafe_reconstruct(axis, SimpleAxis(inds); keys=keys)
    else
        axs = (unsafe_reconstruct(axis, SimpleAxis(eachindex(inds))),)
        return AxisArray{eltype(inds),ndims(inds),typeof(inds),typeof(axs)}(inds, axs)
    end
end
function maybe_unsafe_reconstruct(axis::AbstractOffsetAxis, inds::AbstractArray; keys=nothing)
    if known_step(inds) === 1 && eltype(inds) <: Integer
        return unsafe_reconstruct(axis, SimpleAxis(inds); keys=keys)
    else
        axs = (unsafe_reconstruct(axis, SimpleAxis(eachindex(inds))),)
        return AxisArray{eltype(inds),ndims(inds),typeof(inds),typeof(axs)}(inds, axs)
    end
end

###
### append
###
# TODO document append_keys!
#= append_keys!(x, y) =#
append_keys!(x::AbstractRange, y) = set_length!(x, length(x) + length(y))
function append_keys!(x, y)
    if eltype(x) <: eltype(y)
        for x_i in x
            if x_i in y
                error("Element $x_i appears in both collections in call to append_axis!(collection1, collection2). All elements must be unique.")
            end
        end
        return append!(x, y)
    else
        return append_axis!(x, promote_axis_collections(y, x))
    end
end

@inline ArrayInterface.offsets(axis::AbstractAxis) = (first(axis),)

Base.show(io::IO, axis::AbstractAxis) = _show(io, axis)
Base.show(io::IO, ::MIME"text/plain", axis::AbstractAxis) = _show(io, axis)
function _show(io::IO, axis::AbstractAxis)
    if haskey(io, :compact)
        ks = keys(axis)
        if known_step(ks) === 1
            # this prevents StaticInt from creating long prints to REPL
            print(io, "$(Int(first(ks))):$(Int(last(ks)))")
        else
            print(io, ks)
        end
    else
        print_axis(io, axis)
    end
end


text:
function getDataset(dataset_name::String)
    pathToScripts = joinpath(Pkg.dir("EconDatasets"), "src/getDataset/")

    ## resolve key names and required package
    fileDict = ["FFF" => :(EconDatasets.getFFF()),
                "UMD" => :(EconDatasets.getUMD()),
                "SP500" => :(include(string($pathToScripts,
                                            "getSP500.jl"))),
                "SP500Ticker" => :(include(string($pathToScripts,
                                                  "getSP500TickerSymbols.jl"))),
                "SP500Industries" => :(include(string($pathToScripts,
                                                      "getSP500TickerSymbols.jl"))),
                "Indices" => :(EconDatasets.getIndices()),
                "Treasuries" => :(EconDatasets.getTreasuries()),
                "DieboldLi" => :(EconDatasets.getDieboldLi())
                ]

    downloadFunc = fileDict[dataset_name]
    eval(downloadFunc)
end


text:
"""
	write_mccs(outfile::Function, MCCs::AbstractDict[, segments])
	write_mccs(
		outfile::AbstractString,
		MCCs::AbstractDict,
		segments = get_segments(MCCs),
	)
	write_mccs(of, MCCs::AbstractArray)

Write MCCs to file.
"""
function write_mccs(
	outfile::AbstractString,
	MCCs::AbstractDict,
	mode="w";
	segments = get_segments(MCCs),
)
	write_mccs((s1,s2) -> outfile, MCCs, mode; segments)
end
function write_mccs(
	outfile::Function,
	MCCs::AbstractDict,
	mode = "w";
	segments = get_segments(MCCs),
)
	for i in 1:length(segments), j in (i+1):length(segments)
		write_mccs(outfile(segments[i], segments[j]), MCCs[segments[i], segments[j]], mode)
	end

	return nothing
end

function write_mccs(of, MCCs::AbstractArray, mode="w")
	open(of, mode) do w
		for (i,m) in enumerate(MCCs)
			for x in m[1:end-1]
				write(w, x * ",")
			end
			if i < length(MCCs)
				write(w, m[end] * "\n")
			else
				write(w, m[end])
			end
		end
	end

	return nothing
end

function get_segments(MCCs)
	segments = []
	for k in keys(MCCs)
		push!(segments, k...)
	end

	return unique(sort(segments))
end

function read_mccs(file)
	map(eachline(file)) do m
       String.(split(m, ','))
    end
end

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
using IntervalEigenvalues

import Documenter
Documenter.makedocs(
    modules = [IntervalEigenvalues],
    format = :html,
    sitename = "IntervalEigenvalues.jl",
    root = joinpath(dirname(dirname(@__FILE__)), "docs"),
    pages = Any["Home" => "index.md"],
    strict = true,
    linkcheck = true,
    checkdocs = :exports,
    authors = "David Sanders"
)

using Base.Test

# write your own tests here
@test 1 == 2


text:
using MOSLab
using Plots
using LaTeXStrings
Vg = 0.1 # Gate Voltage
Vd = 1.0 # Dain Voltage
Nx = 101 # Number of grid points on the x direction
Ny = 101 # Number of grid points on the y direction
N_d = 6e16 # Dopants density at the drain and source side  in cm ^-3
T = 300.0 # Temperature in Kelvin
Nb = 5.7e17 # Dopants ensity at the bulk in cm ^-3
l = .180e-4 # Channel length
h = 0.5e-4 # Silicon depth
t_ox = 4e-7 # Oxide thickens in cm
Ms = MOSFETInputDeck(N_d,Nb,T,l,h,h/2.0,l,t_ox,:NMOS) # Create the input deck for an NMOS transistor
Vs = ψs_PSP(Vg,Ms.Gate) # Calculate the surface potential at the gate
G,psi,nn,pp = MOSFETSimulation(Ms,Vs,Vd,Nx,Ny;verbose=false,ξ₀=0.01) # run 2D MOSFET simulation
gr()
X = unique(G[1,:])
Y = unique(G[2,:])
contour(X,Y,psi,fill=true)
surface(G[1,:],G[2,:],psi,fill=true)
xlabel!(L"x \ [\mu m]")
ylabel!(L"y \ [\mu m]")
plot!(zlabel = L"\psi \ [V]")
savefig("Psi_Surface.png")
yvals = unique(G[2,:])
ySel = G[2,:] .== yvals[end-1]
sum(ySel)
minimum(psi[ySel])
plot(G[1,ySel],psi[ySel],label="\$ E_C (x,$(yvals[end-1]) ) \$",legend=:bottomleft)
savefig("Psi_channel.png")
plot(G[1,ySel],nn[ySel],label="\$ E_C (x,$(yvals[end-1]) ) \$",legend=:topright)
savefig("nn_channel.png")


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const SDL_IMAGE_MAJOR_VERSION = 2
const SDL_IMAGE_MINOR_VERSION = 0
const SDL_IMAGE_PATCHLEVEL = 4

# Skipping MacroDefinition: SDL_IMAGE_VERSION ( X ) \
#{ ( X ) -> major = SDL_IMAGE_MAJOR_VERSION ; ( X ) -> minor = SDL_IMAGE_MINOR_VERSION ; ( X ) -> patch = SDL_IMAGE_PATCHLEVEL ; \
#}
# Skipping MacroDefinition: SDL_IMAGE_COMPILEDVERSION SDL_VERSIONNUM ( SDL_IMAGE_MAJOR_VERSION , SDL_IMAGE_MINOR_VERSION , SDL_IMAGE_PATCHLEVEL )
# Skipping MacroDefinition: SDL_IMAGE_VERSION_ATLEAST ( X , Y , Z ) ( SDL_IMAGE_COMPILEDVERSION >= SDL_VERSIONNUM ( X , Y , Z ) )

#const IMG_SetError = SDL_SetError
#const IMG_GetError = SDL_GetError

# begin enum ANONYMOUS_1
const ANONYMOUS_1 = UInt32
const SDL_FALSE = (UInt32)(0)
const SDL_TRUE = (UInt32)(1)
const SDL_bool = Bool

# end enum ANONYMOUS_1

# begin enum ANONYMOUS_49
const ANONYMOUS_49 = UInt32
const IMG_INIT_JPG = (UInt32)(1)
const IMG_INIT_PNG = (UInt32)(2)
const IMG_INIT_TIF = (UInt32)(4)
const IMG_INIT_WEBP = (UInt32)(8)
# end enum ANONYMOUS_49

#$(Expr(:typealias, :IMG_InitFlags, :Void))

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
using LinearAlgebra
using Random
using StatsFuns: softplus
using Roots # for inverse

################################################################################
#                            Planar and Radial Flows                           #
#             Ref: Variational Inference with Normalizing Flows,               #
#               D. Rezende, S. Mohamed(2015) arXiv:1505.05770                  #
################################################################################

###############
# RadialLayer #
###############

# FIXME: using `TrackedArray` for the parameters, we end up with
# nested tracked structures; don't want this.
mutable struct RadialLayer{T1,T2} <: Bijector{1}
    α_::T1
    β::T1
    z_0::T2
end

function RadialLayer(dims::Int, container=Array)
    α_ = container(randn(1))
    β = container(randn(1))
    z_0 = container(randn(dims, 1))
    return RadialLayer(α_, β, z_0)
end

h(α, r) = 1 ./ (α .+ r)     # for radial flow from eq(14)
dh(α, r) = - h(α, r) .^ 2   # for radial flow; derivative of h()

# An internal version of transform that returns intermediate variables
function _transform(flow::RadialLayer, z)
    α = softplus(flow.α_[1])            # from A.2
    β_hat = -α + softplus(flow.β[1])    # from A.2
    r = sqrt.(sum((z .- flow.z_0).^2; dims = 1))
    transformed = z + β_hat .* h(α, r) .* (z .- flow.z_0)   # from eq(14)
    return (transformed=transformed, α=α, β_hat=β_hat, r=r)
end

(b::RadialLayer)(z::AbstractMatrix{<:Real}) = _transform(b, z).transformed
(b::RadialLayer)(z::AbstractVector{<:Real}) = vec(_transform(b, z).transformed)

function _forward(flow::RadialLayer, z)
    transformed, α, β_hat, r = _transform(flow, z)
    # Compute log_det_jacobian
    d = size(flow.z_0, 1)
    h_ = h(α, r)
    log_det_jacobian = @. (
        (d - 1) * log(1.0 + β_hat * h_)
        + log(1.0 +  β_hat * h_ + β_hat * (- h_ ^ 2) * r)
    )   # from eq(14)
    return (rv=transformed, logabsdetjac=vec(log_det_jacobian))
end

forward(flow::RadialLayer, z) = _forward(flow, z)

function forward(flow::RadialLayer, z::AbstractVector{<:Real})
    res = _forward(flow, z)
    return (rv=vec(res.rv), logabsdetjac=res.logabsdetjac[1])
end

function (ib::Inversed{<:RadialLayer})(y)
    flow = ib.orig
    α = softplus(flow.α_[1])            # from A.2
    β_hat = - α + softplus(flow.β[1])   # from A.2
    # Define the objective functional
    f(y) = r -> norm(y - flow.z_0, 2) - r * (1 + β_hat / (α + r))   # from eq(26)
    # Run solver
    rs = [find_zero(f(y[:,i:i]), 0.0, Order16()) for i in 1:size(y, 2)]'    # from A.2
    z_hat = (y .- flow.z_0) ./ (rs .* (1 .+ β_hat ./ (α .+ rs)))            # from eq(25)
    z = flow.z_0 .+ rs .* z_hat # from A.2
    return z
end

function (ib::Inversed{<: RadialLayer})(y::AbstractVector{<:Real})
    return vec(ib(reshape(y, (length(y), 1))))
end

logabsdetjac(flow::RadialLayer, x) = forward(flow, x).logabsdetjac
isclosedform(b::Inversed{<:RadialLayer}) = false


text:
"""
    _guesskey(ks, possibilities)

Unfortunately, many model types that contain dictionares do not have
standardized field names, so we need to try a few possibilities and guess the
best one. The keys used to look for valid field names should be ideally
specified as constants in `src/base/constants.jl`.
"""
function _guesskey(avail, possibilities)
    x = intersect(possibilities, avail)

    if isempty(x)
        @debug "could not find any of keys: $possibilities"
        return nothing
    end

    if length(x) > 1
        @debug "Possible ambiguity between keys: $x"
    end
    return x[1]
end


text:
# testing FP growth in Julia
using Gallium

n1 = FPNode()
n2 = FPNode("a", 3, n1)
@code_warntype FPNode()
@code_warntype FPNode("a", 3, n1)

transacts = [["a", "b", "c"],
             ["a", "c", "d", "e"],
             ["b", "e"],
             ["b", "a", "c"],
             ["c", "d"],
             ["a", "b"],
             ["d", "c", "e"],
             ["a", "b", "c", "d"],
             ["c", "d", "e"]]

find_frequent_patterns(transacts, 0.001)


text:
using LinearAlgebra

mutable struct EpicomposeGramDiagonal{M, P, R} <: Epicompose
    L::M
    f::P
    mu::R
    function EpicomposeGramDiagonal{M, P, R}(L::M, f::P, mu::R) where {M, P, R}
        if mu <= 0
            error("mu must be positive")
        end
        new(L, f, mu)
    end
end

EpicomposeGramDiagonal(L, f, mu) = EpicomposeGramDiagonal{typeof(L), typeof(f), typeof(mu)}(L, f, mu)

function prox!(y::AbstractArray{C}, g::EpicomposeGramDiagonal, x::AbstractArray{C}, gamma::R=one(R)) where {
    R <: Real, C <: Union{R, Complex{R}}
}
    z = (g.L'*x)/g.mu
    p, v = prox(g.f, z, gamma/g.mu)
    mul!(y, g.L, p)
    return v
end

function prox_naive(g::EpicomposeGramDiagonal, x::AbstractArray{C}, gamma::R=one(R)) where {
    R <: Real, C <: Union{R, Complex{R}}
}
    z = (g.L'*x)/g.mu
    p, v = prox_naive(g.f, z, gamma/g.mu)
    y = g.L*p
    return y, v
end


text:
using DSGE, ModelConstructors, Dates, OrderedCollections, Test, CSV, DataFrames, Random, FileIO

regenerate_output = false

## Regime switching and full-distribution
# Initialize model objects
Random.seed!(1793 * 10)
m = Model1002("ss60"; custom_settings = Dict{Symbol, Setting}(:flexible_ait_policy_change =>
                                                              Setting(:flexible_ait_policy_change,
                                                                      false),
                                                              :add_pgap => Setting(:add_pgap, false),
                                                              :add_ygap => Setting(:add_ygap, false))) # Set to false unless you want to re-generate any saved output
df_full = CSV.read(joinpath(dirname(@__FILE__), "../reference/uncertain_altpolicy_data.csv"), DataFrame)
m <= Setting(:forecast_horizons, 20)
m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, # Have to add anticipated rates to conditional data
                                :obs_nominalrate, :obs_longrate,
                                :obs_nominalrate1, :obs_nominalrate2, :obs_nominalrate3,
                                :obs_nominalrate4, :obs_nominalrate5, :obs_nominalrate6])
m <= Setting(:cond_semi_names, [:obs_spread,
                                :obs_nominalrate, :obs_longrate,
                                :obs_nominalrate1, :obs_nominalrate2, :obs_nominalrate3,
                                :obs_nominalrate4, :obs_nominalrate5, :obs_nominalrate6])
m <= Setting(:date_forecast_start, Date(2020, 6, 30))
m <= Setting(:date_conditional_end, Date(2020, 6, 30))
m <= Setting(:forecast_ndraws, 5)
m <= Setting(:forecast_block_size, get_setting(m, :forecast_ndraws))
m <= Setting(:forecast_jstep, 1)
m <= Setting(:use_parallel_workers, false)
setup_regime_switching_inds!(m; cond_type = :full)

# For parameters, use modal ones, but also draw from prior for the COVID-19 shocks
mparas = repeat(ModelConstructors.get_values(m.parameters)', get_setting(m, :forecast_ndraws), 1)
j = length(m.parameters) # To help index for regime-switching parameters
biidc_keys = Dict()
ziid_keys = Dict()
φ_keys = Dict()
for p in m.parameters # Find the keys of the COVID-19 shocks
    if haskey(p.regimes, :value)
        for k in keys(p.regimes[:value])
            if k > 1
                global j += 1
                if p.key == :σ_φ
                    φ_keys[k] = j
                elseif p.key == :σ_biidc
                    biidc_keys[k] = j
                elseif p.key == :σ_ziid
                    ziid_keys[k] = j
                end
            end
        end
    end
end

for i in 1:size(mparas, 1)
    mparas[i, φ_keys[2]] = rand(get(m.parameters[m.keys[:σ_φ]].prior)) # draw from prior for second-regime value
    mparas[i, biidc_keys[2]] = rand(get(m.parameters[m.keys[:σ_biidc]].prior)) # draw from prior for second-regime value
    mparas[i, ziid_keys[2]] = rand(get(m.parameters[m.keys[:σ_ziid]].prior)) # draw from prior for second-regime value
    mparas[i, φ_keys[3]] = rand(get(m.parameters[m.keys[:σ_φ]].prior)) # draw from prior for third-regime value
    mparas[i, biidc_keys[3]] = rand(get(m.parameters[m.keys[:σ_biidc]].prior)) # draw from prior for third-regime value
    mparas[i, ziid_keys[3]] = rand(get(m.parameters[m.keys[:σ_ziid]].prior)) # draw from prior for third-regime value
end

# Function for adding new regimes to parameters
function covid_set_regime_vals(m::AbstractDSGEModel, n::Int)
    if n > 4
        for p in m.parameters
            if haskey(p.regimes, :value)
                for i in 5:n
                    ModelConstructors.set_regime_val!(p, i, ModelConstructors.regime_val(p, 4))
                end
            end
        end
    end
end

# Run forecast!
output_vars = [:forecastobs, :bddforecastobs]
forecast_one(m, :full, :full, output_vars, verbose = :none, params = mparas, df = df_full,
             zlb_method = :temporary_altpolicy, set_regime_vals_altpolicy = covid_set_regime_vals)

# Either test against saved output or re-generate output
output_files = get_forecast_output_files(m, :full, :full, output_vars)
if regenerate_output
    using JLD2, FileIO
    if (VERSION >= v"1.5")
        JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist_v1p5.jld2"),
                     true, true, true, IOStream) do file
            for (k, v) in output_files
                write(file, string(k), load(v, "arr"))
            end
        end
    else
        JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist.jld2"),
                     true, true, true, IOStream) do file
            for (k, v) in output_files
                write(file, string(k), load(v, "arr"))
            end
        end
    end
else
    refdata = (VERSION >= v"1.5") ?
        load(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist_v1p5.jld2")) :
        load(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist.jld2"))
    @testset "Automatic enforcement of ZLB as a temporary alternative policy during full-distribution forecast" begin
        for (k, v) in output_files # Test variables which don't have forward-looking measurement equations

            # measurement eqns for forward looking variables e.g. :obs_longrate are still being updated.
            if k == :bddforecastobs
                @test_matrix_approx_eq refdata[string(k)][:, # this test seems to work but perhaps due to changes in the way the automatic ZLB is
                                                                vcat(1:9, 12:13), :] load(v, "arr")[:, vcat(1:9, 12:13), :] # inferred, it no longer matches reference output. Leaving this test as "broken" for now.
                @test all(refdata[string(k)][:, m.observables[:obs_nominalrate], :] .> -1e-14)
            else
                @test @test_matrix_approx_eq refdata[string(k)][:,
                                                                vcat(1:9, 12:13), :] load(v, "arr")[:, vcat(1:9, 12:13), :]
                @test !all(refdata[string(k)][:, m.observables[:obs_nominalrate], :] .> -1e-14)
            end
        end
    end
    for v in values(output_files)
        rm(v)
    end
end


text:
version = v"1.7.9"
api_version = v"1.3.7"

# Collection of sources required to build this package
#
# The level zero repository contains both the API headers and a loader,
# which are versioned independently.
sources = BinaryBuilder.AbstractSource[
    GitSource("https://github.com/oneapi-src/level-zero.git",
              "78f08cc338d70ac60a8f61084ad194fa0dbc90b0"),
]


text:
using ProjectMoana

sphere = Sphere()
ray1 = Ray(Point(0,0,2), Vec(0,0,-1))
intersection1 = ray_intersection(sphere,ray1)
ray2 = Ray(Point(2,0,0), Vec(-1,0,0))
intersection2 = ray_intersection(sphere,ray2)
ray3 = Ray(Point(0,0,0), Vec(0,1,0))
intersection3 = ray_intersection(sphere,ray3)

sphere2 = Sphere(translation(Vec(10,0,0)))
ray4 = Ray(Point(10,0,2), Vec(0,0,-1))
intersection4 = ray_intersection(sphere2,ray4)
ray5 = Ray(Point(13,0,0), Vec(-1,0,0))
intersection5 = ray_intersection(sphere2,ray5)
ray6 = Ray(Point(0,0,2), Vec(0,0,1))
intersection6 = ray_intersection(sphere,ray6)
ray7 = Ray(Point(-10,0,0), Vec(0,0,-1))
intersection7 = ray_intersection(sphere,ray7)

@testset "Shape: Sphere" begin
	@test intersection1 ≈ HitRecord(Point(0,0,1), Normal(0,0,1), Vec2D(0,0), 1.0, ray1, sphere )
	@test intersection2 ≈ HitRecord(Point(1,0,0), Normal(1,0,0), Vec2D(0,0.5), 1.0, ray2 , sphere)
	@test intersection3 ≈ HitRecord(Point(0,1,0), Normal(0,-1,0), Vec2D(0.25,0.5), 1.0, ray3 , sphere)

	@test intersection4 ≈ HitRecord(Point(10,0,1), Normal(0,0,1), Vec2D(0,0), 1.0, ray4, sphere2 )

	@test intersection5 ≈ HitRecord(Point(11,0,0), Normal(1,0,0), Vec2D(0,0.5), 2.0, ray5, sphere2 )
	@test nothing == intersection6
	@test nothing == intersection7
end

plane = Plane(translation(Vec(0.,0.,0.)), Material(DiffuseBRDF(CheckeredPigment(RGB(0.2, 0.7, 0.3), RGB(0.3, 0.2, 0.7), 8))))
ray1 = Ray(Point(11.5,0.3,10), Vec(0,0,-1))
intersection1 = ray_intersection(plane,ray1)
ray2 = Ray(Point(11.5,0.3,10), Vec(0,0,1))
intersection2 = ray_intersection(plane,ray2)
plane2 = Plane(rotation_x(pi/2.f0)) #plane XZ
ray3 = Ray(Point(1.1,3.2,4.0), Vec(0,-1,0))
intersection3 = ray_intersection(plane2, ray3)

@testset "Shape: Plane" begin
	@test intersection1 ≈ HitRecord(Point(11.5,0.3,0), Normal(0,0,1), Vec2D(0.5,0.3), 10., ray1, plane )
	@test nothing == intersection2
	@test intersection3 ≈ HitRecord(Point(1.1,0.,4.0), Normal(0,1,0), Vec2D(0.1,0.), 3.2, ray3, plane2 )
end

@testset "Shape: AAB  " begin
	cube = AAB()
	ray1 = Ray(Point(0.5, 0.5, 0.5), Vec(1.,0.,0.))
	intersection1 = ray_intersection(cube, ray1)
	@test intersection1 ≈ HitRecord(Point(1., 0.5, 0.5), Normal(-1.0, 0.0, 0.0), Vec2D(0.625, 0.5), 0.5, ray1, cube	)

	ray2 = Ray(Point(0.1, 2.0, 0.1), Vec(0.,-1,0.) )
	intersection2 = ray_intersection(cube, ray2)
	@test intersection2 ≈ HitRecord(Point(0.1, 1., 0.1), Normal(0.0, 1.0, 0.0), Vec2D(1.1/4., 0.1/3.), 1., ray2, cube )

	ray1 = Ray(Point(0, 0, 1.5), Vec(0,1,0))
	intersection1 = ray_intersection(cube, ray1)
	@test isnothing(intersection1)

	ray2 = Ray(Point(1.5, 1, 0), Vec(1,1,0))
	intersection2 = ray_intersection(cube, ray2)
	@test isnothing(intersection2)

	# ray3 = Ray(Point(0, 0, -1), Vec(0,0,-1))
	# intersection3 = ray_intersection(cube, ray3)
	# @test isnothing(intersection3)

	# ray3 = Ray(Point(1.62,1.33,0.29), Vec(-0.62,-1.83,0.21))
	# intersection3 = ray_intersection(cube, ray3)
	# @test intersection3 ≈ HitRecord(Point(1.0, -0.5, 0.5), Normal(1.0, 0.0, 0.0), Vec2D(1.5/4., 0.5/3.), 2.16, ray2, cube )

	# ray4 = Ray(Point(-1.76,1.41,1.75), Vec(1.76,-0.65,-1.02))
	# intersection1 = ray_intersection(cube,ray4)
	# @test intersection1 ≈ HitRecord(Point(0.,0.76,0.74), Normal(-1,0,0), Vec2D(0.74/4.,1.76/3.), 2.13, ray4, cube )

	# ray13 = Ray(Point(1.95,1.18,1.62), Vec(-1.45,-0.68,-0.62))
	# intersection2 = ray_intersection(cube, ray13)
	# @test intersection2 ≈ HitRecord(Point(0.5,0.5,1.), Normal(0,0,1), Vec2D(1.5/4.,1.5/3.), 1.72, ray12, cube )

	# cube1 = AAB( translation(Vec(10., 0.0, 0.0)))
	# ray1 = Ray(Point(0., 0.1 , 0.1), Vec(0,1,0))
	# intersection1 = ray_intersection(cube1, ray1)
	# @test intersection1 ≈ HitRecord(Point(10.0, .1, 0.1), Normal(-1.0, 0.0, 0.0), Vec2D(0.5/4., 1.5/3.), 10.0, ray1, cube1)

	# cube2 = AAB(translation(Vec(1., 0., 0.))*rotation_x(π/2.0))
	# ray2 = Ray(Point(0.5, 0.0, 0.5), Vec(0., 1., 0.))
	# intersection2 = ray_intersection(cube2, ray2)
	# @test intersection2 ≈ HitRecord( Point(1., 0.5, 0.5), Normal(0.0, -1., 0.0), Vec2D(1.5/4., 2.5/3.), 1.0, ray2, cube2 ) ≈ intersection2
end


text:
# exported logarithmic functions

const FP_ILOGB0   = typemin(Int)
const FP_ILOGBNAN = typemin(Int)
const INT_MAX     = typemax(Int)

"""
    ilogb(x)

Returns the integral part of the logarithm of `abs(x)`, using base 2 for the
logarithm. In other words, this computes the binary exponent of `x` such that

    x = significand × 2^exponent,

where `significand ∈ [1, 2)`.

* Exceptional cases (where `Int` is the machine wordsize)
    * `x = 0`    returns `FP_ILOGB0`
    * `x = ±Inf`  returns `INT_MAX`
    * `x = NaN`  returns `FP_ILOGBNAN`
"""
function ilogb(x::T) where {T<:Union{Float32,Float64}}
    e = ilogbk(abs(x))
    x == 0 && (e = FP_ILOGB0)
    isnan(x) && (e = FP_ILOGBNAN)
    isinf(x) && (e = INT_MAX)
    return e
end

"""
    log10(x)

Returns the base `10` logarithm of `x`.
"""
function log10(a::T) where {T<:Union{Float32,Float64}}
    x = T(dmul(logk(a), MDLN10E(T)))

    isinf(a) && (x = T(Inf))
    (a < 0 || isnan(a)) && (x = T(NaN))
    a == 0 && (x = T(-Inf))

    return x
end

"""
    log2(x)

Returns the base `2` logarithm of `x`.
"""
function log2(a::T) where {T<:Union{Float32,Float64}}
    u = T(dmul(logk(a), MDLN2E(T)))

    isinf(a) && (u = T(Inf))
    (a < 0 || isnan(a)) && (u = T(NaN))
    a == 0 && (u = T(-Inf))

    return u
end

const over_log1p(::Type{Float64}) = 1e307
const over_log1p(::Type{Float32}) = 1f38

"""
    log1p(x)

Accurately compute the natural logarithm of 1+x.
"""
function log1p(a::T) where {T<:Union{Float32,Float64}}
    x = T(logk2(dadd2(a, T(1.0))))

    a > over_log1p(T) && (x = T(Inf))
    a < -1 && (x = T(NaN))
    a == -1 && (x = T(-Inf))
    isnegzero(a) && (x = T(-0.0))

    return x
end

@inline function log_kernel(x::Float64)
    c7 = 0.1532076988502701353
    c6 = 0.1525629051003428716
    c5 = 0.1818605932937785996
    c4 = 0.2222214519839380009
    c3 = 0.2857142932794299317
    c2 = 0.3999999999635251990
    c1 = 0.6666666666667333541
    return @horner x c1 c2 c3 c4 c5 c6 c7
end

@inline function log_kernel(x::Float32)
    c3 = 0.3027294874f0
    c2 = 0.3996108174f0
    c1 = 0.6666694880f0
    return @horner x c1 c2 c3
end

"""
    log(x)

Compute the natural logarithm of `x`. The inverse of the natural logarithm is
the natural expoenential function `exp(x)`
"""
function log(d::T) where {T<:Union{Float32,Float64}}
    o = d < floatmin(T)
    o && (d *= T(Int64(1) << 32) * T(Int64(1) << 32))

    e = ilogb2k(d * T(1.0/0.75))
    m = ldexp3k(d, -e)
    o && (e -= 64)

    x  = ddiv(dadd2(T(-1.0), m), dadd2(T(1.0), m))
    x2 = x.hi*x.hi

    t = log_kernel(x2)

    s = dmul(MDLN2(T), T(e))
    s = dadd(s, scale(x, T(2.0)))
    s = dadd(s, x2*x.hi*t)
    r = T(s)

    isinf(d) && (r = T(Inf))
    (d < 0 || isnan(d)) && (r = T(NaN))
    d == 0 && (r = -T(Inf))

    return r
end

# First we split the argument to its mantissa `m` and integer exponent `e` so
# that `d = m \times 2^e`, where `m \in [0.5, 1)` then we apply the polynomial
# approximant on this reduced argument `m` before putting back the exponent
# in. This first part is done with the help of the private function
# `ilogbk(x)` and we put the exponent back using

#     `\log(m \times 2^e) = \log(m) + \log 2^e =  \log(m) + e\times MLN2

# The polynomial we evaluate is based on coefficients from

#     `log_2(x) = 2\sum_{n=0}^\infty \frac{1}{2n+1} \bigl(\frac{x-1}{x+1}^{2n+1}\bigr)`

# That being said, since this converges faster when the argument is close to
# 1, we multiply  `m` by `2` and subtract 1 for the exponent `e` when `m` is
# less than `sqrt(2)/2`

@inline function log_fast_kernel(x::Float64)
    c8 = 0.153487338491425068243146
    c7 = 0.152519917006351951593857
    c6 = 0.181863266251982985677316
    c5 = 0.222221366518767365905163
    c4 = 0.285714294746548025383248
    c3 = 0.399999999950799600689777
    c2 = 0.6666666666667778740063
    c1 = 2.0
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8
end

@inline function log_fast_kernel(x::Float32)
    c5 = 0.2392828464508056640625f0
    c4 = 0.28518211841583251953125f0
    c3 = 0.400005877017974853515625f0
    c2 = 0.666666686534881591796875f0
    c1 = 2f0
    return @horner x c1 c2 c3 c4 c5
end

"""
    log_fast(x)

Compute the natural logarithm of `x`. The inverse of the natural logarithm is
the natural expoenential function `exp(x)`
"""
function log_fast(d::T) where {T<:Union{Float32,Float64}}
    o = d < floatmin(T)
    o && (d *= T(Int64(1) << 32) * T(Int64(1) << 32))

    e = ilogb2k(d * T(1.0/0.75))
    m = ldexp3k(d, -e)
    o && (e -= 64)

    x  = (m - 1) / (m + 1)
    x2 = x * x

    t = log_fast_kernel(x2)

    x = x * t + T(MLN2) * e

    isinf(d) && (x = T(Inf))
    (d < 0 || isnan(d)) && (x = T(NaN))
    d == 0 && (x = -T(Inf))

    return x
end


text:
module ApproxInferenceProblems
using ApproxInferenceBase, Distributions

export ApproxInferenceProblem, BlowFly

struct ApproxInferenceProblem{M,D,P,T}
    model::M
    data::D
    prior::P
    target::T
end
function Base.show(io::IO, mime::MIME"text/plain", p::ApproxInferenceProblem{M}) where {M}
    println(io, "ApproxInferenceProblem{M}")
end

include("blowfly.jl")

end # module


text:
import HTTP
import EzXML
using Suppressor

function get_complete_link(root_url, current_url, link)
    # Complete link, good to go
    if occursin(root_url, link)
        return link
    end
    # Relative link
    if startswith(link, "/")
        return "$current_url$link"
    end
    # External link, discard
    return nothing
end

function grab_links(root_url, current_url, current_depth, max_depth)
    if current_depth > max_depth
        return true
    end
    request = nothing
    try
        request = HTTP.request("GET", "$root_url$current_url")
    catch e
        # Skip, link not acessible
        return true
    end
    body = request.body
    html_obj = EzXML.parsehtml(body)
    nodes = EzXML.findall("//a/@href", html_obj)
    if length(nodes) === 0
        return true
    end
    for node in nodes
        link = EzXML.nodecontent(node)
        complete_link = get_complete_link(root_url, current_url, link)
        if complete_link !== nothing
            if (current_depth > 1)
                pad = repeat(" ", current_depth)
                print("$pad └ ")
            end
            print("$link\n")
            grab_links(root_url, complete_link, current_depth + 1, max_depth)
        end
    end
end

function scan_links(url, max_depth=2)
    print("\nGetting links from $url, with depth $max_depth:\n")
    @suppress_err grab_links(url, "", 1, max_depth)
end

max_stars_repo_name: thecrypters/nis

max_issues_repo_name: thecrypters/nis

max_forks_repo_name: thecrypters/julia-host-scanner


text:
#ДАНО: Робот - в произвольной клетке поля (без внутренних перегородок и маркеров)
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы

function kray(r,side) #функция перемещает робота к краю
    while isborder(r,side) == false
        move!(r,side)
    end
end

function mramka(r,side) #обходит по периметру и ставит марки
    while isborder(r,side) == false
        move!(r,side)
        putmarker!(r)
    end
end

function ramka(r)
    kray(r,Sud)
    kray(r,Ost) #ведет робот в угол
    for side in (HorizonSide(i) for i=0:3)
        mramka(r,side)
    end
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
#=using Revise
using LinearAlgebra, CUDA, BenchmarkTools

##WARNING- comment out when live
include("../Annihilator.jl")=#

#creates crosssections of the data and working arrays
#note this is now only used in testing- use the annihilator matrix instead of the below
function dataxsection(D::NamedTuple, ts::Vector{Int},
  ::Type{TM}, ::Type{TV}, ::Type{T}=eltype(TM),
  ) where {
    T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}

  #technically we don't need this, but if its not sorted the process if much harder to reverse
  @assert issorted(ts)
  @assert sum(ts.==1) == 0
  allts = sort!(unique([1; ts...;]))
  dims= (T = length(allts), Nexpanded=length(ts))
  @assert collect(1:dims.T) == allts #validate integrity over time

  Dkeys = keys(D)
  Dtypes = (typeof).(values(D))

  xsections = Vector{NamedTuple}(undef, dims.T-1)
  allinds = collect(1:dims.Nexpanded)
  for t ∈ 1:(dims.T-1)
    #we build the named tuple by first creating a dictionary, then converting to a named tuple
    xinds = allinds[ts.==t+1] #cache this since we will call it a lot

    xsections[t] = (; broadcast(zip(Dkeys,Dtypes)) do (k,DT) #form a named tuple for each time t
      if DT <: TM
        k=>TM(D[k][xinds,:])
      elseif DT <: TV
        k=>TV(D[k][xinds])
      else
        @assert false
      end
    end...)
  end

  #now make a tuple of the data fields, where each field is a vector for all t
  xsection = (; [k=>(x->x[k]).(xsections) for k ∈ Dkeys]...)

  return xsection
end

abstract type AbstractXYIter{TM<:AbstractMatrix, TV<:AbstractVector, T<:Real} end
abstract type AbstractXYIterControl{TM<:AbstractMatrix, TV<:AbstractVector,
   T<:Real, TxM<:Union{Nothing, Vector{<:AbstractAnnhilator}}} <: AbstractXYIter{TM,TV,T} end

const XYIterControlLevel = AbstractXYIterControl{
  <:AbstractMatrix, <:AbstractVector, <:Real, <:Nothing}
const XYIterControlDemean= AbstractXYIterControl{
  <:AbstractMatrix, <:AbstractVector, <:Real, <:AbstractVector}

function findnonsingularcols(Mraw::AbstractMatrix)
  dims = (rows = size(Mraw,1), cols = size(Mraw,2))
  rank(Mraw) == dims.cols && throw("No singular cols to drop!")

  Mcols = [true; falses(dims.cols-1)]
  for c ∈ 2:(dims.cols)
    Mcols[c] = true
    sM = view(Mraw, :, Mcols)
    Mcols[c] = rank(sM) == sum(Mcols)
  end
  sM = view(Mraw, :, Mcols)
  @assert rank(sM) == size(sM,2) == rank(Mraw)
  return Mcols
end

rankrt(M::AbstractMatrix{T}) where T = rank(M, eps(T)^0.5)

function findnonsingularcols2(Mraw::AbstractMatrix)
  dims = (rows = size(Mraw,1), cols = size(Mraw,2))
  isposdef!(Symmetric(Mraw' * Mraw)) &&  (rankrt(Mraw) == dims.cols) && throw(
    "No singular cols to drop!")

  Mcols = [true; falses(dims.cols-1)]
  for c ∈ 2:(dims.cols)
    Mcols[c] = true
    M = Mraw[:, Mcols]
    Mcols[c] = (rankrt(M) == sum(Mcols)) && isposdef!(Symmetric(M'*M))
  end
  M = Mraw[:, Mcols]
  @assert rankrt(M) == size(M,2)
  @assert isposdef!(Symmetric(M'*M))
  return Mcols
end

#this function groups the controls as well as checking the matrices for degeneracy
function groupcontrolsbyt(Wgroup::TM, ts, ::Type{T}=eltype(TM))where {TM, T}
  xWgroupraw = dataxsection((xWgroup=Wgroup |> Matrix,), ts, Matrix{Float64}, Vector{Float64}, T)

  #checking for control independence
  xWgroupcpu = similar(xWgroupraw.xWgroup)
  KW::Int = size(Wgroup,2)
  droppedcols::Vector{Int} = zeros(Int, KW)
  print("Checking for control column independence t∈2:$(length(xWgroupcpu)).
    Singular cols found, if any, at t=[")
  Wgroupcolsused = trues((maximum(ts)-1)*KW) #this is an index of the control cols
  for (tM1, Wraw) ∈ enumerate(xWgroupraw.xWgroup)
    if (rankrt(Wraw) == size(Wraw,2)) && isposdef!(Symmetric(Wraw' * Wraw))
      xWgroupcpu[tM1] = Wraw
      @assert isposdef!(Symmetric(xWgroupcpu[tM1]' * xWgroupcpu[tM1])) "
        matrix: $(Symmetric(xWgroupcpu[tM1]' * xWgroupcpu[tM1]))"
    else
      nonsingular = view(Wgroupcolsused,((tM1 - 1)*KW+1):(tM1*KW))
      @assert (nonsingular .=== true) |> all
      nonsingular .= findnonsingularcols2(Wraw)
      print("$(tM1+1), ")
      droppedcols .+= (!).(nonsingular)
      #@info "Singular control columns found in group at t=$(tM1+1)." *
      #  "Dropped cols $((1:(length(nonsingular)))[(!).(nonsingular)])"
      xWgroupcpu[tM1] = Wraw[:, nonsingular]
      @assert isposdef!(Symmetric(xWgroupcpu[tM1]' * xWgroupcpu[tM1])) "
        nonsingular: $nonsingular
        matrix: $(Symmetric(xWgroupcpu[tM1]' * xWgroupcpu[tM1]))"

    end

    @assert isposdef!(Symmetric(xWgroupcpu[tM1]' * xWgroupcpu[tM1])) "
      matrix: $(Symmetric(xWgroupcpu[tM1]' * xWgroupcpu[tM1]))"
  end
  println("]", sum(droppedcols) > 0 ?
    reduce(*, ["\ncol c dropped (ntimes): ";
      ((c)->"$(c[1]) ($(c[2])), ").(enumerate(droppedcols));]) : "")

  #convert the type if necessary
  xWgroup = (xWgroup = TM === Matrix{Float64} ? xWgroupcpu : TM.(xWgroupcpu),)

  return (;xWgroup, Wgroupcolsused)
end

#this creates the expanded control matrix
function expandcontrolsbyt(Wgroup::TM, ts, ::Type{T}=eltype(TM)) where {TM,T}

  xWgroup, Wgroupcolsused = groupcontrolsbyt(Wgroup, ts) |> (nt)->(nt.xWgroup, nt.Wgroupcolsused)
  K = size(Wgroup,2)
  N = length(ts)

  expanded = similar(Wgroup, N, K * length(ts |> unique))
  #expanded = similar(Wgroup, N, sum(Wgroupcolsused))
  expanded .= 0.0

  rctr = 1
  for (tM1, xWgroupt) ∈ enumerate(xWgroup.xWgroup)
    Nt = size(xWgroupt,1)

    #be careful here to not capture non-singular columns
    #proceed in two steps
    submatrix = view(expanded, rctr:(rctr + Nt-1), (K*(tM1 - 1)+1):(tM1*K))
    subcolumnsused = view(Wgroupcolsused, (K*(tM1 - 1)+1):(tM1*K))
    submatrix = view(submatrix, :, subcolumnsused)

    @assert length(subcolumnsused) == K
    @assert all(submatrix .== 0.0)
    submatrix .= xWgroupt
    rctr += Nt
  end

  #now run some checks
  @assert rctr == N+1
  #@assert (sum(expanded, dims=2) .≈ sum(Wgroup, dims=2)) |> all
  nonsingularexpanded = expanded[:, Wgroupcolsused]
  checkcolsums = reduce(vcat, (xWgroupt -> sum(xWgroupt,dims=1) |> vec).(xWgroup.xWgroup))
  @assert (checkcolsums .≈ vec(sum(nonsingularexpanded, dims=1))) |> all
  checkrowsums = reduce(vcat, (xWgroupt -> sum(xWgroupt,dims=2) |> vec).(xWgroup.xWgroup))
  @assert (checkrowsums .≈ vec(sum(nonsingularexpanded, dims=2))) |> all

  return (;expanded=nonsingularexpanded, Wgroupcolsused)
end

#holds all unchanging data as well as helper pre-allocations
#this version also includes the projection matrix for demeaning
struct XYIterControl{TM<:AbstractMatrix,
    TV<:AbstractVector,
    T<:Real,
    TxM,
    Txsection} <: AbstractXYIterControl{TM,TV,T, TxM}
  ws::TM
  Rtws::TM
  RLws::TM
  v::TV
  ṽ::TV #ṽ = xM * W

  W::TM
  W̃::TM #W̃ = xM * W

  xsection::Txsection

  #constructor from minimum componeents assuming a Wgroup cross-secitonal component
  function XYIterControl(ws::TM, Rtws::TM, RLws::TM,  v::TV, ts::Tts, W::TM,  Wgroup::TM,
      ::Type{T}=eltype(TM)) where {TM<:AbstractMatrix,TV, Tts<:AbstractVector, T}

    issorted(ts) || error("ts needs to be sorted, otherwise row orderwing will be misalligned
      due to the reduce->vcat")

    #no need for anything fancy if there are no cross-secitonal controls
    if size(Wgroup,2) == 0
      smallxsection = (;xM=nothing, dataxsection((;xws=ws, xRtws=Rtws,xRLws=RLws), ts, TM, TV, T)...)
      return new{TM,TV, T, Nothing, typeof(smallxsection)}(ws, Rtws, RLws, v, v, W, W, smallxsection)
    end

    #(size(Wgroup,2) == 0) && return new{TM,TV, T, Nothing}(ws, RLws, v, v, W, W,nothing)

    local W̃::TM

    #construct a compact form of the annihilator matrix
    #only pre-multiplication of the matrix on a matrix or vector is supported
    #start by sectioning off components. These will be discarded

    #Form the grouped projection matrices
    xv = dataxsection((xv=v,), ts, TM, TV, T)
    xWgroup = groupcontrolsbyt(Wgroup, ts).xWgroup
    xM = (xWgroupt->DenseAnnihilator(xWgroupt)).(xWgroup.xWgroup)

    #project W onto the columnspace of Ws to form W̃. Similarly project v onto Ws
    ṽ = reduce(vcat, ((vt,Mt)->Mt*vt).(xv.xv, xM))

    #if size(W,2) > 0
    xW = dataxsection((xW=W,), ts, TM, TV, T)
    W̃ = reduce(vcat, ((Wt,Mt)->Mt*Wt).(xW.xW, xM))
    #=else
      W̃ = W
    end=#

    #now bundle the relevant components together for easy group calculations
    D = (xws=ws, xRtws=Rtws, xRLws=RLws)
    xsection = (;xM=xM, dataxsection(D, ts, TM, TV, T)...)

    #a quick integrity check
    xsectioncheck= dataxsection((;xṽ=ṽ, xW̃=W̃), ts, TM, TV, T)
    for (xvt, xṽt, xWt, xW̃t, xMt) ∈ zip(
        xv.xv, xsectioncheck.xṽ, xW.xW, xsectioncheck.xW̃, xsection.xM)

      #=if !((xṽt |> Vector) ≈ (
        xvt - xMt.W * (cholesky(Symmetric(xMt.W' * xMt.W))\(xMt.W' * xvt))) |> Vector)
        println("ERROR !((xṽt |> Vector) ≈ (
          xvt - xMt.W * (cholesky(Symmetric(xMt.W' * xMt.W))\\(xMt.W' * xvt))) |> Vector)")
        println("xvt:")
        printmln(xṽt)
        println("xvt - (...* xvt)")
        printmln((
          xvt - xMt.W * (cholesky(Symmetric(xMt.W' * xMt.W))\(xMt.W' * xvt))) |> Vector)
        printmln("xvt: ")
        printmln(xvt)
        error("stop")
      end=#

      heuristicatol = mean(abs.(xvt)*1e-4)
      heuristicisapprox(x,y) = isapprox(x,y,atol=heuristicatol)
      @assert heuristicisapprox(xṽt |> Vector,
        (xvt - xMt.W * (cholesky(Symmetric(xMt.W' * xMt.W))\(xMt.W' * xvt))) |> Vector)# |> all
      @assert heuristicisapprox(xW̃t |> Matrix,
        (xWt - xMt.W * (cholesky(Symmetric(xMt.W' * xMt.W))\(xMt.W' * xWt))) |> Matrix)# |> all
    end

    return new{TM,TV, T, typeof(xM), typeof(xsection)}(ws, Rtws, RLws, v, ṽ, W, W̃, xsection)
  end
end

XYIterControl(ws::TM, Rtws::TM, RLws::TM,  v::TV, ts::Tts, ::Type{T}=eltype(TM);
    W::TM = similar(ws, size(ws,1),0),
    Wgroup::TM = similar(ws, size(ws,1),0), #empty matrix
    ) where {TM,TV, Tts, T} = XYIterControl(ws, Rtws, RLws,  v, ts, W, Wgroup, T)

#entry point for default constructor of XYIter object
XYIter(ws::TM, Rtws::TM, RLws::TM,  v::TV, ts::Tts,;
  PredictionType::Val=Val(PARAM[:predictiontype])
  ) where {TM, TV, Tts} = XYIter(ws, Rtws, RLws, v, ts, PredictionType)

XYIter(ws::TM, Rtws::TM, RLws::TM,  v::TV, Z::TM, ts::Tts, ;
  PredictionType::Val=Val(PARAM[:predictiontype])
  ) where {TM, TV, Tts} = XYIter(ws, Rtws, RLws, v, Z, ts, PredictionType)

XYIter(ws::TM, Rtws::TM, RLws::TM,  v::TV, ts::Tts, ::Val{:leveldemean}
  ) where {TM, TV, Tts} = XYIterDemean(ws, Rtws, RLws, v, ts)

XYIter(ws::TM, Rtws::TM, RLws::TM,  v::TV, ts::Tts, ::Val{:level}
  ) where {TM, TV, Tts} = XYIterLevel(ws, Rtws, RLws, v, ts)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#
# JSON Server with julia
#
# author: Atsushi Sakai
#

using HttpServer
using JSON

function response_json(req)
	println("json requested!!")

	# Parse input json
    jreq=JSON.parse(String(req.data))
	println("request:")
	println(jreq)

    jres = Dict("title"=>"Julia JSON server", "PI" => pi)
    jres["time"] = Dates.format(Dates.now(), "yyyymmddHHMMSS")
	println("response:")
	println(jres)

    return JSON.json(jres)
end

function main()
    println(PROGRAM_FILE," start!!")

    http = HttpHandler() do req::Request, res::Response
        Response(response_json(req))
    end

    server = Server(http)
    host = IPv4(127,0,0,1) # localhost
    port = 8000
    run(server, host=host, port=port)
    println("serving at port", port)

    println(PROGRAM_FILE," Done!!")
end

main()

max_stars_repo_name: AtsushiSakai/JSONServerAndClient

max_issues_repo_name: AtsushiSakai/JSONServerAndClient

max_forks_repo_name: AtsushiSakai/JSONServerAndClient


text:
function household(TM::TaxModel)
    @unpack na,ns,nty,J,κ0,κ1,κ2,τk,a̲  = TM

#INPUTS NEEDED: agrid,Tr
Vfun = zeros(na,ns,nty,J)
lfun = zeros(na,ns,nty,J)
cfun = zeros(na,ns,nty,J)
afun = zeros(na,ns,nty,J)
vpfun = zeros(na,ns,nty,J)

# Last generation
for a in 1:na
    # Calculating earnings and taxes
    earnl = 0.0
    earncap = r*(agrid[a] + Tr)
    lab_tax = tax_gs(κ0,κ1,κ2,earnl)
    cap_tax = τk*earncap # DIFF THAN REAL CODE
    mtaxes = marginal_tax_gs(κ0,κ1,κ2,earnl) + τk*earncap

    # Consume all assets and supply no labor
    l = 0.0
    c = (SS+(1.0+r(1-τk))*(agrid[a]+Tr))/(1.0+τc)

    # Setting function values
	Vfun[a,1:ns,1:nty,J] = U(c,l)
	lfun[a,1:ns,1:nty,J] = l
	cfun[a,1:ns,1:nty,J] = c
	afun[a,1:ns,1:nty,J] = 0.0
	vpfun[a,1:ns,1:nty,J]=(1.0+r*(1.0-τk))*marginal_utility(con,lab)
end

# Rest of time period
for j in 1:(J-1)
    for ty in 1:nty
        for s in 1:ns
            for a in 1:na
                # Retirement periods
                if j >= jr
                    # Don't work in retirement
                    lfun[a,s,ty,j] = 0.0
                    # Bounds on capital
                    lbound = a̲
                    ubound = agrid[na] # assets will be decreasing
                    # earnings
                    earnl = 0.0
                    earncap = r*(agrid[a]+Tr)
                    # Taxes
                    lab_tax = tax_gs(κ0,κ1,κ2,earnl)
                    cap_tax = τk*earncap # DIFF THAN REAL CODE
                    mtaxes = marginal_tax_gs(κ0,κ1,κ2,earnl) + τk*earncap

                    lab=0.0
			        conl = (SS+(1.0+r(1-τk))*(agrid[a]+Tr) - lbound)/(1.0+τc)
			        conh = (SS+(1.0+r(1-τk))*(agrid[a]+Tr) - ubound)/(1.0+τc)


text:
# ===== Showcasing the main functionalities of WaveletsExt.jl =====

using Plots,
      Wavelets,
      WaveletsExt

# 1. Additional Redundant Wavelet Transform Methods ----------------------------------------
#    On top of maximal overlap transform implemented in Wavelets.jl, WaveletsExt.jl has
#    additional redundant transform functions through Autocorrelation Wavelet Transforms
#    (Beylkin, Saito), Stationary Wavelet Transforms (Nason, Silverman), and Shift Invariant
#    Wavelet Transforms (Cohen, Raz, Malah).

x = zeros(1<<8)
x[128] = 1
wt = wavelet(WT.db4)                    # Construct Daubechies 4-tap wavelet filter

# ----- Autocorrelation Wavelet Transforms -----
y = acdwt(x, wt)
p1 = wiggle(y) |> p -> plot!(p, yticks=1:9, title="Autocorrelation WT")

# ----- Stationary Wavelet Transforms -----
y = sdwt(x, wt)
p2 = wiggle(y) |> p -> plot!(p, yticks=1:9, title="Stationary WT")

# Combine and save plot
p = plot(p1, p2, layout=(1,2))
savefig(p, "transforms.png")

# 2. Best Basis Algorithms -----------------------------------------------------------------
#    Compared to Wavelets.jl, WaveletsExt.jl has an extended ability of catering toward
#    multiple signals at once via Joint Best Basis (JBB) and Least Statistically Dependent
#    Basis (LSDB). One may then use the `plot_tfbdry` function implemented in WaveletsExt.jl
#    to analyze the best basis subspace.

# Generate 100 noisy heavisine signals of length 2⁸
x = generatesignals(:heavisine, 8) |> x -> duplicatesignals(x, 100, 2, true, 0.5)
# Wavelet packet decomposition of all signals
xw = wpdall(x, wt, 6)

# ----- Joint Best Basis (JBB)
tree = bestbasistree(xw, JBB())
p1 = plot_tfbdry(tree, 6, nd_col=:green, ln_col=:black, bg_col=:white) |>
     p -> plot!(p, title="JBB")

# ----- Least Statistically Dependent Basis (LSDB)
tree = bestbasistree(xw, LSDB())
p2 = plot_tfbdry(tree, 6, nd_col=:green, ln_col=:black, bg_col=:white) |>
     p -> plot!(p, title="LSDB")

# Combine and save plot
p = plot(p1, p2, layout=(1,2))
savefig(p, "paper/bestbasis.png")

# 3. Denoising Algorithms ------------------------------------------------------------------
#    There are two functions available: `denoise` and `denoiseall`. The former denoises one
#    signal whereas the latter denoises multiple signals at once.

# Generate 6 circularly shifted original heavisine signals
x₀ = generatesignals(:heavisine, 8) |> x -> duplicatesignals(x, 6, 2, false)
# Generate 6 noisy versions of the original signals
x = generatesignals(:heavisine, 8) |> x -> duplicatesignals(x, 6, 2, true, 0.8)

# Decompose each noisy signal
xw = wpdall(x, wt)

# Get best basis tree from the decomposition of signals
bt = bestbasistree(xw, JBB())
# Get best basis coefficients based on best basis tree
y = getbasiscoefall(xw, bt)

# Denoise all signals based on computed best basis tree
x̂ = denoiseall(y, :wpt, wt, tree=bt)

# Plot results
xs = repeat([0,256],6) |> x -> reshape(x, (2,6))
ys = repeat(1:6, inner=2) |> x -> reshape(x, (2,6))
p1 = plot(title="Noisy Signals")
wiggle!(x₀, sc=0.7, FaceColor=:white, ZDir=:reverse)
wiggle!(x, sc=0.7, EdgeColor=:red, FaceColor=:white, ZDir=:reverse)
plot!(p1, xs, ys, lc=:black)

p2 = plot(title="Denoised Signals")
wiggle!(x₀, sc=0.7, FaceColor=:white, ZDir=:reverse)
wiggle!(x̂, sc=0.7, EdgeColor=:blue, FaceColor=:white, ZDir=:reverse)
plot!(p2, xs, ys, lc=:black)

# Combine and save plot
p = plot(p1, p2, layout=(1,2))
savefig(p, "paper/denoising.png")

# 4. Feature Extraction Algorithm For Signals ----------------------------------------------
#    The feature extraction algorithm implemented in WaveletsExt.jl is the Local
#    Discriminant Basis (LDB) (Saito, Coifman).

# Generate 100 signals for each class of cylinder-bell-funnel
X, y = generateclassdata(ClassData(:cbf, 100, 100, 100))
# View sample signals and how each class differs from one another
cylinder = wiggle(X[:,1:5], sc=0.3)
plot!(cylinder, yticks=1:5, title="Cylinder signals")
bell = wiggle(X[:,101:105], sc=0.3)
plot!(bell, yticks=1:5, title="Bell signals")
funnel = wiggle(X[:,201:205], sc=0.3)
plot!(funnel, yticks=1:5, title="Funnel signals")
p1 = plot(cylinder, bell, funnel, layout=(3,1))

# Define Local Discriminant Basis object (`n_features` can be tweaked depending on the
# number of desired features to be used as input into classification model)
wt = wavelet(WT.coif4)
ldb = LocalDiscriminantBasis(wt=wt,
                             max_dec_level=6,
                             dm=SymmetricRelativeEntropy(),
                             en=TimeFrequency(),
                             dp=BasisDiscriminantMeasure(),
                             top_k=10,
                             n_features=10)

# Fit and transform signals into features using LDB
X̂ = fit_transform(ldb, X, y)

# Plot the best basis for feature extraction
p2 = plot_tfbdry(ldb.tree, 6, nd_col=:green, ln_col=:black, bg_col=:white)
plot!(p2, title="Basis Selection using LDB")

p = plot(p1, p2, size=(600,300))
savefig(p, "paper/ldb.png")


text:
abstract type AbstractLayer end

mutable struct LinearLayer <: AbstractLayer
    W::AbstractArray{Float64,2}
    b::AbstractArray{Float64,1}
    input::AbstractArray{Float64,2}
    output::AbstractArray{Float64,2}
    dim_input::Int64
    dim_output::Int64
    # for Momentum SGD, RMSProp, and Adam
    W_momentum::AbstractArray{Float64,2}
    W_momentum_aux::AbstractArray{Float64,2}
    b_momentum::AbstractArray{Float64,1}
    b_momentum_aux::AbstractArray{Float64,1}
    param::Float64
    param_aux::Float64
    function LinearLayer(
        dim_input::Int64,
        dim_output::Int64,
        std::Float64 = 1.0,
        bias::Float64 = 0.0,
        param::Float64 = 0.9,
        param_aux::Float64 = 0.9,
    )
        trunc = truncated(Normal(0.0, std), -2.0 * std, 2.0 * std)
        W = collect(reshape(rand(trunc, dim_output, dim_input), dim_output, dim_input))
        b = ones(dim_output) * bias
        new(
            W,
            b,
            zeros(1, 1),
            zeros(1, 1),
            dim_input,
            dim_output,
            zeros(Float64, size(W)),
            zeros(Float64, size(W)),
            zeros(Float64, size(b)),
            zeros(Float64, size(b)),
            param,
            param_aux,
        )
    end
end

mutable struct SigmoidLayer <: AbstractLayer
    W::AbstractArray{Float64,2}
    b::AbstractArray{Float64,1}
    input::AbstractArray{Float64,2}
    output::AbstractArray{Float64,2}
    dim_input::Int64
    dim_output::Int64
    # for Momentum SGD, RMSProp, and Adam
    W_momentum::AbstractArray{Float64,2}
    W_momentum_aux::AbstractArray{Float64,2}
    b_momentum::AbstractArray{Float64,1}
    b_momentum_aux::AbstractArray{Float64,1}
    param::Float64
    param_aux::Float64
    function SigmoidLayer(
        dim_input::Int64,
        dim_output::Int64,
        std::Float64 = 1.0,
        bias::Float64 = 0.0,
        param::Float64 = 0.9,
        param_aux::Float64 = 0.9,
    )
        trunc = truncated(Normal(0.0, std), -2.0 * std, 2.0 * std)
        W = collect(reshape(rand(trunc, dim_output, dim_input), dim_output, dim_input))
        b = ones(dim_output) * bias
        new(
            W,
            b,
            zeros(1, 1),
            zeros(1, 1),
            dim_input,
            dim_output,
            zeros(Float64, size(W)),
            zeros(Float64, size(W)),
            zeros(Float64, size(b)),
            zeros(Float64, size(b)),
            param,
            param_aux,
        )
    end
end

mutable struct TanhLayer <: AbstractLayer
    W::AbstractArray{Float64,2}
    b::AbstractArray{Float64,1}
    input::AbstractArray{Float64,2}
    output::AbstractArray{Float64,2}
    dim_input::Int64
    dim_output::Int64
    # for Momentum SGD, RMSProp, and Adam
    W_momentum::AbstractArray{Float64,2}
    W_momentum_aux::AbstractArray{Float64,2}
    b_momentum::AbstractArray{Float64,1}
    b_momentum_aux::AbstractArray{Float64,1}
    param::Float64
    param_aux::Float64
    function TanhLayer(
        dim_input::Int64,
        dim_output::Int64,
        std::Float64 = 1.0,
        bias::Float64 = 0.0,
        param::Float64 = 0.9,
        param_aux::Float64 = 0.9,
    )
        trunc = truncated(Normal(0.0, std), -2.0 * std, 2.0 * std)
        W = collect(reshape(rand(trunc, dim_output, dim_input), dim_output, dim_input))
        b = ones(dim_output) * bias
        new(
            W,
            b,
            zeros(1, 1),
            zeros(1, 1),
            dim_input,
            dim_output,
            zeros(Float64, size(W)),
            zeros(Float64, size(W)),
            zeros(Float64, size(b)),
            zeros(Float64, size(b)),
            param,
            param_aux,
        )
    end
end

mutable struct ReLULayer <: AbstractLayer
    W::AbstractArray{Float64,2}
    b::AbstractArray{Float64,1}
    input::AbstractArray{Float64,2}
    output::AbstractArray{Float64,2}
    dim_input::Int64
    dim_output::Int64
    # for Momentum SGD, RMSProp, and Adam
    W_momentum::AbstractArray{Float64,2}
    W_momentum_aux::AbstractArray{Float64,2}
    b_momentum::AbstractArray{Float64,1}
    b_momentum_aux::AbstractArray{Float64,1}
    param::Float64
    param_aux::Float64
    function ReLULayer(
        dim_input::Int64,
        dim_output::Int64,
        std::Float64 = 1.0,
        bias::Float64 = 0.0,
        param::Float64 = 0.9,
        param_aux::Float64 = 0.9,
    )
        trunc = truncated(Normal(0.0, std), -2.0 * std, 2.0 * std)
        W = collect(reshape(rand(trunc, dim_output, dim_input), dim_output, dim_input))
        b = ones(dim_output) * bias
        new(
            W,
            b,
            zeros(1, 1),
            zeros(1, 1),
            dim_input,
            dim_output,
            zeros(Float64, size(W)),
            zeros(Float64, size(W)),
            zeros(Float64, size(b)),
            zeros(Float64, size(b)),
            param,
            param_aux,
        )
    end
end

function forward_propagation(layer::AbstractLayer, x_::AbstractArray{Float64,1})
    @assert layer.dim_input == size(x_)[1]
    x = collect(reshape(x_, size(x_)[1], 1))
    return forward_propagation_impl(layer, x)
end

function forward_propagation(layer::AbstractLayer, X::AbstractArray{Float64,2})
    @assert layer.dim_input == size(X)[1]
    return forward_propagation_impl(layer, X)
end

function forward_propagation_impl(layer::LinearLayer, X::AbstractArray{Float64,2})
    # X is the array of [x1 x2 ... x_N], x1 is a vector of feature_size
    layer.input = copy(X)
    return (layer.W) * X .+ (layer.b)
end

function forward_propagation_impl(layer::SigmoidLayer, X::AbstractArray{Float64,2})
    # X is the array of [x1 x2 ... x_N], x1 is a vector of feature_size
    layer.input = copy(X)
    activation = (layer.W) * X .+ (layer.b)
    layer.output = 1.0 ./ (1.0 .+ exp.(-activation))
    return layer.output
end

function forward_propagation_impl(layer::TanhLayer, X::AbstractArray{Float64,2})
    # X is the array of [x1 x2 ... x_N], x1 is a vector of feature_size
    layer.input = copy(X)
    activation = (layer.W) * X .+ (layer.b)
    layer.output = tanh.(activation)
    return layer.output
end

function forward_propagation_impl(layer::ReLULayer, X::AbstractArray{Float64,2})
    # X is the array of [x1 x2 ... x_N], x1 is a vector of feature_size
    layer.input = copy(X)
    activation = (layer.W) * X .+ (layer.b)
    # layer.output = activation.clip(min=0)
    layer.output = (activation + broadcast(abs, activation)) / 2.0
    return layer.output
end

function activate_derivative(layer::AbstractLayer)
    return activate_derivative_impl(layer)
end

function activate_derivative_impl(layer::LinearLayer)
    # [f'(z_1), f'(z_2),,,, f'(z_N)]
    return ones(layer.dim_output, size(layer.output)[2]) * 1.0
end

function activate_derivative_impl(layer::SigmoidLayer)
    # [f'(z_1), f'(z_2),,,, f'(z_N)]
    return (layer.output) .* (1.0 .- layer.output)
end

function activate_derivative_impl(layer::TanhLayer)
    # [f'(z_1), f'(z_2),,,, f'(z_N)]
    return 1.0 .- (layer.output) .* (layer.output)
end

function activate_derivative_impl(layer::ReLULayer)
    # [f'(z_1), f'(z_2),,,, f'(z_N)]
    return (layer.output .> 0.0) * 1.0
end

function backward_propagation(
    layer::AbstractLayer,
    delta_output::AbstractArray{Float64,2},
    learning_rate::Float64,
    policy::String = "Adam",
)
    @assert size(delta_output)[1] == layer.dim_output
    if policy == "SGD"
        return backward_propagation_impl_sgd(layer, delta_output, learning_rate)
    elseif policy == "MomentumSGD"
        return backward_propagation_impl_msgd(layer, delta_output, learning_rate)
    elseif policy == "RMSProp"
        return backward_propagation_impl_rmsprop(layer, delta_output, learning_rate)
    elseif policy == "Adam"
        return backward_propagation_impl_adam(layer, delta_output, learning_rate)
    else
        error("Policy must be (SGD|MomentumSGD|RMSProp|Adam)")
    end
end

function backward_propagation_impl_sgd(
    layer::AbstractLayer,
    delta_output::AbstractArray{Float64,2},
    learning_rate::Float64,
)
    delta = delta_output .* activate_derivative_impl(layer)
    W = copy(layer.W)
    layer.W -= learning_rate .* (delta * transpose(layer.input))
    layer.b -= learning_rate .* reshape(sum(delta, dims = 2), layer.dim_output)

    return transpose(W) * delta
end

function backward_propagation_impl_msgd(
    layer::AbstractLayer,
    delta_output::AbstractArray{Float64,2},
    learning_rate::Float64,
)
    delta = delta_output .* activate_derivative_impl(layer)
    W = copy(layer.W)
    W_grad = (delta * transpose(layer.input))
    b_grad = reshape(sum(delta, dims = 2), layer.dim_output)

    mu = layer.param
    W_momentum = mu * layer.W_momentum - (1.0 - mu) * learning_rate * W_grad
    b_momentum = mu * layer.b_momentum - (1.0 - mu) * learning_rate * b_grad

    layer.W_momentum = W_momentum
    layer.b_momentum = b_momentum

    layer.W += W_momentum
    layer.b += b_momentum

    return transpose(W) * delta
end

function backward_propagation_impl_rmsprop(
    layer::AbstractLayer,
    delta_output::AbstractArray{Float64,2},
    learning_rate::Float64,
)
    delta = delta_output .* activate_derivative_impl(layer)
    W = copy(layer.W)
    W_grad = (delta * transpose(layer.input))
    b_grad = reshape(sum(delta, dims = 2), layer.dim_output)

    rho = layer.param
    W_momentum = rho * layer.W_momentum + (1.0 - rho) * (W_grad .^ 2)
    b_momentum = rho * layer.b_momentum + (1.0 - rho) * (b_grad .^ 2)

    layer.W_momentum = W_momentum
    layer.b_momentum = b_momentum

    layer.W -= learning_rate * W_grad ./ sqrt.(W_momentum .+ 1e-8)
    layer.b -= learning_rate * b_grad ./ sqrt.(b_momentum .+ 1e-8)

    return transpose(W) * delta
end

function backward_propagation_impl_adam(
    layer::AbstractLayer,
    delta_output::AbstractArray{Float64,2},
    learning_rate::Float64,
)
    delta = delta_output .* activate_derivative_impl(layer)
    W = copy(layer.W)
    W_grad = (delta * transpose(layer.input))
    b_grad = reshape(sum(delta, dims = 2), layer.dim_output)

    rho1 = layer.param
    rho2 = layer.param_aux

    W_momentum = rho1 * layer.W_momentum + (1.0 - rho1) * W_grad # m
    W_momentum_aux = rho2 * layer.W_momentum_aux + (1.0 - rho2) * (W_grad) .^ 2 # v
    W_momentum_hat = W_momentum / (1.0 - rho1) # \hat{m}
    W_momentum_aux_hat = W_momentum_aux / (1.0 - rho2) # \hat{v}
    b_momentum = rho1 * layer.b_momentum + (1.0 - rho1) * b_grad # m
    b_momentum_aux = rho2 * layer.b_momentum_aux + (1.0 - rho2) * (b_grad) .^ 2 # v
    b_momentum_hat = b_momentum / (1.0 - rho1) # \hat{m}
    b_momentum_aux_hat = b_momentum_aux / (1.0 - rho2) # \hat{v}

    layer.W_momentum = W_momentum
    layer.W_momentum_aux = W_momentum_aux
    layer.b_momentum = b_momentum
    layer.b_momentum_aux = b_momentum_aux

    layer.W -= learning_rate * W_momentum_hat ./ sqrt.(W_momentum_aux_hat .+ 1e-8)
    layer.b -= learning_rate * b_momentum_hat ./ sqrt.(b_momentum_aux_hat .+ 1e-8)

    return transpose(W) * delta
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
# Iteration
# =========
#
# Types and methods for iterating over biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Ambiguous nucleotides iterator
# ------------------------------

immutable AmbiguousNucleicAcidIterator{A<:Union{DNAAlphabet,RNAAlphabet}}
    seq::BioSequence{A}
end

ambiguous_positions(seq::BioSequence) = AmbiguousNucleicAcidIterator(seq)

Base.start(it::AmbiguousNucleicAcidIterator) = find_next_ambiguous(it.seq, 1)
Base.done(it::AmbiguousNucleicAcidIterator, nextpos) = nextpos == 0
function Base.next(it::AmbiguousNucleicAcidIterator, nextpos)
    return nextpos, find_next_ambiguous(it.seq, nextpos + 1)
end

Base.iteratorsize(::AmbiguousNucleicAcidIterator) = Base.SizeUnknown()

function find_next_ambiguous{A<:Union{DNAAlphabet{2},RNAAlphabet{2}}}(
        seq::BioSequence{A}, i::Integer)
    # no ambiguity
    return 0
end

function find_next_ambiguous{A<:Union{DNAAlphabet{4},RNAAlphabet{4}}}(
        seq::BioSequence{A}, from::Integer)
    for i in max(from, 1):endof(seq)
        nt = inbounds_getindex(seq, i)
        if isambiguous(nt)
            return i
        end
    end
    # no ambiguity
    return 0
end


text:
# we need to parametrise by type so we know which library to call
abstract type Viewer{T} end

const CPetscViewer         = Ptr{Cvoid}

Base.cconvert(::Type{CPetscViewer}, obj::Viewer) = obj.ptr
Base.unsafe_convert(::Type{Ptr{CPetscViewer}}, obj::Viewer) =
    convert(Ptr{CPetscViewer}, pointer_from_objref(obj))

mutable struct ViewerStdout{T} <: Viewer{T}
    ptr::CPetscViewer
    comm::MPI.Comm
end

@for_libpetsc begin
    function ViewerStdout{$PetscScalar}(comm::MPI.Comm=MPI.COMM_SELF)
        ptr = ccall((:PETSC_VIEWER_STDOUT_, $libpetsc), CPetscViewer, (MPI.MPI_Comm,), comm)
        return ViewerStdout{$PetscScalar}(ptr, comm)
    end
    function Base.push!(viewer::Viewer{$PetscScalar}, format::PetscViewerFormat)
        @chk ccall((:PetscViewerPushFormat, $libpetsc), PetscErrorCode, (CPetscViewer,PetscViewerFormat), viewer, format)
        return nothing
    end
    function Base.pop!(viewer::Viewer{$PetscScalar})
        @chk ccall((:PetscViewerPopFormat, $libpetsc), PetscErrorCode, (CPetscViewer,), viewer)
        return nothing
    end

end

function with(f, viewer::Viewer, format::PetscViewerFormat)
    push!(viewer, format)
    try
        f()
    finally
        pop!(viewer)
    end
end

# ideally we would capture the output directly, but this looks difficult
# easiest option is to redirect stdout
# based on suggestion from https://github.com/JuliaLang/julia/issues/32567
function _show(io::IO, obj)
    old_stdout = stdout
    try
        rd, = redirect_stdout()
        view(obj)
        Libc.flush_cstdio()
        flush(stdout)
        write(io, readavailable(rd))
    finally
        redirect_stdout(old_stdout)
    end
    return nothing
end

#=
# PETSc_jll isn't built with X support
mutable struct ViewerDraw <: Viewer
    ptr::CPetscViewer
    comm::MPI.Comm
end
function ViewerDraw(comm::MPI.Comm)
    ptr = ccall((:PETSC_VIEWER_DRAW_, libpetsc), CPetscViewer, (MPI.MPI_Comm,), comm)
    return ViewerDraw(ptr, comm)
end
=#