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

using Markdown
using InteractiveUtils

# ╔═╡ ad4957b2-5292-11eb-1f64-c96793eb50a3
using DrWatson

# ╔═╡ b993ce0a-52a4-11eb-325d-9762114e31e8
using BenchmarkTools, DynamicalSystems, Neighborhood, Entropies, DataFrames, CSV

# ╔═╡ bbc53e10-53f9-11eb-1166-17295832370f
using Plots

# ╔═╡ 3660b24a-52a5-11eb-06a7-5b4bf9b5f359
md"# Setup"

# ╔═╡ 02804130-5295-11eb-2076-9dd3d942b26c
@quickactivate "FractalDimension"

# ╔═╡ f08543ee-52a4-11eb-1928-89a52d10ffdf
md"Import all needed functions for benchmarks from packages and by direct definition."

# ╔═╡ 00c03bce-52a5-11eb-1fa0-19ffb17b23dd
function _correlationsum_new(tree, X, ε, w, qf)
    theiler = Neighborhood.Theiler(w)
    searchtype = Neighborhood.WithinRange(maximum(ε))
    N = length(X)
    Cq = similar(ε, Float64)
    @inbounds for i in 1:N
        idxs = Neighborhood.isearch(tree, X[i], searchtype, theiler(i))
        isempty(idxs) && continue
        dist = [Neighborhood.evaluate(tree.metric, X[i], X[j]) for j in idxs]
        sort!(dist)
        for k in 1:length(ε)-1
            j = searchsortedfirst(dist, ε[k])
            Cq[k] += qf(j-1)
        end
        Cq[end] += qf(length(dist))
    end
    return Cq ./ (N*qf(N-w-1))
end

# ╔═╡ fa4e2bfc-52a4-11eb-088f-9bd8377cc691
function correlationsum_new(X, ε; q = 2, metric = Euclidean(), w = 0)
    ntype = WithinRange(maximum(ε))
    tree = KDTree(X)
    qf = q == 2 ? identity : x -> float(x)^(q-1)
    _correlationsum_new(tree, X, ε, w, qf)
end

# ╔═╡ 07c6fc50-52a5-11eb-37d6-e546a66a5c28
# Define a function to access data sets.
"""
	sample_data(data, N_sample)
Sample `N_sample` unique random points out of `data`.
"""
function sample_data(data, N_sample)
    N_data = length(data)
    possible_indices = 1:N_data
    indices = unique(rand(possible_indices, 2N_sample))[1:N_sample]
    data[indices] |> Dataset
end

# ╔═╡ 3f6eb602-52a5-11eb-1735-f9d1f28fbcd8
md"# Single Benchmark Functions
## Classic Correlation"

# ╔═╡ 11e7cb6a-52a5-11eb-16de-cd668a8473ff
function benchmark_correlation(data)
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = minimum(maxima($data) - minima($data))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        [correlationsum($data, ε) for ε in εs]
    end
    minimum(bm.times)
end

# ╔═╡ 5a1a167a-52a5-11eb-11a9-958f6cefa7b3
md"## New Correlation"

# ╔═╡ 162b224e-52a5-11eb-0351-955a17fd1524
function benchmark_new_correlation(data)
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = minimum(maxima($data) - minima($data))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        correlationsum_new($data, εs)
    end
    minimum(bm.times)
end

# ╔═╡ 68a515f2-52a5-11eb-29ec-310b27c74895
md"## Boxed Correlation"

# ╔═╡ 1a849636-52a5-11eb-11fe-61afc4a4b453
function benchmark_boxed_correlation(data::Dataset{D,T}, prismdim = D) where {D,T}
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = estimate_r0_buenoorovio($data, $prismdim)
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        boxed_correlationsum($data, εs, ε1, M = $prismdim)
    end
    minimum(bm.times)
end

# ╔═╡ 7da26766-52a5-11eb-2318-593d397abb4c
md"## Moltenos Generalized Entropy"

# ╔═╡ 1f5dbdb8-52a5-11eb-3814-33ca2ad86123
function benchmark_molteno(data)
    bm = @benchmark molteno_boxing($data)
    minimum(bm.times)
end

# ╔═╡ 8a690dba-52a5-11eb-0703-932f4a1ab648
md"## Standard Generalized Entropy"

# ╔═╡ 2384e768-52a5-11eb-38cb-15b49cb5fdcb
function benchmark_genentropy(data)
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = minimum(maxima($data) - minima($data))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
		[Entropies.genentropy($data, ε) for ε in εs]
	end
	minimum(bm.times)
end

# ╔═╡ 9663d82a-52a5-11eb-3fe8-3f753118ab60
md"# Actual Experiments
## Benchmark of Size
The values are spread logarithmically"

# ╔═╡ 289e88c6-52a5-11eb-113a-b78dd25b0154
"""
	benchmark_N(
		lower_number,
		upper_number,
		number_of_numbers;
		base = 10,
		feedback = false
	)
Benchmarks `number_of_numbers` values in the interval `base^lower_number`, `base^higher_number`. These values are distributed logarithmically. If `feedback` is `true` there is a text feedback.
"""
function benchmark_N(
	lower_number,
	upper_number,
	number_of_numbers;
	base = 10,
	feedback = false
)
	result = DataFrame(N = Int[], Model = [], Method = [], Time = [])
	numbers = floor.(Int, base .^ range(lower_number, upper_number, length = number_of_numbers))

	# Define Systems.
	henon_sys = Systems.henon(rand(2))
	lorenz_sys = Systems.lorenz(randn(3))
	N_traj = 5 * floor(base ^ upper_number)

	# Integrate trajectories.
	henon = trajectory(henon_sys, N_traj, Ttr = 100)
	lorenz = trajectory(lorenz_sys, N_traj * 0.1, dt = 0.1, Ttr = 10)

	trajectories = [(:henon, henon), (:lorenz, lorenz)]

	for (name, traj) in trajectories
		feedback && println("Now using the $name system.")
		for N in numbers
			feedback && println("Currently calculating with $N random data points.")
			data = sample_data(traj, N)
			feedback && println("Benchmarking the classic correlation…")
			push!(
				result,
				(N, name, "Correlation",benchmark_correlation(data))
			)
			feedback && println("Benchmarking the new correlation…")
			push!(
				result,
				(N, name, "New Correlation", benchmark_new_correlation(data))
			)
			feedback && println("Benchmarking the boxed correlation…")
			push!(
				result,
				(N, name, "Boxed Correlation", benchmark_boxed_correlation(data))
			)
			feedback && println("Benchmarking the generalized entropy with factor 2 boxing…")
			push!(
				result,
				(N, name, "Molteno", benchmark_molteno(data))
			)
			feedback && println("Benchmarking the generalized entropy with standard boxing…")
			push!(
				result,
				(N, name, "Generalized Dimension", benchmark_genentropy(data))
			)
		end
	end
	wsave(
		datadir(
			"benchmark",
			savename(
				"benchmark_number",
				Dict(
					:lowerN => lower_number,
					:upperN => upper_number,
					:NofN => number_of_numbers,
					:base => base,
				),
				"csv",
			),
		),
		result,
	)
	result
end

# ╔═╡ 0ab12ebc-52a6-11eb-300d-1f864d2fcdf3
md"## Benchmark of dimension"

# ╔═╡ 2f2d2cec-52a5-11eb-0748-0dc568cf81aa
"""
	benchmark_dimension(
		lower_dimension::Int = 4,
		upper_dimension::Int = 16,
		distance_dimension::Int = 1;
		N_sample = 10000,
		feedback = false,
	)
Benchmarks different correlation algorithms for dimensions with distance `distance_dimension` between `lower_dimension` and `upper_dimension`. These algorithms are compared for different distance_dimensions. If `feedback` is `true` a status is printed after each benchmark step.
"""
function benchmark_dimension(
	lower_dimension::Int = 4,
	upper_dimension::Int = 16,
	distance_dimension::Int = 1;
	N_sample = 10000,
	feedback = false,
)
	result = DataFrame(Dimension = Int[], Method = [], Time = [])
	dimensions = lower_dimension:distance_dimension:upper_dimension
	for dimension in dimensions
		feedback && println("The current dimension is $dimension.")
		# The current trajectroy is longer than the data sample taken from the trajectory.
		traj = trajectory(
			Systems.lorenz96(dimension, rand(dimension), F = 8.0),
			N_sample,
			Ttr = 100,
			dt = 0.1,
		)
		# Only 10000 unique data points are sampled from the 100000 trajectory points
		data = sample_data(traj, N_sample)
		# Benchmark the classic correlation.
		feedback && println("Now benchmarking the classic correlation…")
		push!(
			result,
			(dimension, "Classic Correlation", benchmark_correlation(data)),
		)
		# Benchmark the new correlation that uses KDTrees.
		feedback && println("Now benchmarking the new correlation…")
		push!(
			result,
			(dimension, "New Correlation", benchmark_new_correlation(data)),
		)
		# Benchmark the Boxed Correlation without accessing the prism property.
		feedback && println("Now benchmarking the boxed correlation…")
		push!(
			result,
			(dimension, "Boxed Correlation", benchmark_boxed_correlation(data)),
		)
		feedback && println("Now benchmarking the prism correlation…")
		# Benchmark the Prism Correlation with Prism Dimension by Theiler.
		embedding_dim =
			0.75log2(N_sample) < dimension ?
			(0.5log2(N_sample) |> ceil |> Int) :
			dimension
		push!(
			result,
			(dimension, "Prism Correlation", benchmark_boxed_correlation(data, embedding_dim)),
		)
		feedback && println("Now benchmarking the generalized entropy by Molteno…")
		# Benchmark the generalized entropy method by Molteno.
		push!(
			result,
			(dimension, "Molteno", benchmark_molteno(data))
		)
		feedback && println("Now benchmarking the classic generalized entropy…")
		# Benchmark the classic generalized entropy.
		push!(
			result,
			(dimension, "Generalized Entropy", benchmark_genentropy(data))
		)
	end
	wsave(
		datadir(
			"benchmark",
			savename(
				"benchmark_dim",
				Dict(
					:N => N_sample,
					:lowerdim => lower_dimension,
					:upperdim => upper_dimension,
					:ddim => distance_dimension,
				),
				"csv",
			),
		),
		result,
	)
	result
end

# ╔═╡ fa38cdc0-52c0-11eb-08ba-e3b2e399dc2f
data = trajectory(Systems.henon(0.1rand(2)), 1000000, Ttr = 10)

# ╔═╡ fa33702a-52c0-11eb-250b-9d7a7a7506ad
benchmark_correlation(sample_data(data, 1000))

# ╔═╡ 2a048274-52c1-11eb-339c-6b10fb972d11
benchmark_new_correlation(sample_data(data, 1000))

# ╔═╡ 33866012-52c1-11eb-329c-eff5928eb12e
benchmark_boxed_correlation(sample_data(data, 1000))

# ╔═╡ 40b6a85a-52c1-11eb-080a-272e1642425d
benchmark_molteno(sample_data(data, 1000))

# ╔═╡ 4c2583a0-52c1-11eb-35ec-49b3c59c471f
benchmark_genentropy(sample_data(data, 1000))

# ╔═╡ 5c2b6d64-52c1-11eb-14e5-bd81d666bb81
benchmark_boxed_correlation(sample_data(data, 1000), 1)

# ╔═╡ 0124b0f2-53f9-11eb-3ebb-f182fd11a129
md"# Tests"

# ╔═╡ 64d806be-53f9-11eb-2ca9-1d2ece3404bd
sdata = sample_data(data, 10000)

# ╔═╡ e7bf3d5c-5414-11eb-0556-452471c2b65b
q = 1.99

# ╔═╡ bd059f48-5414-11eb-08c2-794f627f15ee

# ╔═╡ bba5b5d8-5415-11eb-35b1-c3e070d1e09d
md"# Results
## `boxed_correlationsum` vs `correlationsum`
It seems like they behave the same for a given ε.

Say `henon` is a dataset with a trajectory of the henon map."

# ╔═╡ c6d939e8-5415-11eb-0c37-c73be604f29c
henon = trajectory(Systems.henon(0.1rand(2)), 10000, Ttr = 22)

# ╔═╡ 0d016116-5416-11eb-3ab5-a534b0d1453a
md"Then we can choose the εs by selecting the upper border `r0` using `estimate_r0_buenoorovio` and the lower border `ε_min` via `min_pairwise_distance` and spread them logarithmically on that interval."

# ╔═╡ 6f5e0968-5416-11eb-0fcb-572aaa902574
r0 = estimate_r0_buenoorovio(henon)

# ╔═╡ a0f41bf2-5416-11eb-2da6-cb4b485084c5
ε_min = min_pairwise_distance(henon)[2]

# ╔═╡ a2ca2444-5416-11eb-0a40-1f65148b816e
εs = 10 .^ range(log10(ε_min), log10(r0), length = 15)

# ╔═╡ 0798e090-53f9-11eb-377a-fdf6c2ade4d1
bcs = boxed_correlationsum(sdata, εs, q = q)

