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

using Markdown
using InteractiveUtils

# ╔═╡ 485292d8-2bde-11eb-0a97-6b44d54596dd
using Plots

# ╔═╡ fc8ff7b8-2be1-11eb-0c9c-b9d2ee7c2880
using Markdown, Images

# ╔═╡ 2411c6ca-2bdd-11eb-050c-0399b3b0d7af
md"""
# Project 1: images and cellular automata

We wrap up the day by exploring convolutions and similar operations on 1-D and 2-D (or even n-D!) matrices. Departing from some basic building blocks we will cover signal processing, image processing and cellular automata!

Learning goals:
- code reuse;
- speed and performance (including optimizations);
- learning to plot and make animations.
"""

# ╔═╡ 3d9c3fea-2bde-11eb-0f09-cf11dcb07c0d
md"""
## 1-D operations
"""

# ╔═╡ e1147d4e-2bee-11eb-0150-d7af1f51f842
md"""
### 1-D convolution

Convolution:

$$y_i = \sum_j x_{j-m} k_{j}$$

Boundary effects.
"""

# ╔═╡ 6115f95c-2bde-11eb-01ce-3d17b78aa8a9
noisy_signal = 0:0.1:10 .|> x -> x * sin(x) + rand()  # use covid cases?

# ╔═╡ 8a996336-2bde-11eb-10a3-cb0046ed5de9
plot(noisy_signal)

# ╔═╡ a1f75f4c-2bde-11eb-37e7-2dc342c7032a
function convolve_1d!(x::Vector, out::Vector, w::Vector)
	n = length(x)
	m = div(length(w), 2)
	@assert length(w) % 2 == 1 "length of `w` has to be odd!"
	@assert n == length(out) "`x` and `out` have to be of the same length!"
	fill!(out, 0.0)
	for (i, xj) in enumerate(x)
		for (j, wj) in enumerate(w)
			k = clamp(i + j - m, 1, n)
			out[i] += w[j] * x[k]
		end
	end
	return out
end

# ╔═╡ 906841ee-2be1-11eb-03f4-49b5a8b47a6b
convolve_1d(x::Vector, w::Vector) = convolve_1d!(x, similar(x), w)

# ╔═╡ ef84b6fe-2bef-11eb-0943-034b8b90c4c4
function uniform_weights(m)
	return ones(2m+1) / (2m+1)
end

# ╔═╡ 294140a4-2bf0-11eb-22f5-858969a4640d
function triangle_weigths(m)
	w = zeros(2m+1)
	for i in 1:(m+1)
		w[i] = i
		w[end-i+1] = i
	end
	w ./= sum(w)
	return w
end

# ╔═╡ fa911a76-2be3-11eb-1c85-3df313eb0fcb
w = triangle_weigths(2)

# ╔═╡ 20c46656-2bf0-11eb-2dc4-af2cc2161f6a
scatter(w, xlabel="j", ylabel="weight")

# ╔═╡ cb5f4a20-2be8-11eb-0000-c59f23a024ef
signal_c = convolve_1d(noisy_signal, w)

# ╔═╡ eb5f8062-2be1-11eb-1e96-1f89be62f3b0
plot(signal_c)

# ╔═╡ af33214a-2be7-11eb-34f5-0d046c9db735
function generalized_convolve_1d!(x::Vector, out::Vector, m, f)
	n = length(x)
	@assert n == length(out) "`x` and `out` have to be of the same length!"
	fill!(out, 0.0)
	for (i, xj) in enumerate(x)
		if  m < i < n - m
			out[i] = f(@view(x[(i-m):(i+m)]))
		else
			out[i] = x[i]
		end
	end
	return out
end

# ╔═╡ 4263e920-2bf1-11eb-1db1-ed644a480a2d

# ╔═╡ 21ce0552-2be9-11eb-2359-152a7c75cc2e
vec_conv(x_s) = sum(x_s .* w)

# ╔═╡ 9bfd0e60-2bf1-11eb-3de4-03308da53e73
mean(x) = sum(x) / length(x)

# ╔═╡ 04a8eb54-2be4-11eb-06af-890cebf38da4
signal_gc = generalized_convolve_1d!(noisy_signal, similar(noisy_signal), 2, minimum)

# ╔═╡ b94aae80-2be8-11eb-0c88-cb3d2820e9ef
plot(signal_gc)

# ╔═╡ b7ba4ed8-2bf1-11eb-24ee-731940d1c29f
md"""
### Sliding window: secundary structure prediction

Chou-Fasman method
"""

# ╔═╡ c924a5f6-2bf1-11eb-3d37-bb63635624e9
spike_sars2 = "MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFSNVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIVNNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLEGKQGNFKNLREFVFKNIDGYFKIYSKHTPINLVRDLPQGFSALEPLVDLPIGINITRFQTLLALHRSYLTPGDSSSGWTAGAAAYYVGYLQPRTFLLKYNENGTITDAVDCALDPLSETKCTLKSFTVEKGIYQTSNFRVQPTESIVRFPNITNLCPFGEVFNATRFASVYAWNRKRISNCVADYSVLYNSASFSTFKCYGVSPTKLNDLCFTNVYADSFVIRGDEVRQIAPGQTGKIADYNYKLPDDFTGCVIAWNSNNLDSKVGGNYNYLYRLFRKSNLKPFERDISTEIYQAGSTPCNGVEGFNCYFPLQSYGFQPTNGVGYQPYRVVVLSFELLHAPATVCGPKKSTNLVKNKCVNFNFNGLTGTGVLTESNKKFLPFQQFGRDIADTTDAVRDPQTLEILDITPCSFGGVSVITPGTNTSNQVAVLYQDVNCTEVPVAIHADQLTPTWRVYSTGSNVFQTRAGCLIGAEHVNNSYECDIPIGAGICASYQTQTNSPRRARSVASQSIIAYTMSLGAENSVAYSNNSIAIPTNFTISVTTEILPVSMTKTSVDCTMYICGDSTECSNLLLQYGSFCTQLNRALTGIAVEQDKNTQEVFAQVKQIYKTPPIKDFGGFNFSQILPDPSKPSKRSFIEDLLFNKVTLADAGFIKQYGDCLGDIAARDLICAQKFNGLTVLPPLLTDEMIAQYTSALLAGTITSGWTFGAGAALQIPFAMQMAYRFNGIGVTQNVLYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVEAEVQIDRLITGRLQSLQTYVTQQLIRAAEIRASANLAATKMSECVLGQSKRVDFCGKGYHLMSFPQSAPHGVVFLHVTYVPAQEKNFTTAPAICHDGKAHFPREGVFVSNGTHWFVTQRNFYEPQIITTDNTFVSGNCDVVIGIVNNTVYDPLQPELDSFKEELDKYFKNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIKWPWYIWLGFIAGLIAIVMVTIMLCCMTSCCSCLKGCCSCGSCCKFDEDDSEPVLKGVKLHYT"

# ╔═╡ 4e6dedf0-2bf2-11eb-0bad-3987f6eb5481
md"""
### Elementary cellular automata
"""

# ╔═╡ 625c4e1e-2bf3-11eb-2c03-193f7d013fbe
function nextstate(l::Bool, s::Bool, r::Bool, rule::UInt8)
	bitstring(rule)[8-(4r+2s+1r)] == '1';
end

# ╔═╡ b03c60f6-2bf3-11eb-117b-0fc2a259ffe6
rule = 110 |> UInt8

# ╔═╡ a9ffc05c-2bf3-11eb-0771-b3732ed7d68d
md"""
Table for rule $rule:

| l      | s     | r    | result |
| :--------: | :-----: | :-----: | :--: |
| 10 | $true | $true | $false |

CAN BE DONE IN TABLE?

"""

# ╔═╡ 5f97da58-2bf4-11eb-26de-8fc5f19f02d2
Dict((l, s, r) => nextstate(l, s, r, rule)
	for l in [true, false]
	for s in [true, false]
	for r in [true, false]
				)

# ╔═╡ 924461c0-2bf3-11eb-2390-71bad2541463
function update1dca!(xnew, x, rule::UInt8)
	n = length(x)
	xnew[1] = nextstate(x[end], x[1], x[2], rule)
	xnew[end] = nextstate(x[end-1], x[end], x[1], rule)
	for i in 2:n-1
		xnew[i] = nextstate(x[i-1], x[i], x[i+1], rule)
	end
	return xnew
end

# ╔═╡ 21440956-2bf5-11eb-0860-11127d727282
next1dca(x, rule::UInt8) = update1dca!(similar(x), x, rule::UInt8)

# ╔═╡ 405a1036-2bf5-11eb-11f9-a1a714dbf7e1
x0_ca = rand(Bool, 200)

# ╔═╡ 6405f574-2bf5-11eb-3656-d7b9c94f145a
next1dca(x0_ca, rule)

# ╔═╡ 756ef0e0-2bf5-11eb-107c-8d1c65eacc45
"""
    simulate(x0, rule; nsteps::Integer=100)

Simulate `nsteps` time steps according to `rule` with `X0` as the initial condition.
Returns a matrix X, where the rows are the state vectors at different time steps.
"""
function simulate(x0, rule::UInt8; nsteps::Integer=100)
	n = length(x0)
    X = zeros(Bool, nsteps+1, n)
	X[1,:] = x0
	for t in 1:nsteps
		x = @view X[t,:]
		xnew = @view X[t+1,:]
		update1dca!(xnew, x, rule)
	end
    return X
end

# ╔═╡ e1dd7abc-2bf5-11eb-1f5a-0f46c7405dd5
X = simulate(x0_ca, rule; nsteps=200)

# ╔═╡ 9ab79ad4-2bf6-11eb-1d0e-bb6434b520a3
heatmap(X)

# ╔═╡ 9dbb9598-2bf6-11eb-2def-0f1ddd1e6b10
show_ca(X) = Gray.(X)

# ╔═╡ fb9a97d2-2bf5-11eb-1b92-ab884f0014a8
show_ca(X)

# ╔═╡ 4810a052-2bf6-11eb-2e9f-7f9389116950
all_cas = Dict(rule=>show_ca(simulate(x0_ca, rule; nsteps=200)) for rule in UInt8(0):UInt8(251))

# ╔═╡ 47e3c2f8-2bf6-11eb-1fbf-d5f3fc899056
all_cas[0]

# ╔═╡ aea84d6c-2be3-11eb-2fe7-63e0d0c28507
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("correct", "Keep working on it!", [text]))

# ╔═╡ dc2230ac-2be3-11eb-24d7-f980be57d697
keep_working()

# ╔═╡ 41e7ef30-2be4-11eb-1b9a-e7f7eae7a115
philip_file = download("https://i.imgur.com/VGPeJ6s.jpg")

# ╔═╡ 473c581c-2be5-11eb-1ddc-2d30a3468c8a
decimate(image, ratio=5) = image[1:ratio:end, 1:ratio:end]

# ╔═╡ 4c9d2e0e-2be4-11eb-2fb3-377f29141076
philip = let
	original = Images.load(philip_file)
	decimate(original, 8)
end

# ╔═╡ 608e930e-2bfc-11eb-09f2-29600cfba1e6
function convolve_2d!(M::AbstractMatrix, out::AbstractMatrix, K::AbstractMatrix)
	n_rows, n_cols = size(M)
	fill!(out, 0.0)
	m = div(size(K, 1), 2)
	for i in 1:n_rows
		for j in 1:n_cols
			for k in -m:m
				for l in -m:m
					out[i,j] += M[clamp(i+k, 1, n_rows), clamp(j+l, 1, n_cols)] * K[k+m+1,l+m+1]
				end
			end
		end
	end
	return out
end

# ╔═╡ 245deb30-2bfe-11eb-3e47-6b7bd5d73ccf
function gaussian_kernel(m; σ)
	K = [exp(-(x^2 + y^2) / σ^2) for x in -m:m, y in -m:m]
	K ./= sum(K)
	return K
end

# ╔═╡ 6278ae28-2bfe-11eb-3c36-9bf8cc703a3a
K = gaussian_kernel(1, σ=10)

# ╔═╡ 6c40b89c-2bfe-11eb-3593-9d9ef6f45b80
heatmap(K)

# ╔═╡ 73c3869a-2bfd-11eb-0597-29cd6c1c90db
convolve_2d(M, K) = convolve_2d!(M, similar(M), K)

# ╔═╡ a34f1e4c-2bfd-11eb-127e-295a7039c301
convolve_2d(randn(1000, 1000), K)

# ╔═╡ c80e4362-2bfe-11eb-2055-5fbf167be551
K_sharpen = [0 -1 0;
			-1 5 -1;
			0 -1 0]

# ╔═╡ 84880dfc-2bfd-11eb-32b7-47750e1c696b
function convolve_image(M, K)
	Mred = convolve_2d(red.(M), K)
	Mgreen = convolve_2d(green.(M), K)
	Mblue = convolve_2d(blue.(M), K)
	return RGB.(Mred, Mgreen, Mblue)
end

# ╔═╡ 8c62855a-2bff-11eb-00f5-3702e142f76f
const Gx = [1 0 -1;
		   2 0 -2;
		   1 0 -1]

# ╔═╡ ae8d6258-2bff-11eb-0182-afc33d9d4efc
const Gy = [1 2 1;
		    0 0 0;
		   -1 -2 -1]

# ╔═╡ 746d3578-2bff-11eb-121a-f3d17397a49e
function edge_detection(M)
	M = Gray.(M) |> Matrix
	return convolve_2d(M, Gx).^2 + convolve_2d(M, Gy).^2
end

# ╔═╡ f41919a4-2bff-11eb-276f-a5edd5b374ac
edge_detection(philip)

# ╔═╡ e1f23724-2bfd-11eb-1063-1d59909bcacc
convolve_image(philip, K)

# ╔═╡ 55bfb46a-2be5-11eb-2be6-1bc153db11ac
RGB.(red.(philip), green.(philip), blue.(philip))

# ╔═╡ 04f3c460-2bed-11eb-1039-238bfb3ec445
images = [rand([RGB(0, 0, 1),RGB(1, 0, 0)], 50, 100) for i in 1:100]