# ╔═╡ 88515afc-5415-11eb-273c-83e703c7392d
cs = correlationsum(sdata, εs, q = q)

# ╔═╡ 9fc5c036-5415-11eb-1d0e-5377b9d92943
bcs ≈ cs

# ╔═╡ be07006c-53f9-11eb-0141-4f9afae62d4b
plot(εs .|> log10, bcs .|> log10)

# ╔═╡ e89063f6-53fa-11eb-2671-5bacb746b401
linear_region(εs .|> log10, bcs .|> log10)[2]*1/(q-1)

# ╔═╡ efe5247c-5416-11eb-235e-6d7ce045b31e
md"Now the two correlation sums can be calculated."

# ╔═╡ fc2bd50a-5416-11eb-38f5-97335565a8eb
correlation_2 = correlationsum(henon, εs)

# ╔═╡ 067a4be0-5417-11eb-2774-915824fb7987
boxed_correlation_2 = boxed_correlationsum(henon, εs)

# ╔═╡ cbcac542-5416-11eb-3b99-791d68ce0745
md"Comparing these two directly works."

# ╔═╡ 1b3d11c0-5417-11eb-0856-b5aa0be46b0b
correlation_2 == boxed_correlation_2

# ╔═╡ 3842c2ec-5417-11eb-0030-83c97d87e157
md"The algorithm differs for `q≠2` so we need to compare this case as well."

# ╔═╡ 4d1e7eea-5417-11eb-1021-d50ee8615148
correlation_q = correlationsum(henon, εs, q = 3.3)

# ╔═╡ 61abb6b6-5417-11eb-32a7-f74a2549ef42
boxed_correlation_q = boxed_correlationsum(henon, εs, q = 3.3)

# ╔═╡ 70a0118a-5417-11eb-12cf-016cdbccb50d
correlation_q == boxed_correlation_q

# ╔═╡ 81510afc-5417-11eb-3111-6dfc809e8f6b
md"This does not come unexpected. The sums are only expected to be approximately equal, since the precision of the types plays a role."

# ╔═╡ c8c092ea-5417-11eb-2175-b3a1c62f1454
correlation_q ≈ boxed_correlation_q

# ╔═╡ 8f48f42e-5419-11eb-1865-c19b60d19b02
begin
	import ChaosTools: boxed_correlationdim
	function boxed_correlationdim(data, εs, r0 = maximum(εs); q = 2, M = size(data, 2))
		@assert M ≤ size(data,2) "Prism dimension has to be lower or equal than " *
		"data dimension."
		if isinf(q)
			dd = boxed_correlationsum(data, εs, r0; q = 3, M = M)
			return 1 / (3 - 1) * linear_region(log.(εs), log.(maximum.(dd)), tol = 0.1)[2]
		else
			dd = boxed_correlationsum(data, εs, r0; q = q, M = M)
			return 1 / (q - 1) * linear_region(log.(εs), log.(dd), tol = 0.1)[2]
		end
	end
end

# ╔═╡ 420e26b2-5418-11eb-0477-5fd4284ff1d0
boxed_correlationdim(henon, q = 3)

# ╔═╡ e3627ac8-5417-11eb-1cf7-e13f24c1b62f
md"## q-order correlation literature
I was not able to find any paper using the explicit q-order correlationsum, yet its definition seems sensible and the plots I retrieve show linear scaling regions for different `q`."

# ╔═╡ d0e25b7c-5424-11eb-04dd-396f8e034b86
begin
	ϵ_min = min_pairwise_distance(henon)[2]
	ϵ_max = maximum(maxima(henon)-minima(henon))
	ϵs = 10 .^ range(log10(ϵ_min), log10(ϵ_max), length = 20)
	plot()
	for q in 1.5:0.5:6
		corsum = correlationsum(henon, ϵs, q = q)
		(p1, p2), dim = linear_region(log10.(ϵs), log10.(corsum))
		dim = 1/(q-1)*dim
		plot!(
			log10.(ϵs),
			log10.(corsum),
			label = "q = $q, Δ = $(round(dim, digits = 3))",
		)
	end
	plot!(
		legend = :bottomright,
		xlims = log10.((ϵs[1], ϵs[20])),
		xguide = "\$ \\log(\\epsilon) \$",
		yguide = "\$ C_q \$",
	)
end

# ╔═╡ Cell order:
# ╟─3660b24a-52a5-11eb-06a7-5b4bf9b5f359
# ╠═ad4957b2-5292-11eb-1f64-c96793eb50a3
# ╠═02804130-5295-11eb-2076-9dd3d942b26c
# ╟─f08543ee-52a4-11eb-1928-89a52d10ffdf
# ╠═b993ce0a-52a4-11eb-325d-9762114e31e8
# ╠═fa4e2bfc-52a4-11eb-088f-9bd8377cc691
# ╠═00c03bce-52a5-11eb-1fa0-19ffb17b23dd
# ╠═07c6fc50-52a5-11eb-37d6-e546a66a5c28
# ╟─3f6eb602-52a5-11eb-1735-f9d1f28fbcd8
# ╠═11e7cb6a-52a5-11eb-16de-cd668a8473ff
# ╟─5a1a167a-52a5-11eb-11a9-958f6cefa7b3
# ╠═162b224e-52a5-11eb-0351-955a17fd1524
# ╟─68a515f2-52a5-11eb-29ec-310b27c74895
# ╠═1a849636-52a5-11eb-11fe-61afc4a4b453
# ╟─7da26766-52a5-11eb-2318-593d397abb4c
# ╠═1f5dbdb8-52a5-11eb-3814-33ca2ad86123
# ╟─8a690dba-52a5-11eb-0703-932f4a1ab648
# ╠═2384e768-52a5-11eb-38cb-15b49cb5fdcb
# ╟─9663d82a-52a5-11eb-3fe8-3f753118ab60
# ╠═289e88c6-52a5-11eb-113a-b78dd25b0154
# ╟─0ab12ebc-52a6-11eb-300d-1f864d2fcdf3
# ╠═2f2d2cec-52a5-11eb-0748-0dc568cf81aa
# ╠═fa38cdc0-52c0-11eb-08ba-e3b2e399dc2f
# ╠═fa33702a-52c0-11eb-250b-9d7a7a7506ad
# ╠═2a048274-52c1-11eb-339c-6b10fb972d11
# ╠═33866012-52c1-11eb-329c-eff5928eb12e
# ╠═40b6a85a-52c1-11eb-080a-272e1642425d
# ╠═4c2583a0-52c1-11eb-35ec-49b3c59c471f
# ╠═5c2b6d64-52c1-11eb-14e5-bd81d666bb81
# ╟─0124b0f2-53f9-11eb-3ebb-f182fd11a129
# ╠═64d806be-53f9-11eb-2ca9-1d2ece3404bd
# ╠═e7bf3d5c-5414-11eb-0556-452471c2b65b
# ╠═0798e090-53f9-11eb-377a-fdf6c2ade4d1
# ╠═88515afc-5415-11eb-273c-83e703c7392d
# ╠═9fc5c036-5415-11eb-1d0e-5377b9d92943
# ╠═bbc53e10-53f9-11eb-1166-17295832370f
# ╠═be07006c-53f9-11eb-0141-4f9afae62d4b
# ╠═e89063f6-53fa-11eb-2671-5bacb746b401
# ╠═bd059f48-5414-11eb-08c2-794f627f15ee
# ╟─bba5b5d8-5415-11eb-35b1-c3e070d1e09d
# ╠═c6d939e8-5415-11eb-0c37-c73be604f29c
# ╟─0d016116-5416-11eb-3ab5-a534b0d1453a
# ╠═6f5e0968-5416-11eb-0fcb-572aaa902574
# ╠═a0f41bf2-5416-11eb-2da6-cb4b485084c5
# ╠═a2ca2444-5416-11eb-0a40-1f65148b816e
# ╟─efe5247c-5416-11eb-235e-6d7ce045b31e
# ╠═fc2bd50a-5416-11eb-38f5-97335565a8eb
# ╠═067a4be0-5417-11eb-2774-915824fb7987
# ╟─cbcac542-5416-11eb-3b99-791d68ce0745
# ╠═1b3d11c0-5417-11eb-0856-b5aa0be46b0b
# ╟─3842c2ec-5417-11eb-0030-83c97d87e157
# ╠═4d1e7eea-5417-11eb-1021-d50ee8615148
# ╠═61abb6b6-5417-11eb-32a7-f74a2549ef42
# ╠═70a0118a-5417-11eb-12cf-016cdbccb50d
# ╟─81510afc-5417-11eb-3111-6dfc809e8f6b
# ╠═c8c092ea-5417-11eb-2175-b3a1c62f1454
# ╠═8f48f42e-5419-11eb-1865-c19b60d19b02
# ╠═420e26b2-5418-11eb-0477-5fd4284ff1d0
# ╟─e3627ac8-5417-11eb-1cf7-e13f24c1b62f
# ╟─d0e25b7c-5424-11eb-04dd-396f8e034b86

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
## Calibrate camera
    # Initial guess at reference camera

    # Richard Shenton email:
    # Camera height: 2.065m
    #
    # Camera offset from centre line: 0.6m
    # Pitch: 19.16 degrees down
    # Yaw: 7.28 degrees to the left
    # Roll: Assumed to be zero

    # The camera frame rate is 25fps for this data. Regarding train details, there are some details here https://en.wikipedia.org/wiki/British_Rail_Class_357 .
    #
    # It looks to me that the carriage with the cab (seems to be called DMOS) is 20.75(L) x 2.8 (W) x 3.78(H) and weighs 40.7 t.

using TrainView

trackprop = TrackProperties(track_gauge = 1.435 + 0.065)

cameraposition_reference = [0.0,0.6,-2.165];
ψθφ_ref = [0.0,-19.1,-7.28] .* (pi/180.0)

camera_reference = VideoCamera(cameraposition_reference;
    focal_length = 5.8e-3,
    pixelspermeter = 1 / 5.5e-6,
    ψθφ = ψθφ_ref
)

using CSV, DataFrames
using Statistics, LinearAlgebra
using Plots
gr(yflip = true)
GR.inline("png")

# args = ["../../data/output_results_centre_line.csv", "output.csv"]

file = "../../data/output_results_centre_line.csv";
# df = CSV.read(file, DataFrame);
df = CSV.File(file;

lines = open(file) do f
    [l for l in eachline(f)]
end;

header = lines[1];
lines = lines[2:end];

max_v = 720
max_u = 1280

uv_data = map(lines) do l
    data = parse.(Int,split(l,',')[2:end-1])

    Llen = data[1]
    Rlen = data[2]

    us = data[3:2:end] .- (max_u / 2.0);
    vs = data[4:2:end] .- (max_v / 2.0);

    Lu = us[1:Llen];
    Lv = vs[1:Llen];

    Ru = us[Llen+1:Rlen+Llen];
    Rv = vs[Llen+1:Rlen+Llen];

    return [Lu,Lv,Ru,Rv]
end;

using Plots
gr(yflip = true)

i = 100
Lu = uv_data[i][1];
Lv = uv_data[i][2];
Ru = uv_data[i][3];
Rv = uv_data[i][4];

Lu0 = left_track_image_u(camera_reference, trackprop, Lv);
Ru0 = right_track_image_u(camera_reference, trackprop, Rv);

plot(Lu,Lv, lab = "")
plot!(Ru,Rv, lab = "")

plot!(Lu0,Lv, lab = "")
plot!(Ru0,Rv, lab = "")

choose_distortions = [:Y,:Z,:θ,:φ];
choose_distortions = [:Y,:Z,:θ,:φ,:α,:β];

skip_distortion = repeat([0.0],length(choose_distortions))

distortions = map(uv_data) do uv

    left_uvs = map(uv[1],uv[2]) do u, v
        [u,v]
    end

    right_uvs = map(uv[3],uv[4]) do u, v
        [u,v]
    end

    distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop;
        choose_distortions = choose_distortions,
        iterations = 4)

    if abs(distortion[:Y]) > 0.1 + 0.5*abs(camera_reference.xyz[2]) || abs(distortion[:Z]) >= 0.5*abs(camera_reference.xyz[3])
        return skip_distortion
    else
        return [distortion[k] for k in choose_distortions]
    end
end;

filter!(d ->  norm(d - skip_distortion) > 0.0,distortions);

mean_distortion = mean(distortions);
std_distortion = std(distortions);

std_theshold = std_distortion * 0.4
calibrate_inds = [isempty(findall(abs.(mean_distortion - d) .> std_theshold)) for d in distortions];
count(calibrate_inds)

mean_distortion2 = mean(distortions[calibrate_inds]);
camera = VideoCamera(camera_reference, Dict(choose_distortions .=> mean_distortion2))

Lu2 = left_track_image_u(camera, trackprop, Lv);
Ru2 = right_track_image_u(camera, trackprop, Rv);

plot(Lu,Lv, lab = "")
plot!(Ru,Rv, lab = "")

plot!(Lu2,Lv,linestyle=:dash, lab = "")
plot!(Ru2,Rv,linestyle=:dash, lab = "")