# ╔═╡ b3eed0e0-2bed-11eb-21d3-19323b600edf
images[1]

# ╔═╡ 2a3b0be4-2bee-11eb-3e26-f1d850c752b6
@gif for i ∈ 1:100
	plot(rand([RGB(0, 0, 1),RGB(1, 0, 0)], 100, 100))
end

# ╔═╡ Cell order:
# ╠═2411c6ca-2bdd-11eb-050c-0399b3b0d7af
# ╠═3d9c3fea-2bde-11eb-0f09-cf11dcb07c0d
# ╠═e1147d4e-2bee-11eb-0150-d7af1f51f842
# ╠═6115f95c-2bde-11eb-01ce-3d17b78aa8a9
# ╠═8a996336-2bde-11eb-10a3-cb0046ed5de9
# ╠═485292d8-2bde-11eb-0a97-6b44d54596dd
# ╠═a1f75f4c-2bde-11eb-37e7-2dc342c7032a
# ╠═fa911a76-2be3-11eb-1c85-3df313eb0fcb
# ╠═20c46656-2bf0-11eb-2dc4-af2cc2161f6a
# ╠═906841ee-2be1-11eb-03f4-49b5a8b47a6b
# ╠═cb5f4a20-2be8-11eb-0000-c59f23a024ef
# ╠═eb5f8062-2be1-11eb-1e96-1f89be62f3b0
# ╠═ef84b6fe-2bef-11eb-0943-034b8b90c4c4
# ╠═294140a4-2bf0-11eb-22f5-858969a4640d
# ╠═af33214a-2be7-11eb-34f5-0d046c9db735
# ╠═4263e920-2bf1-11eb-1db1-ed644a480a2d
# ╠═21ce0552-2be9-11eb-2359-152a7c75cc2e
# ╠═9bfd0e60-2bf1-11eb-3de4-03308da53e73
# ╠═04a8eb54-2be4-11eb-06af-890cebf38da4
# ╠═b94aae80-2be8-11eb-0c88-cb3d2820e9ef
# ╠═b7ba4ed8-2bf1-11eb-24ee-731940d1c29f
# ╠═c924a5f6-2bf1-11eb-3d37-bb63635624e9
# ╠═4e6dedf0-2bf2-11eb-0bad-3987f6eb5481
# ╠═625c4e1e-2bf3-11eb-2c03-193f7d013fbe
# ╠═b03c60f6-2bf3-11eb-117b-0fc2a259ffe6
# ╠═a9ffc05c-2bf3-11eb-0771-b3732ed7d68d
# ╠═5f97da58-2bf4-11eb-26de-8fc5f19f02d2
# ╠═924461c0-2bf3-11eb-2390-71bad2541463
# ╠═21440956-2bf5-11eb-0860-11127d727282
# ╠═405a1036-2bf5-11eb-11f9-a1a714dbf7e1
# ╠═6405f574-2bf5-11eb-3656-d7b9c94f145a
# ╠═756ef0e0-2bf5-11eb-107c-8d1c65eacc45
# ╠═e1dd7abc-2bf5-11eb-1f5a-0f46c7405dd5
# ╠═9ab79ad4-2bf6-11eb-1d0e-bb6434b520a3
# ╠═9dbb9598-2bf6-11eb-2def-0f1ddd1e6b10
# ╠═fb9a97d2-2bf5-11eb-1b92-ab884f0014a8
# ╠═4810a052-2bf6-11eb-2e9f-7f9389116950
# ╠═47e3c2f8-2bf6-11eb-1fbf-d5f3fc899056
# ╠═fc8ff7b8-2be1-11eb-0c9c-b9d2ee7c2880
# ╠═aea84d6c-2be3-11eb-2fe7-63e0d0c28507
# ╠═dc2230ac-2be3-11eb-24d7-f980be57d697
# ╠═41e7ef30-2be4-11eb-1b9a-e7f7eae7a115
# ╠═4c9d2e0e-2be4-11eb-2fb3-377f29141076
# ╠═473c581c-2be5-11eb-1ddc-2d30a3468c8a
# ╠═608e930e-2bfc-11eb-09f2-29600cfba1e6
# ╠═245deb30-2bfe-11eb-3e47-6b7bd5d73ccf
# ╠═6278ae28-2bfe-11eb-3c36-9bf8cc703a3a
# ╠═6c40b89c-2bfe-11eb-3593-9d9ef6f45b80
# ╠═73c3869a-2bfd-11eb-0597-29cd6c1c90db
# ╠═a34f1e4c-2bfd-11eb-127e-295a7039c301
# ╠═c80e4362-2bfe-11eb-2055-5fbf167be551
# ╠═84880dfc-2bfd-11eb-32b7-47750e1c696b
# ╠═8c62855a-2bff-11eb-00f5-3702e142f76f
# ╠═ae8d6258-2bff-11eb-0182-afc33d9d4efc
# ╠═746d3578-2bff-11eb-121a-f3d17397a49e
# ╠═f41919a4-2bff-11eb-276f-a5edd5b374ac
# ╠═e1f23724-2bfd-11eb-1063-1d59909bcacc
# ╠═55bfb46a-2be5-11eb-2be6-1bc153db11ac
# ╠═04f3c460-2bed-11eb-1039-238bfb3ec445
# ╠═b3eed0e0-2bed-11eb-21d3-19323b600edf
# ╠═2a3b0be4-2bee-11eb-3e26-f1d850c752b6


text:
function get_qzss_orbital_elements(obj::QZSSPropagator)
    return jcall(obj, "getQZSSOrbitalElements", QZSSOrbitalElements, ())
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PersistentVolumeStatus is the current status of a persistent volume.

    IoK8sApiCoreV1PersistentVolumeStatus(;
        message=nothing,
        phase=nothing,
        reason=nothing,
    )

    - message::String : A human-readable message indicating details about why the volume is in this state.
    - phase::String : Phase indicates if a volume is available, bound to a claim, or released by a claim. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#phase
    - reason::String : Reason is a brief CamelCase string that describes any failure and is meant for machine parsing and tidy display in the CLI.
"""
mutable struct IoK8sApiCoreV1PersistentVolumeStatus <: SwaggerModel
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason

    function IoK8sApiCoreV1PersistentVolumeStatus(;message=nothing, phase=nothing, reason=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeStatus, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1PersistentVolumeStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        validate_property(IoK8sApiCoreV1PersistentVolumeStatus, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeStatus

const _property_map_IoK8sApiCoreV1PersistentVolumeStatus = Dict{Symbol,Symbol}(Symbol("message")=>Symbol("message"), Symbol("phase")=>Symbol("phase"), Symbol("reason")=>Symbol("reason"))
const _property_types_IoK8sApiCoreV1PersistentVolumeStatus = Dict{Symbol,String}(Symbol("message")=>"String", Symbol("phase")=>"String", Symbol("reason")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeStatus }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeStatus[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolumeStatus }, name::Symbol, val)
end


text:
function generate(obj::ConvexHullGenerator2D, arg0::Collection)
    return jcall(obj, "generate", ConvexHull2D, (Collection,), arg0)
end


text:
using StokesDiffEq
using Base.Test

##Stokes
println("Stokes Tests")
@time @test include("stokes/stokes_tests.jl")
println("DGS Internals Test")
@time @test include("stokes/stokes_dgs_tests.jl")


text:
@testset "TAMU Parsing " begin
    # Added due to inconsistencies in the incoming data. If the data is fixed, this test_logs call will fail
    sys = PSB.build_system(PSB.PSYTestSystems, "tamu_ACTIVSg2000_sys"; force_build = true)
    @test isa(sys, PSY.System)
    # Test bad input
    @test_throws PowerSystems.DataFormatError TamuSystem("")
end


text:
function plot_price_ma_trade_data(index::Int64, duration::Int64, window::Int64)
    # Retrieve data from various helper functions
    Price_df, _ , Vol_df = get_price_data_single(currencies[index])

    # Make sure that duration does not exceed the number of rows - max(windows) in the DataFrame
    # This allows calculation of MA for the longest duration
    if duration > size(Price_df)[1]-maximum(windows)
        duration = size(Price_df)[1]-maximum(windows)
    end

    ################# Daily average data #################
    trace1 = PlotlyJS.scatter(;x = Price_df[1:duration,:Date], y = Price_df[1:duration,2], mode="markers+lines", name = "$(currencies[index]) price")

    ################# Moving averages data #################
    sort!(Price_df, :Date) # oldest date first, newest at the bottom
    Price_SMA, Price_WMA, Price_EMA = moving_averages(Price_df, duration, window)

    trace2 = PlotlyJS.scatter(;x = Price_df[!,:Date][end-length(Price_SMA)+1:end], y = Price_SMA, mode="lines", name = "$(names(Price_df)[2]) SMA over $(window) days")
    trace3 = PlotlyJS.scatter(;x = Price_df[!,:Date][end-length(Price_WMA)+1:end], y = Price_WMA, mode="lines", name = "$(names(Price_df)[2]) WMA over $(window) days")
    trace4 = PlotlyJS.scatter(;x = Price_df[!,:Date][end-length(Price_EMA)+1:end], y = Price_EMA, mode="lines", name = "$(names(Price_df)[2]) EMA over $(window) days")

    ################# Daily trade data #################
    sort!(Price_df, :Date, rev = true) # newest date first, oldest at the bottom
    trace5 = PlotlyJS.bar(;x = Price_df[1:duration,:Date], y = Price_df[1:duration,2].*Vol_df[1:duration,2] , mode="markers+lines", name = "$(currencies[index]) daily trade")

    return trace1, trace2, trace3, trace4, trace5
end

function plot_candle_vol_data(index::Int64, duration::Int64)

    # Retrieve data from various helper functions
    _ , Candle_df, Vol_df = get_price_data_single(currencies[index])

    # Make sure that duration does not exceed the number of rows in the DataFrame
    if duration > size(Candle_df)[1]
        duration = size(Candle_df)[1]
    end

    ################# Daily candlestick data #################
    open_col, high_col, low_col, close_col = [Float64[] for i = 1:4]

    for i = 1:duration
        push!(open_col, Candle_df[!,2][i][1])
        push!(high_col, Candle_df[!,2][i][2])
        push!(low_col, Candle_df[!,2][i][3])
        push!(close_col, Candle_df[!,2][i][4])
    end

    trace1 = PlotlyJS.candlestick(; x = Candle_df[1:duration,:Date], open = open_col, high = high_col, low = low_col,
                                   close = close_col, name = "$(currencies[index])")

    ################# Daily volume data #################
    trace2 = PlotlyJS.bar(;x = Vol_df[1:duration,:Date], y = Vol_df[1:duration,2], name = "$(currencies[index]) volume")

    return trace1, trace2
end

function plot_cumul_daily_return_hist(index::Int64, duration::Int64)

    # Retrieve data from various helper functions
    Price_df , _ , _ = get_price_data_single(currencies[index])

    # Make sure that duration does not exceed the number of rows in the DataFrame
    if duration > size(Price_df)[1]
        duration = size(Price_df)[1]
    end

    # Reverse the order (oldest date first, newest at the bottom)
    sort!(Price_df, :Date)

    ################# Cumulative return #################

    trace1 = PlotlyJS.scatter(; x = Price_df[end-duration+2:end,:Date], y = (cumsum(diff(Price_df[end-duration+1:end,2])) ./ Price_df[end-duration+1, 2]) .* 100,
                                mode="markers+lines", name = "$(currencies[index]) cumulative return")

    ################# Daily return #################

    X = Price_df[end-duration+2:end,:Date]
    Y = (diff(Price_df[end-duration+1:end,2]) ./ Price_df[end-duration+1:end-1,2]) .* 100

    # Split into two datasets (green: positive change, red: negative change)
    green_Y = Y[Y .≥ 0.0]
    green_X = X[Y .≥ 0.0]

    red_Y = Y[Y .< 0.0]
    red_X = X[Y .< 0.0]

    green_share = round((length(green_Y)/length(Y))*100, digits = 2)
    red_share = 100.0 - green_share

    trace2_green = PlotlyJS.bar(; x = green_X, y = green_Y, marker_color = "green", name = "$(currencies[index]) increase, share = $(green_share) %")
    trace2_red   = PlotlyJS.bar(; x = red_X, y = red_Y, marker_color = "red", name = "$(currencies[index]) decrease, share = $(red_share) %")

    ################# Daily change histogram / Daily volatility #################

    σ = round(std(Y), digits = 2)
    trace3 = plotly_hist(Y, 75; normalize = false)

    return trace1, trace2_green, trace2_red, trace3, σ, duration
end

function plotly_hist(x::AbstractVector{T}, nbins::Integer; normalize::Bool = true) where T <: Number

    # use StatsBase to create a histogram object
    hist = StatsBase.fit(Histogram, x, nbins = nbins, closed = :left)

    # obtain bar positions -> center of each interval
    bins = similar(x, length(hist.edges[1]) - 1)
    edges = hist.edges[1]

    for k in eachindex(bins)
        bins[k] = (edges[k] + edges[k + 1]) * 0.5
    end

    if normalize
        y = hist.weights ./ length(x) # we need a new array
    else
        y = hist.weights
    end

    trace = bar(x = bins, y = y)

    return trace
end

function plot_fcas_data(index::Int64)

    us, fs, ds, ms, fr = get_ratings_data(currencies[index])

    ################# FCAS metrics data #################
    trace1 = PlotlyJS.bar(;x = ["Utility", "FCAS", "Developer", "Market maturity"], y = [us, fs, ds, ms], width = 0.25)
    return trace1, fr
end

function plot_macd_signal(index::Int64, duration::Int64)

    # Retrieve data from various helper functions
    Price_df, _ , _ = get_price_data_single(currencies[index])

    # Make sure that duration does not exceed the number of rows in the DataFrame
    if duration > size(Price_df)[1]
        duration = size(Price_df)[1]
    end

    sort!(Price_df, :Date)                            # oldest date first, newest at the bottom

    Price_df = Price_df[end-duration+1-26-9+1:end, :] # filter based on selected duration and effective
                                                      # window size of 26+9

    df_ema_all = calculate_macd(Price_df)             # get EMA and MACD data into a DataFrame

    ################# Daily average, EMA-12 and EMA-26 data #################
    trace1 = PlotlyJS.scatter(;x = df_ema_all[!,:date], y = df_ema_all[!,:Raw],
                              mode="markers+lines", name = "Daily average")

	trace2 = PlotlyJS.scatter(;x = df_ema_all[!,:date], y = df_ema_all[!,:EMA_long],
                              mode="markers+lines", name = "EMA-26")

	trace3 = PlotlyJS.scatter(;x = df_ema_all[!,:date], y = df_ema_all[!,:EMA_short],
                              mode="markers+lines", name = "EMA-12")

    ################# MACD and its signal data #################
    trace4 = PlotlyJS.scatter(;x = df_ema_all[!,:date], y = df_ema_all[!,:MACD],
                              mode="markers+lines", name = "MACD")

	trace5 = PlotlyJS.scatter(;x = df_ema_all[!,:date], y = df_ema_all[!,:Signal],
                              mode="markers+lines", name = "Signal (EMA-9)")

    ################# Distance between MACD and its signal line #################

    X = df_ema_all[!,:date]
	Y = df_ema_all[!,:MACD] - df_ema_all[!,:Signal]

	# Split into two datasets (green: positive change, red: negative change)
	green_Y = Y[Y .≥ 0.0]
    green_X = X[Y .≥ 0.0]

    red_Y = Y[Y .< 0.0]
    red_X = X[Y .< 0.0]

	trace6_green = PlotlyJS.bar(; x = green_X, y = green_Y, marker_color = "green",
                                  name = "MACD above signal ")

    trace6_red   = PlotlyJS.bar(; x = red_X, y = red_Y, marker_color = "red",
                                  name = "MACD below signal")

    return trace1, trace2, trace3, trace4, trace5, trace6_green, trace6_red

end

function plot_linear_regression(index::Int64, duration::Int64)

    # Retrieve data from various helper functions
    Price_df, _ , _ = get_price_data_single(currencies[index])

    # Make sure that duration does not exceed the number of rows in the DataFrame
    if duration > size(Price_df)[1]
        duration = size(Price_df)[1]
    end

    # Filter on duration and create index column (model fit does not work with dates)
	df_fit = sort!(Price_df[1:duration, :])
	df_fit.Index = 1:size(df_fit)[1]

	# Rename column to price to make df_fit generic
	rename!(df_fit, Dict(Symbol(names(df_fit)[2]) => "Price"))

	# Get model parameters and R²
	model_params = lm(@formula(Price ~ Index), df_fit)
    R² = round(r2(model_params), digits = 2)

	# Predicted data
	price_from_model = Float64.(predict(model_params, df_fit))

	# Get standard deviation σ
	σ = std(price_from_model)

    ################# Raw price data #################

    trace1 = PlotlyJS.scatter(;x = df_fit[!,:Date], y = df_fit[!,:Price],
                              mode="markers+lines", name = "Actual price")

    ################# Linear regression channel #################

	trace2 = PlotlyJS.scatter(;x = df_fit[!,:Date], y = price_from_model,
                              mode="lines", name = "Linear regression line")

	trace3 = PlotlyJS.scatter(;x = df_fit[!,:Date], y = price_from_model .+ 2*σ,
                              mode="markers", name = "Upper channel (+2σ)")

	trace4 = PlotlyJS.scatter(;x = df_fit[!,:Date], y = price_from_model .- 2*σ,
                              mode="markers", name = "Lower channel (-2σ)")

    return trace1, trace2, trace3, trace4, R²
end

function plot_dev_comm_data(index::Int64)

    # Convert currency symbol to lowercase and fetch data from CoinGecko
    df_dev, df_comm = get_dev_comm_data(lowercase(currencies[index]))

    ################# Developer data #################
    trace1 = PlotlyJS.bar(; x = df_dev[!, :Metric], y = df_dev[!, :Value],
                            name = "Developer data")

    ################# Community data #################
    trace2 = PlotlyJS.bar(; x = df_comm[!, :Metric], y = df_comm[!, :Value],
                            name = "Community data")

    return trace1, trace2
end

function plot_exchange_vol_data(index::Int64, num_exchanges::Int64 = 10)

    # Convert currency symbol to lowercase and fetch data from CoinGecko
    df_ex_vol = get_exchange_vol_data(lowercase(currencies[index]), num_exchanges)

    ################# Coin volume data #################
    trace1 = PlotlyJS.bar(; x = df_ex_vol[!, :Name], y = df_ex_vol[!, :Coin_volume],
                            name = "Volume data in coins")

    ################# USD volume data #################
    trace2 = PlotlyJS.bar(; x = df_ex_vol[!, :Name], y = df_ex_vol[!, :USD_volume],
                            name = "Volume data in USD")

    return trace1, trace2
end

function plot_overall_vol_data(duration::Int64, num_exchanges::Int64 = 10)

    # Collect all traces for all exchanges
    all_traces = GenericTrace{Dict{Symbol, Any}}[]

    # Fetch overall volume data from CoinGecko for given historical duration
    df_ex_vol = get_overall_vol_data(duration, num_exchanges)

    if ~isempty(df_ex_vol)

        # First column is for duration
        exchanges = names(df_ex_vol)[2:end]

        for i = 1:length(exchanges)

            trace = PlotlyJS.bar(;x = df_ex_vol[!,:Time], y = df_ex_vol[!,i+1],
            mode="markers+lines", name = "$(exchanges[i])")

            push!(all_traces, trace)
        end
    end

    return all_traces
end


text:
# composite type for Bordered arrays
import Base: copy, copyto!, eltype, zero
import LinearAlgebra: norm, dot, length, similar, axpy!, axpby!, rmul!, mul!

"""
	x = BorderedArray(vec1, vec2)

This defines an array (although not `<: AbstractArray`) to hold two arrays or an array and a scalar. This is useful when one wants to add constraints (phase, ...) to a functional for example. It is used throughout the package for the Pseudo Arc Length Continuation, for the continuation of Fold / Hopf points, for periodic orbits... It is also used to define periodic orbits as (orbit, period). As such, it is a convenient alternative to `cat`, `vcat` and friends. We chose not make it a subtype of AbstractArray as we wish to apply the current package to general "arrays", see [Requested methods for Custom State](@ref). Finally, it proves useful for the GPU where the operation `x[end]` can be slow.
"""
mutable struct BorderedArray{vectype1, vectype2}
	u::vectype1
	p::vectype2
end

eltype(b::BorderedArray{vectype, T}) where {T, vectype} = eltype(b.p)
similar(b::BorderedArray{vectype, T}, ::Type{S} = eltype(b)) where {S, T, vectype} = BorderedArray(similar(b.u, S), similar(b.p, S))
similar(b::BorderedArray{vectype, T}, ::Type{S} = eltype(b)) where {S, T <: Real, vectype} = BorderedArray(similar(b.u, S), S(0))

# a version of copy which cope with our requirements concerning the methods
# available for
_copy(b) = copyto!(similar(b), b)
copy(b::BorderedArray) = copyto!(similar(b), b)# BorderedArray(copy(b.u), copy(b.p))

copyto!(dest::BorderedArray{vectype, T1}, src::BorderedArray{vectype, T2}) where {vectype, T1, T2 } = (copyto!(dest.u, src.u); copyto!(dest.p, src.p); dest)
copyto!(dest::BorderedArray{vectype, T1}, src::BorderedArray{vectype, T2}) where {vectype, T1 <: Number, T2 <: Number} = (copyto!(dest.u, src.u); dest.p = src.p;dest)

length(b::BorderedArray{vectype, T}) where {vectype, T} = length(b.u) + length(b.p)
length(b::BorderedArray{vectype, T}) where {vectype, T <: Number} = length(b.u) + 1

dot(a::BorderedArray{Tv1, Tp1}, b::BorderedArray{Tv2, Tp2}) where {Tv1, Tv2, Tp1, Tp2} = dot(a.u, b.u) + dot(a.p, b.p)
norm(b::BorderedArray{vectype, T}, p::Real) where {vectype, T} = max(norm(b.u, p), norm(b.p, p))
zero(b::BorderedArray{vectype, T}) where {vectype, T} = BorderedArray(zero(b.u), zero(b.p))
################################################################################
function rmul!(A::BorderedArray{vectype, Tv}, a::T, b::T) where {vectype, T <: Number, Tv}
	# Scale an array A by a scalar b overwriting A in-place
	rmul!(A.u, a)
	rmul!(A.p, b)
	return A
end

function rmul!(A::BorderedArray{vectype, Tv}, a::T, b::T) where {vectype, Tv <: Number, T <: Number }
	# Scale an array A by a scalar b overwriting A in-place
	rmul!(A.u, a)
	A.p = A.p * b
	return A
end

rmul!(A::BorderedArray{vectype, Tv}, a::T) where {vectype, T <: Number, Tv} = rmul!(A, a, a)
################################################################################
function mul!(A::BorderedArray{Tv1, Tp1}, B::BorderedArray{Tv2, Tp2}, α::T) where {Tv1, Tv2, Tp1, Tp2, T <: Number}
	mul!(A.u, B.u, α)
	mul!(A.p, B.p, α)
	return A
end

function mul!(A::BorderedArray{Tv1, Tp1}, B::BorderedArray{Tv2, Tp2}, α::T) where {Tv1, Tv2, Tp1 <: Number, Tp2 <: Number, T <: Number}
	mul!(A.u, B.u, α)
	A.p = B.p * α
	return A
end

mul!(A::BorderedArray{Tv1, Tp1}, α::T, B::BorderedArray{Tv2, Tp2}) where {Tv1, Tv2, Tp1, Tp2, T} = mul!(A, B, α)
################################################################################
function axpy!(a::T, X::BorderedArray{Tv1, Tp1}, Y::BorderedArray{Tv2, Tp2}) where {Tv1, Tv2, T <: Number, Tp1, Tp2}
	# Overwrite Y with a*X + Y, where a is scalar
	axpy!(a, X.u, Y.u)
	axpy!(a, X.p, Y.p)
	return Y
end

function axpy!(a::T, X::BorderedArray{Tv1, Tp1}, Y::BorderedArray{Tv2, Tp2}) where {Tv1, Tv2, T <: Number, Tp1 <: Number, Tp2 <: Number}
	# Overwrite Y with a*X + Y, where a is scalar
	axpy!(a, X.u, Y.u)
	Y.p = a * X.p + Y.p
	return Y
end
################################################################################
function axpby!(a::T, X::BorderedArray{vectype, Tv1}, b::T, Y::BorderedArray{vectype, Tv2}) where {vectype, T <: Number, Tv1, Tv2}
	# Overwrite Y with a * X + b * Y, where a, b are scalar
	axpby!(a, X.u, b, Y.u)
	axpby!(a, X.p, b, Y.p)
	return Y
end

function axpby!(a::T, X::BorderedArray{vectype, Tv1}, b::T, Y::BorderedArray{vectype, Tv2}) where {vectype, Tv1 <: Number, Tv2 <: Number, T <: Number}
	# Overwrite Y with a * X + b * Y, where a is a scalar
	axpby!(a, X.u, b, Y.u)
	Y.p = a * X.p + b * Y.p
	return Y
end
################################################################################
# this function is actually axpy!(-1, y, x)
#
# 	`minus!(x, y)`
#
# computes x-y into x and returns x
minus!(x, y) = axpy!(convert(eltype(x), -1), y, x)
minus!(x::vec, y::vec) where {vec <: AbstractArray} = (x .= x .- y)
minus!(x::T, y::T) where {T <: Number} = (x = x - y)
minus!(x::BorderedArray{vectype, T}, y::BorderedArray{vectype, T}) where {vectype, T} = (minus!(x.u, y.u); minus!(x.p, y.p))
function minus!(x::BorderedArray{vectype, T}, y::BorderedArray{vectype, T}) where {vectype, T <: Number}
	minus!(x.u, y.u)
	# Carefull here. If I use the line below, then x.p will be left unaffected
	# minus_!(x.p, y.p)
	x.p = x.p - y.p
	return x
end
################################################################################
#
#	`minus(x,y)`
#
# returns x - y
minus(x, y) = (x1 = copyto!(similar(x), x); minus!(x1, y); return x1)
# minus(x, y) = (x - y)
minus(x::vec, y::vec) where {vec <: AbstractArray} = (return x .- y)
minus(x::T, y::T) where {T <:Real} = (return x - y)
minus(x::BorderedArray{vectype, T}, y::BorderedArray{vectype, T}) where {vectype, T} = (return BorderedArray(minus(x.u, y.u), minus(x.p, y.p)))
minus(x::BorderedArray{vectype, T}, y::BorderedArray{vectype, T}) where {vectype, T <: Number} = (return BorderedArray(minus(x.u, y.u), x.p - y.p))


text:
"""
    iswrsparse(S)

Returns `true` if `S` is sparse, and `false` otherwise.