# Create a gif
    n = 1;
    uv = uv_data[n]
    ulims = 1.5 .* (min(minimum(uv[1]),minimum(uv[3])),max(maximum(uv[1]),maximum(uv[3])))
    vlims = 1.5 .* (min(minimum(uv[2]),minimum(uv[4])),max(maximum(uv[2]),maximum(uv[4])))

    Lv = uv[2]; Rv = uv[4];
    Lu = left_track_image_u(camera, trackprop, Lv);
    Ru = right_track_image_u(camera, trackprop, Rv);

    plot(Lu,Lv,linestyle=:dash, lab = "", color = :blue)
    plot!(Ru,Rv,linestyle=:dash, lab = "", color = :blue)

    maxframes = length(uv_data)

    # anim = @animate for n in Int.(round.(LinRange(1,maxframes,600)))
    #     uv = uv_data[n];
    #
    #     plot(uv[1],uv[2], lab = "", color = :red )
    #     plot!(uv[3],uv[4], lab = "", color = :red )
    #
    #     plot!(Lu,Lv,linestyle=:dash, lab = "", color = :blue)
    #     plot!(Ru,Rv,linestyle=:dash, lab = "", color = :blue)
    #
    #     plot!(camera, xlims = ulims, ylims = vlims, xguide = "", yguide = "",axis = false)
    #
    # end
    #
    # # gif(anim, "../../images/videos/only-curvature.gif", fps = 7)
    # gif(anim, "../images/videos/compare-reference-to-track-center-line.gif", fps = 7)

# Calculate train car movement from frames
    choose_distortions = [:Y,:Z,:θT,:φT,:α,:β];

    distortions = map(uv_data) do uv

        left_uvs = map(uv[1],uv[2]) do u, v
            [u,v]
        end

        right_uvs = map(uv[3],uv[4]) do u, v
            [u,v]
        end

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera, trackprop;
            choose_distortions = choose_distortions,
            iterations = 1)

        return [distortion[k] for k in choose_distortions]
    end;

# create gif of distortions
    n = 1;
    uv = uv_data[n]
    ulims = 1.9 .* (min(minimum(uv[1]),minimum(uv[3])),max(maximum(uv[1]),maximum(uv[3])))
    vlims = 1.2 .* (min(minimum(uv[2]),minimum(uv[4])),max(maximum(uv[2]),maximum(uv[4])))

    Lv = uv[2]; Rv = uv[4];
    Lu = left_track_image_u(camera, trackprop, Lv);
    Ru = right_track_image_u(camera, trackprop, Rv);

    plot(Lu,Lv,linestyle=:dash, lab = "", color = :blue)
    plot!(Ru,Rv,linestyle=:dash, lab = "", color = :blue)

    plot!(camera, xlims = ulims, ylims = vlims, xguide = "", yguide = "",axis = false)

    Lu_ref = left_track_image_u(camera, trackprop, Lv);
    Ru_ref = right_track_image_u(camera, trackprop, Rv);

    ΔuL = left_track_Δu(camera, trackprop, Lv; choose_distortions=choose_distortions);
    ΔuR = right_track_Δu(camera, trackprop, Rv; choose_distortions=choose_distortions);

    maxframes = length(uv_data)

    anim = @animate for n in Int.(round.(LinRange(1,maxframes,500)))
        uv = uv_data[n];

        Lu = Lu_ref + ΔuL * distortions[n];
        Ru = Ru_ref + ΔuR * distortions[n];

        plot(uv[1],uv[2], lab = "", color = :red )
        plot!(uv[3],uv[4], lab = "recorded data", color = :red )

        plot!(Lu,Lv,linestyle=:dash, lab = "", color = :blue)
        plot!(Ru,Rv,linestyle=:dash, lab = "modeled track", color = :blue)

        p = plot!(camera, xlims = ulims, ylims = vlims,
            xguide = "", yguide = "", axis = nothing,
            xlab = "pixel u", ylab = "pixel v",
            legend = :topleft)
    end

    gif(anim, "../../images/videos/model-track-center-line.gif", fps = 7)

## gif that tracks distortions

    curve_distortions = [:α,:β]

    ΔuL = left_track_Δu(camera, trackprop, Lv; choose_distortions=curve_distortions);
    ΔuR = right_track_Δu(camera, trackprop, Rv; choose_distortions=curve_distortions);

    indαβ = [findfirst(s .== choose_distortions) for s in curve_distortions]
    indYZ = [findfirst(s .== choose_distortions) for s in [:Y,:Z]]
    indθT = [findfirst(s .== choose_distortions) for s in [:θT,:φT]]

    maxYZ = 1.5 .* mean([abs.(d[indYZ]) for d in distortions])
    maxθT = 5 .* mean([abs.(d[indθT]) for d in distortions])

    anim = @animate for n in Int.(round.(LinRange(1,maxframes,500)))
        uv = uv_data[n];

        Lu = Lu_ref + ΔuL * distortions[n][indαβ];
        Ru = Ru_ref + ΔuR * distortions[n][indαβ];

        plot(uv[1],uv[2], lab = "", color = :red )
        plot!(uv[3],uv[4], lab = "recorded data", color = :red )

        plot!(Lu,Lv,linestyle=:dash, lab = "", color = :blue)
        plot!(Ru,Rv,linestyle=:dash, lab = "model curvature", color = :blue)

        p = plot!(camera, xlims = ulims, ylims = vlims,
            xguide = "", yguide = "", axis = nothing,
            xlab = "pixel u", ylab = "pixel v",
            legend = :topleft)

        dYdZ = distortions[n][indYZ]
        dθdφ = distortions[n][indθT]

        pYZ = scatter([-dYdZ[1]],[-dYdZ[2]],
            label = "", title = "camera position",
            markersize = 6, markercolor = :green,
            xlims = (-maxYZ[1],maxYZ[1]),
            ylims = (-maxYZ[2],maxYZ[2]),
            axis = nothing,
            xlab = "y", ylab = "z"
        )
        ptheta = scatter([dθdφ[1]],[dθdφ[2]],
            label = "", title = "cabin orientation",
            markersize = 6, markercolor = :purple,
            xlims = (-maxθT[1],maxθT[1]),
            ylims = (-maxθT[2],maxθT[2]),
            axis = nothing,
            xlab = "θ", ylab = "φ"
        )

        l = @layout [
            a{0.6w,1.0h} [b{0.5h}
                     c{0.5h}  ]
        ]
        plot(p,pYZ,ptheta, layout = l)
    end

    gif(anim, "../../images/videos/track-distortions-track-center-line.gif", fps = 5)

    # p1 = plot(x, y) # Make a line plot
    # p2 = scatter(x, y) # Make a scatter plot
    # p3 = plot(x, y, xlabel = "This one is labelled", lw = 3, title = "Subtitle")
    # p4 = histogram(x, y) # Four histograms each with 10 points? Why not!
    # plot(p1, p2, p3, p4, layout = (2, 2), legend = false)

ext: jl


text:
module CureMIDI

import FluidSynth
import SampledSignals
import MIDI

export synth
export tick_to_frame, tick_to_ms, ms_to_tick

"""
Synthesize a block of 16 bit audio samples to audio buffers.
# Arguments
- `track::MIDI.MIDITrack`
- `tpq::Int16`
- 'bpm::Real'
- 'sample_rate::Real'
- 'sound_font::AbstractString'
"""
function synth(track::MIDI.MIDITrack, tpq::Int16, bpm::Real, sample_rate::Real, sound_font::AbstractString)

    settings = FluidSynth.Settings()
    synth = FluidSynth.Synth(settings)
    # settings
    FluidSynth.set_sample_rate(synth, Float32(sample_rate))

    # load sound font
    sfont_id = FluidSynth.sfload(synth, sound_font)
    FluidSynth.program_select(synth, Int32(0), sfont_id, Int32(0), Int32(0))

    # allocate full buffer
    full_ticks = 0
    for note in  MIDI.getnotes(track)
        full_ticks = max(full_ticks, note.position + note.duration)
    end

    full_frames = tick_to_frame(full_ticks, tpq, bpm, sample_rate)
    left_full_buf = zeros(Float32, full_frames)
    right_full_buf = zeros(Float32, full_frames)

    left_note_buf = zeros(Float32, full_frames)
    right_note_buf = zeros(Float32, full_frames)
    for note in MIDI.getnotes(track)
        note_frames = tick_to_frame(note.position, tpq, bpm, sample_rate)
        fill!(left_note_buf, 0)
        fill!(right_note_buf, 0)

        channel = Int32(note.channel)
        pitch = Int32(note.pitch)
        velocity = Int32(note.velocity)
        FluidSynth.noteon(synth, channel, pitch, velocity)
        start_frames = tick_to_frame(note.position, tpq, bpm, sample_rate)
        duration_frames = tick_to_frame(note.duration, tpq, bpm, sample_rate)

        FluidSynth.write_float(synth, Int32(duration_frames), left_note_buf, Int32(start_frames) , Int32(1), right_note_buf, Int32(start_frames), Int32(1))
        try
            FluidSynth.noteoff(synth, channel, pitch)
        catch
            # sometimes, noteoff fail.
            # But not affect playing. so ignore error.
        end

        left_full_buf += left_note_buf
        right_full_buf += right_note_buf
    end

    buf = hcat(left_full_buf, right_full_buf)
    sampled_buf = SampledSignals.SampleBuf(buf, sample_rate)
    return sampled_buf
end

# support functions
"""
Convert from tick to synthsized result buffer index.
# Arguments
- 'tick::Uint' : ticks.
- `tpq::Int16` : Ticks per quarter note.
- `bpm::Real` : Beat per min. Beat means quater notes in this library.
- `sample_rate::Real` : Beat per min. Beat means quater notes in this library.
"""
function tick_to_frame(tick::UInt, tpq::Int16, bpm::Real, sample_rate::Real)
    ms_per_tick = MIDI.ms_per_tick(tpq, bpm)
    full_ms = ms_per_tick * tick
    frames = sample_rate * full_ms / 1000
    return UInt(ceil(frames))
end

"""
Convert from tick to millisecond.
# Arguments
- 'tick::Uint' : ticks.
- `tpq::Int16` : Ticks per quarter note.
- `bpm::Real` : Beat per min. Beat means quater notes in this library.
"""
function tick_to_ms(tick::UInt, tpq::Int16, bpm::Real)
    ms_per_tick = MIDI.ms_per_tick(tpq, bpm)
    return  tick * ms_per_tick
end

"""
Convert from millisecond to tick.
# Arguments
- 'ms::Real' : milliseconds.
- `tpq::Int16` : Ticks per quarter note.
- `bpm::Real` : Beat per min. Beat means quater notes in this library.
"""
function ms_to_tick(ms::Real, tpq::Int16, bpm::Real)
    ms_per_tick = MIDI.ms_per_tick(tpq, bpm)
    return  ms / ms_per_tick
end

end #module


text:
module builder

export Builder

include("load_kwant.jl")
import .. AbstractKwantObject

struct Builder <:AbstractKwantObject
    o::PyObject
    Builder() = new(kwant.Builder())
    Builder(x;kwargs...) = new(kwant.Builder(x;kwargs...))
end

Base.setproperty!(b::Builder, name::Symbol, x::AbstractKwantObject) =
    setproperty!(b.o,name,x.o)

Base.setindex!(b::Builder,val,keys...) = set!(b.o,keys,val)
Base.setindex!(b::Builder,val,key) = set!(b.o,key,val)
Base.setindex!(b::Builder,val,key::AbstractKwantObject) = set!(b.o,key.o,val)

struct HoppingKind <: AbstractKwantObject
    o::PyObject
    (hk::HoppingKind)(syst::Builder) = hk.o(syst)
    function HoppingKind(a::Tuple,b::AbstractKwantObject,c::AbstractKwantObject)
        new(pycall(kwant.builder.HoppingKind,PyObject,a,b.o,c.o))
    end
    # function HoppingKind(a::Tuple,b::PyObject,c::PyObject)
        # new(pycall(kwant.builder.HoppingKind,PyObject,a,b,c))
    # end
    HoppingKind(a::Tuple) = HoppingKind(a...)
end

function Base.deleteat!(b::Builder,o::PyObject)
    py"""
    del $(b.o)[$o]
    """
    return nothing
end

function Base.setindex!(b::Builder,val,keys::Array{HoppingKind})
    for k ∈ keys
        b[k]=val
    end
    return nothing
end

end # module

# struct FiniteSystem <:AbstractKwantObject o::PyObject end


text:
# Test codegencontext
@testset "CodeGenContext" begin
    @test_throws ArgumentError Automa.CodeGenContext(generator=:fdjfhkdj)
    @test_throws ArgumentError Automa.CodeGenContext(generator=:simd)
    @test_throws ArgumentError Automa.CodeGenContext(generator=:simd, checkbounds=false, loopunroll=2)
    @test_throws ArgumentError Automa.CodeGenContext(generator=:simd, checkbounds=false, getbyte=identity)
end

import Automa
const re = Automa.RegExp
import Automa.RegExp: @re_str