# Examples
```jldoctest
julia> sv = sparsevec([1, 4], [2.3, 2.2], 10)
10-element SparseVector{Float64,Int64} with 2 stored entries:
  [1 ]  =  2.3
  [4 ]  =  2.2

julia> iswrsparse(sv)
true

julia> iswrsparse(Array(sv))
false

```
"""
iswrsparse(::T) where T<:AbstractArray = iswrsparse(T)
iswrsparse(T::Type) = walk_wrapper(_iswrsparse, T)
_iswrsparse(::Type) = false
_iswrsparse(::Type{<:AbstractSparseArray}) = true

indextype(::T) where T<:AbstractArray = indextype(T)
indextype(T::Type) = walk_wrapper(_indextype, T)
_indextype(::Type{<:AbstractSparseArray{<:Any,Ti}}) where Ti = Ti
_indextype(::Type) = Int

nnz_estimation(A::AbstractMatrix) = walk_wrapper(_nnz_estimation, A)
_nnz_estimation(::Any) = missing
_nnz_estimation(A::SparseMatrixCSC) = nnz(A)
_nnz_estimation(A::SparseVector) = nnz(A)
_nnz_estimation(A::Array) = length(A)

depth(::T) where T<:AbstractArray = depth(T)
depth(T::Type) = walk_wrapper(_depth, T, x->x+1)
_depth(::Type) = 0

array_storage(::T) where T<:AbstractArray = array_storage(T)
array_storage(T::Type) = walk_wrapper(_array_storage, T)
_array_storage(T::Type{<:AbstractArray}) = T
_array_storage(T::Type) = Nothing

import LinearAlgebra: Symmetric, Hermitian, LowerTriangular, UnitLowerTriangular
import LinearAlgebra: UpperTriangular, UnitUpperTriangular, Transpose, Adjoint

isupper(A::AbstractArray) = walk_wrapper(_isupper, A)
_isupper(::Any) = missing
_isupper(::Union{UpperTriangular,UnitUpperTriangular}) = true
_isupper(A::Union{Symmetric,Hermitian}) = A.uplo == 'U'
_islower(::Union{LowerTriangular,UnitLowerTriangular}) = true
_islower(A::Union{Symmetric,Hermitian}) = A.uplo == 'L'

walk_wrapper(f::Function, x::Any, g::Function=identity) = f(x)
for wr in (Symmetric, Hermitian,
           LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular,
           Transpose, Adjoint,
           SubArray, Conjugate,
           Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal, HermiteTridiagonal )

    pl = wr === SubArray ? :($wr{<:Any,<:Any,T}) : :($wr{<:Any,T})
    @eval function walk_wrapper(f::Function, ::Type{<:$pl}, g::Function=identity) where T
        g(walk_wrapper(f, T, g))
    end
    @eval function walk_wrapper(f::Function, A::$pl) where T
        res = f(A)
        ismissing(res) && ( res = walk_wrapper(f, parent(A)) )
        ifelse(ismissing(res), false, res)
    end
end

"""
    inflate(A::AbstractMatrix)

Reduce depth in the case of nested wrappers of an abstract array.
If `dept(A) <= 1` convert to SparseMatrixCSC or Array.
Otherwise convert the deepest parents of the nesting and leave the rest as is.
"""
inflate(A::AbstractArray) = _inflate(A)
_inflate(A::AbstractArray) = iswrsparse(A) ? sparsecsc(A) : A isa Array ? A : Array(A)
for ty in (Symmetric, Hermitian)
    @eval function inflate(A::$ty)
        depth(A) == 1 ? _inflate(A) : $ty(inflate(parent(A)), up(A))
    end
end
function inflate(A::SubArray)
    depth(A) == 1 ? _inflate(A) : SubArray(inflate(parent(A)), A.indices)
end
for ty in ( LowerTriangular, UnitLowerTriangular,
            UpperTriangular, UnitUpperTriangular,
            Conjugate, Transpose, Adjoint)

    @eval function inflate(A::$ty)
        depth(A) == 1 ? _inflate(A) : $ty(inflate(parent(A)))
    end
end
for ty in (Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal, HermiteTridiagonal)
    @eval inflate(A::$ty) = dropzeros!(sparse(A))
end

import SparseArrays.SparseMatrixCSC
for wr in (Symmetric, Hermitian,
           UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular,
           Transpose, Adjoint, SubArray, Conjugate,
           Bidiagonal, Tridiagonal, SymTridiagonal, HermiteTridiagonal)

    @eval SparseMatrixCSC(A::$wr) = sparsecsc(A)
    @eval SparseMatrixCSC{Tv}(A::$wr{Tv}) where Tv = sparsecsc(A)
    @eval SparseMatrixCSC{Tv}(A::$wr) where Tv = SparseMatrixCSC{Tv}(sparsecsc(A))
    @eval SparseMatrixCSC{Tv,Ti}(A::$wr) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(sparsecsc(A))
end

"""
    unwrap(A::AbstractMatrix)

In case A is a wrapper type (`SubArray, Symmetric, Adjoint, SubArray, Triangular, Tridiagonal`, etc.)
convert to `Matrix` or `SparseMatrixCSC`, depending on final storage type of A.
For other types return A itself.
"""
unwrap(A::AbstractArray) = iswrsparse(A) ? convert(SparseMatrixCSC, A) : convert(Array, A)
unwrap(A::Any) = A

import Base.copy
copy(A::SubArray) = getindex(unwrap(parent(A)), A.indices...)

"""
    sparsecsc(A::AbstractArray)