@testset "SIMD generator" begin
    machine = let
        seq = re"[A-Z]+"
        name = re"[a-z]+"
        rec = re">"  * name * re"\n" * seq
        Automa.compile(re.opt(rec) * re.rep(re"\n" * rec))
    end

    context = Automa.CodeGenContext(generator=:simd, checkbounds=false)

    @eval function is_valid_fasta(data::String)
        $(Automa.generate_init_code(context, machine))
        p_end = p_eof = ncodeunits(data)
        $(Automa.generate_exec_code(context, machine, nothing))
        return p == ncodeunits(data) + 1
    end

    s1 = ">seq\nTAGGCTA\n>hello\nAJKGMP"
    s2 = ">seq1\nTAGGC"
    s3 = ">verylongsequencewherethesimdkicksin\nQ"

    for (seq, isvalid) in [(s1, true), (s2, false), (s3, true)]
        @test is_valid_fasta(seq) == isvalid
    end
end


text:
# # Dynamic Von Mises Truss

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/dynamic_von_mises_truss.ipynb)
#

# ![](../assets/diagramDynVonMises.png)

using StructuralDynamicsODESolvers, Plots


text:
function Ground(;name)
    @named port = PositiveMagneticPort()
    eqs = [port.V_m ~ 0]
    ODESystem(eqs, t, [], [], systems=[port], name=name)
end

function Idle(;name)
    @named two_port = TwoPort()
    @unpack Phi = two_port
    eqs = [
        Phi ~ 0,
    ]
    extend(ODESystem(eqs, t, [Phi], [], systems=[], name=name), two_port)
end

function Short(;name)
    port_p = PositiveMagneticPort()
    port_n = NegativeMagneticPort()
    eqs = [
        connect(port_p, port_n),
    ]
    ODESystem(eqs, t, [], [], systems=[port_p, port_n], name=name)
end

function Crossing(;name)
    port_p1 = PositiveMagneticPort()
    port_p2 = PositiveMagneticPort()
    port_n1 = NegativeMagneticPort()
    port_n2 = NegativeMagneticPort()
    eqs = [
        connect(port_p1, port_p2),
        connect(port_n1, port_n2),
    ]
    ODESystem(eqs, t, [], [], systems=[port_p1, port_p2, port_n1, port_n2], name=name)
end

function ConstantPermeance(;name, G_m=1.0)
    val = G_m
    @named two_port = TwoPort()
    @unpack V_m, Phi = two_port
    @variables G_m(t)
    eqs = [
        Phi ~ G_m * V_m,
    ]
    extend(ODESystem(eqs, t, [V_m, Phi, G_m], [], systems=[], defaults=Dict(G_m => val), name=name), two_port)
end

function ConstantReluctance(;name, R_m=1.0)
    val = R_m
    @named two_port = TwoPort()
    @unpack V_m, Phi = two_port
    @variables R_m(t)
    eqs = [
        V_m ~ Phi * R_m,
    ]
    extend(ODESystem(eqs, t, [V_m, Phi, R_m], [], systems=[], defaults=Dict(R_m => val), name=name), two_port)
end


text:
#=##############################################################################
# DESCRIPTION
    Vehicle with vortex-lattice method (VLM) and VPM-based propeller models
    shedding VPM wakes.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

################################################################################
# UNSTEADY VLM VEHICLE TYPE
################################################################################
"""
    `UVLMVehicle(system; optargs...)`

Type handling all geometries and subsystems that define a flight vehicle made
out of VLM (Wing, WingSystem, Rotor) components.

# ARGUMENTS
* `system::vlm.WingSystem`:        System of all FLOWVLM objects. This system
                                    is considered as the entire vehicle. Not all
                                    components in this system will be solved,
                                    but they will all be rotated and translated
                                    during maneuver.
# OPTIONAL ARGUMENTS
* `tilting_systems::Tuple(vlm.WingSystem, ...)`:   Tuple of all FLOWVLM
                                    tilting objects, where `tilting_systems[i]`
                                    contains the i-th FLOWVLM system of lifting
                                    surfaces and rotors that tilt together.
* `rotors_systems::Tuple(Array{vlm.Rotor,1}, ...)`:   Tuple of groups of Rotors
                                    that share a common RPM.
* `vlm_system::vlm.WingSystem`:    System of all FLOWVLM objects to be solved
                                    through the VLM solver.
* `wake_system::vlm.WingSystem`:   System of all FLOWVLM objects that will
                                    shed a VPM wake.
* `grids::Array{gt.GridTypes, 1}`: Array of grids that will be translated and
                                    rotated along with `system`.
"""
struct UVLMVehicle{N, M, R} <: AbstractVLMVehicle{N, M, R}

    # Required inputs
    system::vlm.WingSystem

    # Optional inputs
    tilting_systems::NTuple{N, vlm.WingSystem}
    rotor_systems::NTuple{M, Array{vlm.Rotor, 1}}
    vlm_system::vlm.WingSystem
    wake_system::vlm.WingSystem
    grids::Array{gt.GridTypes, 1}

    # Internal properties
    V::Array{R, 1}                          # Current vehicle velocity
    W::Array{R, 1}                          # Current vehicle angular velocity
    prev_data::Array{Any, 1}                # Information about previous step
    grid_O::Array{Array{R, 1}, 1}           # Origin of every grid

    UVLMVehicle{N, M, R}(
                    system;
                    tilting_systems=NTuple{0, vlm.WingSystem}(),
                    rotor_systems=NTuple{0, Array{vlm.Rotor, 1}}(),
                    vlm_system=vlm.WingSystem(),
                    wake_system=vlm.WingSystem(),
                    grids=Array{gt.GridTypes, 1}(),
                    V=zeros(3), W=zeros(3),
                    prev_data=[deepcopy(vlm_system), deepcopy(wake_system),
                                                    deepcopy(rotor_systems)],
                    grid_O=Array{Array{Float64, 1}, 1}(),
                ) where {N, M, R} = new(
                    system,
                    tilting_systems,
                    rotor_systems,
                    vlm_system,
                    wake_system,
                    grids,
                    V, W,
                    prev_data,
                    grid_O,
                )
end

# Implicit N and M constructor
UVLMVehicle(system::vlm.WingSystem;
        V::Array{R, 1}=zeros(3), W::Array{R, 1}=zeros(3),
        tilting_systems::NTuple{N, vlm.WingSystem}=NTuple{0, vlm.WingSystem}(),
        rotor_systems::NTuple{M, Array{vlm.Rotor, 1}}=NTuple{0, Array{vlm.Rotor, 1}}(),
        grids=Array{gt.GridTypes, 1}(),
        optargs...
        ) where {N, M, R} = UVLMVehicle{N, M, R}( system;
                                V=V, W=W,
                                tilting_systems=tilting_systems,
                                rotor_systems=rotor_systems,
                                grids=Array{gt.GridTypes, 1}(grids),
                                grid_O=[zeros(R, 3) for i in 1:length(grids)],
                                optargs...)

# Alias
const VLMVehicle = UVLMVehicle

##### FUNCTIONS  ###############################################################
function shed_wake(self::VLMVehicle, Vinf::Function,
                            pfield::vpm.ParticleField, dt::Real, nt::Int; t=0.0,
                            unsteady_shedcrit=-1.0, p_per_step=1,
                            sigmafactor=1.0, overwrite_sigma=nothing)
    if nt!=0
        VLM2VPM(self.wake_system, pfield, dt, Vinf; t=t,
                    prev_system=_get_prev_wake_system(self),
                    unsteady_shedcrit=unsteady_shedcrit,
                    p_per_step=p_per_step, sigmafactor=sigmafactor,
                    overwrite_sigma=overwrite_sigma, check=false)
    end
end

function generate_static_particle_fun(pfield::vpm.ParticleField,
                                                self::VLMVehicle, sigma::Real)

    if sigma<=0
        error("Invalid smoothing radius $sigma.")
    end

    function static_particles_function(pfield, args...)

        # Particles from vlm system
        _static_particles(pfield, self.vlm_system, sigma)

        # Particles from rotor systems
        for rotors in self.rotor_systems
            for rotor in rotors
                _static_particles(pfield, rotor, sigma)
            end
        end

        return nothing
    end

    return static_particles_function
end

save_vtk(self::VLMVehicle, args...;
                        optargs...) = save_vtk_base(self, args...; optargs...)

##### INTERNAL FUNCTIONS  ######################################################
function _static_particles(pfield::vpm.ParticleField,
                            system::Union{vlm.Wing, vlm.WingSystem, vlm.Rotor},
                            sigma::Real)

    # Adds a particle for every bound vortex of the VLM
    for i in 1:vlm.get_m(system)
        (Ap, A, B, Bp, _, _, _, Gamma) = vlm.getHorseshoe(system, i)
        for (i,(x1, x2)) in enumerate(((Ap,A), (A,B), (B,Bp)))
            vpm.add_particle(pfield, (x1+x2)/2, Gamma*(x2-x1), sigma; vol=0)
        end
    end
end
##### END OF VEHICLE ###########################################################

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
using Turing
#using Plots; gr()

# NOTE: chiange-point model is taken from http://www.mit.edu/~ilkery/papers/GibbsSampling.pdf
@model change_point(x) = begin
  N = length(x)
  n ~ Uniform(1, N)
  lam1 ~ Gamma(2, 1)
  lam2 ~ Gamma(2, 1)
  for i = 1:floor(Int, n)
    x[i] ~ Poisson(lam1)
  end
  for i = floor(Int, n)+1:N
    x[i] ~ Poisson(lam2)
  end
end

_N = 50
_n = floor(Int, rand(Uniform(1, _N)))
_lam1 = rand(Gamma(2, 1))
_lam2 = rand(Gamma(2, 1))
x = Vector{Int}(undef, _N)
for i = 1:_n
  x[i] = rand(Poisson(_lam1))
end
for i = _n+1:_N
  x[i] = rand(Poisson(_lam2))
end
println("N = $_N, n = $_n, lam1 = $_lam1, lam2 = $_lam2")

#plot(x, linetype=[:scatter], lab="data point")
#plot!([_n], linetype=:vline, lab="change point")

modelf = change_point(x)

chn = sample(modelf, SMC(2000))

lam1_ = mean(chn[:lam1]); lam2_ = mean(chn[:lam2]); n_ = mean(chn[:n])

#plot!([i <= n_ ? lam1_ : lam2_ for i = 1:_N], linetype=:steppre, lab="inference result")
#title!("Change-point model (raw data & inference results)")
#xlabel!("n"); ylabel!("counts")


text:
function solve()
    h, w = [parse(Int, x) for x in split(readline())]
    s = []
    for i=1:h
        push!(s, readline())
    end
    ret = 0
    for i=1:h-1, j=1:w-1
        c = length([x for x in [s[i][j], s[i+1][j], s[i][j+1], s[i+1][j+1]] if x == '#'])
        if c == 1 || c == 3
            ret += 1
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module StringParserPEG

export StandardCache, Node, transform, Grammar, Rule, MatchRule
export parseGrammar, parseDefinition
export standardrules

import Base: show, parse, ==

include("rules.jl")
include("grammar.jl")
include("comparison.jl")
include("standardactions.jl")
include("Node.jl")
include("parse.jl")
include("transform.jl")
include("grammarparsing.jl")
include("standardrules.jl")

end


text:
using Documenter, UtilityFunctionsLH, FilesLH

makedocs(
    modules = [UtilityFunctionsLH],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    checkdocs = :exports,
    sitename = "UtilityFunctionsLH.jl",
    pages = Any["index.md"]
)

pkgDir = rstrip(normpath(@__DIR__, ".."), '/');
@assert endswith(pkgDir, "UtilityFunctionsLH")
deploy_docs(pkgDir; trialRun = false);

# deploydocs(
#     repo = "github.com/hendri54/UtilityFunctionsLH.jl.git",
# )

# ----------


text:
struct SGD_prox_DNN_iterable{R<:Real,Tg,Tf,Tp}
    F::Union{Array{Tf},Tf}            # smooth term
    g::Tg                   # nonsmooth term
    opt_params::Tp                  # initial point
    N::Int64                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    μ::Maybe{Union{Array{R},R}}  # convexity moduli of the gradients
    γ::Maybe{R}             # stepsize
    plus::Bool              # plus version (diminishing stepsize)
    η0::R
    η_tilde::R
    data::Tuple
    DNN_config!
end

mutable struct SGD_prox_DNN_state{R<:Real,Tx,Mx}
    γ::R                    # stepsize
    z::Tx
    cind::Int               # current interation index
    idxr::Int               # current index
    # some extra placeholders
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx
    temp_x::Mx
    temp_y::Tx
end

function SGD_prox_DNN_state(γ::R, z::Tx, cind, temp_x::Mx, temp_y) where {R,Tx,Mx}
    return SGD_prox_DNN_state{R,Tx,Mx}(γ, z, cind, Int(0), copy(z), copy(z), temp_x, temp_y)
end

function Base.iterate(iter::SGD_prox_DNN_iterable{R}) where {R}
    println("started")
    N = iter.N
    ind = collect(1:N)
    # updating the stepsize
    if iter.γ === nothing && !iter.plus
        if iter.L === nothing
            @warn "smoothness or convexity parameter absent"
            return nothing
        else
            L_M = maximum(iter.L)
            γ = 1/ (2*L_M)
        end
    else
        γ = iter.γ # provided γ
    end
    if iter.plus
        println("diminishing stepsize version")
        γ = 0.0
    end
    # initializing
    cind = 0
    z = iter.DNN_config!()
    temp_x = zeros(size(iter.opt_params[1],2))
    temp_y = zeros(size(iter.opt_params[length(iter.opt_params)],1))
    state = SGD_prox_DNN_state(γ, z, cind, temp_x, temp_y)
    return state, state
end

batchmemaybe(x) = tuple(x)
batchmemaybe(x::Tuple) = x

function Base.iterate(iter::SGD_prox_DNN_iterable{R}, state::SGD_prox_DNN_state{R}) where {R}
    # The inner cycle
    for i=1:iter.N # just for speed in implementation
        state.cind += 1
        if iter.plus
            state.γ = iter.η0/(1 + iter.η_tilde * floor(state.cind/iter.N))
        end
        # if mod(state.cind,100) == 0
        #     println("iter: $(i)")
        #     flush(stdout)
        # end

        # state.tupl .= (iter.data[1][:,i],iter.data[2][:,i])
        state.temp_x .= iter.data[1][:,i]
        state.temp_y .= iter.data[2][:,i]
        # gradient!(state.∇f_temp, iter.F[i], state.z)
        # state.∇f_temp .= (Flux.gradient(() -> iter.F[i](state.z), params(state.z)))[state.z][:,1]

        gs = Flux.gradient(iter.opt_params) do # sampled gradient
        iter.F(state.temp_x,state.temp_y)
        # iter.F(batchmemaybe(tupl)...)
        end
        state.∇f_temp .= iter.DNN_config!(gs=gs)
        state.z .= iter.DNN_config!()

        state.∇f_temp .*= - state.γ
        state.∇f_temp ./= iter.N
        state.∇f_temp .+= state.z

        CIAOAlgorithms.prox!(state.z, iter.g, state.∇f_temp, state.γ)
        iter.DNN_config!(gs=state.z)
    end

    return state, state
end

solution(state::SGD_prox_DNN_state) = state.z


text:
module DataREPL

using Markdown

_data_repl_help = md"""
## DataSets Data REPL

Press `>` to enter the data repl. Press TAB to complete commands.

|   Command    | Alias     | Action      |
|:----------   |:--------- | :---------- |
| `help`       | `?`       | Show this message |
| `list`       | `ls`      | List all datasets by name |
| `show $name` |           | Preview the content of dataset `$name` |
| `stack` | `st`           | Manipulate the global data search stack |
| `stack list` | `st ls`   | List all projects in the global data search stack |
| `stack push $path` | `st push` | Add data project `$path` to front of the search stack |
| `stack pop`  | `st pop`  | Remove data project from front of the search stack    |

"""

using ..DataSets

using ResourceContexts

using REPL
using ReplMaker
# using URIs

#-------------------------------------------------------------------------------
# Utilities for browsing dataset content
function show_dataset(name)
    out_stream = stdout
    @context begin
        data = @! open(dataset(name))
        _show_dataset(out_stream, data)
    end
end

# hex dump in xxd format
function hexdump(out_stream, buf; groups_per_line=8, group_size=2, max_lines=typemax(Int))
    linesize = groups_per_line*group_size
    for line = 1:min(max_lines, div(length(buf), linesize, RoundUp))
        linebuf = buf[(line-1)*linesize+1 : min(line*linesize,end)]
        address = (line-1)*linesize
        print(out_stream, string(address, base=16, pad=4), ": ")
        for group = 1:groups_per_line
            for i=1:group_size
                j = (group-1)*group_size+i
                if j <= length(linebuf)
                    print(out_stream, string(linebuf[j], base=16, pad=2))
                else
                    print(out_stream, "  ")
                end
            end
            print(out_stream, ' ')
        end
        print(out_stream, ' ')
        for j = 1:linesize
            c = Char(j <= length(linebuf) ? linebuf[j] : ' ')
            print(out_stream, isprint(c) ? c : '.')
        end
        print(out_stream, '\n')
    end
end

function _show_dataset(out_stream::IO, blob::Blob)
    @context begin
        io = @! open(IO, blob)
        N = 1024
        buf = read(io, N)
        str = String(copy(buf))
        n_textlike = count(str) do c
            isvalid(c) || return false
            isprint(c) || c in ('\n', '\r', '\t')
        end
        display_lines, _ = displaysize(out_stream)
        max_lines = max(5, display_lines ÷ 2)
        if n_textlike / length(str) > 0.95
            # It's approximately UTF-8 encoded text data - print as text
            lines = split(str, '\n', keepempty=true)
            nlines = min(lastindex(lines), max_lines)
            print(out_stream, join(lines[1:nlines], '\n'))
            println(out_stream)
            if !eof(io) || nlines < length(lines)
                println(out_stream, "⋮")
            end
        else
            # It's something else, perhaps binary or another text
            # encoding. Do a hex dump instead.
            println(out_stream, "Binary data:")
            hexdump(out_stream, buf; max_lines=max_lines)
            if !eof(io)
                println(out_stream, "⋮")
            end
        end
    end
end

function _show_dataset(out_stream::IO, tree::BlobTree)
    show(out_stream, MIME("text/plain"), tree)
end

function _show_dataset(out_stream::IO, x)
    show(out_stream, MIME("text/plain"), x)
end

#-------------------------------------------------------------------------------
# REPL command handling and completions

# function split_command(str)
#     token_ranges = []
#     i = 1
#     while true
#         rng = findnext(r"[^\s]+", full, i)
#         !isnothing(rng) || break
#         push!(token_ranges, rng)
#         i = last(rng)+1
#     end
#     tokens = getindex.(full, token_ranges)
#     token_ranges, tokens
# end

function complete_command_list(cmd_prefix, commands)
    # Completions for basic commands
    completions = String[]
    for cmdset in commands
        for cmd in cmdset
            if cmd == cmd_prefix
                # Space after full length command
                return ([" "], "", true)
            end
            if startswith(cmd, cmd_prefix)
                push!(completions, cmd*" ")
                break
            end
        end
    end
    return completions
end

function path_str(path_completion)
    path = REPL.REPLCompletions.completion_text(path_completion)
    if Sys.iswindows()
        # On windows, REPLCompletions.complete_path() adds extra escapes for
        # use within a normal string in the Juila REPL but we don't need those.
        path = replace(path, "\\\\"=>'\\')
    end
    return path
end

function complete(str_to_complete)
    tokens = split(str_to_complete, r" +", keepempty=true)
    cmd = popfirst!(tokens)
    if isempty(tokens)
        # Completions for basic commands
        completions = complete_command_list(cmd, [
            ("list","ls"),
            ("show",),
            ("stack",),
            ("help","?")
        ])
        # Empty completion => return anyway to show user their prefix is wrong
        return (completions, cmd, !isempty(completions))
    end
    if cmd == "show"
        if length(tokens) <= 1
            name_prefix = isempty(tokens) ? "" : tokens[1]
            completions = String[]
            ks = sort!(collect(keys(DataSets.PROJECT)))
            for k in ks
                if startswith(k, name_prefix) && k != name_prefix
                    push!(completions, k)
                end
            end
            return (completions, name_prefix, !isempty(completions))
        end
    elseif cmd == "stack"
        if length(tokens) <= 1
            subcmd_prefix = isempty(tokens) ? "" : tokens[1]
            # Completions for project stack subcommands
            completions = complete_command_list(subcmd_prefix, [
                ("push",),
                ("pop",),
                ("list","ls",)
            ])
            return (completions, tokens[1], !isempty(completions))
        elseif length(tokens) == 2
            subcmd = popfirst!(tokens)
            if subcmd == "push"
                path_prefix = isempty(tokens) ? "" : tokens[1]
                (path_completions, range, should_complete) =
                    REPL.REPLCompletions.complete_path(path_prefix, length(path_prefix))
                completions = [path_str(c) for c in path_completions]
                return (completions, path_prefix[range], should_complete)
            end
        end
    end
    return ([], "", false)
end

# Translate `data>` REPL syntax into an Expr to be evaluated in the REPL
# backend.
function parse_data_repl_cmd(cmdstr)
    # Use shell tokenization rules for familiarity
    tokens = Base.shell_split(cmdstr)
    cmd = tokens[1]
    popfirst!(tokens)
    if cmd in ("list", "ls")
        return quote
            $DataSets.DataProject($DataSets.PROJECT)
        end
    elseif cmd == "stack" && length(tokens) >= 1
        subcmd = popfirst!(tokens)
        if subcmd == "push"
            path = popfirst!(tokens)
            return quote
                proj = $DataSets.data_project_from_path($path)
                stack = $DataSets.PROJECT
                pushfirst!(stack, proj)
                stack
            end
        elseif subcmd == "pop"
            return quote
                stack = $DataSets.PROJECT
                popfirst!(stack)
                stack
            end
        elseif subcmd in ("list", "ls")
            return quote
                $DataSets.PROJECT
            end
        end
    elseif cmd == "show"
        name = tokens[1]
        return quote
            $DataREPL.show_dataset($name)
        end
    elseif cmd in ("help", "?")
        return _data_repl_help
    end
    error("Invalid data REPL syntax: \"$cmdstr\"")
end

#-------------------------------------------------------------------------------
# Integration with REPL / ReplMaker
struct DataCompletionProvider <: REPL.LineEdit.CompletionProvider
end

function REPL.complete_line(provider::DataCompletionProvider,
                            state::REPL.LineEdit.PromptState)::
                            Tuple{Vector{String},String,Bool}
    # See REPL.jl complete_line(c::REPLCompletionProvider, s::PromptState)
    partial = REPL.beforecursor(state.input_buffer)
    full = REPL.LineEdit.input_string(state)
    if partial != full
        # For now, only complete at end of line
        return ([], "", false)
    end
    complete(full)
end

function init_data_repl(; start_key = ">")
    ReplMaker.initrepl(parse_data_repl_cmd,
                       repl = Base.active_repl,
                       # valid_input_checker = Support for multiple lines syntax?
                       prompt_text = "data> ",
                       prompt_color = :red,
                       start_key = start_key,
                       sticky_mode = true,
                       mode_name = "DataSets",
                       completion_provider = DataCompletionProvider(),
                       startup_text = false)
    nothing
end

function __init__()
    isinteractive() && init_data_repl()
end

end


text:
using MultivariateStats
using LinearAlgebra
using Test
using StatsBase
import Statistics: mean, cov
import Random