Return `A` if it is a `SparseMatrixCSC` or `SparseVector`, otherwise convert to that type
in an efficient manner.
"""
function sparsecsc(@nospecialize A::AbstractArray)
    if iswrsparse(A)
        if depth(A) >= 1
            sparsecsc(inflate(A))
        else
            A
        end
    else
        Tv = eltype(A)
        invoke(SparseMatrixCSC{Tv,Int}, Tuple{AbstractMatrix}, A)
    end
end

sparsecsc(A::SparseVector) = A
sparsecsc(A::UpperTriangular{T,<:AbstractSparseMatrix}) where T = triu(A.data)
sparsecsc(A::LowerTriangular{T,<:AbstractSparseMatrix}) where T = tril(A.data)
function sparsecsc(A::Transpose{<:Any,<:AbstractSparseVector})
    B = parent(A);
    copy(reshape(transpose.(B), 1, size(B,1)))
end
function sparsecsc(A::Adjoint{<:Any,<:AbstractSparseVector})
    B = parent(A);
    copy(reshape(adjoint.(B), 1, size(B,1)))
end

sparsecsc(A::Transpose{<:Any,<:UpperTriangularPlain}) = _sparse(nzrangeup, transpose, A)
sparsecsc(A::Transpose{<:Any,<:LowerTriangularPlain}) = _sparse(nzrangelo, transpose, A)
sparsecsc(A::Adjoint{<:Any,<:UpperTriangularPlain}) = _sparse(nzrangeup, adjoint, A)
sparsecsc(A::Adjoint{<:Any,<:LowerTriangularPlain}) = _sparse(nzrangelo, adjoint, A)

function sparsecsc(A::AbstractTriangular{<:Any,<:SparseMatrixCSC})
    _sparse(nzrangeup, A)
end
function sparsecsc(A::Symmetric{<:Any,<:SparseMatrixCSC})
    _sparse(A.uplo == 'U' ? nzrangeup : nzrangelo, transpose, A)
end
function sparsecsc(A::Hermitian{<:Any,<:SparseMatrixCSC})
    _sparse(A.uplo == 'U' ? nzrangeup : nzrangelo, adjoint, A)
end
function sparsecsc(S::SubArray{<:Any,2,<:SparseMatrixCSC})
    getindex(S.parent,S.indices...)
end
function sparsecsc(S::Conjugate{<:Any,<:SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti}
    A = parent(S)
    SparseMatrixCSC{Tv,Ti}(A.m, A.n, A.colptr, A.rowval, conj.(A.nzval))
end

# 2 cases: Unit(Upper|Lower)Triangular{Tv,SparseMatrixCSC}
function _sparse(fnzrange::Function, A::AbstractTriangular{Tv}) where {Tv}
    S = A.data
    rowval = rowvals(S)
    nzval = nonzeros(S)
    m, n = size(S)
    Ti = eltype(rowval)
    unit = sA isa Union{UnitUpperTriangular,UnitLowerTriangular}
    newcolptr = Vector{Ti}(undef, n+1)
    newrowval = Vector{Ti}(undef, nnz(S))
    newnzval = Vector{Tv}(undef, nnz(S))
    newcolptr[1] = 1
    uplo = fnzrange == nzrangeup
    newk = 1
    @inbounds for j = 1:n
        newkk = newk
        if unit
            newk += !uplo
        end
        r = fnzrange(S, j); r1 = r.start; r2 = r.stop
        for k = r1:r2
            i = rowval[k]
            if i != j || i == j && !unit
                newrowval[newk] = i
                newnzval[newk] = nzval[k]
                newk += 1
            end
        end
        if unit
            uplo && (newkk = newk)
            newrowval[newkk] = j
            newnzval[newkk] = one(Tv)
            newk += uplo
        end
        newcolptr[j+1] = newk
    end
    nz = newcolptr[n+1] - 1
    resize!(newrowval, nz)
    resize!(newnzval, nz)
    SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval)
end

# 4 cases: (Symmetric|Hermitian) variants (:U|:L)
function _sparse(fnzrange::Function, fadj::Function, sA::SparseMatrixCSCSymmHerm{Tv}) where {Tv}
    A = sA.data
    rowval = rowvals(A)
    nzval = nonzeros(A)
    m, n = size(A)
    Ti = eltype(rowval)
    newcolptr = Vector{Ti}(undef, n+1)
    diagmap = fadj == transpose ? identity : real

    newcolptr[1] = 1
    colrange = fnzrange === nzrangeup ? (1:n) : (n:-1:1)
    @inbounds for j = colrange
        r = fnzrange(A, j); r1 = r.start; r2 = r.stop
        newcolptr[j+1] = r2 - r1 + 1
        for k = r1:r2
            row = rowval[k]
            if row != j
                newcolptr[row+1] += 1
            end
        end
    end
    cumsum!(newcolptr, newcolptr)
    nz = newcolptr[n+1] - 1
    newrowval = Vector{Ti}(undef, nz)
    newnzval = Vector{Tv}(undef, nz)
    @inbounds for j = 1:n
        newk = newcolptr[j]
        for k = fnzrange(A, j)
            i = rowval[k]
            nzv = nzval[k]
            if i != j
                newrowval[newk] = i
                newnzval[newk] = nzv
                newk += 1
                ni = newcolptr[i]
                newrowval[ni] = j
                newnzval[ni] = fadj(nzv)
                newcolptr[i] = ni + 1
            else
                newrowval[newk] = i
                newnzval[newk] = diagmap(nzv)
                newk += 1
            end
        end
        newcolptr[j] = newk
    end
    _sparse_gen(m, n, newcolptr, newrowval, newnzval)
end

# 8 cases: (Transpose|Adjoint){Tv,[Unit](Upper|Lower)Triangular}
function _sparse(fnzrange::Function, fadj::Function, taA::Union{Transpose{Tv,<:AbstractTriangular},Adjoint{Tv,<:AbstractTriangular}}) where {Tv}

    sA = taA.parent
    A = sA.data
    rowval = rowvals(A)
    nzval = nonzeros(A)
    m, n = size(A)
    Ti = eltype(rowval)
    newcolptr = Vector{Ti}(undef, n+1)
    unit = sA isa Union{UnitUpperTriangular,UnitLowerTriangular}
    uplo = A isa Union{UpperTriangular,UnitUpperTriangular}

    fill!(newcolptr, 1unit)
    newcolptr[1] = 1
    @inbounds for j = 1:n
        for k = fnzrange(A, j)
            i = rowval[k]
            if i != j || i == j && !unit
                newcolptr[i+1] += 1
            end
        end
    end
    cumsum!(newcolptr, newcolptr)
    nz = newcolptr[n+1] - 1
    newrowval = Vector{Ti}(undef, nz)
    newnzval = Vector{Tv}(undef, nz)

    @inbounds for j = 1:n
        if !uplo && unit
            ni = newcolptr[j]
            newrowval[ni] = j
            newnzval[ni] = fadj(one(Tv))
            newcolptr[j] = ni + 1
        end
        for k = fnzrange(A, j)
            i = rowval[k]
            nzv = nzval[k]
            if i != j || i == j && !unit
                ni = newcolptr[i]
                newrowval[ni] = j
                newnzval[ni] = fadj(nzv)
                newcolptr[i] = ni + 1
            end
        end
        if uplo && unit
            ni = newcolptr[j]
            newrowval[ni] = j
            newnzval[ni] = fadj(one(Tv))
            newcolptr[j] = ni + 1
        end
    end
    _sparse_gen(n, m, newcolptr, newrowval, newnzval)
end

# not used
# serves as a blueprint for accessing all nonzeros of an AbstractSparseMatrix
function sparsecopy(A::T) where T<:AbstractMatrix{Tv} where Tv
    iswrsparse(A) ? sparsecopy_sparse(A) : SparseMatrixCSC{Tv,indextype(A)}(A)
end

function sparsecopy_sparse(A::AbstractMatrix{Tv}) where {Tv}
    sA = sparseaccess(A)
    rowval = rowvals(sA)
    nzval = nonzeros(sA)
    m, n = size(A)
    nz = nnz_estimation(A)
    Ti = indextype(A)
    newcolptr = Vector{Ti}(undef, n+1)
    newrowval = Vector{Ti}(undef, nz)
    newnzval = Vector{Tv}(undef, nz)
    ip = 1
    for j = 1:n
        newcolptr[j] = ip
        if ip + m - 1 > nz
            nz += max(nz, m)
            resize!(newrowval, nz)
            resize!(newnzval, nz)
        end
        for k in nzrange(sA, j)
            newrowval[ip] = rowval[k]
            newnzval[ip] = nzval[k]
            ip += 1
        end
    end
    newcolptr[n+1] = ip
    ip -= 1
    resize!(newrowval, ip)
    resize!(newnzval, ip)
    SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval)
end

function _sparse_gen(m, n, newcolptr, newrowval, newnzval)
    @inbounds for j = n:-1:1
        newcolptr[j+1] = newcolptr[j]
    end
    newcolptr[1] = 1
    SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval)
end


text:
using DataStructures
using Base.Test

# construction

@test typeof(OrderedDict()) == OrderedDict{Any,Any}
@test typeof(OrderedDict('a',1)) == OrderedDict{Char,Int}
@test typeof(OrderedDict([("a",1),("b",2)])) == OrderedDict{ASCIIString,Int}

# empty dictionary
d = OrderedDict(Char, Int)
@test length(d) == 0
@test isempty(d)
@test_throws KeyError d['c'] == 1
d['c'] = 1
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification

for c in 'a':'z'
    d[c] = c-'a'+1
end

@test (d['a'] += 1) == 2
@test 'a' in keys(d)
@test haskey(d, 'a')
@test get(d, 'B', 0) == 0
@test !('B' in keys(d))
@test !haskey(d, 'B')
@test pop!(d, 'a') == 2

@test collect(keys(d)) == ['b':'z']
@test collect(values(d)) == [2:26]
@test collect(d) == [(a,i) for (a,i) in zip('b':'z', 2:26)]


text:
using JuMP
using Ipopt
using Hiop
using Printf
using DelimitedFiles
using acopf

export solve, model, initialPt_IPOPT, outputAll, computeAdmitances
function solve(opfmodel, opf_data)

#
# Initial point - needed especially for pegase cases
#
Pg0,Qg0,Vm0,Va0 = initialPt_IPOPT(opf_data)
opfmodel[:Pg] = Pg0
opfmodel[:Qg] = Qg0
opfmodel[:Vm] = Vm0
opfmodel[:Va] = Va0

optimize!(opfmodel)
status = termination_status(opfmodel)
if status != MOI.LOCALLY_SOLVED
    println("Could not solve the model to optimality.")
end
return opfmodel,status
end

function model(opf_data; max_iter=100, solver="Ipopt")
#shortcuts for compactness
lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

nbus  = length(buses); nline = length(lines); ngen  = length(generators)

#branch admitances
YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = acopf.computeAdmitances(lines, buses, baseMVA)

#
# JuMP model now
#
if solver == "Hiop"
  opfmodel = Model(optimizer_with_attributes(Hiop.Optimizer))
else
  opfmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter" => max_iter))
end

@variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax)
@variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax)

@variable(opfmodel, Va[1:nbus])
@variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus] <= buses[i].Vmax)
# fix the voltage angle at the reference bus
set_lower_bound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)
set_upper_bound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)

# minimize active power
coeff0 = Vector{Float64}(undef, ngen)
for (i,v) in enumerate(generators)
    coeff0[i] = v.coeff[v.n]
end
coeff1 = Vector{Float64}(undef, ngen)
for (i,v) in enumerate(generators)
    coeff1[i] = v.coeff[v.n - 1]
end
coeff2 = Vector{Float64}(undef, ngen)
for (i,v) in enumerate(generators)
    coeff2[i] = v.coeff[v.n - 2]
end

@NLobjective(opfmodel, Min, sum( coeff2[i]*(baseMVA*Pg[i])^2
                                + coeff1[i]*(baseMVA*Pg[i])
                                + coeff0[i] for i=1:ngen))

#
# power flow balance
#

for b in 1:nbus
    # real part
    @NLconstraint(
    opfmodel,
    ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b]^2
    + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *( YftR[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftI[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
    + sum( Vm[b]*Vm[busIdx[lines[l].from]]*( YtfR[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfI[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
    - ( sum(baseMVA*Pg[g] for g in BusGeners[b]) - buses[b].Pd ) / baseMVA      # Sbus part
    ==0)
    # imaginary part
end
for b in 1:nbus
    @NLconstraint(
    opfmodel,
    ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b]^2
    + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *(-YftI[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftR[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
    + sum( Vm[b]*Vm[busIdx[lines[l].from]]*(-YtfI[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfR[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
    - ( sum(baseMVA*Qg[g] for g in BusGeners[b]) - buses[b].Qd ) / baseMVA      #Sbus part
    ==0)
end
#
# branch/lines flow limits
#
nlinelim=0
for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
    nlinelim += 1
    flowmax=(lines[l].rateA/baseMVA)^2

    #branch apparent power limits (from bus)
    Yff_abs2=YffR[l]^2+YffI[l]^2;        Yft_abs2=YftR[l]^2+YftI[l]^2
    Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
    @NLconstraint(
        opfmodel,
        Vm[busIdx[lines[l].from]]^2 *
        ( Yff_abs2*Vm[busIdx[lines[l].from]]^2 + Yft_abs2*Vm[busIdx[lines[l].to]]^2
        + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]
            *(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
            -Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
            )
        )
        - flowmax <=0)

    #branch apparent power limits (to bus)
    Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
    Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
    @NLconstraint(
        opfmodel,
        Vm[busIdx[lines[l].to]]^2 *
        ( Ytf_abs2*Vm[busIdx[lines[l].from]]^2 + Ytt_abs2*Vm[busIdx[lines[l].to]]^2
        + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]
            *(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
            -Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
            )
        )
        - flowmax <=0)
    end
end

@printf("Buses: %d  Lines: %d  Generators: %d\n", nbus, nline, ngen)
println("Lines with limits  ", nlinelim)

return opfmodel, Pg, Qg, Va, Vm
end
# Compute initial point for IPOPT based on the values provided in the case data
function initialPt_IPOPT(opfdata)
Pg=zeros(length(opfdata.generators)); Qg=zeros(length(opfdata.generators)); i=1
for g in opfdata.generators
    # set the power levels in in between the bounds as suggested by matpower
    # (case data also contains initial values in .Pg and .Qg - not used with IPOPT)
    Pg[i]=0.5*(g.Pmax+g.Pmin)
    Qg[i]=0.5*(g.Qmax+g.Qmin)
    i=i+1
end
@assert i-1==length(opfdata.generators)

Vm=zeros(length(opfdata.buses)); i=1;
for b in opfdata.buses
    # set the ini val for voltage magnitude in between the bounds
    # (case data contains initials values in Vm - not used with IPOPT)
    Vm[i]=0.5*(b.Vmax+b.Vmin);
    i=i+1
end
@assert i-1==length(opfdata.buses)

# set all angles to the angle of the reference bus
Va = opfdata.buses[opfdata.bus_ref].Va * ones(length(opfdata.buses))

return Pg,Qg,Vm,Va
end

function outputAll(opfmodel, opf_data, Pg, Qg, Va, Vm)
  #shortcuts for compactness
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  # OUTPUTING
  println("Objective value: ", objective_value(opfmodel), "USD/hr")
  VM = Array{Float64}(undef, size(Vm, 1))
  VM=value.(Vm); VA=value.(Va)
  PG=value.(Pg); QG=value.(Qg)

  println("============================= BUSES ==================================")
  println("  BUS    Vm     Va   |   Pg (MW)    Qg(MVAr) ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")

  println("                     |     (generation)      ")
  println("----------------------------------------------------------------------")
  for i in 1:nbus
    @printf("%4d | %6.2f  %6.2f | %s  | \n",
	    buses[i].bus_i, VM[i], VA[i]*180/pi,
	    length(BusGeners[i])==0 ? "   --          --  " : @sprintf("%7.2f     %7.2f", baseMVA*PG[BusGeners[i][1]], baseMVA*QG[BusGeners[i][1]]))
  end
  println("\n")

  within=20 # percentage close to the limits

  nflowlim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nflowlim += 1
    end
  end

  if nflowlim>0
    println("Number of lines with flow limits: ", nflowlim)

    optvec=zeros(2*nbus+2*ngen)
    optvec[1:ngen]=PG
    optvec[ngen+1:2*ngen]=QG
    optvec[2*ngen+1:2*ngen+nbus]=VM
    optvec[2*ngen+nbus+1:2*ngen+2*nbus]=VA

    d = JuMP.NLPEvaluator(opfmodel)
    MOI.initialize(d, [:Jac])

    consRhs = zeros(2*nbus+2*nflowlim)
    MOI.eval_constraint(d, consRhs, optvec)

    #println(consRhs)

    @printf("================ Lines within %d %s of flow capacity ===================\n", within, "%")
    println("Line   From Bus    To Bus    At capacity")

    nlim=1
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        flowmax=(lines[l].rateA/baseMVA)^2
        idx = 2*nbus+nlim

        if( (consRhs[idx]+flowmax)  >= (1-within/100)^2*flowmax )
          @printf("%3d      %3d      %3d        %5.3f%s\n", l, lines[l].from, lines[l].to, 100*sqrt((consRhs[idx]+flowmax)/flowmax), "%" )
          # @printf("%7.4f   %7.4f    %7.4f \n", consRhs[idx], consRhs[idx]+flowmax,  flowmax)
        end
        nlim += 1
      end
    end
  end

  #println(getvalue(Vm))
  #println(getvalue(Va)*180/pi)

  #println(getvalue(Pg))
  #println(getvalue(Qg))

  return
end


text:
# #Metadata wrangling

# ## Get sample-specific metadata from airtable database

using Resonance

samplemeta = airtable_metadata() # having set ENV["AIRTABLE_KEY"]
# Now to add important metadata to samples

@chain samplemeta begin
    groupby([:subject, :timepoint])
    @transform!(:has_metabolomics = any(!ismissing, :Metabolomics_batch))
    sort!([:subject, :timepoint])
end

fmp_samples = DataFrame(XLSX.readtable("data/resonance_fmp/Fecal_All_033022.xlsx", "Sheet1", infer_eltypes=true)...)
rename!(fmp_samples, Dict(:studyID=>:subject, :collectionNum=> :timepoint))
@rsubset! fmp_samples begin
    :subject in samplemeta.subject
    !ismissing(:timepoint)
end

samplemeta = leftjoin(samplemeta, select(fmp_samples, [:subject, :timepoint, :collectionDate]), on=[:subject, :timepoint])
unique!(samplemeta)

##

# ## Metadata stored in filemaker pro database
#
# Read in the different tables as `DataFrame`s,
# then normalize certain columns.
# First, subject-specific data

fmp_subject = DataFrame(XLSX.readtable("data/resonance_fmp/Subject_Centric_060222.xlsx", "Sheet1", infer_eltypes=true)...)
rename!(fmp_subject, Dict(:studyID=>:subject))
@rsubset! fmp_subject :subject in samplemeta.subject

# Then, timepoint-specific data

fmp_timepoint = DataFrame(XLSX.readtable("data/resonance_fmp/Timepoint_Centric_033122.xlsx", "Sheet1", infer_eltypes=true)...)
rename!(fmp_timepoint, Dict(:studyID=>:subject))

# # and COVID-specific samples

fmp_covid = DataFrame(XLSX.readtable("data/resonance_fmp/COVID_Fecal_040122.xlsx", "Sheet1", infer_eltypes=true)...)
rename!(fmp_covid, Dict(:studyID=>:subject))
fmp_covid = leftjoin(fmp_covid, fmp_subject, on=:subject)

# ## Getting data joined together
@rsubset!(fmp_timepoint, !ismissing(:subject))
fmp_alltp = leftjoin(fmp_timepoint, fmp_subject, on=[:subject])

codebreastfeeding!(fmp_alltp)

count(!ismissing, fmp_alltp.bfcalculated)

# Add info about brain data (just if it's there)

brain = let
    fcleft = brain_ingest("data/brain/freesurfer_curvature_leftHemi_oct2021.csv"; label="curvature_leftHemi")
    fcright = brain_ingest("data/brain/freesurfer_curvature_rightHemi_oct2021.csv"; label="curvature_rightHemi")
    ftleft = brain_ingest("data/brain/freesurfer_thickness_leftHemi_oct2021.csv"; label="thickness_leftHemi")
    ftright = brain_ingest("data/brain/freesurfer_thickness_rightHemi_oct2021.csv"; label="thickness_rightHemi")
    seg = brain_ingest("data/brain/segmentationVolumeMeasurements_oct2021.csv"; label="segmentation")

    outerjoin(fcleft, fcright, ftleft, ftright, seg, on=[:subject, :timepoint], makeunique=true)
end
brain.has_segmentation .= true

## Validation

for row in eachrow(brain)
    all(h-> !ismissing(h) && h, row[r"has_"]) || @info row[r"has_"]
end

fmp_alltp = leftjoin(fmp_alltp, brain, on=[:subject, :timepoint])
fmp_alltp.has_segmentation = [ismissing(x) ? false : x for x in fmp_alltp.has_segmentation]

# transform!(fmp_alltp, :sid_old => ByRow(id-> ismissing(id) ? id : replace(id, r"_(\d+)F_"=>s"_\1E_")) => :sid_old_etoh)
# etoh_map = Dict((old=>new for (old, new) in zip(samplemeta.sid_old, samplemeta.sample)))
# transform!(fmp_alltp, :sid_old_etoh => ByRow(id-> (ismissing(id) || !haskey(etoh_map, id)) ? missing : etoh_map[id]) => :sample_etoh)

sort!(fmp_alltp, [:subject, :timepoint])

# ## Getting values for presence of data types
#
# We want to plot set intersections for having various kinds of data.
# In some cases, additional wrangling is necessary

fmp_alltp.race = map(fmp_alltp."Merge_Dem_Child_Race") do r
    ismissing(r) && return missing
    r == "Unknown" && return missing
    r == "Decline to Answer" && return missing
    contains(r, "\n") && return "Mixed"
    r ∈ ("Mixed", "Mixed Race") && return "Mixed"
    r ∈ ("Other Asian", "Asian ") && return "Asian"
    return r
end

# unique!(fmp_alltp, [:subject, :timepoint])

fmp_alltp.has_race = .!ismissing.(fmp_alltp."race")

##

subj = groupby(fmp_alltp, :subject)
transform!(subj, nrow => :n_samples; ungroup=false)

transform!(subj, AsTable(r"subject|breast"i) => (s -> begin
    if any(!ismissing, s.breastFedPercent)
        return fill(true, length(s[1]))
    else
        return fill(false, length(s[1]))
    end
    fill(0, length(s[1]))
end) => :has_bfperc_subj; ungroup=false)

transform!(subj, AsTable(r"subject|breast"i) => (s -> begin
    if any(!ismissing, s.breastFedPercent)
        return fill(true, length(s[1]))
    else
        return fill(false, length(s[1]))
    end
    fill(0, length(s[1]))
end) => :has_bfperc_subj; ungroup=false)

# fmp_alltp = DataFrames.transform(subj, AsTable([:timepoint, :sample]) => (s -> begin
#     has_stool = .!ismissing.(s.sample)
#     has_prevstool = fill(false, length(s[1]))
#     for i in 1:length(s[1])
#         i == 1 && continue
#         any(!ismissing, s.sample[1:i-1]) && (has_prevstool[i] = true)
#     end
#     (; has_stool, has_prevstool)
# end)=> [:has_stool, :has_prevstool])

fmp_alltp.has_everbreast = .!ismissing.(fmp_alltp."everBreastFed")
fmp_alltp.has_bfperc = .!ismissing.(fmp_alltp."breastFedPercent")

# ## Dealing with ages

fmp_alltp.ageMonths = map(eachrow(fmp_alltp)) do row
    if ismissing(row.scanAgeMonths)
        (row.assessmentAgeDays ./ 365 .* 12) .+ row.assessmentAgeMonths
    else
        (row.scanAgeDays ./ 365 .* 12) .+ row.scanAgeMonths
    end
end
count(row-> !ismissing(row."Blood Pressure::Diastolic") && !ismissing(row.ageMonths),
                eachrow(fmp_alltp))

fmp_alltp.age0to3mo   = map(a-> !ismissing(a) && a < 3,       fmp_alltp.ageMonths)
fmp_alltp.age3to6mo   = map(a-> !ismissing(a) && 3 <= a < 6,  fmp_alltp.ageMonths)
fmp_alltp.age6to12mo  = map(a-> !ismissing(a) && 6 <= a < 12, fmp_alltp.ageMonths)
fmp_alltp.age12moplus = map(a-> !ismissing(a) && 12 <= a,     fmp_alltp.ageMonths)

for n in names(fmp_alltp)
    fmp_alltp[!, n] = map(fmp_alltp[!, n]) do v
        v isa AbstractString ? replace(v, r"\n"=>"___") : v
    end
end

CSV.write("data/wrangled/timepoints.csv", fmp_alltp)
CSV.write("data/wrangled/omnisamples.csv", @rsubset(samplemeta, :Fecal_EtOH == "F"))
CSV.write("data/wrangled/etohsamples.csv", @rsubset(samplemeta, :Fecal_EtOH == "E"))
CSV.write("data/wrangled/covid.csv", fmp_covid)

##

using Airtable

base = AirBase("appSWOVVdqAi5aT5u")
tab = AirTable("Samples", base)
joinedsamples = leftjoin(select(samplemeta, [:airtable_id, :subject, :timepoint, :sample, :childAgeMonths]),
                         select(fmp_alltp, [:subject, :timepoint, :ageMonths, :has_segmentation]),
                         on=[:subject, :timepoint]
)

agepatches = AirRecord[]
brainpatches = AirRecord[]

for row in eachrow(joinedsamples)
    if !ismissing(row[:childAgeMonths]) && ismissing(row[:ageMonths])
        push!(agepatches, AirRecord(row[:airtable_id], tab, (; childAgeMonths=row[:ageMonths])))
    end
    if coalesce(row[:has_segmentation], false)
        push!(brainpatches, AirRecord(row[:airtable_id], tab, (; concurrent_scan=row[:has_segmentation])))
    end
end

!isempty(agepatches) && Airtable.patch!(tab, agepatches)

unique!(brainpatches)
!isempty(brainpatches) && Airtable.patch!(tab, brainpatches)

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
test = ["8A004A801A8002F478",
        "620080001611562C8802118E34",
        "C0015000016115A2E0802F182340",
        "A0016C880162017C3686B18A3D4780"]

@enum PacketType begin
    SUM          = 0
    PROD         = 1
    MINIMUM      = 2
    MAXIMUM      = 3
    LITERAL      = 4
    GREATER_THAN = 5
    LESS_THAN    = 6
    EQUAL        = 7
end

function parse_packet(packet::AbstractString)
    (packet, rest) = all(∈(['0', '1']), packet) ?
        parse_bin_packet(packet) :
        parse_hexa_packet(packet)
    @assert all(==('0'), rest)
    return packet
end

function parse_hexa_packet(hexa::AbstractString)
    chars = split(hexa, "")
    binary = join(string.(parse.(Int, chars; base = 16); base = 2, pad = 4))
    return parse_bin_packet(binary)
end

function parse_bin_packet(bin::AbstractString)
    version = unpack(bin[1:3])
    type = PacketType(unpack(bin[4:6]))
    (payload, rest) = parse_bin_packet(Val(type), bin[7:end])
    return ((; version, type = type, payload), rest)
end

function parse_bin_packet(::Val{LITERAL}, bin)
    stop = 5 * findfirst(==('0'), bin[1:5:end])
    bits = reshape(split(bin, "")[1:stop], 5, :)
    @assert all(==("1"), bits[1, 1:end - 1]) && bits[1, end] == "0"
    return (unpack(join(bits[2:end, 1:end])), bin[stop + 1:end])
end

function parse_bin_packet(::Val, bin)
    subs = []
    if bin[1] == '0'
        len = unpack(bin[2:16])
        rest = bin[17:17 + len - 1]
        while !isempty(rest)
            (sub, rest) = parse_bin_packet(rest)
            push!(subs, sub)
        end
        return (subs, bin[17 + len:end])
    else
        rest = bin[13:end]
        for i in 1:unpack(bin[2:12])
            (sub, rest) = parse_bin_packet(rest)
            push!(subs, sub)
        end
    end
    return (subs, rest)
end

unpack(bits::AbstractString) = parse(Int, bits; base = 2)

function q1(packet)
    return packet.version +
        (packet.type != LITERAL ? sum(q1, packet.payload) : 0)
end

evaluate(packet) = evaluate(Val(packet.type), packet.payload)
evaluate(::Val{LITERAL},      payload) = payload
evaluate(::Val{SUM},          payload) = sum(evaluate, payload)
evaluate(::Val{PROD},         payload) = prod(evaluate, payload)
evaluate(::Val{MINIMUM},      payload) = payload .|> evaluate |> minimum
evaluate(::Val{MAXIMUM},      payload) = payload .|> evaluate |> maximum
evaluate(::Val{LESS_THAN},    payload) = <(evaluate.(payload)...)
evaluate(::Val{GREATER_THAN}, payload) = >(evaluate.(payload)...)
evaluate(::Val{EQUAL},        payload) = ==(evaluate.(payload)...)

q2(packet) = evaluate(packet)

@assert parse_packet("D2FE28") ==
    (version = 6, type = LITERAL, payload = 2021)

@assert parse_packet("38006F45291200") ==
    (version = 1,
     type = LESS_THAN,
     payload = [(version = 6, type = LITERAL, payload = 10),
                (version = 2, type = LITERAL, payload = 20)])

@assert parse_packet("EE00D40C823060") ==
    (version = 7,
     type = MAXIMUM,
     payload = [(version = 2, type = LITERAL, payload = 1),
                (version = 4, type = LITERAL, payload = 2),
                (version = 1, type = LITERAL, payload = 3)])

@assert sum(q1, parse_packet.(test)) == 16 + 12 + 23 + 31

packet = parse_packet(strip(read("day16.in", String)))
println("Q1: ", q1(packet))
println("Q2: ", q2(packet))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
module Husk

include("HtmlAssembler.jl")
include("HtmlParser.jl")
include("HtmlError.jl")

end

max_stars_repo_name: haverzard/Husk

max_issues_repo_name: haverzard/Husk

max_forks_repo_name: haverzard/Husk


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ContainerStateTerminated <: SwaggerModel
    containerID::Any # spec type: Union{ Nothing, String } # spec name: containerID
    exitCode::Any # spec type: Union{ Nothing, Int32 } # spec name: exitCode
    finishedAt::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: finishedAt
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    signal::Any # spec type: Union{ Nothing, Int32 } # spec name: signal
    startedAt::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: startedAt

    function IoK8sApiCoreV1ContainerStateTerminated(;containerID=nothing, exitCode=nothing, finishedAt=nothing, message=nothing, reason=nothing, signal=nothing, startedAt=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("containerID"), containerID)
        setfield!(o, Symbol("containerID"), containerID)
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("exitCode"), exitCode)
        setfield!(o, Symbol("exitCode"), exitCode)
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("finishedAt"), finishedAt)
        setfield!(o, Symbol("finishedAt"), finishedAt)
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("signal"), signal)
        setfield!(o, Symbol("signal"), signal)
        validate_property(IoK8sApiCoreV1ContainerStateTerminated, Symbol("startedAt"), startedAt)
        setfield!(o, Symbol("startedAt"), startedAt)
        o
    end
end # type IoK8sApiCoreV1ContainerStateTerminated

const _property_map_IoK8sApiCoreV1ContainerStateTerminated = Dict{Symbol,Symbol}(Symbol("containerID")=>Symbol("containerID"), Symbol("exitCode")=>Symbol("exitCode"), Symbol("finishedAt")=>Symbol("finishedAt"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"), Symbol("signal")=>Symbol("signal"), Symbol("startedAt")=>Symbol("startedAt"))
const _property_types_IoK8sApiCoreV1ContainerStateTerminated = Dict{Symbol,String}(Symbol("containerID")=>"String", Symbol("exitCode")=>"Int32", Symbol("finishedAt")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("message")=>"String", Symbol("reason")=>"String", Symbol("signal")=>"Int32", Symbol("startedAt")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiCoreV1ContainerStateTerminated }) = collect(keys(_property_map_IoK8sApiCoreV1ContainerStateTerminated))
Swagger.property_type(::Type{ IoK8sApiCoreV1ContainerStateTerminated }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ContainerStateTerminated[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ContainerStateTerminated }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ContainerStateTerminated[property_name]

function check_required(o::IoK8sApiCoreV1ContainerStateTerminated)
    (getproperty(o, Symbol("exitCode")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ContainerStateTerminated }, name::Symbol, val)
end


text:
using Random

"""
Generate indices of train-test split given training indices
and percent training data
"""
function sample(inds, percent)
    N = length(inds)
    num = (Int64)(round(percent * N))
    #train = sort(rand(1:490, (num,)))
    ind_train = sort(randperm(N)[1:num])
    ind_test = copy(inds)
    deleteat!(ind_test, ind_train)
    #work with indices up until this point
    x_train = xs[ind_train]
    x_test = xs[ind_test]
    return ind_train, x_train, ind_test, x_test
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@testset "TAMU Parsing " begin
    files = readdir(TAMU_DIR)
    if length(files) == 0
        error("No test files in the folder")
    end

    sys = PSB.build_system(PSB.PSYTestSystems, "tamu_ACTIVSg2000_sys")
    @info "Successfully parsed $TAMU_DIR to System struct"

    # Test bad input
    @test_throws PowerSystems.DataFormatError TamuSystem("")
end


text:
module GriddedTests

include("gridded.jl")
include("mixed.jl")

end


text:
function selectionsort!(x::AbstractArray)
    n = length(x)
    for i ∈ 1:n
        # mi = minimum index
        mi = i
        for j ∈ (i+1):n
            @inbounds if x[j] ≤ x[mi]
                mi = j
            end
        end
        @inbounds tmp   = x[i]
        @inbounds x[i]  = x[mi]
        @inbounds x[mi] = tmp
    end
    return nothing
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct FullMachine <: Machine
        R::Float64
        R_f::Float64
        R_1d::Float64
        R_1q::Float64
        L_d::Float64
        L_q::Float64
        L_ad::Float64
        L_aq::Float64
        L_f1d::Float64
        L_ff::Float64
        L_1d::Float64
        L_1q::Float64
        ext::Dict{String, Any}
        inv_d_fluxlink::Array{Float64,2}
        inv_q_fluxlink::Array{Float64,2}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameter of a full order flux stator-rotor model without zero sequence flux in the stator.
 The derivative of stator fluxes (ψd and ψq) is NOT neglected. Only one q-axis damping circuit is considered. All parameters are in machine per unit.
 Refer to Chapter 3 of Power System Stability and Control by P. Kundur or Chapter 11 of Power System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby, for more details.
 Note that the models are somewhat different (but equivalent) due to the different Park Transformation used in both books.

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: `(0, nothing)`
- `R_f::Float64`: Field rotor winding resistance in per unit, validation range: `(0, nothing)`
- `R_1d::Float64`:  Damping rotor winding resistance on d-axis in per unit. This value is denoted as RD in Machowski., validation range: `(0, nothing)`
- `R_1q::Float64`: Damping rotor winding resistance on q-axis in per unit. This value is denoted as RQ in Machowski., validation range: `(0, nothing)`
- `L_d::Float64`: Inductance of fictitious damping that represent the effect of the three-phase stator winding in the d-axis of the rotor, in per unit. This value is denoted as L_ad + L_l in Kundur (and Ld in Machowski)., validation range: `(0, nothing)`
- `L_q::Float64`: Inductance of fictitious damping that represent the effect of the three-phase stator winding in the q-axis of the rotor, in per unit. This value is denoted as L_aq + L_l in Kundur., validation range: `(0, nothing)`
- `L_ad::Float64`: Mutual inductance between stator winding and rotor field (and damping) winding inductance on d-axis, in per unit, validation range: `(0, nothing)`
- `L_aq::Float64`: Mutual inductance between stator winding and rotor damping winding inductance on q-axis, in per unit, validation range: `(0, nothing)`
- `L_f1d::Float64`: Mutual inductance between rotor field winding and rotor damping winding inductance on d-axis, in per unit, validation range: `(0, nothing)`
- `L_ff::Float64`: Field rotor winding inductance, in per unit, validation range: `(0, nothing)`
- `L_1d::Float64`: Inductance of the d-axis rotor damping circuit, in per unit, validation range: `(0, nothing)`
- `L_1q::Float64`: Inductance of the q-axis rotor damping circuit, in per unit, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `inv_d_fluxlink::Array{Float64,2}`: Equations 3.127, 3.130, 3.131 From Kundur
- `inv_q_fluxlink::Array{Float64,2}`: Equations 3.128, 3.132 From Kundur
- `states::Vector{Symbol}`: The states are:
	ψd: d-axis stator flux,
	ψq: q-axis stator flux,
	ψf: field rotor flux,
	ψ1d: d-axis rotor damping flux,
	ψ1q: q-axis rotor damping flux
- `n_states::Int`: FullMachine has 5 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct FullMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Field rotor winding resistance in per unit"
    R_f::Float64
    " Damping rotor winding resistance on d-axis in per unit. This value is denoted as RD in Machowski."
    R_1d::Float64
    "Damping rotor winding resistance on q-axis in per unit. This value is denoted as RQ in Machowski."
    R_1q::Float64
    "Inductance of fictitious damping that represent the effect of the three-phase stator winding in the d-axis of the rotor, in per unit. This value is denoted as L_ad + L_l in Kundur (and Ld in Machowski)."
    L_d::Float64
    "Inductance of fictitious damping that represent the effect of the three-phase stator winding in the q-axis of the rotor, in per unit. This value is denoted as L_aq + L_l in Kundur."
    L_q::Float64
    "Mutual inductance between stator winding and rotor field (and damping) winding inductance on d-axis, in per unit"
    L_ad::Float64
    "Mutual inductance between stator winding and rotor damping winding inductance on q-axis, in per unit"
    L_aq::Float64
    "Mutual inductance between rotor field winding and rotor damping winding inductance on d-axis, in per unit"
    L_f1d::Float64
    "Field rotor winding inductance, in per unit"
    L_ff::Float64
    "Inductance of the d-axis rotor damping circuit, in per unit"
    L_1d::Float64
    "Inductance of the q-axis rotor damping circuit, in per unit"
    L_1q::Float64
    ext::Dict{String, Any}
    "Equations 3.127, 3.130, 3.131 From Kundur"
    inv_d_fluxlink::Array{Float64,2}
    "Equations 3.128, 3.132 From Kundur"
    inv_q_fluxlink::Array{Float64,2}
    "The states are:
	ψd: d-axis stator flux,
	ψq: q-axis stator flux,
	ψf: field rotor flux,
	ψ1d: d-axis rotor damping flux,
	ψ1q: q-axis rotor damping flux"
    states::Vector{Symbol}
    "FullMachine has 5 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function FullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, ext=Dict{String, Any}(), )
    FullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, ext, inv([[-L_d L_ad L_ad]; [-L_ad L_ff L_f1d]; [-L_ad L_f1d L_1d]]), inv([[-L_q L_aq]; [-L_aq L_1q]]), [:ψd, :ψq, :ψf, :ψ1d, :ψ1q], 5, InfrastructureSystemsInternal(), )
end

function FullMachine(; R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, ext=Dict{String, Any}(), inv_d_fluxlink=inv([[-L_d L_ad L_ad]; [-L_ad L_ff L_f1d]; [-L_ad L_f1d L_1d]]), inv_q_fluxlink=inv([[-L_q L_aq]; [-L_aq L_1q]]), states=[:ψd, :ψq, :ψf, :ψ1d, :ψ1q], n_states=5, internal=InfrastructureSystemsInternal(), )
    FullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, ext, inv_d_fluxlink, inv_q_fluxlink, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function FullMachine(::Nothing)
    FullMachine(;
        R=0,
        R_f=0,
        R_1d=0,
        R_1q=0,
        L_d=1,
        L_q=1,
        L_ad=2,
        L_aq=2,
        L_f1d=1,
        L_ff=2,
        L_1d=1,
        L_1q=1,
        ext=Dict{String, Any}(),
    )
end

"""Get [`FullMachine`](@ref) `R`."""
get_R(value::FullMachine) = value.R
"""Get [`FullMachine`](@ref) `R_f`."""
get_R_f(value::FullMachine) = value.R_f
"""Get [`FullMachine`](@ref) `R_1d`."""
get_R_1d(value::FullMachine) = value.R_1d
"""Get [`FullMachine`](@ref) `R_1q`."""
get_R_1q(value::FullMachine) = value.R_1q
"""Get [`FullMachine`](@ref) `L_d`."""
get_L_d(value::FullMachine) = value.L_d
"""Get [`FullMachine`](@ref) `L_q`."""
get_L_q(value::FullMachine) = value.L_q
"""Get [`FullMachine`](@ref) `L_ad`."""
get_L_ad(value::FullMachine) = value.L_ad
"""Get [`FullMachine`](@ref) `L_aq`."""
get_L_aq(value::FullMachine) = value.L_aq
"""Get [`FullMachine`](@ref) `L_f1d`."""
get_L_f1d(value::FullMachine) = value.L_f1d
"""Get [`FullMachine`](@ref) `L_ff`."""
get_L_ff(value::FullMachine) = value.L_ff
"""Get [`FullMachine`](@ref) `L_1d`."""
get_L_1d(value::FullMachine) = value.L_1d
"""Get [`FullMachine`](@ref) `L_1q`."""
get_L_1q(value::FullMachine) = value.L_1q
"""Get [`FullMachine`](@ref) `ext`."""
get_ext(value::FullMachine) = value.ext
"""Get [`FullMachine`](@ref) `inv_d_fluxlink`."""
get_inv_d_fluxlink(value::FullMachine) = value.inv_d_fluxlink
"""Get [`FullMachine`](@ref) `inv_q_fluxlink`."""
get_inv_q_fluxlink(value::FullMachine) = value.inv_q_fluxlink
"""Get [`FullMachine`](@ref) `states`."""
get_states(value::FullMachine) = value.states
"""Get [`FullMachine`](@ref) `n_states`."""
get_n_states(value::FullMachine) = value.n_states
"""Get [`FullMachine`](@ref) `internal`."""
get_internal(value::FullMachine) = value.internal

"""Set [`FullMachine`](@ref) `R`."""
set_R!(value::FullMachine, val) = value.R = val
"""Set [`FullMachine`](@ref) `R_f`."""
set_R_f!(value::FullMachine, val) = value.R_f = val
"""Set [`FullMachine`](@ref) `R_1d`."""
set_R_1d!(value::FullMachine, val) = value.R_1d = val
"""Set [`FullMachine`](@ref) `R_1q`."""
set_R_1q!(value::FullMachine, val) = value.R_1q = val
"""Set [`FullMachine`](@ref) `L_d`."""
set_L_d!(value::FullMachine, val) = value.L_d = val
"""Set [`FullMachine`](@ref) `L_q`."""
set_L_q!(value::FullMachine, val) = value.L_q = val
"""Set [`FullMachine`](@ref) `L_ad`."""
set_L_ad!(value::FullMachine, val) = value.L_ad = val
"""Set [`FullMachine`](@ref) `L_aq`."""
set_L_aq!(value::FullMachine, val) = value.L_aq = val
"""Set [`FullMachine`](@ref) `L_f1d`."""
set_L_f1d!(value::FullMachine, val) = value.L_f1d = val
"""Set [`FullMachine`](@ref) `L_ff`."""
set_L_ff!(value::FullMachine, val) = value.L_ff = val
"""Set [`FullMachine`](@ref) `L_1d`."""
set_L_1d!(value::FullMachine, val) = value.L_1d = val
"""Set [`FullMachine`](@ref) `L_1q`."""
set_L_1q!(value::FullMachine, val) = value.L_1q = val
"""Set [`FullMachine`](@ref) `ext`."""
set_ext!(value::FullMachine, val) = value.ext = val
"""Set [`FullMachine`](@ref) `inv_d_fluxlink`."""
set_inv_d_fluxlink!(value::FullMachine, val) = value.inv_d_fluxlink = val
"""Set [`FullMachine`](@ref) `inv_q_fluxlink`."""
set_inv_q_fluxlink!(value::FullMachine, val) = value.inv_q_fluxlink = val


text:
using Test

using AlgebraicPetri
using AlgebraicPetri.Epidemiology
using AlgebraicPetri.BilayerNetworks
using LabelledArrays

using Catlab
using Catlab.CategoricalAlgebra
using Catlab.Graphics
import Catlab.CategoricalAlgebra: migrate!
using Catlab.WiringDiagrams
using Catlab.Programs.RelationalPrograms
display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

# #### SIR Model:

# define model
sir = @relation (s,i,r) begin
    infection(s,i)
    recovery(i,r)
end
display_uwd(sir)
#-
#

# Extract the Petri network form of the model.

psir = apex(oapply_epi(sir))
psir
Graph(psir)

# Convert the Petri network into a bilayer network and draw it.
# This model uses a computation graph to express the computation of the vector field of the Petri net.

bnsir = BilayerNetwork()
lab_bnsir = LabelledBilayerNetwork()
migrate!(bnsir, psir)
migrate!(lab_bnsir, psir)
bnsir
@test typeof(to_graphviz(bnsir)) == Graph
@test typeof(to_graphviz(lab_bnsir)) == Graph

bnsir_test = @acset BilayerNetwork begin
    Qin = 3
    Qout = 3
    Win = 3
    Wa = 3
    Wn = 3
    Box = 2
    arg = [1,2,2]
    call = [1,1,2]
    efflux = [1,1,2]
    effusion = [1,2,2]
    influx = [1,1,2]
    infusion = [2,2,3]
end

@test bnsir == bnsir_test

######################
# Round Trip testing #
######################
function roundtrip(pn::AbstractPetriNet, bn::AbstractBilayerNetwork)
    roundtrippetri = PetriNet()
    migrate!(roundtrippetri, bn)
    pn_structure = PetriNet()
    copy_parts!(pn_structure, pn)
    return roundtrippetri, pn_structure
end

function roundtrip(pn::AbstractLabelledPetriNet, bn::AbstractLabelledBilayerNetwork)
    roundtrippetri = LabelledPetriNet()
    migrate!(roundtrippetri, bn)
    return roundtrippetri, pn
end

function test_roundtrip(pn::AbstractPetriNet, bn::AbstractBilayerNetwork)
    roundtrippetri, pn_structure = roundtrip(pn, bn)
    @test roundtrippetri == pn_structure
end

bnsir = BilayerNetwork()
migrate!(bnsir, psir)
test_roundtrip(psir, bnsir)

seir = @relation (s,e,i,r) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
end

pseir = apex(oapply_epi(seir))
bnseir = BilayerNetwork()
migrate!(bnseir, pseir)
test_roundtrip(pseir, bnseir)

lbnseir = LabelledBilayerNetwork()
migrate!(lbnseir, pseir)
test_roundtrip(pseir, lbnseir)

########################
# Compile and Evaluate #
########################

du = LVector(S=0.0,I=0.0,R=0.0)
u = LVector(S=10.0,I=1.0,R=0.0)
params = LVector(inf=0.1,rec=0.05)

bn_du =  AlgebraicPetri.BilayerNetworks.evaluate(lab_bnsir, u; inf=0.1,rec=0.05)
vectorfield(psir)(du, u, params, 0)

@test all(abs.(bn_du .- du[[:S,:I,:R]]) .< 1e-9)

eval(AlgebraicPetri.BilayerNetworks.compile(lab_bnsir, :du, :ϕ, :u; inf=0.1, rec=0.05))
comp_du = [0.0,0.0,0.0]
comp_ϕ = [0.0,0.0,0.0]

f!(comp_du, comp_ϕ, u, 0)

@test all(abs.(comp_du .- du[[:S,:I,:R]]) .< 1e-9)

p, params = params, nothing
eval(AlgebraicPetri.BilayerNetworks.compile(lab_bnsir, :du, :ϕ, :u, :p))
comp_du = [0.0,0.0,0.0]
comp_ϕ = [0.0,0.0,0.0]

f!(comp_du, comp_ϕ, u, p, 0)

@test all(abs.(comp_du .- du[[:S,:I,:R]]) .< 1e-9)

# Test for arbitrary variable names
eval(AlgebraicPetri.BilayerNetworks.compile(lab_bnsir, :a, :b, :c, :d))
comp_du = [0.0,0.0,0.0]
comp_ϕ = [0.0,0.0,0.0]

f!(comp_du, comp_ϕ, u, p, 0)

@test all(abs.(comp_du .- du[[:S,:I,:R]]) .< 1e-9)

##############
# Edge cases #
##############

function test_sir_equivalence(bn, pn)
  du = LVector(S=0.0,I=0.0,R=0.0)
  u = LVector(S=10.0,I=1.0,R=0.0)
  params = LVector(inf=0.1,rec=0.05)

  bn_du =  AlgebraicPetri.BilayerNetworks.evaluate(bn, u; inf=0.1,rec=0.05)
  vectorfield(pn)(du, u, params, 0)

  @test all(abs.(bn_du .- du[[:S,:I,:R]]) .< 1e-9)
end

# Ensure that the migration function properly balances the BLN
bnsir_edge = @acset LabelledBilayerNetwork begin
    Qin = 3
    Qout = 3
    Win = 3
    Wa = 2
    Wn = 2
    Box = 2
    arg = [1,2,2]
    call = [1,1,2]
    efflux = [1,2]
    effusion = [1,2]
    influx = [1,2]
    infusion = [2,3]
    parameter=[:inf, :rec]
    variable=[:S, :I, :R]
    tanvar=[:S, :I, :R]
end

test_sir_equivalence(bnsir_edge, psir)

edge_sir = LabelledPetriNet()
migrate!(edge_sir, bnsir_edge)
test_sir_equivalence(bnsir_edge, edge_sir)

bnsir_edge_add = @acset LabelledBilayerNetwork begin
    Qin = 3
    Qout = 3
    Win = 3
    Wa = 4
    Wn = 4
    Box = 2
    arg = [1,2,2]
    call = [1,1,2]
    efflux = [1,1,1,2]
    effusion = [1,2,2,2]
    influx = [1,1,1,2]
    infusion = [2,2,2,3]
    parameter=[:inf, :rec]
    variable=[:S, :I, :R]
    tanvar=[:S, :I, :R]
end

test_sir_equivalence(bnsir_edge_add, psir)

edge_sir = LabelledPetriNet()
migrate!(edge_sir, bnsir_edge_add)
test_sir_equivalence(bnsir_edge, edge_sir)

bnsir_migrate_break = @acset LabelledBilayerNetwork begin
    Qin = 3
    Qout = 3
    Win = 2
    Wa = 2
    Wn = 2
    Box = 2
    arg = [1,2]
    call = [1,1]
    efflux = [1,2]
    effusion = [1,2]
    influx = [1,2]
    infusion = [2,3]
    parameter=[:inf, :rec]
    variable=[:S, :I, :R]
    tanvar=[:S, :I, :R]
end
bn_error = "Mass action does not allow species $(bnsir_migrate_break[2, :variable]) to be "*
           "removed without contributing to the rate."
@test_throws ErrorException(bn_error) migrate!(edge_sir, bnsir_migrate_break)

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
using Hungarian
using Random
include("hungarian.jl")

function benching()
    seed = Random.seed!(7)

    for n in [800,1000,2000,4000,8000,16000,20000,30000,40000]
        Base.GC.gc()

        A = rand(seed, UInt16, n, n)
        A[A .> 50000] .= 50000
        B = copy(A)
        C = copy(A)
        println("n: ", n)
        println("Lib: ")
        @time matching_lib =  Hungarian.munkres!(A)
        matching_lib = [findfirst(matching_lib[:,i].==Hungarian.STAR) for i = 1:n]
        println("Own: ")
        @time matching = run_hungarian!(B)

        sum_lib = 0.0
        c = 1
        for r in matching_lib
            sum_lib += C[r,c]
            c += 1
        end
        sum_own = 0.0
        c = 1
        for r in matching
            sum_own += C[r,c]
            c += 1
        end
        println("Sum lib: ", sum_lib)
        println("Sum own: ", sum_own)
        @assert sum_lib == sum_own
    end
end

benching()

max_stars_repo_name: opensourcesblog/hungarian

max_issues_repo_name: opensourcesblog/hungarian

max_forks_repo_name: opensourcesblog/hungarian


text:
@testset "sender_connect.jl" begin
    @testset "SYN sent test" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_state(sender, JLSponge.SYN_SENT)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        @test sender.bytes_in_flight == 1
    end

    @testset "SYN acked test" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_state(sender, JLSponge.SYN_SENT)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        @test sender.bytes_in_flight == 1
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        @test sender.bytes_in_flight == 0
        expect_no_seg(sender)
    end

    @testset "SYN -> wrong ack test" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_state(sender, JLSponge.SYN_SENT)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        @test sender.bytes_in_flight == 1
        ack_received_test(sender, WrappingInt32(isn))
        expect_state(sender, JLSponge.SYN_SENT)
        @test sender.bytes_in_flight == 1
        expect_no_seg(sender)
    end

    @testset "SYN acked, data" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_state(sender, JLSponge.SYN_SENT)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        @test sender.bytes_in_flight == 1
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        @test sender.bytes_in_flight == 0
        expect_no_seg(sender)
        write_bytes!(sender, "abcdefgh")
        tick!(sender, 1)
        expect_state(sender, JLSponge.SYN_ACKED)
        expect_seg(sender, data="abcdefgh", seqno=WrappingInt32(isn + 1))
        @test sender.bytes_in_flight == 8
        ack_received_test(sender, WrappingInt32(isn + 9))
        expect_state(sender, JLSponge.SYN_ACKED)
        @test sender.bytes_in_flight == 0
        expect_no_seg(sender)
        @test next_seqno(sender) == WrappingInt32(isn + 9)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
"""
    State{T}

    state information in maximal coordinates for Body at time steps: 1, 2, 3.
    information at time step 3 is recovered using configurations at time step 2 and velocities at time step 2.5.

    x1: position at previous time step
    q1: orientation (Quaternion) at previous time step
    v15: linear velocity at time step 1.5 (midpoint)
    ϕ15: angular velocity at time step 1.5 (midpoint)

    x2: position at current time step
    q2: orientation (Quaternion) at current time step
    F2: linear impulse (force * time step) applied at current time step
    τ2: angular impulse (torque * timestep) applied at current time step

    vsol: linear velocity at time step 2.5 (midpoint); contains current value (index 1) and candidate value (index 2)
    ϕsol: angular velocity at time step 2.5 (midpoint); contains current value (index 1) and candidate value (index 2)

    d: implicit dynamics evaluator (zero vector indicates physics are satisfied)
    D: Jacobian of implicit dynamics