@testset "Multi-class LDA" begin

    #Test equivalence of eigenvectors/singular vectors taking into account possible
    #phase (sign) differences
    #This may end up in Base.Test; see JuliaLang/julia#10651
    function test_approx_eq_vecs(a::StridedVecOrMat{S},
        b::StridedVecOrMat{T}, error=nothing) where {S<:Real, T<:Real}
        n = size(a, 2)
        @test n==size(b, 2) && size(a, 1)==size(b, 1)
        error==nothing && (error=n^3*(eps(S)+eps(T)))
        for i=1:n
            ev1, ev2 = a[:, i], b[:, i]
            deviation = min(abs(norm(ev1-ev2)), abs(norm(ev1+ev2)))
            if !isnan(deviation)
                @test isapprox.(deviation, 0.0, atol=error)
            end
        end
    end

    Random.seed!(34568)

    ## prepare data
    d = 5
    ns = [10, 15, 20]
    nc = length(ns)
    n = sum(ns)
    Xs = Matrix{Float64}[]
    ys = Vector{Int}[]
    Ss = Matrix{Float64}[]
    cmeans = zeros(d, nc)

    for k = 1:nc
        R = qr(randn(d, d)).Q
        nk = ns[k]

        Xk = R * Diagonal(2 * rand(d) .+ 0.5) * randn(d, nk) .+ randn(d)
        yk = fill(k, nk)
        uk = vec(mean(Xk, dims=2))
        Zk = Xk .- uk
        Sk = Zk * Zk'

        push!(Xs, Xk)
        push!(ys, yk)
        push!(Ss, Sk)
        cmeans[:,k] .= uk
    end

    X = hcat(Xs...)
    y = vcat(ys...)
    mv = vec(mean(X, dims=2))
    Sw = zeros(d, d)
    for k = 1:nc
        Sw += Ss[k]
    end

    Sb = zeros(d, d)
    for k = 1:nc
        dv = cmeans[:,k] - mv
        Sb += ns[k] * (dv * dv')
    end

    Z = X .- mv
    Sa = Z * Z'
    @test Sw + Sb ≈ Sa

    @assert size(X) == (d, n)
    @assert size(y) == (n,)

    ## Stats

    S = multiclass_lda_stats(nc, X, y)

    @test S.dim == d
    @test S.nclasses == nc

    @test classweights(S) == ns
    @test classmeans(S) ≈ cmeans
    @test mean(S) ≈ mv

    @test withclass_scatter(S) ≈ Sw
    @test betweenclass_scatter(S) ≈ Sb

    covestimator = SimpleCovariance(corrected=true)
    Sce = multiclass_lda_stats(nc, X, y; covestimator_between=covestimator, covestimator_within=covestimator)

    Swcorr = Sw * sum(ns)/(sum(ns)-1)
    Sbcorr = Sb * nc/(nc-1)

    @test withclass_scatter(Sce) ≈ Swcorr
    @test betweenclass_scatter(Sce) ≈ Sbcorr

    ## Solve

    emax = maximum(eigvals(Sw))
    Sw_c = copy(Sw)
    Sb_c = copy(Sb)

    lambda = 1.0e-3
    Sw_r = Sw + (lambda * emax) * Matrix(I, d, d)

    P1 = mclda_solve(Sb, Sw, :gevd, nc-1, lambda)
    @test P1' * Sw_r * P1 ≈ Matrix(I, nc-1, nc-1)
    U = Sb * P1
    V = Sw_r * P1
    # test whether U is proportional to V,
    # which indicates P is the generalized eigenvectors
    @test U ≈ V*Diagonal(vec(mean(U./V, dims=1)))

    P2 = mclda_solve(Sb, Sw, :whiten, nc-1, lambda)
    @test P2' * Sw_r * P2 ≈ Matrix(I, nc-1, nc-1)

    test_approx_eq_vecs(P1, P2)

    ## MC-LDA

    for T in (Float32, Float64)
        M = fit(MulticlassLDA, nc, convert(Matrix{T}, X), y; method=:gevd, regcoef=convert(T, lambda))
        @test size(M) == (d, nc - 1)
        @test projection(M) ≈ P1
        @test M.pmeans ≈ M.proj'cmeans
        @test predict(M, X) ≈ M.proj'X

        M = fit(MulticlassLDA, nc, convert(Matrix{T}, X), y; method=:whiten, regcoef=convert(T, lambda))
        @test size(M) == (d, nc - 1)
        # @test projection(M) P2  # signs may change
        @test M.pmeans ≈ M.proj'cmeans
        @test predict(M, X) ≈ M.proj'X
    end

    ## High-dimensional LDA (subspace LDA)

    # Low-dimensional case (no singularities)

    centers = [zeros(5) [10.0;zeros(4)] [0.0;10.0;zeros(3)]]

    # Case 1: 3 groups of 500
    dX = randn(5,1500);
    for i = 0:500:1000
        dX[:,(1:500).+i] .= dX[:,(1:500).+i] .- mean(dX[:,(1:500).+i], dims=2)  # make the mean of each 0
    end
    dX1 = dX
    X1 = [dX[:,1:500].+centers[:,1] dX[:,501:1000].+centers[:,2] dX[:,1001:1500].+centers[:,3]]
    label1 = [fill(1,500); fill(2,500); fill(3,500)]
    # Case 2: 3 groups, one with 1000, one with 100, and one with 10
    dX = randn(5,1110);
    dX[:,   1:1000] .= dX[:,   1:1000] .- mean(dX[:,   1:1000], dims=2)
    dX[:,1001:1100] .= dX[:,1001:1100] .- mean(dX[:,1001:1100], dims=2)
    dX[:,1101:1110] .= dX[:,1101:1110] .- mean(dX[:,1101:1110], dims=2)
    dX2 = dX
    X2 = [dX[:,1:1000].+centers[:,1] dX[:,1001:1100].+centers[:,2] dX[:,1101:1110].+centers[:,3]]
    label2 = [fill(1,1000); fill(2,100); fill(3,10)]

    for (X, dX, label) in ((X1, dX1, label1), (X2, dX2, label2))
        w = [length(findall(label.==1)), length(findall(label.==2)), length(findall(label.==3))]
        M = fit(SubspaceLDA, X, label)
        @test indim(M) == 5
        @test outdim(M) == 2
        totcenter = vec(sum(centers.*w', dims=2)./sum(w))
        @test mean(M) ≈ totcenter
        @test classmeans(M) ≈ centers
        @test classweights(M) == w
        x = rand(5)
        @test predict(M, x) ≈ projection(M)'*x
        dcenters = centers .- totcenter
        Hb = dcenters.*sqrt.(w)'
        Sb = Hb*Hb'
        Sw = dX*dX'
        # When there are no singularities, we need to solve the "regular" LDA eigenvalue equation
        proj = projection(M)
        @test Sb*proj ≈ Sw*proj*Diagonal(M.λ)
        @test proj'*Sw*proj ≈ Matrix(I,2,2)

        # also check that this is consistent with the conventional algorithm
        Mld = fit(MulticlassLDA, 3, X, label)
        for i = 1:2
            @test abs(dot(normalize(proj[:,i]), normalize(Mld.proj[:,i]))) ≈ 1
        end
    end

    # High-dimensional case (undersampled => singularities)
    X = randn(10^6, 9)
    label = rand(1:3, 9); label[1:3] = 1:3
    M = fit(SubspaceLDA, X, label)
    centers = M.cmeans
    for i = 1:3
        flag = label.==i
        @test centers[:,i] ≈ mean(X[:,flag], dims=2)
        @test M.cweights[i] == sum(flag)
    end
    dcenters = centers .- mean(X, dims=2)
    Hb = dcenters*Diagonal(sqrt.(M.cweights))
    Hw = X - centers[:,label]
    @test (M.projw'*Hb)*(Hb'*M.projw)*M.projLDA ≈ (M.projw'*Hw)*(Hw'*M.projw)*M.projLDA*Diagonal(M.λ)
    Gw = projection(M)'*Hw
    @test Gw*Gw' ≈ Matrix(I,2,2)

    # Test that nothing breaks if one class has no members
    label[1:4] .= 1
    label[5:9] .= 3  # no data points have label == 2
    M = fit(SubspaceLDA, X, label)
    centers = M.cmeans
    dcenters = centers .- mean(X, dims=2)
    Hb = dcenters*Diagonal(sqrt.(M.cweights))
    Hw = X - centers[:,label]
    @test (M.projw'*Hb)*(Hb'*M.projw)*M.projLDA ≈ (M.projw'*Hw)*(Hw'*M.projw)*M.projLDA*Diagonal(M.λ)

    # Test normalized LDA
    function gen_ldadata_2(centers, n1, n2)
        d = size(centers, 1)
        X = randn(d, n1+n2)
        X[:,1:n1]       .-= vec(mean(X[:,1:n1], dims=2))
        X[:,n1+1:n1+n2] .-= vec(mean(X[:,n1+1:n1+n2], dims=2))
        dX = copy(X)
        X[:,1:n1]       .+= centers[:,1]
        X[:,n1+1:n1+n2] .+= centers[:,2]
        label = [fill(1,n1); fill(2,n2)]
        X, dX, label
    end
    centers = zeros(2,2); centers[1,2] = 10
    n1 = 100
    for n2 in (100, 1000, 10000)
        X, dX, label = gen_ldadata_2(centers, n1, n2)
        M = fit(SubspaceLDA, X, label; normalize=true)
        proj = projection(M)
        Hb = centers .- mean(centers, dims=2)   # not weighted by number in each cluster
        Sb = Hb*Hb'
        Sw = dX*dX'/(n1+n2)
        @test Sb*proj ≈ Sw*proj*Diagonal(M.λ)
    end

    # Test various input data types
    for (T, nrm) in Iterators.product((Float64, Float32), (false, true))
        n2 = 100
        X, dX, label = gen_ldadata_2(centers, n1, n2)
        M = fit(SubspaceLDA, convert(Matrix{T}, X), label; normalize=nrm)
        proj = projection(M)
        @test eltype(proj) === T
    end

end


text:
abstract type NamingConvention end

abstract type CamelCase <: NamingConvention end
abstract type SnakeCase <: NamingConvention end

struct SnakeCaseLower <: SnakeCase
    value::Any
    SnakeCaseLower(value) =
        is_snake_case(value) && lowercase(value) == value ? new(value) : error("Invalid string format $value")
end

struct SnakeCaseUpper <: SnakeCase
    value::Any
    SnakeCaseUpper(value) =
        is_snake_case(value) && uppercase(value) == value ? new(value) : error("Invalid string format $value")
end

struct CamelCaseLower <: CamelCase
    value::Any
    CamelCaseLower(value) =
        is_camel_case(value) && uppercase(value[1]) != value[1] ? new(value) : error("Invalid string format $value")
end

struct CamelCaseUpper <: CamelCase
    value::Any
    CamelCaseUpper(value) =
        is_camel_case(value) && uppercase(value[1]) == value[1] ? new(value) : error("Invalid string format $value")
end

Base.split(str::SnakeCase) = split(str.value, "_")

function Base.split(str::CamelCase)
    strval = str.value
    lowercase(strval) == strval && return [strval]
    length(strval) > 1 && lowercase(strval[2:end]) == strval[2:end] && return [strval]
    reg_upper = r"(([A-Z]+|\d+)(?=(([A-Z]+|\d+)|$))|([A-Z]{1})[a-z]*?(?=($|([A-Z]|\d))))"
    if uppercase(strval[1]) == strval[1] # CamelCaseUpper
        matches = getproperty.(collect(eachmatch(reg_upper, strval)), :match)
    else
        first = match(r"[a-z]+(?=([A-Z]|\d))", strval).match
        matches = [first, getproperty.(collect(eachmatch(reg_upper, strval)), :match)...]
    end
    matches
end

SnakeCaseLower(parts::AbstractArray) = SnakeCaseLower(lowercase(snake_case(parts)))
SnakeCaseUpper(parts::AbstractArray) = SnakeCaseUpper(uppercase(snake_case(parts)))

CamelCaseLower(parts::AbstractArray) = CamelCaseLower(camel_case(lowercase.(parts)))
CamelCaseUpper(parts::AbstractArray) = CamelCaseUpper(uppercasefirst(camel_case(lowercase.(parts))))

snake_case(parts::AbstractArray) = join(parts, "_")
camel_case(parts::AbstractArray) = length(parts) == 1 ? parts[1] : join([parts[1], uppercasefirst.(parts[2:end])...])

Base.convert(T::Type{SnakeCaseLower}, str::SnakeCaseUpper) = T(lowercase(str.value))
Base.convert(T::Type{SnakeCaseUpper}, str::SnakeCaseLower) = T(uppercase(str.value))
Base.convert(T::Type{CamelCaseLower}, str::CamelCaseUpper) = T(lowercase(str.value[1]) * str.value[2:end])
Base.convert(T::Type{CamelCaseUpper}, str::CamelCaseLower) = T(uppercasefirst(str.value))
Base.convert(T::Type{<:CamelCase}, str::SnakeCase) = T(split(str))
Base.convert(T::Type{<:SnakeCase}, str::CamelCase) = T(split(str))
nc_convert(T::Type{<:NamingConvention}, str::AbstractString) =
    Base.convert(T, (detect_convention(str, instance = true))).value
nc_convert(T::Type{<:NamingConvention}, sym::Symbol) = Symbol(nc_convert(T, string(sym)))

is_camel_case(str) = !occursin("_", str)
is_snake_case(str) = lowercase(str) == str || uppercase(str) == str

is_camel_case(str::NamingConvention) = is_camel_case(str.value)
is_snake_case(str::NamingConvention) = is_snake_case(str.value)

function remove_parts(str::T, discarded_parts) where {T<:NamingConvention}
    splitted_str = split(str)
    parts_to_keep = 1:length(splitted_str) |> x -> filter(y -> y ∉ discarded_parts, x) |> collect
    kept_parts = getindex.(Ref(splitted_str), parts_to_keep)
    T(kept_parts)
end

remove_parts(str; discarded_parts = [1]) = remove_parts(detect_convention(str, instance = true), discarded_parts)

remove_prefix(name::T) where {T<:NamingConvention} = T(split(name)[2:end])
remove_prefix(str) = remove_prefix(detect_convention(str, instance = true)).value

function detect_convention(str; instance = false)
    instanced(T, x) = instance ? T(x) : T
    is_camel_case(str) && lowercase(str)[1] == str[1] && return instanced(CamelCaseLower, str)
    is_camel_case(str) && uppercase(str)[1] == str[1] && return instanced(CamelCaseUpper, str)
    is_snake_case(str) && lowercase(str) == str && return instanced(SnakeCaseLower, str)
    is_snake_case(str) && uppercase(str) == str && return instanced(SnakeCaseUpper, str)
    error("No convention detected for string $str")
end


text:
function colorize(ano, default::Colorant, id_color::Pair...)
    dest = similar(ano, Int32)
    colorize!(dest, ano, default, id_color...)
end

function colorize(ano::AxisArray, default::Colorant, id_color::Pair...)
    dest = similar(parent(ano), Int32)
    cdest = colorize!(dest, ano, default, id_color...)
    AxisArray(cdest, axes(ano))
end

colorize!(ano, default::Colorant, id_color::Pair...) =
    colorize!(ano, ano, default, id_color...)

function colorize!(ano::AxisArray, default::Colorant, id_color::Pair...)
    dest = colorize!(parent(ano), parent(ano), default, id_color...)
    AxisArray(dest, axes(ano))
end

function colorize!(dest::AbstractArray{T}, ano, default::Colorant, id_color::Pair...) where T
    colorlist = [RGBA{Float32}(0,0,0,0)]   # background (not inside brain)
    relabeldict = Dict{eltype(ano),eltype(dest)}(0=>1)
    # Add the default
    push!(colorlist, default)
    for (ks, v) in id_color
        push!(colorlist, v)
        for k in ks
            relabeldict[k] = length(colorlist)
        end
    end
    colorize!(dest, ano, relabeldict, colorlist)
end

function colorize!(dest::AbstractArray{To}, ano::AbstractArray{Ti},
                   relabeldict::Dict{Ti,To},
                   colorlist::AbstractVector{C}) where {To<:Integer, Ti<:Integer,C<:Colorant}
    indices(dest) == indices(ano) || throw(DimensionMismatch("inputs must match, got $(indices(dest)) and $(indices(ano))"))
    @inbounds for i in CartesianRange(indices(dest))
        dest[i] = get(relabeldict, ano[i], To(2))
    end
    IndirectArray(dest, colorlist)
end


text:
data() = readlines("day05.txt")

# This is really just the binary representation of an integer in disguise!
row(s) = parse(Int, replace(replace(s[1:7], "F"=>"0"), "B"=>"1"); base = 2)
col(s) = parse(Int, replace(replace(s[8:10], "L"=>"0"), "R"=>"1"); base = 2)
id(s) = row(s) * 8 + col(s)

# quick test
# row("BFFFBBFRRR")
# col("BFFFBBFRRR")
# id("BFFFBBFRRR")

# id("FFFBBBFRRR") == 119
# id("BBFFBBFRLL") == 820

# read data
boarding_passes = id.(data())

# part 1
part1() = maximum(boarding_passes)  # 980

# part 2
function part2()
    sorted = sort(boarding_passes)
    idx = argmax([sorted[i] - sorted[i-1] > 1 for i in 2:length(sorted)])  # 567
    return sorted[idx] + 1
end

# using BenchmarkTools

# From Aaron on Discord. The mapfoldl function is awesome in terms of performance!
seatnumber(line) = mapfoldl(c->c in "BR", (x, y)->2x+y, line, init=0)
seatnumber("BFFFBBFRRR")


text:
using VortexDistributions

# simple test field
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi = Torus(copy(psi0),x,y)
xp = 1.133
yp = 2.787
vp = PointVortex(xp,yp,1)

## Test differnt calls to ScalarVortex()
sp = ScalarVortex(vp)
sp1 = ScalarVortex(1.0, vp)
@test sp1[1].core === sp.core
@test sp1[1].vort === sp.vort

sp2 = ScalarVortex(0.8, vp)
@test !(sp2[1].core===sp.core)
@test sp2[1].vort === sp.vort

sp3 = ScalarVortex([1.0, 1.0], [vp])
@test sp3[1].core === sp3[2].core
@test sp3[1].vort === sp3[2].vort

sp4 = ScalarVortex([1.0, 0.8], [vp])
@test !(sp4[1].core === sp4[2].core)
@test (sp4[1].vort === sp4[2].vort)

sp5 = ScalarVortex([vp])
@test sp5[1].core === sp1[1].core
@test sp5[1].vort === sp1[1].vort

rs1 = rand_scalarvortex()
@test typeof(rs1) <: ScalarVortex{Exact}

rs2 = rand_scalarvortex(10)
@test length(rs2) == 10
@test typeof(rs2[1]) <: ScalarVortex{Exact}

rs3 = rand_scalarvortex(10, psi)
@test typeof(rs3) == typeof(rs2)
@test length(rs3) == 10

rs4 = rand_scalarvortex(psi)
@test typeof(rs4) <: ScalarVortex{Exact}

rv1 = rand_vortex()
@test typeof(rv1) <: ScalarVortex{Exact}

rv2 = rand_vortex(10)
@test length(rv2) == 10
@test typeof(rv2[1]) <: ScalarVortex{Exact}

rv3 = rand_scalarvortex(10, psi)
@test typeof(rv3) == typeof(rv2)
@test length(rv3) == 10


text:
@testset "msgpass" begin
    T = Float32
    in_channel = 10
    out_channel = 5
    num_V = 6
    num_E = 7

    adj = T[0. 1. 0. 0. 0. 0.;
        1. 0. 0. 1. 1. 1.;
        0. 0. 0. 0. 0. 1.;
        0. 1. 0. 0. 1. 0.;
        0. 1. 0. 1. 0. 1.;
        0. 1. 1. 0. 1. 0.]

    struct NewLayer <: MessagePassing
        weight
    end
    NewLayer(m, n) = NewLayer(randn(T, m,n))

    function (l::NewLayer)(fg::FeaturedGraph, X::AbstractMatrix)
        _, x, _ = GeometricFlux.propagate(l, fg, edge_feature(fg), X, global_feature(fg), +)
        x
    end

    (l::NewLayer)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

    X = Array{T}(reshape(1:num_V*in_channel, in_channel, num_V))
    fg = FeaturedGraph(adj, nf=X, ef=Fill(zero(T), 0, num_E))

    l = NewLayer(out_channel, in_channel)

    @testset "no message or update" begin
        fg_ = l(fg)

        @test GraphSignals.adjacency_matrix(fg_) == adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (0, num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.message(l::NewLayer, x_i, x_j, e_ij) = l.weight * x_j
    @testset "message function" begin
        fg_ = l(fg)

        @test GraphSignals.adjacency_matrix(fg_) == adj
        @test size(node_feature(fg_)) == (out_channel, num_V)
        @test size(edge_feature(fg_)) == (0, num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.update(l::NewLayer, m, x) = l.weight * x + m
    @testset "message and update" begin
        fg_ = l(fg)

        @test GraphSignals.adjacency_matrix(fg_) == adj
        @test size(node_feature(fg_)) == (out_channel, num_V)
        @test size(edge_feature(fg_)) == (0, num_E)
        @test size(global_feature(fg_)) == (0,)
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.7.0-beta3
#     language: julia
#     name: julia-1.7
# ---

# # Remote Access To Climate Model Output
#
# - Climate model archive ([CMIP6](https://esgf-node.llnl.gov/search/cmip6/)) is accessed using [AWS.jl](https://github.com/JuliaCloud/AWS.jl) and [Zarr.jl](https://github.com/meggart/Zarr.jl) via [ClimateModels.jl](https://github.com/gaelforget/ClimateModels.jl)
# - Choose `institution_id`, `source_id`, `variable_id` (inside `CMIP6.jl`)
# - Compute and plot (1) time mean global map and (2) time evolving global mean (inside `CMIP6.jl`)

using ClimateModels
p=dirname(pathof(ClimateModels))
include(joinpath(p,"../examples/CMIP6.jl"))

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
module SudokuSolver

export solvesudoku

include("heap.jl")
include("sudoku.jl")
include("puzzleexamples.jl")
end


text:
import jInv.LinearSolvers.copySolver;
import jInv.Utils.clear!
import jInv.LinearSolvers.AbstractSolver
import jInv.LinearSolvers.AbstractDirectSolver
import jInv.LinearSolvers.factorLinearSystem!
import jInv.LinearSolvers.solveLinearSystem

export MUMPSsolver, getMUMPSsolver,copySolver

"""
mutable struct MUMPSsolver

Fields:

	Ainv    - holds MUMPSfactorization
	doClear - flag to clear factorization
	ooc     - flag for out-of-core option
	sym     - 0=unsymmetric, 1=symm. pos def, 2=general symmetric
	nFac    - number of factorizations performed
	facTime - cumulative time for factorizations
	nSolve  - number of solves
	solveTime - cumnulative time for solves

Example:

	Ainv = getMUMPSsolver()
"""
mutable struct MUMPSsolver<: AbstractDirectSolver
	Ainv
	doClear::Int
	ooc::Int
	sym::Int
	nFac::Int
	facTime::Real
	nSolve::Int
	solveTime::Real
end

# See if MUMPS solver is installed and include code for it.

if hasMUMPS

	function getMUMPSsolver(Ainv=[],doClear=1,ooc=0,sym=0)
		return MUMPSsolver(Ainv,doClear,ooc,sym,0,0.0,0,0.)
	end

    function factorLinearSystem!(A,param::MUMPSsolver)
		if param.doClear == 1
			clear!(param)
		end

		param.facTime+= @elapsed begin
			param.Ainv = factorMUMPS(A, param.sym, param.ooc)
		end
		param.nFac+=1
	end

	function solveLinearSystem!(A,B,X,param::MUMPSsolver,doTranspose::Int=0)
		if param.doClear == 1
			clear!(param)
		end

		if param.Ainv==[]
			param.facTime+= @elapsed begin
			param.Ainv = factorMUMPS(A, param.sym, param.ooc)
			end
			param.nFac+=1
		end

		param.solveTime+= @elapsed begin
			applyMUMPS!(param.Ainv, B,X,doTranspose)
		end
		param.nSolve+=1

		return X, param
	end # function solveLinearSystem MUMPSsolver

	function clear!(M::MUMPSsolver)
		if isa(M.Ainv,MUMPSfactorization)
			destroyMUMPS(M.Ainv)
			M.Ainv = []
		end
	end

	function copySolver(Ainv::MUMPSsolver)
		return MUMPSsolver([],Ainv.doClear,Ainv.ooc,Ainv.sym,0,0.0,0,0.);
	end

	hasMUMPS=true
end


text:
using AxisKeys
using InteractiveUtils
using NamedDims
using Statistics
using Printf
using Tullio

export loo, psis_loo, loo_from_psis, PsisLoo

#####################
###### STRUCTS ######
#####################

# """
#     PsisLooMethod

# Use Pareto-smoothed importance sampling together with leave-one-out cross validation to
# estimate the out-of-sample predictive accuracy.
# """
# struct PsisLooMethod <: AbstractCVMethod end

"""
    PsisLoo <: AbstractCV

A struct containing the results of leave-one-out cross validation computed with Pareto
smoothed importance sampling.

$CV_DESC

See also: [`loo`]@ref, [`bayes_cv`]@ref, [`psis_loo`]@ref, [`Psis`]@ref
"""
struct PsisLoo{
    RealType <: Real,
    ArrayType <: AbstractArray{RealType},
    VectorType <: AbstractVector{RealType},
} <: AbstractCV
    estimates::KeyedArray{RealType, 2, <:NamedDimsArray, <:Any}
    pointwise::KeyedArray{RealType, 2, <:NamedDimsArray, <:Any}
    psis_object::Psis{RealType, ArrayType, VectorType}
    gmpd::RealType
    mcse::RealType
end

function Base.show(io::IO, ::MIME"text/plain", loo_object::PsisLoo)
    table = loo_object.estimates
    _throw_pareto_k_warning(loo_object.pointwise(:pareto_k))
    post_samples = loo_object.psis_object.posterior_sample_size
    data_size = loo_object.psis_object.data_size
    println(
        "Results of PSIS-LOO-CV with $post_samples Monte Carlo samples and $data_size " *
        Printf.@sprintf("data points. Total Monte Carlo SE of %.2g.", loo_object.mcse),
    )
    return pretty_table(
        table;
        compact_printing=false,
        header=table.column,
        row_names=table.statistic,
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end

#####################
#### LOO METHODS ####
#####################

"""
    function loo(args...; kwargs...) -> PsisLoo

Compute an approximate leave-one-out cross-validation score.

Currently, this function only serves to call `psis_loo`, but this could change in the
future. The default methods or return type may change without warning, so we recommend
using `psis_loo` instead if reproducibility is required.

See also: [`psis_loo`](@ref), [`PsisLoo`](@ref).
"""
function loo(args...; kwargs...)
    return psis_loo(args...; kwargs...)
end

"""
    function psis_loo(
        log_likelihood::AbstractArray{<:Real} [, args...];
        [, chain_index::Vector{Int}, kwargs...]
    ) -> PsisLoo

Use Pareto-Smoothed Importance Sampling to calculate the leave-one-out cross validation
score.

# Arguments

  - $LOG_LIK_ARR
  - $ARGS [`psis`](@ref).
  - $CHAIN_INDEX_DOC
  - $KWARGS [`psis`](@ref).

See also: [`psis`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).
"""
function psis_loo(log_likelihood::AbstractArray{<:Real, 3}, args...; kwargs...)
    psis_object = psis(-log_likelihood, args...; kwargs...)
    return loo_from_psis(log_likelihood, psis_object)
end

function psis_loo(
    log_likelihood::AbstractMatrix{<:Real},
    args...;
    chain_index::AbstractVector=_assume_one_chain(log_likelihood),
    kwargs...,
)
    chain_index = Int.(chain_index)
    new_log_ratios = _convert_to_array(log_likelihood, chain_index)
    return psis_loo(new_log_ratios, args...; kwargs...)
end

"""
    loo_from_psis(
        log_likelihood::AbstractArray{<:Real}, psis_object::Psis;
        chain_index::Vector{<:Integer}
    )

Use a precalculated `Psis` object to estimate the leave-one-out cross validation score.

# Arguments

  - $LOG_LIK_ARR
  - `psis_object`: A precomputed `Psis` object used to estimate the LOO-CV score.
  - $CHAIN_INDEX_DOC

See also: [`psis`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).

"""
function loo_from_psis(log_likelihood::AbstractArray{<:Real, 3}, psis_object::Psis)
    dims = size(log_likelihood)
    data_size = dims[1]
    mcmc_count = dims[2] * dims[3]  # total number of samples from posterior
    log_count = log(mcmc_count)

    # TODO: Add a way of using score functions other than ELPD
    # log_likelihood::ArrayType = similar(log_likelihood)
    # log_likelihood .= score(log_likelihood)

    weights = psis_object.weights
    ξ = psis_object.pareto_k
    r_eff = psis_object.r_eff

    @tullio pointwise_loo[i] := weights[i, j, k] * exp(log_likelihood[i, j, k]) |> log
    @tullio pointwise_naive[i] := exp(log_likelihood[i, j, k] - log_count) |> log
    pointwise_p_eff = pointwise_naive - pointwise_loo
    pointwise_mcse = _calc_mcse(weights, log_likelihood, pointwise_loo, r_eff)

    pointwise = KeyedArray(
        hcat(pointwise_loo, pointwise_naive, pointwise_p_eff, pointwise_mcse, ξ);
        data=1:length(pointwise_loo),
        statistic=[:cv_elpd, :naive_lpd, :p_eff, :mcse, :pareto_k],
    )

    table = _generate_loo_table(pointwise)

    gmpd = exp.(table(column=:mean, statistic=:cv_elpd))

    mcse = sum(abs2, pointwise_mcse) |> sqrt
    return PsisLoo(table, pointwise, psis_object, gmpd, mcse)
end

function loo_from_psis(
    log_likelihood::AbstractMatrix{<:Real}, psis_object::Psis, args...;
    chain_index::AbstractVector=_assume_one_chain(log_likelihood), kwargs...
)
    chain_index = Int.(chain_index)
    new_log_ratios = _convert_to_array(log_likelihood, chain_index)
    return loo_from_psis(new_log_ratios, psis_object, args...; kwargs...)
end

function _generate_loo_table(pointwise::AbstractMatrix{<:Real})

    data_size = size(pointwise, :data)
    # create table with the right labels
    table = KeyedArray(
        similar(NamedDims.unname(pointwise), 3, 4);
        statistic=[:cv_elpd, :naive_lpd, :p_eff],
        column=[:total, :se_total, :mean, :se_mean],
    )

    # calculate the sample expectation for the total score
    to_sum = pointwise([:cv_elpd, :naive_lpd, :p_eff])
    @tullio avgs[statistic] := to_sum[data, statistic] |> _ / data_size
    avgs = reshape(avgs, 3)
    table(:, :mean) .= avgs

    # calculate the sample expectation for the average score
    table(:, :total) .= table(:, :mean) .* data_size

    # calculate the sample expectation for the standard error in the totals
    se_mean = std(to_sum; mean=avgs', dims=1) / sqrt(data_size)
    se_mean = reshape(se_mean, 3)
    table(:, :se_mean) .= se_mean

    # calculate the sample expectation for the standard error in averages
    table(:, :se_total) .= se_mean * data_size

    if table(:p_eff, :total) ≤ 0
        @warn "The calculated effective number of parameters is negative, which should " *
        "not be possible. PSIS has failed to approximate the target distribution."
    end

    return table
end

function _calc_mcse(weights, log_likelihood, pointwise_loo, r_eff)
    pointwise_gmpd = exp.(pointwise_loo)
    @tullio pointwise_var[i] :=
        (weights[i, j, k] * (exp(log_likelihood[i, j, k]) - pointwise_gmpd[i]))^2
    # If MCMC draws follow a log-normal distribution, then their log has this std. error:
    @. pointwise_var = log1p(pointwise_var / pointwise_gmpd^2)
    # (google "log-normal method of moments" for a proof)
    # apply MCMC correlation correction:
    return @. sqrt(pointwise_var / r_eff)
end

max_stars_repo_name: ParadaCarleton/JuLOOa

max_issues_repo_name: ParadaCarleton/JuLOOa

max_forks_repo_name: ParadaCarleton/JuLOOa


text:
##: ------------------------------------------------
#: general tree
##: ------------------------------------------------

abstract type AbstractTree end

function children(t::T)::Vector{T} where {T<:AbstractTree}
    return t.C
end

label(t::AbstractTree) = t.R

#: my favourite
function DFS(t::T, f::Function, g::Function, h::Function)::Any where {T <: AbstractTree}
    f(t)
    V = Any[DFS(c, f, g, h) for c ∈ children(t)]
    g(t)
    return h((V,t))
end

@inline DFS(t::T, f::Function, g::Function) where {T <: AbstractTree} = DFS(t, f, g, identity)

@inline DFS(t::T, f::Function) where {T <: AbstractTree} = DFS(t, f, identity, identity)

function num_nodes(t::T) where {T <: AbstractTree}
    tot = 0
    DFS(t, x->(tot+=1; 0))
    return tot
    #return 1 + sum(Int[num_nodes(i) for i in children(t)])
end

function max_depth(t::T) where {T <: AbstractTree}
    return DFS(t, x->nothing, x->nothing, x->(length(x[1])==0 ? 0 : maximum(x[1])+1))
end

function min_depth(t::T) where {T <: AbstractTree}
    return DFS(t, x->nothing, x->nothing, x->(length(x[1])==0 ? 0 : minimum(x[1])+1))
end

function collect_action_dfs(t::T, action::Function) where {T <: AbstractTree}
    return DFS( t,  x->nothing,
                    x->nothing,
                    x->vcat(x[1]...,Any[action(x[2]),]) )
end

function collect_action_bfs(t::T, action::Function) where {T <: AbstractTree}
    return DFS( t,  x->nothing,
                    x->nothing,
                    x->vcat(Any[action(x[2]),], x[1]...) )
end

function tree_print(
    t::T;
    propfunc=label,
    colfunc=x->(255,255,255),
    header=y->repeat("|--",max(0,y)),
    offset=0
    ) where {T <: AbstractTree}

    println_rgb(rgb_t) = (print("\e[1m\e[38;2;91;91;91;249m",rgb_t[2]);
                          println("\e[1m\e[38;2;$(rgb_t[1][1]);$(rgb_t[1][2]);$(rgb_t[1][3]);249m",rgb_t[3]);)
    nl = []
    level = [offset,]
    @inline s2s(x) = ((x isa Symbol) ? ":"*string(x) : string(x))
    f(x) = (push!(nl, [colfunc(x), header(level[end]), s2s(propfunc(x))]); push!(level,level[end]+1); 0)
    g(x) = (pop!(level); 0)

    DFS(t, f, g)

    println_rgb.(nl)
    return nothing
end

##: ------------------------------------------------
#: specially designed tree
##: ------------------------------------------------

import Base.copy
import Base.isequal

#: ------------  Block  -------------
mutable struct Block{TR} <: AbstractTree
    # number of repetitions
    n::Int
    # x
    # full range of the block
    x::IntRange
    # p
    # pattern for signle-line
    # undef for multi-line
    p::TPattern
    # R
    # identifier, can be pattern hash, or DFS traverse data
    R::TR
    # C
    # children
    C::Vector{Block{TR}}
    # DATA
    # for single line, store extracted data,
    # for multi line, empty
    # DATA is organized by parsee_File!() as follows
    # range => Vector{Pair{String,Any}}[...]
    # each v in [...] is a line of data from parsing single line
    DATA::DATATYPE
end

copy(s::Block{TR}) where TR = Block{TR}(s.n, s.x, copy(s.p), copy(s.R), copy.(s.C), copy(s.DATA))

isequal(s1::Block{TR},s2::Block{TR}) where TR = (s1.n==s2.n && s1.R==s2.R)

is_multi(s::Block{TR}) where TR = length(s.C)>0

is_single(s::Block{TR}) where TR = length(s.C)==0

@inline patt_dfs(t::Block) = collect_action_dfs(t, x->(is_single(x) ? x.p : __M1_PATT__))

@inline patt_bfs(t::Block) = collect_action_bfs(t, x->(is_single(x) ? x.p : __M1_PATT__))

#: ----------- hash -----------

#NOTE : the second parameter in hash(. , .) is crucial.
#> it guarantees that the folded blocks have exactly the same structure.
#> otherwise, when folding several blocks of similar structures with different repetition b.n
#> this information is lost after folding
khash(b::Block)::UInt64 = (is_single(b) ? hash(b.p,UInt(b.n)) : hash(khash.(b.C),UInt(b.n)))

is_valid_C(C::Vector{Block{TR}}) where TR = (length(C)>0 ? mapreduce(c->length(getfield(c,:x)), +, C)==length(concat0(C[1].x,C[end].x)) : true)

is_valid_x(M) = (is_single(M) ? length(M.x)==M.n : length(M.x)==M.n*sum([length(z.x) for z ∈ children(M)]))

function verify_block(b::Block{TR})::Bool where TR
    function h(x::Tuple{Vector,Block{TR}})
        if ! is_valid_x(x[2])
            @show x[2]
            return false
        else
            return all(x[1])
        end
    end
    DFS(b, x->nothing, x->nothing, h)
end

is_valid_block(b::Block{TR}) where TR = is_valid_C(children(b)) && verify_block(b)

#+ ========= compute_label ==========

compute_label(b::Block) = (khash(b), patt_dfs(b))

#+ ==================================

function Block(patt::TPattern, rg::IntRange)::Block{RTYPE}
    b = Block{RTYPE}(length(rg), rg, patt, __DEFAULT__R__, Block{RTYPE}[], [])
    b.R = compute_label(b)
    return b
end

Block(patt::TPattern, i::Int)::Block{RTYPE} = Block(patt, i:i)

Block()::Block{RTYPE} = Block(__DEFAULT_PATT__, 0)

function Block(C::Vector{Block{RT}})::Block{RT} where RT
    if length(C)==1  return C[1]  end
    rconcat = concat0(C[1].x, C[end].x)
    b = Block{RT}(1, rconcat, __DEFAULT_PATT__, __DEFAULT__R__, C, [])
    b.R = compute_label(b)
    return b
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
@testset "helper functions" begin
    @testset "_checked_rrule" begin
        try
            @eval cool(x,y) = x + y
            # We defined a 2-arg method for `cool` but no `rrule`
            ChainRules._checked_rrule(cool, 1.0, 2.0)
        catch e
            @test e isa ArgumentError
            @test e.msg == "can't differentiate `cool(::Float64, ::Float64)`; no " *
                "matching `rrule` is defined"
        end
    end
end


text:
function compile_clearswi(path;
        app_name="clearswi",
        filter_stdlibs=true,
        precompile_execution_file=abspath(joinpath(@__DIR__, "..", "test", "clearswi_test.jl")),
        clean=true,
        kw...)
    if !isdir(pathof("CLEARSWI"))
        download_pkg("CLEARSWI")
    end
    clearswiapppath = pathof("ClearswiApp")
    if !isdir(clearswiapppath)
        download_pkg("ClearswiApp")
    end
    create_app(clearswiapppath, path; app_name, filter_stdlibs, precompile_execution_file, kw...)
    test_clearswi(path, app_name)
    if clean
        clean_app(path, app_name) # remove unneccesary artifacts dir (600MB)
    end
    copy_matlab(path)
    printstyled("Success! CLEARSWI compiled and tested!\n"; color=:green)
    @warn("Relocatability has to be tested manually!")
end

function test_clearswi(path, app_name)
    file = tempname()
    phasefile = abspath(joinpath(@__DIR__, "..", "test", "data", "small", "Phase.nii"))
    magfile = abspath(joinpath(@__DIR__, "..", "test", "data", "small", "Mag.nii"))
    args = ["-p", phasefile, "-m", magfile, "-o", file, "-t", "1:3"]
    name = app_name * (Sys.iswindows() ? ".exe" : "")
    clearswifile = joinpath(path, "bin", name)
    @assert isfile(clearswifile)
    cmd = `$clearswifile $args`
    @assert success(run(cmd))
end

function update_clearswi()
    try
        rm(pathof("ClearswiApp"); force=true, recursive=true)
    catch
        @warn "Couldn't remove the old ClearswiApp folder! ($(pathof("ClearswiApp"))) Maybe it is opened in another App"
    end
    download_pkg("ClearswiApp")
end


text:
function circle!(c,r)
    θ = LinRange(0,2π,180)
    plot!(r*cos.(θ) .+ c[1], r*sin.(θ) .+ c[2], label = :none)
end

@userplot Circle
@recipe function f(cc::Circle)
    c = cc.args[1]
    r = cc.args[2]

    θ = LinRange(0,2π,180)
    x = r*cos.(θ) .+ c[1]
    y = r*sin.(θ) .+ c[2]
    @series begin
        x, y
    end
end