"""
mutable struct State{T}
    # previous state
    x1::SVector{3,T}
    q1::Quaternion{T}
    v15::SVector{3,T}
    ϕ15::SVector{3,T}

    # current state
    x2::SVector{3,T}
    q2::Quaternion{T}
    F2::SVector{3,T}
    τ2::SVector{3,T}

    # solution estimate [before step; after step]
    vsol::Vector{SVector{3,T}}
    ϕsol::Vector{SVector{3,T}}

    # dynamics and Jacobian evaluation
    d::SVector{6,T}
    D::SMatrix{6,6,T,36}

    function State{T}() where T
        x1 = szeros(T, 3)
        q1 = one(Quaternion{T})
        v15 = szeros(T, 3)
        ϕ15 = szeros(T, 3)

        x2 = szeros(T, 3)
        q2 = one(Quaternion{T})
        F2 = szeros(T, 3)
        τ2 = szeros(T, 3)

        vsol = [szeros(T, 3) for i=1:2]
        ϕsol = [szeros(T, 3) for i=1:2]

        d = szeros(T, 6)
        D = szeros(T, 6, 6)

        new{T}(x1, q1, v15, ϕ15, x2, q2, F2, τ2, vsol, ϕsol, d, D)
    end
end


text:
import SVMLightLoader: load_svmlight_file
import SoftConfidenceWeighted: init, fit!, predict

# The code is from
# http://thirld.com/blog/2015/05/30/julia-profiling-cheat-sheet/
function benchmark(func, args...; filename="trace.def")
    # Any setup code goes here.

    # Run once, to force compilation.
    println("======================= First run:")
    @time func(args...)

    # Run a second time, with profiling.
    println("======================= Second run:")

    Profile.init()
    Profile.clear()
    Profile.clear_malloc_data()

    @profile @time func(args...)

    # Write profile results to profile.bin.
    f = open(filename, "w")
    Profile.print(f)
    close(f)
end

function calc_accuracy(results, answers)
    n_correct_answers = 0
    for (result, answer) in zip(results, answers)
        if result == answer
            n_correct_answers += 1
        end
    end
    accuracy = n_correct_answers / length(y)
end

X, y = load_svmlight_file("data.txt")

model = init(10.0, 10.0)

benchmark(fit!, model, X, y, filename="trace.fit.txt")

model = init(10.0, 10.0)
model = fit!(model, X, y)
results = predict(model, X)

accuracy = calc_accuracy(results, y)
println("Accuracy: $accuracy")


text:
# https://www.hackerrank.com/challenges/30-review-loop/problem

function get_input(n::Int)
    s = String[]
    for _ in 1:n
        append!(s, readlines())
    end
    return s
end

function print_sol(sa::Array{String})
    for i in eachindex(sa)
        s1 = sa[i][1:2:end]
        s2 = sa[i][2:2:end]
        println("$s1 $s2")
    end
end

n = parse(Int, readline())
s = get_input(n)
print_sol(s)

max_stars_repo_path: Practice/Tutorials/30 Days of Code/Day 06 Let's Review.jl

max_stars_repo_name: mtaufanr/hackerrank-code-julia

max_issues_repo_name: mtaufanr/hackerrank-code-julia

max_forks_repo_name: mtaufanr/hackerrank-code-julia


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE

using MortarContact2DAD, Test
using MortarContact2DAD: project_from_master_to_slave_ad,
                         project_from_slave_to_master_ad
using ForwardDiff

function get_func_1(max_iterations)
    function func_1(u)
        x1 = [0.0, 0.0] + u[1:2]
        x2 = [0.0, 2.0] + u[3:4]
        x3 = [0.5, 0.5] + u[5:6]
        t1_ = (x2-x1) / norm(x2-x1)
        n1 = n2 = [0.0 1.0; -1.0 0.0] * t1_
        x1_ = DVTI(x1, x2)
        n1_ = DVTI(n1, n2)
        slave = Element(Seg2, [1, 2])
        xi1 = project_from_master_to_slave_ad(slave, x1_, n1_, x3;
                                              debug=true, max_iterations=max_iterations)
        return [xi1]
    end
    return func_1
end

J = ForwardDiff.jacobian(get_func_1(10), zeros(6))
@test isapprox(J, [-0.25 -0.75 0.25 -0.25 0.0 1.0])
@test_throws Exception ForwardDiff.jacobian(get_func_1(0), zeros(6))

function get_func_2(max_iterations)
    function func_2(u)
        x1 = [0.0, 0.0] + u[1:2]
        x2 = [0.0, 2.0] + u[3:4]
        x3 = [0.5, -1.0] + u[5:6]
        x4 = [0.5,  1.0] + u[7:8]
        t1 = (x2-x1) / norm(x2-x1)
        n1 = n2 = [0.0 1.0; -1.0 0.0] * t1
        x1_ = DVTI(x1, x2)
        n1_ = DVTI(n1, n2)
        x2_ = DVTI(x3, x4)
        master = Element(Seg2, [3, 4])
        xi2 = project_from_slave_to_master_ad(master, x1, n1, x2_;
                                              max_iterations=max_iterations)
        return [xi2]
    end
    return func_2
end

J = ForwardDiff.jacobian(get_func_2(10), zeros(8))
@test isapprox(J, [0.25  1.0  -0.25  0.0  0.0  -0.5  0.0  -0.5])
@test_throws Exception ForwardDiff.jacobian(get_func_2(0), zeros(8))


text:
immutable HLineGeometry <: Gadfly.GeometryElement
    color::Union(ColorValue, Nothing)
    size::Union(Measure, Nothing)

    function HLineGeometry(; color=nothing,
                           size::Union(Measure, Nothing)=nothing)
        new(color === nothing ? nothing : Color.color(color),
            size)
    end
end

const hline = HLineGeometry

function element_aesthetics(::HLineGeometry)
    [:yintercept]
end

# Generate a form for the hline geometry
function render(geom::HLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.hline", aes, :yintercept)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size

    return compose!(
        context(),
        Compose.line([[(0w, y), (1w, y)] for y in aes.yintercept]),
        stroke(color),
        linewidth(size),
        svgclass("xfixed"))
end


text:
"""
    struct ScaledMeasure{R,M,X} <: AbstractMeasure{X}
        logscale :: R
        base :: M
    end

"""
struct ScaledMeasure{R,M,X} <: AbstractMeasure{X}
    logscale :: R
    base :: M
end

function logdensity(sm::ScaledMeasure{R,M,X}, x::X) where {X, R <: Real, M <: AbstractMeasure{X}}
    logdensity(sm.base, x) + sm.logscale
end

function Base.:*(k::Number, m::AbstractMeasure{X}) where {X}
    if iszero(k)
        return TrivialMeasure(X)
    else
        return ScaledMeasure{typeof(k), typeof(m),X}(log(k),m)
    end
end

Base.:*(m::AbstractMeasure{X}, k::Real) where {X} = k * m

≪(::M, ::ScaledMeasure{R,M,X}) where {R,M,X} = true
≪(::ScaledMeasure{R,M,X}, ::M) where {R,M,X} = true


text:
## A small driver file to run the 4
# synthetic data tests in parallel
# data calibrated to the RossmanKaggleData
#Run like this
#julia -p 4 -L syntheticDataHarness.jl test_SyntheticData.jl numRuns outPathStub d usePoisson
#passed arguments
	#ARGS[1] numRuns
	#ARGs[2] partial path for output.
	#ARGS[3] d
	#ARGS[4] bool for using the poisson amount of data assumption.
	#ARGS[5] Optional: bool for only using ShrunkenPolicies
using Distributed

const numRuns = parse(Int, ARGS[1])
const spath = ARGS[2]
const param_path = "../RossmanKaggleData/CleanedData/"
const d = parse(Int, ARGS[3])
const usePoisson = parse(Bool, ARGS[4])
const onlyShrunken = length(ARGS) >= 5 ? parse(Bool, ARGS[5]) : false

const s = .95
K_grid = vcat(round.(Int, 2 .^(3:.25:10)), 1115)

#Used for most experiments.  This flag is a cludge since larger grid only used for some things
if !onlyShrunken
	N_grid = [10]
else
	N_grid = [10, 15, 20, 30, 40]
end

outPath = "$(spath)_syntheticRossman_$(s)_$(usePoisson)_$(4*numRuns)"
#old naming specification.  keep for a bit.
#outPath = "$(spath)_Ross_$(maximum(K_grid))_$(d)_$(N)_$(s)_$(usePoisson)_$(4*numRuns)"

#First read in the data and parse it appropriately
#Do this on a single processor bc it should be fast.
ps_full = readdlm("../RossmanKaggleData/CleanedData/ps_full$(d).csv", ',')
supp_full = readdlm("../RossmanKaggleData/CleanedData/support$(d).csv", ',')

start_time = time_ns()
file_a = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_a_", N_grid, s, seed=8675309, usePoisson=usePoisson, onlyShrunken=onlyShrunken)
file_b = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_b_", N_grid, s, seed=5167462266, usePoisson=usePoisson, onlyShrunken=onlyShrunken)
file_c = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_c_", N_grid, s, seed=5164174290, usePoisson=usePoisson, onlyShrunken=onlyShrunken)
file_d = @spawn convInKtest(numRuns, K_grid, supp_full, ps_full, "$(outPath)_d_", N_grid, s, seed=112456059, usePoisson=usePoisson, onlyShrunken=onlyShrunken)

# ######
file_a = fetch(file_a)
file_b = fetch(file_b)
file_c = fetch(file_c)
file_d = fetch(file_d)

time_stamp = (time_ns() - start_time) * 1e-9

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)
data_t[:, 1] .+= numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] .+= 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] .+= 3numRuns
data = vcat(data, data_t)

f = open("$(outPath).csv", "w")
writedlm(f, header, ',')
writedlm(f, data, ',')
close(f)

println("Num Paths: \t $(4*numRuns) \t Time:", time_stamp )


text:
# Simulation Name
simname = basename(pwd())

# NOTE: write to scratchdir
function outdir_suffix(scale, kmcmc, seed_data)
  return "seed_$(seed_data)/scale_$(scale)/Kmcmc_$(kmcmc)"
end

# NOTE: modify
settings = [Dict(:simname => simname,
                 :repfam_dist_scale => scale,
                 :N => [300, 300],
                 :Z_idx => 4,  # J = 20
                 :thin_samps => 2,
                 :nburn => 10000,
                 :nsamps => 2000,
                 # :nburn => 10, # NOTE: test
                 # :nsamps => 20, # NOTE: test
                 :Lmcmc => Dict(0 => 2, 1 => 2),
                 :Kmcmc => kmcmc,
                 :seed_data => seed_data,
                 :seed_mcmc => 0,
                 :outdir_suffix => outdir_suffix(scale, kmcmc, seed_data))
            for scale in [0, 10],
            kmcmc in [3, 4, 5, 6, 7, 15],
            seed_data in 1:5]

settings = vec(settings)