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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ fb195cf0-8062-11eb-3c99-fde4a4241d8b
begin
	using JLD2
	using LinearAlgebra
	using StatsBase, Statistics
	using Plots, JSServe, PlutoUI

	import WGLMakie, PyPlot
end

# ╔═╡ 949bde72-81d3-11eb-3e45-dd042b57b2bb
using LaTeXStrings

# ╔═╡ 0de243c4-8063-11eb-3580-136f64776460
Page()

# ╔═╡ 19d832ec-8063-11eb-39cf-b3f130287cc4
function makie(f::Function)
    scene = WGLMakie.Scene(resolution = (600, 400))
	f(scene)

    return scene
end

# ╔═╡ 9a89abe2-8062-11eb-2614-f7cb7a8fed05
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ c96f10f0-8062-11eb-3e4e-2fb11d63ba86
DataIO = ingredients("../src/io.jl").DataIO

# ╔═╡ dedf4344-8137-11eb-3dd6-33edcac405d9
DataFilter = ingredients("../src/filter.jl").DataFilter

# ╔═╡ 0c8e9d0a-81cc-11eb-304f-8d4b03ed97b1
PointCloud = ingredients("../src/geo.jl").PointCloud

# ╔═╡ d3ca4d32-811e-11eb-2993-8bc55dd24808
FIGS = "/home/nolln/root/data/seqspace/eyal/figs"

# ╔═╡ 93bf354e-8124-11eb-3b5a-414a3603360f
SAVE = false

# ╔═╡ 6aa479a6-813a-11eb-0c97-353a3f53b9b2
gr(html_output_format=:png)

# ╔═╡ 1d7393d0-811f-11eb-349a-bba38050379e
SAMPLE = "sample_4"

# ╔═╡ 6fda1e64-8138-11eb-13ed-977e26b0e0cc
DATA_DIR = "/home/nolln/data/star/$SAMPLE"

# ╔═╡ ec2e7784-8062-11eb-16b4-3b009b7c4f60
RAW_INPUT = "$DATA_DIR/Solo.out/Gene/filtered/matrix.mtx"

# ╔═╡ 82119d96-8138-11eb-0138-1b6e3b959ae2
JLD_INPUT = "$DATA_DIR/glm_fit.jld2"

# ╔═╡ de1e1208-8062-11eb-0ee0-a9db551ffcd0
scrna = open(RAW_INPUT) do io
	DataIO.read_mtx(io)
end;

# ╔═╡ 84185814-806d-11eb-28f7-fd158f7638bb
μ = vec(mean(scrna,dims=2)); n = vec(sum(scrna.>0,dims=2))

# ╔═╡ d208a304-806d-11eb-15d0-ab7366d88511
begin
	PyPlot.clf()
	PyPlot.hexbin(log10.(μ), log10.(n.+1),
		gridsize=(30,30),
		cmap="inferno"
	)
	PyPlot.xlabel("log avg count/cell")
	PyPlot.ylabel("log number of expressing cells")
	if SAVE
		PyPlot.savefig("$FIGS/qc/$SAMPLE/count_vs_expressing_cells.png")
	end
	PyPlot.gcf()
end

# ╔═╡ 33a625ba-806f-11eb-1a68-bfadc0cc660b
begin
	h = histogram(log10.(μ),bins=100,legend=false)
	vline!([-3])
	xlabel!("avg log count/per cell")
	if SAVE
		savefig("$FIGS/qc/$SAMPLE/count_per_cell_histogram.png")
	end
	h
end

# ╔═╡ 597a17d0-8063-11eb-3738-3793619ecd8c
data = DataFilter.genes(scrna);

# ╔═╡ 86f746e0-8074-11eb-23b3-972e81ea370a
md"""
## Data filtering
We need to scale the data so so that highly expressed genes and lowly expressed genes are equally weighted - i.e. focus on the variation between them not overall scale.
This is clear after looking at the principal components: dominated by inter-gene variation.

Define ``x_{in}`` to be the expression count for gene i and cell n

###### Idea 1:
Naively, we just want to fit each gene to a given parameterized distribution independently.
Genes can then be re-expressed
However, we want to regress against our degrees of freedom, the total expression for each cell ``g_n \equiv \sum_i x_{in}``, to deal with systematic variations in sequencing depth.
We assume that the expression of gene i takes the form:

`` x_{in} = \text{exp}\left[\alpha_i + \beta_i \text{log}(g_n)\right] + \epsilon_{in} ``

Where ``\epsilon_{in}`` is a noise parameter with zero mean. Then, given a distribution ``\rho`` parameterized by it's mean we could then write the log likelihood as

`` \ell(\vec{\beta}_i) = \sum_n \rho(x_{in}|\mu_i = \text{exp}\left[\alpha_{i} + \beta_{i} \text{log}(g_n)\right]) ``

Thus for example, if we think the data is Poisson distributed, the above relation reduces to

`` \ell(\alpha_i, \beta_i) = \sum_n x_{in}\left(\alpha_i + \beta_i \text{log}(g_n)\right) - \text{exp}\left[\alpha_i + \beta_i \text{log}(g_n) \right] ``

Conversely, if we want to model an overdispersed counting distribution, we can turn to the negative binomial (with mean``=\mu`` and variance``=\mu + \mu^2 / \Theta``):

`` \ell(\alpha_i, \beta_i, \gamma_i) = \sum_n \text{log}(\Gamma(x_{in} +  \gamma_i)) - \text{log}(\Gamma(x_{in})) - \text{log}(\Gamma( \gamma_i)) +  \gamma_i\text{log}\left(\frac{\mu}{\mu +  \gamma_i}\right) + x_{in}\text{log}\left(\frac{\gamma_i}{\mu + \gamma_i}\right) ``

This procedure introduces ``3`` parameters per gene and thus will have to fit ``3n_g`` number of parameters in total.
Hafemeister, C. & Satija, R. (2019) claims this overfits and thus put an heurestic prior on the parameters.
Let's test this
"""

# ╔═╡ 91c38718-8138-11eb-0f18-85916a6aa8ca
@load JLD_INPUT fits

# ╔═╡ 24263ca8-80fe-11eb-0692-07698ae10d02

# ╔═╡ 21c52cf2-813b-11eb-1453-557d5b272ca1
M = vec(maximum(data,dims=2))

# ╔═╡ ba8c6388-8116-11eb-35d2-df884093be9f
begin
	E = map(fits) do r
		r.likelihood
	end
	α = map(fits) do r
		r.parameters[1]
	end
	β = map(fits) do r
		r.parameters[2]
	end
	γ = map(fits) do r
		r.parameters[3]
	end

	δα = map(fits) do r
		sqrt(abs(r.uncertainty[1]))
	end
	δβ = map(fits) do r
		sqrt(abs(r.uncertainty[2]))
	end
	δγ = map(fits) do r
		sqrt(abs(r.uncertainty[3]))
	end
	χ = vec(mean(data,dims=2))[1:length(E)]
end

# ╔═╡ d7070fe2-8116-11eb-35c0-496f1cad3008
begin
	l1 = scatter(χ,E,
		yscale=:log,
		xscale=:log,
		alpha=0.1,
		legend=false,
		xlabel="average count/cell",
		ylabel="average negative log likelihood /cell"
	)
	if SAVE
		savefig("$FIGS/qc/$SAMPLE/count_per_cell_vs_likelihood.png")
	end
	l1
end

# ╔═╡ 2f377f1c-8119-11eb-10a8-71fc4394927c
begin
	l2 = scatter(χ, α, marker_z=log10.(δα),
		xscale=:log10,
		alpha=0.1,
		label="",
		xlabel="average count/cell",
		ylabel="estimated α",
		colorbar_title="log uncertainty"
	)
	if SAVE
		savefig("$FIGS/qc/$SAMPLE/count_per_cell_vs_α.png")
	end
	l2
end

# ╔═╡ f9824c96-8119-11eb-3adb-c3cbc2986947
begin
	l3 = scatter(χ, β,
		marker_z=log10.(δβ),
		alpha=0.1,
		xscale=:log10,
		label="",
		ylims=(-5,5),
		xlabel="average count/cell",
		ylabel="estimated β",
		colorbar_title="log uncertainty"
	)
	if SAVE
		savefig("$FIGS/qc/$SAMPLE/count_per_cell_vs_β.png")
	end
	l3
end

# ╔═╡ 08a8018e-811a-11eb-0ce4-6711063fce7d
begin
	l4 = scatter(χ, γ,
		marker_z=log10.(δγ),
		label="",
		alpha=0.1,
		yscale=:log10,
		xscale=:log10,
		xlabel="average count/cell",
		ylabel="estimated γ",
		colorbar_title="log uncertainty"
	)
	if SAVE
		savefig("$FIGS/qc/$SAMPLE/count_per_cell_vs_γ.png")
	end
	l4
end

# ╔═╡ c81b6f9a-813a-11eb-00d1-1d92878cb882
begin
	PyPlot.clf()
	PyPlot.hexbin(log10.(χ), log10.(γ),
		gridsize=(40,40),
		cmap="inferno"
	)
	PyPlot.xlabel("log avg count/cell")
	PyPlot.ylabel("log estimated γ")
	if SAVE
		PyPlot.savefig("$FIGS/qc/$SAMPLE/count_per_cell_vs_γ_hexbin.png")
	end
	PyPlot.gcf()
end

# ╔═╡ 0e619b6e-813b-11eb-188c-b7ad7b1c6e3a
begin
	l5 = scatter(χ, M,
		marker_z=log10.(γ),
		label="",
		alpha=0.1,
		yscale=:log10,
		xscale=:log10,
		xlabel="average count/cell",
		ylabel="max count/cell",
		colorbar_title="log uncertainty"
	)
	if SAVE
		savefig("$FIGS/qc/$SAMPLE/max_count_per_cell_vs_γ.png")
	end
	l5
end

# ╔═╡ dfe3851e-8064-11eb-049a-4dfed6e9ff3c
function null(data)
	ñ = zeros(size(data))
	for i in 1:size(data,1)
		sample!(data[i,:], view(ñ,i,:))
	end
	return ñ
end

# ╔═╡ 2d358d88-814b-11eb-0c11-23089a86b856
norm = vcat((fit.residuals' for (i,fit) ∈ enumerate(fits) if M[i] > 1 )...)

# ╔═╡ 18a61724-8156-11eb-07db-3de6b3104a1c
function gene_names()
	names = open("$DATA_DIR/Solo.out/Gene/filtered/features.tsv") do io
		[split(line)[2] for line in eachline(io)]
	end
	# NOTE: many details here depend strongly on how we filter...
	#		should just make scrna matrix into a struct so names are always bound
	@assert length(names) == size(scrna,1)

	μ 	  = vec(mean(scrna,dims=2))
	names = names[μ .> 1e-3]
	names = names[M .> 1]

	@assert length(names) == size(norm,1)

	return names
end

# ╔═╡ 8d4fa794-814b-11eb-2a06-1ba01e048084
F = svd(norm);

# ╔═╡ ded15dc2-814b-11eb-35d5-63bf4d616d34
F̃ = svd(null(norm));

# ╔═╡ d7a04a06-814b-11eb-0f00-49c8745fc905
begin
	p1 = plot(F.S/sum(F.S), xscale=:log10, yscale=:log10)
	plot!(F̃.S/sum(F̃.S), xscale=:log10, yscale=:log10)
	xlabel!("principal component")
	ylabel!("singular value")
	if SAVE
		savefig("$FIGS/geo/$SAMPLE/residual_matrix_singular_values.png")
	end
	p1
end

# ╔═╡ 8b400948-814c-11eb-059d-f1fa06164139
ipr(ψ) = sum(ψ.^4, dims=1) ./ (sum(ψ.^2, dims=1)).^2

# ╔═╡ cb23cd18-814c-11eb-3314-674bcab45578
ψ = ipr(F.U)

# ╔═╡ 075494be-81d6-11eb-123e-d3f005c9cb0d
localized_state, localized_size = getindex(argmin(1 ./ ψ),2), minimum(1 ./ ψ)

# ╔═╡ d50cb378-814c-11eb-1863-db97f0895d64
begin
	p2 = plot(1 ./ ψ', yscale=:log10, xscale=:log10)
	xlabel!("component")
	ylabel!("participation ratio")
	if SAVE
		savefig("$FIGS/geo/$SAMPLE/residual_matrix_participation_ratio.png")
	end
	p2
end

# ╔═╡ eb4592e0-8156-11eb-2564-49ede281d83a
names = gene_names()

# ╔═╡ 36a6f78a-81dc-11eb-187f-e91cdaa1f611
localized_genes = names[sortperm(F.U[:,localized_state].^2; rev=true)]

# ╔═╡ fff6ffa0-8156-11eb-0a6b-195479c7b49b
names[argmax(abs.(F.U[:,22]))]

# ╔═╡ 0a30bc3e-814d-11eb-3a4c-77164eb08043
begin
	plot(abs.(F.U[:,1]),yscale=:log10,alpha=.5)
	plot!(abs.(F.U[:,20]),alpha=.5)
end

# ╔═╡ 22a82bd4-81cb-11eb-0c85-af09449b8b9e
dₑ = 50; Ñ = F.U[:,1:dₑ] * Diagonal(F.S[1:dₑ]) * F.Vt[1:dₑ,:]

# ╔═╡ acc97a86-81df-11eb-39bd-6b7c3e599e1c
dd = 500; approx = F.U[:,1:dd] * Diagonal(F.S[1:dd]) * F.Vt[1:dd,:]

# ╔═╡ 6e7c5dc4-81de-11eb-1e29-6d94b51d76b8
let p
	p = plot(sqrt.(cumsum(F.S.^2)/sum(F.S.^2)), linewidth=2)
	xlabel!("rank approximation")
	ylabel!("correlation with original")

	if SAVE
		savefig("$FIGS/geo/$SAMPLE/svd_reduction_accuracy.png")
	end
	p
end

# ╔═╡ 40885c20-81dd-11eb-1b6d-733ca67aa312
ind_sample = sample(CartesianIndices(norm), 2000)

# ╔═╡ 574a4dd8-81dd-11eb-16a9-d9066f374f9d
scatter(approx[ind_sample], norm[ind_sample])

# ╔═╡ 0823bfd8-81cd-11eb-144d-9bb76ef0ea65
g̃ = PointCloud.geodesics(Ñ,12)

# ╔═╡ 69e3f8f8-81f7-11eb-3cc7-a7ac4ca669a5
if SAVE
	@save "$DATA_DIR/geodesic_k12.jld2" Ñ g̃
end

# ╔═╡ 4d0509e8-81cf-11eb-1e7b-8d8ad2ecf672
PyPlot.clf(); PyPlot.matshow(g̃); PyPlot.gcf()

# ╔═╡ ff7617a2-81cf-11eb-1c8f-c5dbacc6d897
ρ, Rs = PointCloud.scaling(g̃, 1000)

# ╔═╡ dce339b0-81d2-11eb-1f94-6906505d3346
begin
	p3 = plot(Rs, ρ', alpha=0.03, color=:cyan3, label="")
	plot!(Rs, mean(ρ, dims=1)', linewidth=5, alpha=1, color=:cyan3, label="mean")
	plot!(Rs, 1e-3*Rs.^3, linestyle=:dashdot, linewidth=2, alpha=1, color=:coral3, label=L"\alpha R^3")
	plot!(Rs, 5e-2*Rs.^2, linestyle=:dashdot, linewidth=2, alpha=1, color=:darkgoldenrod1, label=L"\alpha R^2")

	xaxis!("radius", :log10, (10, 400))
	yaxis!("number of points", :log10, (10, 3200))

	if SAVE
		savefig("$FIGS/geo/$SAMPLE/ball_scaling.png")
	end
	p3
end

# ╔═╡ 4e9e490e-81d4-11eb-2948-6f1b312907fc
ξ = PointCloud.mds(g̃.^2, 10)

# ╔═╡ 708ef216-81d4-11eb-3673-e394d0fa118f
begin

residual = zeros(size(ξ,2))
for d ∈ 1:size(ξ,2)
	D̄ = .√PointCloud.distance²(ξ[:,1:d]')

	residual[d] = cor(
		PointCloud.upper_tri(D̄),
		PointCloud.upper_tri(g̃)
	)
end

end;

# ╔═╡ 97052e74-81d4-11eb-115f-0b6b19bea33a
begin
	p4 = plot(residual, linewidth=2, label="")
	xlabel!("dimension")
	ylabel!("correlation of distance matrices")
	if SAVE
		savefig("$FIGS/geo/$SAMPLE/isomap_distance_correlation.png")
	end
	p4
end

# ╔═╡ eb79c294-81d4-11eb-2cd6-633742a3ca0e
gene_index = Dict(g => i for (i,g) in enumerate(names))

# ╔═╡ 91e387a0-81da-11eb-1c9f-857e7079fe59
@bind THETA PlutoUI.Slider(range(0,180,length=100))

# ╔═╡ 9bb40ce6-81da-11eb-3814-670fc5029879
@bind PHI PlutoUI.Slider(range(0,360,length=100))

# ╔═╡ 449c5422-81e9-11eb-293d-2b3f8db9e449
GENE = "SOX11"

# ╔═╡ a9989792-81d9-11eb-2123-632f4c14bcd0
begin
	PyPlot.clf()
	PyPlot.scatter3D(ξ[:,1], ξ[:,2], ξ[:,3],
		c=vec(Ñ[gene_index[GENE],:]),
		cmap="inferno",
	)
	PyPlot.view_init(THETA,PHI)
	if false
		PyPlot.savefig("$FIGS/geo/$SAMPLE/isomap_coordinates_$GENE.png")
	end
	PyPlot.gcf()
end

# ╔═╡ 736505e4-81d6-11eb-0ed5-db36e1f6a6ec
#=WGLMakie.scatter(ξ[:,1], ξ[:,2], ξ[:,3],
	color=vec(Ñ[gene_index["PAX6"],:]),
	markersize=5
)=#

# ╔═╡ Cell order:
# ╠═fb195cf0-8062-11eb-3c99-fde4a4241d8b
# ╠═949bde72-81d3-11eb-3e45-dd042b57b2bb
# ╠═0de243c4-8063-11eb-3580-136f64776460
# ╠═19d832ec-8063-11eb-39cf-b3f130287cc4
# ╠═9a89abe2-8062-11eb-2614-f7cb7a8fed05
# ╠═c96f10f0-8062-11eb-3e4e-2fb11d63ba86
# ╠═dedf4344-8137-11eb-3dd6-33edcac405d9
# ╠═0c8e9d0a-81cc-11eb-304f-8d4b03ed97b1
# ╠═d3ca4d32-811e-11eb-2993-8bc55dd24808
# ╠═93bf354e-8124-11eb-3b5a-414a3603360f
# ╠═6aa479a6-813a-11eb-0c97-353a3f53b9b2
# ╠═1d7393d0-811f-11eb-349a-bba38050379e
# ╟─6fda1e64-8138-11eb-13ed-977e26b0e0cc
# ╟─ec2e7784-8062-11eb-16b4-3b009b7c4f60
# ╟─82119d96-8138-11eb-0138-1b6e3b959ae2
# ╟─de1e1208-8062-11eb-0ee0-a9db551ffcd0
# ╠═84185814-806d-11eb-28f7-fd158f7638bb
# ╟─d208a304-806d-11eb-15d0-ab7366d88511
# ╟─33a625ba-806f-11eb-1a68-bfadc0cc660b
# ╠═597a17d0-8063-11eb-3738-3793619ecd8c
# ╠═18a61724-8156-11eb-07db-3de6b3104a1c
# ╟─86f746e0-8074-11eb-23b3-972e81ea370a
# ╠═91c38718-8138-11eb-0f18-85916a6aa8ca
# ╟─24263ca8-80fe-11eb-0692-07698ae10d02
# ╟─21c52cf2-813b-11eb-1453-557d5b272ca1
# ╟─ba8c6388-8116-11eb-35d2-df884093be9f
# ╟─d7070fe2-8116-11eb-35c0-496f1cad3008
# ╟─2f377f1c-8119-11eb-10a8-71fc4394927c
# ╟─f9824c96-8119-11eb-3adb-c3cbc2986947
# ╟─08a8018e-811a-11eb-0ce4-6711063fce7d
# ╟─c81b6f9a-813a-11eb-00d1-1d92878cb882
# ╟─0e619b6e-813b-11eb-188c-b7ad7b1c6e3a
# ╠═dfe3851e-8064-11eb-049a-4dfed6e9ff3c
# ╠═2d358d88-814b-11eb-0c11-23089a86b856
# ╠═8d4fa794-814b-11eb-2a06-1ba01e048084
# ╠═ded15dc2-814b-11eb-35d5-63bf4d616d34
# ╠═d7a04a06-814b-11eb-0f00-49c8745fc905
# ╠═8b400948-814c-11eb-059d-f1fa06164139
# ╠═cb23cd18-814c-11eb-3314-674bcab45578
# ╠═075494be-81d6-11eb-123e-d3f005c9cb0d
# ╠═36a6f78a-81dc-11eb-187f-e91cdaa1f611
# ╠═d50cb378-814c-11eb-1863-db97f0895d64
# ╟─eb4592e0-8156-11eb-2564-49ede281d83a
# ╠═fff6ffa0-8156-11eb-0a6b-195479c7b49b
# ╠═0a30bc3e-814d-11eb-3a4c-77164eb08043
# ╠═22a82bd4-81cb-11eb-0c85-af09449b8b9e
# ╠═acc97a86-81df-11eb-39bd-6b7c3e599e1c
# ╠═6e7c5dc4-81de-11eb-1e29-6d94b51d76b8
# ╠═40885c20-81dd-11eb-1b6d-733ca67aa312
# ╠═574a4dd8-81dd-11eb-16a9-d9066f374f9d
# ╠═0823bfd8-81cd-11eb-144d-9bb76ef0ea65
# ╠═69e3f8f8-81f7-11eb-3cc7-a7ac4ca669a5
# ╠═4d0509e8-81cf-11eb-1e7b-8d8ad2ecf672
# ╠═ff7617a2-81cf-11eb-1c8f-c5dbacc6d897
# ╠═dce339b0-81d2-11eb-1f94-6906505d3346
# ╠═4e9e490e-81d4-11eb-2948-6f1b312907fc
# ╠═708ef216-81d4-11eb-3673-e394d0fa118f
# ╠═97052e74-81d4-11eb-115f-0b6b19bea33a
# ╠═eb79c294-81d4-11eb-2cd6-633742a3ca0e
# ╠═91e387a0-81da-11eb-1c9f-857e7079fe59
# ╠═9bb40ce6-81da-11eb-3814-670fc5029879
# ╠═449c5422-81e9-11eb-293d-2b3f8db9e449
# ╠═a9989792-81d9-11eb-2123-632f4c14bcd0
# ╠═736505e4-81d6-11eb-0ed5-db36e1f6a6ec

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
abstract type TwuAlphaModel <: AlphaModel end

struct TwuAlphaParam <: EoSParam
    M::SingleParam{Float64}
    N::SingleParam{Float64}
    L::SingleParam{Float64}
end

@newmodelsimple TwuAlpha TwuAlphaModel TwuAlphaParam

export TwuAlpha
function TwuAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv","alpha/Twu/Twu_like.csv"]; userlocations=userlocations, verbose=verbose)
    M = params["M"]
    N = params["N"]
    L = params["L"]
    packagedparams = TwuAlphaParam(M,N,L)
    model = TwuAlpha(packagedparams, verbose=verbose)
    return model
end

function α_function(model::CubicModel,V,T,z,alpha_model::TwuAlphaModel)
    Tc = model.params.Tc.values
    M  = alpha_model.params.M.values
    N  = alpha_model.params.N.values
    L  = alpha_model.params.L.values
    Tr = T ./ Tc
    α = @. Tr^(N*(M-1))*exp(L*(1-Tr^(N*M)))
    return α
end

ext: jl


text:
using Taro
using DataFrames
using Requests
using TextAnalysis
using Languages
#

# initialize JVM
try
  # fix compatibility issue
  if is_apple()
      # dont initialize coz apple initializes jvm by default
      # JVM.init()
      Taro.init() # init once
  else
      println("initializing JVM and Taro...")
      Taro.init() # init once
  end

catch
  # something is not right
end

#function to clean article
function preCleaning(pdfText)
  sd = StringDocument(pdfText) # convert to SD
  remove_corrupt_utf8!(sd) # remove currupt utf8
  remove_html_tags!(sd) # remove html
  prepare!(sd, strip_punctuation) # remove punctuation
  #prepare!(sd,strip_punctuation | strip_numbers | strip_case | strip_whitespace)
  return(sd)
end

# function to fetch the entire full articles from https://www.ncbi.nlm.nih.gov/pmc
function fetchFullArticleFromPmc(dataFrame,pubYear=2018)
  pmcids = dataFrame[:pmcid]
  tempDir = mktempdir() # create temp directory
  nrows =nrow(dataFrame)
#  tracker = Progress(nrows, 1)
  # define the new columns
  dataFrame[:pmcUrl] = ""
  dataFrame[:fullText] = ""

  for row in eachrow(dataFrame)
    try
      pmcid = row[:pmcid]
      pmcUrl = "https://www.ncbi.nlm.nih.gov/pmc/articles/$pmcid/pdf";
      println("downloading: $pmcid")
      pdfTempPath = Requests.save(Requests.get(pmcUrl), joinpath(tempDir,"$pmcid.pdf"))
      pdfMeta, pdfText = Taro.extract(pdfTempPath)

      # append row
      row[:pmcUrl] = pmcUrl
      sd = preCleaning(pdfText) # clean article
      pubYear =string(pubYear)
      timestamp!(sd, pubYear) # add pub year

      row[:fullText] =  TextAnalysis.text(sd)
    catch err
      # catch errors
      #warn(err)
    end
      #next!(tracker)
  end
  return(dataFrame)
end

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
module RayTracing
    export Vec
    export Sphere, Rectangle, Sketch
    export ConstantTexture, ImageTexture, Lambertian, Metal, Dielectric, DiffuseLight
    export Object, Scene, Camera
    export x, y, z, addObject, render

    using StaticArrays
    using LinearAlgebra
    using Images
    using Random
    using Printf

    include("vector.jl")
    include("ray.jl")
    include("surface.jl")
    include("rectangle.jl")
    include("sphere.jl")
    include("sketch.jl")
    include("hitpoint.jl")
    include("texture.jl")
    include("constant_texture.jl")
    include("image_texture.jl")
    include("material.jl")
    include("lambertian.jl")
    include("metal.jl")
    include("dielectric.jl")
    include("diffuse_light.jl")
    include("object.jl")
    include("scene.jl")
    include("camera.jl")

    function color(scene::Scene, ray::Ray, bg_color::Union{AbstractVec, Function}, depth::Int = 0)::AbstractVec
        local succeeded::Bool
        succeeded, hp::Union{Nothing, HitPoint}, material::Union{Nothing, Material} = hit(scene, ray, (min = 0.001, max = Inf))
        if succeeded
            emitted::AbstractVec = emit(material, hp.u, hp.v, hp.normal)
            if depth < 50
                succeeded, attenuation::Union{Nothing, AbstractVec}, scattered::Union{Nothing, Ray} = scatter(material, ray, hp)
                if succeeded
                    return emitted + attenuation .* color(scene, scattered, bg_color, depth + 1)
                else
                    return emitted
                end
            end
        end

        if bg_color isa AbstractVec
            bg_color
        else
            bg_color(normalize(ray.direction))
        end
    end

    function render(scene::Scene, camera::Camera, size::NamedTuple{(:width, :height), Tuple{T, T}} where T <: Integer, samples::Integer, bg_color::Union{AbstractVec, Function})::Array{RGB, 2}
        image::Array{RGB, 2} = zeros(RGB, size.height, size.width)

        function render_row(j::Integer)::Nothing
            for i = 1:size.width
                col::AbstractVec = Vec(0.0, 0.0, 0.0)
                for s = 1:samples
                    u::AbstractFloat = (i - rand()) / size.width
                    v::AbstractFloat = (j - rand()) / size.height
                    ray::Ray = get_ray(camera, u, v)
                    col += color(scene, ray, bg_color)
                end

                col /= samples
                image[size.height + 1 - j, i] = RGB(sqrt(x(col)), sqrt(y(col)), sqrt(z(col)))
            end
        end

        rows::Vector{Integer} = collect(1:size.height)
        shuffle!(rows)

        m::Threads.Mutex = Threads.Mutex()
        count::Integer = 0

        Threads.@threads for row in rows
            render_row(row)
            lock(m)
                count += 1
                @printf("\r%.0f%%", 100.0 * count / size.height)
            unlock(m)
        end

        image
    end

end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
@testset "Quadratic equation" begin
    Base.info("Testing quadratic equation...")
    solver = Solver()
    @test isapprox(solver(1.0,-4.0,3.0), [1.0,3.0])
    @test isapprox(solver(1.0,-2.0,1.0), ones(2))
    @test isapprox(solver(1.0,0.0,1.0), [-im,im])
    @test_throws ErrorException solver(0.0, 1.0, 2.0)
    Base.info("Done testing quadratic equation.")
end


text:
using Kafka

# create KafkaClient using single bootstrap broker
kc = KafkaClient("127.0.0.1", 9092)

# get supported API versions
vers = api_versions(kc)

# get metadata about a topic(s)
md_channel = metadata(kc, ["test"])
md = take!(md_channel)
# if you prefer synchronous logic, use one-linear
take!(metadata(kc, ["test"]))
take!(metadata(kc))

# get earliest and latest available offsets for topic "test" and partition 0
take!(earliest_offset(kc, "test", 0))
take!(latest_offset(kc, "test", 0))

# produce new messages
keys = [convert(Vector{UInt8}, key) for key in ["1", "2", "3"]]
values = [convert(Vector{UInt8}, value) for value in ["feel", "good", "inc."]]
messages = collect(zip(keys, values))
offset = take!(produce(kc, "test", 0, messages))

# fetch messages
start_offset = 0
offset_messages = take!(fetch(kc, "test", 0, start_offset))


text:
import jcontrol: FULLTRACK

import JSONTables: objecttable, jsontable
using DataFrames

dd = Dict(
    "X" => FULLTRACK.X,
    "Y" => FULLTRACK.Y,
    "curvature" => FULLTRACK.curvature,   # κ value at each waypoint
    "S" => FULLTRACK.S,
    "θ" => FULLTRACK.θ,
    "width" => FULLTRACK.width.(FULLTRACK.S)
)

df = DataFrame(dd)
# @info df dd
open("./track.json", "w") do f
    write(f, objecttable(df))
end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
using BinDeps

@BinDeps.setup

libglut = library_dependency("libglut", aliases = ["libglut", "freeglut"])

@windows_only begin
	using WinRPM
	provides(WinRPM.RPM, "freeglut", libglut, os = :Windows)
end

@linux_only begin
    provides(AptGet, {"freeglut3-dev" => libglut})
    provides(Yum, {"freeglut-devel" => libglut})
end

@BinDeps.install [:libglut => :libglut]


text:
module ApiImpl

using Swagger

include("api/Kubernetes.jl")
using .Kubernetes
import .Kubernetes: getAPIVersions

import Base: convert, get, put!, delete!, show
import Swagger: SwaggerModel

include("typealiases.jl")
include("apialiases.jl")

end # module


text:
##############################################################################
##
## setindex!()
##
##############################################################################

# Will automatically add a new column if needed

# Modify Dataset
function insert_single_column!(ds::Dataset, v::AbstractVector, col_ind::ColumnIndex)
    if ncol(ds) != 0 && nrow(ds) != length(v)
        throw(ArgumentError("New columns must have the same length as old columns"))
    end
    # if isa(v, AbstractRange)
    #     dv = collect(allowmissing(v))
    # elseif isa(v, BitVector)
    #     dv = convert(Vector{Union{Bool, Missing}}, v)
    # elseif
    #
    # dv = isa(v, AbstractRange) ?  : (isa(v, BitVector) ? convert(Vector{Union{Bool, Missing}}, v) : allowmissing(v))
    dv = _preprocess_column(v, nrow(ds), false)
    firstindex(dv) != 1 && _onebased_check_error()

    if haskey(index(ds), col_ind)
        j = index(ds)[col_ind]
        # if the modified column is with in sorting columns empty the gattributes
        # TODO we can be more clever about this
        if j ∈ index(ds).sortedcols
          _reset_grouping_info!(index(ds))
        end
        # change the type of column j
        _columns(ds)[j] = dv
        removeformat!(ds, j)
        _modified(_attributes(ds))
    else
        if col_ind isa SymbolOrString
            push!(index(ds), Symbol(col_ind))
            push!(_columns(ds), dv)
            _modified(_attributes(ds))
        else
            throw(ArgumentError("Cannot assign to non-existent column: $col_ind"))
        end
    end
    return ds
end

# Modify Dataset
function insert_single_entry!(ds::Dataset, v::Any, row_ind::Integer, col_ind::ColumnIndex)
    if haskey(index(ds), col_ind)
      # single entry doesn't remove format
        colidx = index(ds)[col_ind]
        _columns(ds)[colidx][row_ind] = v
        # if the modified column is with in sorting columns empty the gattributes
       if colidx ∈ index(ds).sortedcols
          _reset_grouping_info!(index(ds))
        end
        _modified(_attributes(ds))
        return v
    else
        throw(ArgumentError("Cannot assign to non-existent column: $col_ind"))
    end
end

# ds[!, SingleColumnIndex] = AbstractVector

# Modify Dataset
function Base.setindex!(ds::Dataset, v::AbstractVector, ::typeof(!), col_ind::ColumnIndex)
    insert_single_column!(ds, v, col_ind)
    return ds
end

function Base.setindex!(ds::Dataset, v::DatasetColumn, ::typeof(!), col_ind::ColumnIndex)
    insert_single_column!(ds, __!(v), col_ind)
    return ds
end

# ds.col = AbstractVector
# separate methods are needed due to dispatch ambiguity

# Modify Dataset
# v must be promotable to ds[!, col_ind], thus, we can keep the format but sorting based on it should be revised
function Base.setproperty!(ds::Dataset, col_ind::Symbol, v::AbstractVector)
    insert_single_column!(ds, v, col_ind)
    return ds
end

# Modify Dataset
function Base.setproperty!(ds::Dataset, col_ind::AbstractString, v::AbstractVector)
    insert_single_column!(ds, v, col_ind)
    return ds
end

function Base.setproperty!(ds::Dataset, col_ind::Symbol, v::DatasetColumn)
    insert_single_column!(ds, __!(v), col_ind)
    return ds
end

function Base.setproperty!(ds::Dataset, col_ind::AbstractString, v::DatasetColumn)
    insert_single_column!(ds, __!(v), col_ind)
    return ds
end

# Modify Dataset
Base.setproperty!(::Dataset, col_ind::Symbol, v::Any) =
    throw(ArgumentError("It is only allowed to pass a vector as a column of a Dataset. " *
                        "Instead use `ds[!, col_ind] .= v` if you want to use broadcasting."))

# Modify Dataset
Base.setproperty!(::Dataset, col_ind::AbstractString, v::Any) =
    throw(ArgumentError("It is only allowed to pass a vector as a column of a Dataset. " *
                        "Instead use `ds[!, col_ind] .= v` if you want to use broadcasting."))

# ds[SingleRowIndex, SingleColumnIndex] = Single Item

# Modify Dataset
function Base.setindex!(ds::Dataset, v::Any, row_ind::Integer, col_ind::ColumnIndex)
    insert_single_entry!(ds, v, row_ind, col_ind)
    return ds
end

# ds[SingleRowIndex, MultiColumnIndex] = value
# the method for value of type DatasetRow, AbstractDict and NamedTuple
# is defined in datasetrow.jl

# Modify Dataset
for T in MULTICOLUMNINDEX_TUPLE
    @eval function Base.setindex!(ds::Dataset,
                                  v::Union{Tuple, AbstractArray},
                                  row_ind::Integer,
                                  col_inds::$T)
        idxs = index(ds)[col_inds]
        if length(v) != length(idxs)
            throw(DimensionMismatch("$(length(idxs)) columns were selected but the assigned " *
                                    "collection contains $(length(v)) elements"))
        end
        for (i, x) in zip(idxs, v)
            ds[row_ind, i] = x
        end
        return ds
    end
end

# ds[MultiRowIndex, SingleColumnIndex] = AbstractVector

# Modify Dataset
for T in (:AbstractVector, :Not, :Colon)
    @eval function Base.setindex!(ds::Dataset,
                                  v::AbstractVector,
                                  row_inds::$T,
                                  col_ind::ColumnIndex)
        if row_inds isa Colon && !haskey(index(ds), col_ind)
            ds[!, col_ind] = copy(v)
            return ds
        end
        colidx = index(ds)[col_ind]
        _columns(ds)[colidx][row_inds] = v
        if col_ind ∈ index(ds).sortedcols
            _reset_grouping_info!(ds)
        end
        _modified(_attributes(ds))
        return ds
    end
end

# ds[MultiRowIndex, MultiColumnIndex] = AbstractDataset

# Modify Dataset
for T1 in (:AbstractVector, :Not, :Colon),
    T2 in MULTICOLUMNINDEX_TUPLE
    @eval function Base.setindex!(ds::Dataset,
                                  new_ds::AbstractDataset,
                                  row_inds::$T1,
                                  col_inds::$T2)
        idxs = index(ds)[col_inds]
        if view(_names(ds), idxs) != _names(new_ds)
            throw(ArgumentError("column names in source and target do not match"))
        end
        for (j, col) in enumerate(idxs)
            ds[row_inds, col] = new_ds[!, j]
        end
        return ds
    end
end

# Modify Dataset
for T in MULTICOLUMNINDEX_TUPLE
    @eval function Base.setindex!(ds::Dataset,
                                  new_ds::AbstractDataset,
                                  row_inds::typeof(!),
                                  col_inds::$T)
        idxs = index(ds)[col_inds]
        if view(_names(ds), idxs) != _names(new_ds)
            throw(ArgumentError("Column names in source and target data sets do not match"))
        end
        for (j, col) in enumerate(idxs)
            # make sure we make a copy on assignment
            ds[!, col] = new_ds[:, j]
        end
        return ds
    end
end

# ds[MultiRowIndex, MultiColumnIndex] = AbstractMatrix

# Modify Dataset
for T1 in (:AbstractVector, :Not, :Colon, :(typeof(!))),
    T2 in MULTICOLUMNINDEX_TUPLE
    @eval function Base.setindex!(ds::Dataset,
                                  mx::AbstractMatrix,
                                  row_inds::$T1,
                                  col_inds::$T2)
        idxs = index(ds)[col_inds]
        if size(mx, 2) != length(idxs)
            throw(DimensionMismatch("number of selected columns ($(length(idxs))) " *
                                    "and number of columns in " *
                                    "matrix ($(size(mx, 2))) do not match"))
        end
        for (j, col) in enumerate(idxs)
            ds[row_inds, col] = (row_inds === !) ? mx[:, j] : view(mx, :, j)
        end
        return ds
    end
end

##############################################################################
##
## Mutating methods
##
##############################################################################

"""
    insertcols!(ds::Dataset[, col], (name=>val)::Pair...;
                makeunique::Bool=false, copycols::Bool=true)

Insert a column into a data set in place. Return the updated `Dataset`.
If `col` is omitted it is set to `ncol(ds)+1`
(the column is inserted as the last column).

# Arguments
- `ds` : the Dataset to which we want to add columns
- `col` : a position at which we want to insert a column, passed as an integer
  or a column name (a string or a `Symbol`); the column selected with `col`
  and columns following it are shifted to the right in `ds` after the operation
- `name` : the name of the new column
- `val` : an `AbstractVector` giving the contents of the new column or a value of any
  type other than `AbstractArray` which will be repeated to fill a new vector;
  As a particular rule a values stored in a `Ref` or a `0`-dimensional `AbstractArray`
  are unwrapped and treated in the same way.
- `makeunique` : Defines what to do if `name` already exists in `ds`;
  if it is `false` an error will be thrown; if it is `true` a new unique name will
  be generated by adding a suffix
- `copycols` : whether vectors passed as columns should be copied

If `val` is an `AbstractRange` then the result of `collect(val)` is inserted.

# Examples
```jldoctest
julia> ds = Dataset(a=1:3)
3×1 Dataset
 Row │ a
     │ Int64
─────┼───────
   1 │     1
   2 │     2
   3 │     3

julia> insertcols!(ds, 1, :b => 'a':'c')
3×2 Dataset
 Row │ b     a
     │ Char  Int64
─────┼─────────────
   1 │ a         1
   2 │ b         2
   3 │ c         3

julia> insertcols!(ds, 2, :c => 2:4, :c => 3:5, makeunique=true)
3×4 Dataset
 Row │ b     c      c_1    a
     │ Char  Int64  Int64  Int64
─────┼───────────────────────────
   1 │ a         2      3      1
   2 │ b         3      4      2
   3 │ c         4      5      3
```
"""
function insertcols!(ds::Dataset, col::ColumnIndex, name_cols::Pair{Symbol, <:Any}...;
                     makeunique::Bool=false, unsupported_copy_cols = true)
    copycols = unsupported_copy_cols
# Modify Dataset
    col_ind = Int(col isa SymbolOrString ? columnindex(ds, col) : col)
    if !(0 < col_ind <= ncol(ds) + 1)
        throw(ArgumentError("attempt to insert a column to a data set with " *
                            "$(ncol(ds)) columns at index $col_ind"))
    end

    if !makeunique
        if !allunique(first.(name_cols))
            throw(ArgumentError("Names of columns to be inserted into a data set " *
                                "must be unique when `makeunique=true`"))
        end
        for (n, _) in name_cols
            if hasproperty(ds, n)
                throw(ArgumentError("Column $n is already present in the data set " *
                                    "which is not allowed when `makeunique=true`"))
            end
        end
    end

    if ncol(ds) == 0
        target_row_count = -1
    else
        target_row_count = nrow(ds)
    end

    for (n, v) in name_cols
        if v isa AbstractVector
            if target_row_count == -1
                target_row_count = length(v)
            elseif length(v) != target_row_count
                if target_row_count == nrow(ds)
                    throw(DimensionMismatch("length of new column $n which is " *
                                            "$(length(v)) must match the number " *
                                            "of rows in data set ($(nrow(ds)))"))
                else
                    throw(DimensionMismatch("all vectors passed to be inserted into " *
                                            "a data set must have the same length"))
                end
            end
        elseif v isa AbstractArray && ndims(v) > 1
            throw(ArgumentError("adding AbstractArray other than AbstractVector as " *
                                "a column of a data set is not allowed"))
        end
    end
    if target_row_count == -1
        target_row_count = 1
    end

    for (name, item) in name_cols
        if !(item isa AbstractVector)
            if item isa Union{AbstractArray{<:Any, 0}, Ref}
                x = item[]
                item_new = fill!(allocatecol(Union{typeof(x), Missing}, target_row_count), x)
            else
                @assert !(item isa AbstractArray)
                item_new = fill!(allocatecol(Union{typeof(item), Missing}, target_row_count), item)
            end
        elseif item isa AbstractRange
            item_new = allowmissing(collect(item))
        elseif copycols
            item_new = copy(allowmissing(item))
        else
            item_new = allowmissing(item)
        end

        firstindex(item_new) != 1 && _onebased_check_error()

        if ncol(ds) == 0
            ds[!, name] = item_new
        else
            if hasproperty(ds, name)
                @assert makeunique
                k = 1
                while true
                    nn = Symbol("$(name)_$k")
                    if !hasproperty(ds, nn)
                        name = nn
                        break
                    end
                    k += 1
                end
            end
            # insert! modifies index, thus it should modifies gattributes
            insert!(index(ds), col_ind, name)
            item_new = _preprocess_column(item_new, nrow(ds), false)
            insert!(_columns(ds), col_ind, item_new)
            _modified(_attributes(ds))
        end
        col_ind += 1
    end
    return ds
end

# Modify Dataset
insertcols!(ds::Dataset, col::ColumnIndex, name_cols::Pair{<:AbstractString, <:Any}...;
                     makeunique::Bool=false) =
    insertcols!(ds, col, (Symbol(n) => v for (n, v) in name_cols)...,
                makeunique=makeunique)

# Modify Dataset
insertcols!(ds::Dataset, name_cols::Pair{Symbol, <:Any}...;
            makeunique::Bool=false) =
    insertcols!(ds, ncol(ds)+1, name_cols..., makeunique=makeunique)

# Modify Dataset
insertcols!(ds::Dataset, name_cols::Pair{<:AbstractString, <:Any}...;
            makeunique::Bool=false) =
    insertcols!(ds, (Symbol(n) => v for (n, v) in name_cols)...,
                makeunique=makeunique)

# Modify Dataset
function insertcols!(ds::Dataset, col::Int=ncol(ds)+1; makeunique::Bool=false, name_cols...)
    if !(0 < col <= ncol(ds) + 1)
        throw(ArgumentError("attempt to insert a column to a data set with " *
                            "$(ncol(ds)) columns at index $col"))
    end
    if !isempty(name_cols)
        # an explicit error is thrown as keyword argument was supported in the past
        throw(ArgumentError("inserting colums using a keyword argument is not supported, " *
                            "pass a Pair as a positional argument instead"))
    end
    return ds
end

###########################
# DatasetColumn
##############################

Base.getindex(col1::DatasetColumn{Dataset,<:AbstractVector}, i::Integer) = getindex(__!(col1), i)
Base.getindex(col1::DatasetColumn{Dataset,<:AbstractVector}, i::CartesianIndex) = getindex(__!(col1), i)
Base.getindex(col1::DatasetColumn{Dataset,<:AbstractVector}, row_inds::AbstractVector) = getindex(__!(col1), row_inds)
Base.getindex(col1::DatasetColumn{Dataset,<:AbstractVector}, row_inds::Colon) = getindex(__!(col1), row_inds)

# changing one observation reserves format but if the column belongs to sorting columns reset all sorting info
function Base.setindex!(col1::DatasetColumn{Dataset,<:AbstractVector}, v, i)
    __!(col1)[i] = v
    _modified(_attributes(col1.ds))
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end
function Base.setindex!(col1::DatasetColumn{Dataset,<:AbstractVector}, v, i::Colon)
    __!(col1)[i] = v
    _modified(_attributes(col1.ds))
    removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end
# changing all values should reset the formatting
function Base.copy!(col1::DatasetColumn{Dataset,<:AbstractVector}, src)
    copy!(__!(col1), src)
    _modified(_attributes(col1.ds))
    removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end

function Base.map!(f, col1::DatasetColumn{Dataset,<:AbstractVector}, args...)
    map!(f, __!(col1), args...)
    _modified(_attributes(col1.ds))
    removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end
function Base.map!(f, col1::DatasetColumn{Dataset,<:AbstractVector}, col2::DatasetColumn{Dataset,<:AbstractVector})
    map!(f, __!(col1), __!(col2))
    _modified(_attributes(col1.ds))
    removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end

# ################
# SubDatasetColumn
# #################

Base.getindex(col1::SubDatasetColumn{Dataset,<:AbstractVector}, i::Integer) = getindex(__!(col1), i)
Base.getindex(col1::SubDatasetColumn{Dataset,<:AbstractVector}, i::CartesianIndex) = getindex(__!(col1), i)
Base.getindex(col1::SubDatasetColumn{Dataset,<:AbstractVector}, row_inds::AbstractVector) = getindex(__!(col1), row_inds)
Base.getindex(col1::SubDatasetColumn{Dataset,<:AbstractVector}, row_inds::Colon) = getindex(__!(col1), row_inds)

# changing one observation reserves format but if the column belongs to sorting columns reset all sorting info
function Base.setindex!(col1::SubDatasetColumn{Dataset,<:AbstractVector}, v, i)
    __!(col1)[i] = v
    _modified(_attributes(col1.ds))
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end
# changing all values should reset the formatting
function Base.copy!(col1::SubDatasetColumn{Dataset,<:AbstractVector}, src)
    copy!(__!(col1), src)
    _modified(_attributes(col1.ds))
    # removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end

function Base.map!(f, col1::SubDatasetColumn{Dataset,<:AbstractVector}, args...)
    map!(f, __!(col1), args...)
    _modified(_attributes(col1.ds))
    # removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end
function Base.map!(f, col1::SubDatasetColumn{Dataset,<:AbstractVector}, col2::SubDatasetColumn{Dataset,<:AbstractVector})
    map!(f, __!(col1), __!(col2))
    _modified(_attributes(col1.ds))
    # removeformat!(col1.ds, col1.col)
    if col1.col ∈ index(col1.ds).sortedcols
        _reset_grouping_info!(col1.ds)
    end
    col1
end


text:
module TestAst
  using Base.Test
  import JSON
  using RethinkDB

  r = RethinkDB

  t1 = r.db("bookstore") |>
    r.to_ast |> JSON.json
  @test t1 == "[1,[14,[\"bookstore\"]]]"

  t1 = r.table("authors") |>
    r.to_ast |> JSON.json
  @test t1 == "[1,[15,[\"authors\"]]]"

  t1 = r.db("bookstore") |>
    d -> r.table(d, "authors") |>
    r.to_ast |> JSON.json
  @test t1 == "[1,[15,[[14,[\"bookstore\"]],\"authors\"]]]"
end

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb


text:
using Pretend
using Test

Pretend.activate()

@testset "Pretend.jl" begin
    include("test_basic.jl")
    include("test_anon.jl")
    include("test_multiple_dispatch.jl")
    include("test_spy.jl")
    include("test_misc.jl")
    include("test_3rdparty.jl")
    include("test_inner_macro.jl")
end


text:
using OrdinaryDiffEq, DiffEqDevTools, Sundials, ParameterizedFunctions, Plots, LSODA
gr() # gr(fmt=:png)
using LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

const k1=.35e0
const k2=.266e2
const k3=.123e5
const k4=.86e-3
const k5=.82e-3
const k6=.15e5
const k7=.13e-3
const k8=.24e5
const k9=.165e5
const k10=.9e4
const k11=.22e-1
const k12=.12e5
const k13=.188e1
const k14=.163e5
const k15=.48e7
const k16=.35e-3
const k17=.175e-1
const k18=.1e9
const k19=.444e12
const k20=.124e4
const k21=.21e1
const k22=.578e1
const k23=.474e-1
const k24=.178e4
const k25=.312e1

function f_pollu(dy,y,p,t)
 r1  = k1 *y[1]
 r2  = k2 *y[2]*y[4]
 r3  = k3 *y[5]*y[2]
 r4  = k4 *y[7]
 r5  = k5 *y[7]
 r6  = k6 *y[7]*y[6]
 r7  = k7 *y[9]
 r8  = k8 *y[9]*y[6]
 r9  = k9 *y[11]*y[2]
 r10 = k10*y[11]*y[1]
 r11 = k11*y[13]
 r12 = k12*y[10]*y[2]
 r13 = k13*y[14]
 r14 = k14*y[1]*y[6]
 r15 = k15*y[3]
 r16 = k16*y[4]
 r17 = k17*y[4]
 r18 = k18*y[16]
 r19 = k19*y[16]
 r20 = k20*y[17]*y[6]
 r21 = k21*y[19]
 r22 = k22*y[19]
 r23 = k23*y[1]*y[4]
 r24 = k24*y[19]*y[1]
 r25 = k25*y[20]

 dy[1]  = -r1-r10-r14-r23-r24+
          r2+r3+r9+r11+r12+r22+r25
 dy[2]  = -r2-r3-r9-r12+r1+r21
 dy[3]  = -r15+r1+r17+r19+r22
 dy[4]  = -r2-r16-r17-r23+r15
 dy[5]  = -r3+r4+r4+r6+r7+r13+r20
 dy[6]  = -r6-r8-r14-r20+r3+r18+r18
 dy[7]  = -r4-r5-r6+r13
 dy[8]  = r4+r5+r6+r7
 dy[9]  = -r7-r8
 dy[10] = -r12+r7+r9
 dy[11] = -r9-r10+r8+r11
 dy[12] = r9
 dy[13] = -r11+r10
 dy[14] = -r13+r12
 dy[15] = r14
 dy[16] = -r18-r19+r16
 dy[17] = -r20
 dy[18] = r20
 dy[19] = -r21-r22-r24+r23+r25
 dy[20] = -r25+r24
end

function fjac_pollu(J,y,p,t)
      J .= 0.0
      J[1,1]   = -k1-k10*y[11]-k14*y[6]-k23*y[4]-k24*y[19]
      J[1,11]  = -k10*y[1]+k9*y[2]
      J[1,6]   = -k14*y[1]
      J[1,4]   = -k23*y[1]+k2*y[2]
      J[1,19]  = -k24*y[1]+k22
      J[1,2]   = k2*y[4]+k9*y[11]+k3*y[5]+k12*y[10]
      J[1,13]  = k11
      J[1,20]  = k25
      J[1,5]   = k3*y[2]
      J[1,10]  = k12*y[2]

      J[2,4]   = -k2*y[2]
      J[2,5]   = -k3*y[2]
      J[2,11]  = -k9*y[2]
      J[2,10]  = -k12*y[2]
      J[2,19]  = k21
      J[2,1]   = k1
      J[2,2]   = -k2*y[4]-k3*y[5]-k9*y[11]-k12*y[10]

      J[3,1]   = k1
      J[3,4]   = k17
      J[3,16]  = k19
      J[3,19]  = k22
      J[3,3]   = -k15

      J[4,4]   = -k2*y[2]-k16-k17-k23*y[1]
      J[4,2]   = -k2*y[4]
      J[4,1]   = -k23*y[4]
      J[4,3]   = k15

      J[5,5]   = -k3*y[2]
      J[5,2]   = -k3*y[5]
      J[5,7]   = 2k4+k6*y[6]
      J[5,6]   = k6*y[7]+k20*y[17]
      J[5,9]   = k7
      J[5,14]  = k13
      J[5,17]  = k20*y[6]

      J[6,6]   = -k6*y[7]-k8*y[9]-k14*y[1]-k20*y[17]
      J[6,7]   = -k6*y[6]
      J[6,9]   = -k8*y[6]
      J[6,1]   = -k14*y[6]
      J[6,17]  = -k20*y[6]
      J[6,2]   = k3*y[5]
      J[6,5]   = k3*y[2]
      J[6,16]  = 2k18

      J[7,7]   = -k4-k5-k6*y[6]
      J[7,6]   = -k6*y[7]
      J[7,14]  = k13

      J[8,7]   = k4+k5+k6*y[6]
      J[8,6]   = k6*y[7]
      J[8,9]   = k7

      J[9,9]   = -k7-k8*y[6]
      J[9,6]   = -k8*y[9]

      J[10,10] = -k12*y[2]
      J[10,2]  = -k12*y[10]+k9*y[11]
      J[10,9]  = k7
      J[10,11] = k9*y[2]

      J[11,11] = -k9*y[2]-k10*y[1]
      J[11,2]  = -k9*y[11]
      J[11,1]  = -k10*y[11]
      J[11,9]  = k8*y[6]
      J[11,6]  = k8*y[9]
      J[11,13] = k11

      J[12,11] = k9*y[2]
      J[12,2]  = k9*y[11]

      J[13,13] = -k11
      J[13,11] = k10*y[1]
      J[13,1]  = k10*y[11]

      J[14,14] = -k13
      J[14,10] = k12*y[2]
      J[14,2]  = k12*y[10]

      J[15,1]  = k14*y[6]
      J[15,6]  = k14*y[1]

      J[16,16] = -k18-k19
      J[16,4]  = k16

      J[17,17] = -k20*y[6]
      J[17,6]  = -k20*y[17]

      J[18,17] = k20*y[6]
      J[18,6]  = k20*y[17]

      J[19,19] = -k21-k22-k24*y[1]
      J[19,1]  = -k24*y[19]+k23*y[4]
      J[19,4]  = k23*y[1]
      J[19,20] = k25

      J[20,20] = -k25
      J[20,1]  = k24*y[19]
      J[20,19] = k24*y[1]

      return
end

u0 = zeros(20)
u0[2]  = 0.2
u0[4]  = 0.04
u0[7]  = 0.1
u0[8]  = 0.3
u0[9]  = 0.01
u0[17] = 0.007
prob = ODEProblem(ODEFunction(f_pollu, jac=fjac_pollu),u0,(0.0,60.0))

sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

# group 1
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>lsoda())]

# High Tolerance
abstols = 1 ./ 10 .^ (5:8)
reltols = 1 ./ 10 .^ (1:4);

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
DIAGRAMS["POLLU-HighTol-Group1"] = plot(wp)

# group 2
setups = [Dict(:alg=>Rodas4()),
          Dict(:alg=>Rodas5()),
          Dict(:alg=>Rodas5P()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>lsoda())]

# Low Tolerance
abstols = 1 ./ 10 .^ (7:13)
reltols = 1 ./ 10 .^ (4:10)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
DIAGRAMS["POLLU-LowTol-Group2"] = plot(wp)


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelayEmbeddings, DynamicalSystemsBase
using DifferentialEquations
using Random
using Test
using DelimitedFiles

include("../../src/pecuzal_method.jl")

## In this script we evaluate the dependence of the returned reconstruction
# parameters on the parameter α, needed for the computation of the
# contiuity-statistic

## Dependence on α

# The time series have been computed in the script `/test/timeseries/produce_timeseries.jl`
s = readdlm("./test/timeseries/lorenz_pecora_uni_x.csv")
s = vec(s[1:5000]) # input timeseries = x component of lorenz
w = estimate_delay(s, "mi_min")
Tmax = 100
KNN = 3
samplesize = 1
Tw = 20
K = 14
alphas = [0.05 0.01]

τ_vals = []
Ls = []
sizes = []
εs = []
for alpha in alphas
    display(alpha)
    YY, τ_valss, _, Lss , ε = pecuzal_embedding_update(s;
                                τs = 0:Tmax , w = w, samplesize = samplesize,
                                K = K, KNN = KNN, Tw = Tw, α = alpha)
    push!(sizes,size(YY,2))
    push!(τ_vals, τ_valss)
    push!(Ls, Lss)
    push!(εs, ε)
end

writedlm("./scripts/computed data/dependence_on_alpha_alphas.csv",alphas)
writedlm("./scripts/computed data/dependence_on_alpha_sizes.csv",sizes)
writedlm("./scripts/computed data/dependence_on_alpha_Ls.csv",Ls)
writedlm("./scripts/computed data/dependence_on_alpha_tau_vals.csv",τ_vals)
writedlm("./scripts/computed data/dependence_on_alpha_epsilons.csv",εs)

## Plot results (see script `plot_results_params_dependencies.jl`)

# alphas = readdlm("./scripts/computed data/dependence_on_alpha_alphas.csv")
# sizes = readdlm("./scripts/computed data/dependence_on_alpha_sizes.csv")
# Ls = readdlm("./scripts/computed data/dependence_on_alpha_Ls.csv")
# τ_vals = readdlm("./scripts/computed data/dependence_on_alpha_tau_vals.csv")
# εs = readdlm("./scripts/computed data/dependence_on_alpha_epsilons.csv")
#
# using PyPlot
# pygui(true)
#
# panelnames = ["A" "B"]
#
# lwg = 3         # linewidth of the graph
# lwa = 2         # linewidth of the axis
# fsp = 20        # Fontsize of panelnames
# fsa = 16        # Fontsize of the axis
# fsl = 12        # Fontsize of the legendentries
# fst = 16        # Fontsize of title
# axislabelsize = 16 # axislabelsize
# ticklabelsize = 12  # labelsize of ticks
# ms = 10         # markersize of maxima
# ms_style = "*"  # markerstyle
#
# taus = vec(0:Tmax)
# fig = figure(figsize=[15,8])
#
# subplots_adjust(left = 0.07,
# right = 0.96,
# bottom = 0.08,
# top = 0.9,
# wspace = 0.25,
# hspace = 0.25)
#
# rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
# font0 = Dict(
#         "font.size" => fsa,
#         "axes.labelweight" => "normal",
#         "axes.labelsize" => axislabelsize,
#         "axes.linewidth" => lwa,
#         "xtick.labelsize" => ticklabelsize,
#         "ytick.labelsize" => ticklabelsize,
#         "legend.fontsize" => fsl)
# merge!(rcParams, font0)
#
# ylims = [0, 1.4]
#
# ax = fig.add_subplot(1, 2, 1)
# plot(taus,εs[1][:,1], label="τs=$(τ_vals[1][1])", linewidth = lwg)
# plot(taus,εs[1][:,2], label="τs=$(τ_vals[1][2])", linewidth = lwg)
# plot(taus,εs[1][:,3], label="τs=$(τ_vals[1][3])", linewidth = lwg)
# plot([taus[τ_vals[1][2]+1]], [εs[1][τ_vals[1][2]+1,1]], ms_style, label="", markersize = ms, color = "black")
# plot([taus[τ_vals[1][3]+1]], [εs[1][τ_vals[1][3]+1,2]], ms_style, label="", markersize = ms, color = "black")
# legend()
# ax.set_ylim(ylims)
# ax.text(-0.15, 1.1, panelnames[1], transform=ax.transAxes,
#      fontsize=fsp, fontweight="bold", va="top")
# title("α=0.05")
# xlabel("τ")
# ylabel("⟨ε★⟩")
# grid()
#
# ax = fig.add_subplot(1, 2, 2)
# plot(taus,εs[2][:,1], label="τs=$(τ_vals[2][1])", linewidth = lwg)
# plot(taus,εs[2][:,2], label="τs=$(τ_vals[2][2])", linewidth = lwg)
# plot(taus,εs[2][:,3], label="τs=$(τ_vals[2][3])", linewidth = lwg)
# plot([taus[τ_vals[2][2]+1]], [εs[2][τ_vals[2][2]+1,1]], ms_style, label="", markersize = ms, color = "black")
# plot([taus[τ_vals[2][3]+1]], [εs[2][τ_vals[2][3]+1,2]], ms_style, label="", markersize = ms, color = "black")
# legend()
# ax.set_ylim(ylims)
# ax.text(-0.15, 1.1, panelnames[2], transform=ax.transAxes,
#      fontsize=fsp, fontweight="bold", va="top")
# title("α=0.01")
# xlabel("τ")
# grid()

max_stars_repo_path: scripts/dependence on input parameters/dependence_on_alpha.jl


text:
module HandCombinations

using LinearAlgebra: dot
using PlayingCards
using Combinatorics

export prime,
    straight_ranks,
    quad_ranks,
    full_house_ranks,
    is_straight,
    flush_ranks,
    trip_ranks,
    two_pair_ranks,
    pair_ranks,
    high_card_ranks

##### Helpers

const primes = (41,2,3,5,7,11,13,17,19,23,29,31,37)
prime(i::Int8) = primes[i]
prime(i::Int) = prime(Int8(i))
prime(card::Card) = primes[rank(card)]

function sort_combinations(comb, sort_by_first_val = false)
    if length(comb[1])==1 || sort_by_first_val
        return sort(comb; by=x->high_value(x[1]), rev=true)
    else
        coeffs = (10 .^ (length(comb[1]):-1:2)..., 1)
        return sort(comb; by=x->dot(high_value.(x), coeffs), rev=true)
    end
end

##### Rows 1:10 (Straight flush) | Rows 1600:1609 (Off-suit straight)

function straight_ranks()
    (vcat(10:13, 1:1),9:(9+4),8:(8+4),7:(7+4),6:(6+4),5:(5+4),4:(4+4),3:(3+4),2:(2+4),1:(1+4))
end

##### Rows 11:166 (4 of a kind)

function ranks_old()
    r = ranks()
    return [r[2:end]..., r[1]]
end

function ranks_high_to_low()
    r = ranks()
    return [r[1], r[end:-1:2]...]
end

function quad_ranks()
    quads_combos = collect(combinations(full_deck(), 4))
    sorted_quads_combos = sort.(quads_combos; by = x->high_value(x), rev=true)
    sorted_quads_combos = filter(sorted_quads_combos) do x
        high_value(x[1])==high_value(x[2])==high_value(x[3])==high_value(x[4])
    end
    sorted_quads_combos = sort_combinations(sorted_quads_combos)
    card_ranks = [begin
        (rank.(quads)..., kicker)
    end for quads in sorted_quads_combos for kicker in ranks_old()[end:-1:1] if rank(quads[1]) ≠ kicker]
    return card_ranks
end

##### Rows 167:322 (full house)

function full_house_ranks()
    three_quarter_deck = filter(x->suit(x) == ♣ || suit(x) == ♡ || suit(x) == ♠, full_deck())
    half_deck = filter(x->suit(x) == ♣ || suit(x) == ♡, full_deck())
    trip_combos = collect(combinations(three_quarter_deck, 3))
    pair_combos = collect(combinations(half_deck, 2))
    trip_combos = sort.(trip_combos; by = x->high_value(x), rev=true)
    pair_combos = sort.(pair_combos; by = x->high_value(x), rev=true)

    sorted_pair_combos = filter(pair_combos) do x
        high_value(x[1])==high_value(x[2])
    end
    sorted_trip_combos = filter(trip_combos) do x
        high_value(x[1])==high_value(x[2])==high_value(x[3])
    end
    sorted_trip_combos = sort_combinations(sorted_trip_combos, true)
    sorted_pair_combos = sort_combinations(sorted_pair_combos, true)
    card_ranks = [begin
        (rank.(tc)..., rank.(pc)...)
    end for tc in sorted_trip_combos for pc in sorted_pair_combos if rank(tc[1]) ≠ rank(pc[1])]
    return card_ranks
end

##### Rows 323:1599 (flush)

consecutive(tup) = all(ntuple(i->tup[i]+1==tup[i+1], 4))
function is_straight(cards)
    ranks = sort(collect(high_value.(cards)))
    ranks_low = sort(collect(low_value.(cards)))
    high_straight = consecutive(ranks)
    low_straight = consecutive(ranks_low)
    return low_straight || high_straight
end

function flush_ranks()
    club_combos = combinations(filter(x->suit(x) == ♣, full_deck()), 5)
    sorted_club_combos = sort.(club_combos; by = x->high_value(x), rev=true)
    sorted_club_combos = sort_combinations(sorted_club_combos)
    card_ranks = [begin
        rank.(cards)
    end for cards in sorted_club_combos if !is_straight(cards)]
end

##### Rows 1610:2467 (trips)

function trip_ranks()
    club_deck = filter(x->suit(x) == ♣, full_deck())
    club_kicker_combos = collect(combinations(club_deck, 2))
    sorted_club_kicker_combos = sort.(club_kicker_combos; by = x->high_value(x), rev=true)
    sorted_club_kicker_combos = sort_combinations(sorted_club_kicker_combos)
    card_ranks = [begin
        (rank_trips,rank_trips,rank_trips,rank.(kickers)...)
    end for rank_trips in ranks_high_to_low() for kickers in sorted_club_kicker_combos
        if rank_trips ≠ rank(kickers[1]) && rank_trips ≠ rank(kickers[2])]
    return card_ranks
end

##### Rows 2468:3325 (two pair)

function two_pair_ranks()
    half_deck = filter(x->suit(x) == ♣ || suit(x) == ♡, full_deck())
    combos = collect(combinations(half_deck, 4))
    combos = sort.(combos; by = x->high_value(x), rev=true)
    two_pair_combos = filter(x->high_value(x[1])==high_value(x[2]) && high_value(x[3])==high_value(x[4]), combos)
    sorted_two_pair_combos = sort_combinations(two_pair_combos)
    card_ranks = [begin
        (rank.(two_pair)...,kicker)
    end for two_pair in sorted_two_pair_combos for kicker in ranks_high_to_low()
        if rank(two_pair[1]) ≠ kicker && rank(two_pair[3]) ≠ kicker]
    return card_ranks
end

##### Rows 3326:6185 (pair)

function pair_ranks()
    three_quarters_deck = filter(x->suit(x) == ♣, full_deck())
    combos = collect(combinations(three_quarters_deck, 3))
    combos = sort.(combos; by = x->high_value(x), rev=true)
    combos = sort_combinations(combos)
    card_ranks = [begin
        (ranks_pair, ranks_pair, rank.(kickers)...)
    end for ranks_pair in ranks_high_to_low() for kickers in combos
        if ranks_pair ≠ rank(kickers[1]) && ranks_pair ≠ rank(kickers[2]) && ranks_pair ≠ rank(kickers[3])]
    return card_ranks
end

##### Rows 6186:7462 (high card)

function high_card_ranks()
    club_deck = filter(x->suit(x) == ♣, full_deck())
    combos = collect(combinations(club_deck, 5))
    combos = sort.(combos; by = x->high_value(x), rev=true)
    combos = sort_combinations(combos)
    card_ranks = [kickers for kickers in combos if !is_straight(kickers)]
    return card_ranks
end

end # module


text:
#===============================================================================

    your_own_analysis.jl

    This is a template for you to create your own analytic function.
    Everything here is commented out, so no code here has any effect.
    You need to uncomment and complete the code to make it work.
    You also need to register your own functionality in "AmpModel.jl":
    see the lines that contain `YourOwnAnalysis`.

===============================================================================#

# function act(
#     ::Type{Val{your_own_analysis}},
#     req ::Associative;
#     # whatever keyword arguments you want to add
#     )
#     ## put your own code here
# end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# ---
# title: 802. Find Eventual Safe States
# id: problem802
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Graph
# link: <https://leetcode.com/problems/find-eventual-safe-states/description/>
# hidden: true
# ---
#
# In a directed graph, we start at some node and every turn, walk along a
# directed edge of the graph.  If we reach a node that is terminal (that is, it
# has no outgoing directed edges), we stop.
#
# Now, say our starting node is _eventually safe  _if and only if we must
# eventually walk to a terminal node.  More specifically, there exists a natural
# number `K` so that for any choice of where to walk, we must have stopped at a
# terminal node in less than `K` steps.
#
# Which nodes are eventually safe?  Return them as an array in sorted order.
#
# The directed graph has `N` nodes with labels `0, 1, ..., N-1`, where `N` is
# the length of `graph`.  The graph is given in the following form: `graph[i]`
# is a list of labels `j` such that `(i, j)` is a directed edge of the graph.
#
#
#
#     **Example:**
#     Input: graph = [[1,2],[2,3],[5],[0],[5],[],[]]
#     Output: [2,4,5,6]
#     Here is a diagram of the above graph.
#
#
#
# ![Illustration of graph](https://s3-lc-
# upload.s3.amazonaws.com/uploads/2018/03/17/picture1.png)
#
# **Note:**
#
#   * `graph` will have length at most `10000`.
#   * The number of edges in the graph will not exceed `32000`.
#   * Each `graph[i]` will be a sorted list of different integers, chosen within the range `[0, graph.length - 1]`.
#
#
## @lc code=start
using LeetCode

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


text:
##### Function to use the photoevaporation boundary in Carrera et al 2018 to separate planets:

function photoevap_boundary_Carrera2018(R::Float64, P::Float64)
    #R is the planet radius in Earth radii, P is the period in days
    #This function returns 1 if the planet is above the boundary, and 0 if the planet is below the boundary as defined by Eq. 5 in Carrera et al 2018
    Rtrans = 2.6*P^(-0.1467)
    if R > Rtrans
        return 1
    elseif R < Rtrans
        return 0
    end
end


text:
using MLDataUtils, LinearAlgebra
using DeepGP
using AugmentedGaussianProcesses
using KernelFunctions
using Flux
using Zygote
using Distributions

X,y = noisy_sin(100,-1,1)
X = reshape(X,:,1)

ngp = 1; ndim = 15
μ = zeros(ndim); Σ = diagm(ones(ndim));
Z = collect(reshape(range(-1,1,length=ndim),:,1))
l = 10.0
l1 = DeepGP.SVGPLayer(ngp,μ,Σ,Z,SqExponentialKernel(l),ZeroMean())
l2 = DeepGP.SVGPLayer(1,μ,Σ,rand(ndim,ngp),SqExponentialKernel(l),ZeroMean())
v1 = DeepGP.propagate(l1,X)
v2 = DeepGP.propagate(l2,v1)
m = DeepGP.DeepGPModel(l1)
DeepGP.expec_log_like(m,X,y)
DeepGP.ELBO(m,X,y) |> display

Zygote.refresh()
Zygote.gradient(x->DeepGP.ELBO(m,x,y),X)
[X for _ in 1:nSamples]
t


text:
int_rules_2_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 2.2*(c+d*x)^m*(F^(g*(e+f*x)))^n*(a+b*(F^(g*(e+f*x)))^n)^p =#
    @apply_utils Antiderivative(((~(c') + ~(d') * ~x) ^ ~(m') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n')) / (~a + ~(b') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n')), ~x) => ((~c + ~d * ~x) ^ ~m / (~b * ~f * ~g * ~n * log(~F))) * log(1 + (~b * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) / ~a) - ((~d * ~m) / (~b * ~f * ~g * ~n * log(~F))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * log(1 + (~b * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) / ~a), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n') * (~(a') + ~(b') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n')) ^ ~(p'), ~x) => ((~c + ~d * ~x) ^ ~m * (~a + ~b * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) ^ (~p + 1)) / (~b * ~f * ~g * ~n * (~p + 1) * log(~F)) - ((~d * ~m) / (~b * ~f * ~g * ~n * (~p + 1) * log(~F))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * (~a + ~b * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && NeQ(~p, -1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n') * (~(a') + ~(b') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n')) ^ ~(p'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n * (~a + ~b * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) ^ ~p, ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(k') * (~G) ^ (~(j') * (~(h') + ~(i') * ~x))) ^ ~(q') * (~(a') + ~(b') * ((~F) ^ (~(g') * (~(e') + ~(f') * ~x))) ^ ~(n')) ^ ~(p'), ~x) => ((~k * (~G) ^ (~j * (~h + ~i * ~x))) ^ ~q / ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) * Antiderivative((~c + ~d * ~x) ^ ~m * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n * (~a + ~b * ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n) ^ ~p, ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~j, ~k, ~m, ~n, ~p, ~q], ~x) && (EqQ(~f * ~g * ~n * log(~F) - ~i * ~j * ~q * log(~G), 0) && NeQ((~k * (~G) ^ (~j * (~h + ~i * ~x))) ^ ~q - ((~F) ^ (~g * (~e + ~f * ~x))) ^ ~n, 0))
end

max_stars_repo_path: src/rules/2 Exponentials/2.2 (c+d x)^m (F^(g (e+f x)))^n (a+b (F^(g (e+f x)))^n)^p.jl


text:
using NCDatasets
using Test
using Dates
using Printf
using BenchmarkTools

sz = (400,400)
data = randn(sz)

filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDataset(filename,"c")

# define the dimension "lon" and "lat"
defDim(ds,"lon",sz[1])
defDim(ds,"lat",sz[2])

T = Float64
v = defVar(ds,"var-$T",T,("lon","lat"))
# write array
v[:,:] = data
close(ds)

function sequential_access!(filename,diskless,d)
    sz = size(d)
    ds = NCDataset(filename,"r", diskless = true)
    v = ds["var-Float64"]

    for j = 1:sz[2]
        for i = 1:sz[1]
            d[i,j] = v[i,j]
        end
    end

    close(ds)
end

function random_access!(filename,diskless,d)
    sz = size(d)
    ds = NCDataset(filename,"r", diskless = diskless)
    v = ds["var-Float64"]

    for I in Random.shuffle(CartesianIndices(sz)[:])
        d[I] = v[I]
    end
    close(ds)
end

d = similar(data)

for diskless in [false,true]
    println("sequential access, diskless=$(diskless)")
    @btime sequential_access!($filename,$diskless,$d)
    @test data == d

    println("random access, diskless=$(diskless)")
    @btime random_access!($filename,$diskless,$d)
    @test data == d
end


text:
#Base
Base.parse(::T, ::Missing, kwargs...) where T = missing

# Dates.jl
Dates.Date(::Missing, kwargs...) = missing
Dates.DateTime(::Missing, kwargs...) = missing
Dates.year(::Missing) = missing
Dates.month(::Missing) = missing
Dates.day(::Missing) = missing
Dates.canonicalize(::Missing) = missing
Dates.yearmonthday(::Missing) = missing
Dates.yearmonth(::Missing) = missing
Dates.monthday(::Missing) = missing
Dates.week(::Missing) = missing
Dates.hour(::Missing) = missing
Dates.minute(::Missing) = missing
Dates.second(::Missing) = missing
Dates.millisecond(::Missing) = missing
Dates.dayofmonth(::Missing) = missing
Dates.microsecond(::Missing) = missing
Dates.nanosecond(::Missing) = missing
Dates.dayofweek(::Missing) = missing
Dates.isleapyear(::Missing) = missing
Dates.daysinmonth(::Missing) = missing
Dates.daysinyear(::Missing) = missing
Dates.dayofyear(::Missing) = missing
Dates.dayname(::Missing) = missing
Dates.dayabbr(::Missing) = missing
Dates.dayofweekofmonth(::Missing) = missing
Dates.daysofweekinmonth(::Missing) = missing
Dates.monthname(::Missing) = missing
Dates.monthabbr(::Missing) = missing
Dates.quarterofyear(::Missing) = missing
Dates.dayofquarter(::Missing) = missing
Dates.unix2datetime(::Missing) = missing
Dates.datetime2unix(::Missing) = missing
Dates.firstdayofweek(::Missing) = missing
Dates.lastdayofweek(::Missing) = missing
Dates.firstdayofmonth(::Missing) = missing
Dates.lastdayofmonth(::Missing) = missing
Dates.firstdayofyear(::Missing) = missing
Dates.lastdayofyear(::Missing) = missing
Dates.firstdayofquarter(::Missing) = missing
Dates.lastdayofquarter(::Missing) = missing


text:
#!/usr/bin/env julia

using HDF5

function matFromCool(coolName,startCoord,endCoord,resolution)
	attrs = h5readattr(coolName,resolution)
	binsize = attrs["bin-size"]

	binStart = startCoord/binsize
	binEnd = endCoord/binsize
	matwidth = floor(Int,binEnd-binStart+1)

	matrix = zeros(1,1)
	h5open(coolName) do cin
		onek=cin[resolution]
		pixels = onek["pixels"]
		bin1_id = read(pixels,"bin1_id")
		bin2_id = read(pixels,"bin2_id")
		count = read(pixels,"count")

		matrix = zeros(matwidth,matwidth)
		nnv = length(count)
		for binIdx in 1:nnv
			row = bin1_id[binIdx]+1
			col = bin2_id[binIdx]+1
			if (row <= matwidth) && (col <= matwidth)
				matrix[row,col] = count[binIdx]
			end
		end
	end
	return matrix
end

coolName = "MB286.mcool"
startCoord=1000
endCoord = 10000000
resolution = "resolutions/1000"

matrix = matFromCool(coolName,startCoord,endCoord,resolution)
println(length(matrix))

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
struct Block <: FEMContainer
    block_id::BlockID  # TODO: maybe change to BlockID so types are more verbose
    num_elem::Int64
    num_nodes_per_elem::Int64
    elem_type::String
    conn::Array{Int64}
    function Block(exo_id::ExoID, block_id::BlockID)
        element_type, num_elem, num_nodes, _, _, _ =
        read_element_block_parameters(exo_id::ExoID, block_id::BlockID)
        conn = read_block_connectivity(exo_id, block_id)
        # conn = reshape(conn, (num_elem, num_nodes))  # for easier access downstream
        return new(block_id, num_elem, num_nodes, element_type, conn)
    end
end
Base.show(io::IO, block::Block) =
print(io, "Block:\n",
          "\tBlock ID           = ", block.block_id, "\n",
          "\tNum elem           = ", block.num_elem, "\n",
          "\tNum nodes per elem = ", block.num_nodes_per_elem, "\n",
          "\tElem type          = ", block.elem_type, "\n")

# more verbose types
#
BlockIDs = Vector{BlockID}
BlockIDsPtr = Ref{BlockID}
Blocks = Vector{Block}

# methods below
#

function read_block_ids(exo_id::ExoID, num_elem_blk::Int64)::BlockIDs
    block_ids = BlockIDs(undef, num_elem_blk)
    error = ccall((:ex_get_ids, libexodus), ExodusError,
                (ExoID, ExodusConstant, BlockIDsPtr),
                exo_id, EX_ELEM_BLOCK, block_ids)
    exodus_error_check(error, "read_block_ids")
    return block_ids
end

function read_element_block_parameters(exo_id::ExoID, block_id::BlockID)
    element_type = Vector{UInt8}(undef, MAX_STR_LENGTH)
    num_elem = Ref{Int64}(0)
    num_nodes = Ref{Int64}(0)
    num_edges = Ref{Int64}(0)
    num_faces = Ref{Int64}(0)
    num_attributes = Ref{Int64}(0)
    error = ccall((:ex_get_block, libexodus), ExodusError,
                  (ExoID, ExodusConstant, BlockID,
                   Ptr{UInt8}, Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{Int64}),
                  exo_id, EX_ELEM_BLOCK, block_id,
                  element_type, num_elem, num_nodes, num_edges, num_faces, num_attributes)

    exodus_error_check(error, "read_element_blocK_parameters")

    element_type = unsafe_string(pointer(element_type))

    return element_type, num_elem[], num_nodes[], num_edges[], num_faces[], num_attributes[]
end

function read_block_connectivity(exo_id::ExoID, block_id::BlockID)
    element_type, num_elem, num_nodes, num_edges, num_faces, num_attributes =
    read_element_block_parameters(exo_id::ExoID, block_id::BlockID)

    conn = Vector{Int32}(undef, num_nodes * num_elem)
    conn_face = Vector{Int32}(undef, num_nodes * num_elem)  # Not using these currently
    conn_edge = Vector{Int32}(undef, num_nodes * num_elem)  # Not using these currently

    # TODO: look into why the connectivity arrays need to be 32 bit.
    #
    error = ccall((:ex_get_conn, libexodus), ExodusError,
                  (ExoID, ExodusConstant, BlockID, Ref{Int32}, Ref{Int32}, Ref{Int32}),
                  exo_id, EX_ELEM_BLOCK, block_id, conn, conn_face, conn_edge)
    exodus_error_check(error, "read_block_connectivity")
    return conn
end

function read_blocks!(exo_id::ExoID, block_ids::BlockIDs, blocks::Blocks)
    for (n, block_id) in enumerate(block_ids)
        blocks[n] = Block(exo_id, block_id)
    end
end

function read_blocks(exo_id::ExoID, block_ids::BlockIDs)
    blocks = Vector{Block}(undef, size(block_ids, 1))
    read_blocks!(exo_id, block_ids, blocks)
    return blocks
end


text:
int_rules_4_3_1_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.1.3*(d*sin)^m*(a+b*tan)^n =#
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~x) ^ ~n) / ((~b) ^ 2 + (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, ~b * tan(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && IntegerQ(~m / 2)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative(expand(sin(~e + ~f * ~x) ^ ~m * (~a + ~b * tan(~e + ~f * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ((~m - 1) / 2) && IGtQ(~n, 0))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative((sin(~e + ~f * ~x) ^ ~m * (~a * cos(~e + ~f * ~x) + ~b * sin(~e + ~f * ~x)) ^ ~n) / cos(~e + ~f * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ((~m - 1) / 2) && (ILtQ(~n, 0) && (LtQ(~m, 5) && GtQ(~n, -4) || EqQ(~m, 5) && EqQ(~n, -1))))
    @apply_utils Antiderivative((~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~d * csc(~e + ~f * ~x)) ^ FracPart(~m) * (sin(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~n / (sin(~e + ~f * ~x) / ~d) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~n], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(m') * sin(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative(cos(~e + ~f * ~x) ^ (~m - ~n) * sin(~e + ~f * ~x) ^ ~p * (~a * cos(~e + ~f * ~x) + ~b * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~p], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * cos(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * cot(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative(sin(~e + ~f * ~x) ^ (~m - ~n) * cos(~e + ~f * ~x) ^ ~p * (~a * sin(~e + ~f * ~x) + ~b * cos(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~p], ~x) && IntegerQ(~n)
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.1.3 (d sin)^m (a+b tan)^n.jl


text:
@block SebastianM_C ["3d", "diffeq"] begin
    @cell "DifferentialEquations path animation" [lines, meshscatter, animated, record, diffeq] begin
        using OrdinaryDiffEq
        using Statistics

        function lorenz(du,u,p,t)
            du[1] = 10.0*(u[2]-u[1])
            du[2] = u[1]*(28.0-u[3]) - u[2]
            du[3] = u[1]*u[2] - (8/3)*u[3]
        end

        u0 = [1.0; 0.0; 0.0]
        tspan = (0.0, 10.0)
        prob = ODEProblem(lorenz, u0, tspan)
        sol = solve(prob, Tsit5())

        t = Node(tspan[1])
        trajectory = lift(t; init = [Point3f0(sol(t[], idxs = [1,2,3]))]) do t
            push!(trajectory[], Point3f0(sol(t, idxs=[1,2,3])))
        end
        endpoint = lift(t->[Point3f0(sol(t, idxs = [1,2,3]))], t)

        xm, xM = extrema(sol[1,:])
        ym, yM = extrema(sol[2,:])
        zm, zM = extrema(sol[3,:])

        limits = FRect3D((xm,ym,zm), (xM-xm,yM-ym,zM-zm))

        scene = lines(trajectory, limits = limits, markersize = 0.7)
        meshscatter!(scene, endpoint, limits = limits, markersize = 0.5, color = :red)
        eyepos = Vec3f0(100, 5.0, 25)
        N = 1000
        record(scene, @replace_with_a_path(mp4), LinRange(tspan..., N)) do tᵢ
            global eyepos
            t[] = tᵢ
            eyepos = Vec3f0(eyepos[1] - (N / 10000), eyepos[2], eyepos[3])
            update_cam!(scene, eyepos, mean(trajectory[]))
            sleep(0.001)
        end
    end

end

@block AnshulSinghvi ["interactive", "diffeq"]  begin

    @cell "Interactive Differential Equation" [lines, slider, animated, record, diffeq, interactive] begin
        using OrdinaryDiffEq, ParameterizedFunctions
        import AbstractPlotting: textslider
        lorenz = @ode_def Lorenz begin           # define the system
            dx = σ * (y - x)
            dy = x * (ρ - z) - y
            dz = x * y - β*z
        end σ ρ β

        u0 = [1.0,0.0,0.0]                       # initial conditions
        tspan0 = (0.0,100.0)                      # initial timespan
        p0 = [10.0,28.0,8/3]                      # initial parameters
        prob = ODEProblem(lorenz, u0, tspan0, p0)  # define the problem

        ## setup sliders and plotting

        "The order of magnitude to range between."
        OME = 8

        sσ, oσ = textslider(exp10.(-OME:0.001:OME), "σ", start = p0[1]);

        sρ, oρ = textslider(exp10.(-OME:0.001:OME), "ρ", start = p0[2]);

        sβ, oβ = textslider(exp10.(-OME:0.001:OME), "β", start = p0[3]);

        st, ot = textslider(exp10.(-OME:0.001:OME), "tₘₐₓ", start = tspan0[end]);

        sr, or = textslider(100:10000, "resolution", start = 2000);

        trange = lift(ot, or) do tmax, resolution
            LinRange(0.0, tmax, resolution)
        end

        data = lift(oσ, oρ, oβ, trange) do σ, ρ, β, ts
            Point3f0.(
                solve(
                    remake(
                        prob;
                        p = [σ, ρ, β],
                        tspan = (ts[1], ts[end])
                    ),
                    Tsit5(),
                )(ts).u
            )  # change to fit the dimensionality - maybe even return 2 arrays, or a set of `Point2`s...
        end
        parent = Scene(resolution = (1000, 500))

        three = lines(
            data, linewidth = 2,
            transparency = true, color = ("#fe4a49", 0.4),
            show_axis = false
        )

        scatter!(
            three, data, markersize = 0.3, color = (:white, 0.3),
            strokecolor = :black, strokewidth = 1
        )
        on(data) do x
            # center camera etc
            update!(three)
        end
        scene = vbox(hbox(sσ, sρ, sβ, st, sr), three, parent = parent)

        # Do not execute beyond this point!

        RecordEvents(scene, @replace_with_a_path)

    end

end


text:
const PLOTBM = true     # save plot to img?
const PROCLEN = 8
const STDNRUNS = 5000000
const STDCP = 2:PROCLEN-1
const STDCFGS = [2,4,8,12,26,50,200,1000,4000,20000,120000,400000,1000000,6000000]

include("bm_setup.jl")
if PLOTBM
    include("bm_plot.jl")
end

# scale, scalar
type Scale_incx1 <: FuncType end
type Scale_incxnu <: FuncType end
type Scale_2d_incx1 <: FuncType end
type Scale_oop_incx1 <: FuncType end
type Scale_oop_incxnu <: FuncType end
# scale, array
type Scalearr_incx1 <: FuncType end
type Scalearr_incxnu <: FuncType end
type Scalearr_oop_incx1 <: FuncType end
type Scalearr_oop_incxnu <: FuncType end
# add, scalar

# add, array
type Addarr_incx1 <: FuncType end
type Addarr_oop_incx1 <: FuncType end
# addscal
type Addscal_incx1 <: FuncType end
type Addscal_oop_incx1 <: FuncType end
# copy
type Copy_incx1 <: FuncType end
type Copy_incxnu <: FuncType end
# fill
type Fill_incx1 <: FuncType end
type Fill_incxnu <: FuncType end

const METHODLIST = (
    (Scale_incx1,       "x = a*x (i=1, inc=1)",     STDNRUNS, STDCP, STDCFGS),
    (Scale_incxnu,      "x = a*x (i=2, inc=40)",    STDNRUNS<<3, STDCP, STDCFGS[6:end]),
    (Scale_2d_incx1,    "x = a*x 2-d (columns)",    STDNRUNS>>5, STDCP, STDCFGS),
    (Scale_oop_incx1,   "x = a*y (i=1, inc=1)",     STDNRUNS, STDCP, STDCFGS),
    (Scale_oop_incxnu,  "x = a*y (i=2, inc=40)",    STDNRUNS<<1, STDCP, STDCFGS[6:end]),
    (Scalearr_incx1,    "x = x.*y (i=1, inc=1)",    STDNRUNS>>3, [2,4,6,7], STDCFGS),
    (Scalearr_incxnu,   "x = x.*y (i=2, inc=40)",   STDNRUNS, [4,6,7], STDCFGS[6:end]),
    (Scalearr_oop_incx1, "x = y.*z (i=1, inc=1)",   STDNRUNS>>3, [2,4,6,7], STDCFGS),
    (Scalearr_oop_incxnu, "x = y.*z (i=2, inc=40)", STDNRUNS, [4,6,7], STDCFGS[6:end]),
    (Addarr_incx1,      "x = x + y (i=1, inc=1)",   STDNRUNS>>3, [4,5,6,7], STDCFGS),
    (Addarr_oop_incx1,  "x = y + z (i=1, inc=1)",   STDNRUNS>>3, [4,5,6,7], STDCFGS),
    (Addscal_incx1,     "x = x + a*y (i=1, inc=1)", STDNRUNS>>3, [4,5,6,7], STDCFGS),
    (Addscal_oop_incx1, "x = y + a*z (i=1, inc=1)", STDNRUNS>>3, [4,5,6,7], STDCFGS),
    (Copy_incx1,        "x = y (i=1, inc=1)",       STDNRUNS, STDCP, STDCFGS),
    (Copy_incxnu,       "x = y (i=2, inc=40)",      STDNRUNS<<3, [2,4,5,6,7], STDCFGS),
    (Fill_incx1,        "x = a (a=0, i=1, inc=1)",  STDNRUNS<<3, [2,4,6,7], STDCFGS),
    (Fill_incxnu,       "x = a (a=0, i=2, inc=40)", STDNRUNS<<3, [2,4,6,7], STDCFGS),
          )
const NM = length(METHODLIST)

global bm_method, bm_title, chooseproc, cfgs, FUNT, Nruns
for mi = NM-1:NM   #NM
    FUNT, bm_title, Nruns, chooseproc, cfgs = METHODLIST[mi]
    bm_method = string(FUNT())
    println(bm_title)
    include(string(bm_method,".jl"))
    t_a = time()
    include("bm_perform.jl")
    t_b = time()
    println("time: ", string(round(t_b-t_a, 1)), "s\n")
    if PLOTBM
        savetableimg(rtable, bm_title, string("./master/", string(FUNT())), FAOproc, chooseproc)
    end
end


text:
# Cassette Based Memoization
Cassette.@context MemoizeCtx

function Cassette.overdub(ctx::MemoizeCtx, x::RandVar, args...)
  if x.id in keys(ctx.metadata)
    ctx.metadata[x.id]
  else
    Cassette.recurse(ctx, x, args...)
  end
end

function Cassette.posthook(ctx::MemoizeCtx, retval, x::RandVar, args...)
  ctx.metadata[x.id] = retval
end

# Ω Based Memoization


text:
"""
  `value(ev::Event) :: Any`

Returns the value of the event.
"""
function value(ev::Event) :: Any
  return ev.value
end

function state(ev::Event) :: EVENT_STATE
  ev.state
end

"""
  `append_callback(ev::Event, cb::Function, args...)` :: Function

Adds a callback function to an event, i.e. a function having as first argument an object of type `Simulation` and as second argument the event. Optional arguments can be specified by `args...`.

If the event is being processed an [`EventProcessing`](@ref) exception is thrown.
"""
function append_callback(ev::Event, cb::Function, args::Any...) :: Function
  if ev.state == processed
    throw(EventProcessed())
  end
  func = (sim::Simulation, ev::Event)->cb(sim, ev, args...)
  ev.callbacks[func] = ev.cid+=one(UInt)
  return func
end

function remove_callback(ev::Event, func::Function)
  dequeue!(ev.callbacks, func)
end

type StateValue
  state :: EVENT_STATE
  value :: Any
  function StateValue(state::EVENT_STATE, value::Any=nothing)
    new(state, value)
  end
end

function Event(eval::Function, fev::Event, events...)
  oper = Event()
  event_state_values = Dict{Event, StateValue}()
  for ev in tuple(fev, events...)
    event_state_values[ev] = StateValue(ev.state)
    append_callback(ev, check, oper, eval, event_state_values)
  end
  return oper
end

function check(sim::Simulation, ev::Event, oper::Event, eval::Function, event_state_values::Dict{Event, StateValue})
  if oper.state == idle
    if isa(ev.value, Exception)
      schedule(sim, oper, value=ev.value)
    else
      event_state_values[ev] = StateValue(ev.state, ev.value)
      if eval(collect(values(event_state_values)))
        schedule(sim, oper, value=event_state_values)
      end
    end
  elseif oper.state == triggered
    if isa(ev.value, Exception)
      schedule!(sim, oper, priority=true, value=ev.value)
    else
      event_state_values[ev] = StateValue(ev.state, ev.value)
    end
  end
end

function eval_and(state_values::Vector{StateValue})
  return all(map((sv)->sv.state == processed, state_values))
end

function eval_or(state_values::Vector{StateValue})
  return any(map((sv)->sv.state == processed, state_values))
end

function (&)(ev1::Event, ev2::Event)
  return Event(eval_and, ev1, ev2)
end

function (|)(ev1::Event, ev2::Event)
  return Event(eval_or, ev1, ev2)
end


text:
export difference

# alias for set difference
import Base: \

const IA = IntervalArithmetic

"""
    \\(X::LazySet, Y::LazySet)

Convenience alias for set difference.

### Input

- `X` -- a set
- `Y` -- another set

### Output

The set difference between `X` and `Y`.

### Notes

If `X` and `Y` are intervals, `X \\ Y` is used in some libraries to denote
the left division, as the example below shows. However, it should not be
confused with the *set difference*. For example,

```jldoctest
julia> X = Interval(0, 2); Y = Interval(1, 4);

julia> X \\ Y   # computing the set difference
Interval{Float64, IntervalArithmetic.Interval{Float64}}([0, 1])

julia> X.dat \\ Y.dat  # computing the left division
[0.5, ∞]
```
"""
\(X::LazySet, Y::LazySet) = difference(X, Y)

# =================================
# Set difference between intervals
# =================================

"""
    difference(I1::IN, I2::IN) where {N, IN<:Interval{N}}

Return the set difference between the given intervals.

The set difference is defined as:

```math
    I₁ \\setminus I₂ = \\{x: x ∈ I₁ \\text{ and } x ∉ I₂ \\}
```

The backslash symbol, `\\`, can be used as an alias.

### Input

- `I1` -- first interval
- `I2` -- second interval

### Output

Depending on the position of the intervals, the output is one of the following:

- An `EmptySet`.
- An `Interval`.
- A `UnionSet` of two `Interval` sets.

### Algorithm

Let ``I₁ = [a, b]`` and ``I₂ = [c, d]`` be intervals. Their set difference is
``I₁ \\setminus I₂ = \\{x: x ∈ I₁ \\text{ and } x ∉ I₂ \\}`` and depending on their
position three different results may occur:

- If ``I₁`` and ``I₂`` do not overlap, i.e. if their intersection is empty, then the
  set difference is just ``I₁``.
- Otherwise, let `I₁₂ = I₁ ∩ I₂` and assume that it is not empty, then either
  ``I₁₂`` splits `I₁` into one interval or into two intervals. The latter case
  happens when the inclusion is strict on both ends of ``I₂``.

To check for strict inclusion, we assume that the inclusion is strict and then
check if the resulting intervals that cover `I₁` (one to its left and one to its
right, let them be `Ileft` and `Iright`), obtained by intersection with `I₂`,
are flat or not. Three cases may arise:

- If both `Ileft` and `Iright` are flat then it means that `I₁ = I₂`, then the
  set difference is the empty set.
- If only `Ileft` is flat, then the remaining interval not covered by `I₂` is
  `Iright`. In a similar manner, if only `Iright` is flat, then `Ileft` is returned.
- Finally, if none of the intervals is flat, then `I₂` is strictly contained in
  `I₁` and the set union of `Ileft` and `Iright` is returned.
"""
function difference(I1::IN, I2::IN) where {N, IN<:Interval{N}}
    I12 = intersection(I1, I2)
    if isempty(I12)
        return I1
    else
        Ileft = Interval(min(I1), min(I12))
        Iright = Interval(max(I12), max(I1))

        flat_left = isflat(Ileft)
        flat_right = isflat(Iright)

        if flat_left && flat_right
            return EmptySet{N}(1)
        elseif flat_left && !flat_right
            return Iright
        elseif !flat_left && flat_right
            return Ileft
        else
            return UnionSet(Ileft, Iright)
        end
    end
end

# ============================================
# Set difference between hyperrectangular sets
# ============================================

"""
    difference(X::AbstractHyperrectangle{N}, Y::AbstractHyperrectangle{N}) where {N}

Return the set difference between the given hyperrectangular sets.

### Input

- `X` -- first hyperrectangular set
- `Y` -- second hyperrectangular set

The set difference is defined as:

```math
    X \\setminus Y = \\{x: x ∈ X \\text{ and } x ∉ Y \\}
```

### Output

A `UnionSetArray` consisting of the union of hyperrectangles. Note that this
union is in general not convex.

### Algorithm

This function calls the implementation in `IntervalArithmetic.setdiff`.

### Notes

The backslash symbol, `\\`, can be used as an alias.
"""
function difference(X::AbstractHyperrectangle{N}, Y::AbstractHyperrectangle{N}) where {N}
    Xib = convert(IA.IntervalBox, X)
    Yib = convert(IA.IntervalBox, Y)
    return UnionSetArray(convert.(Hyperrectangle, IA.setdiff(Xib, Yib)))
end


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK

const ArrayType = CLIMA.array_type()

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

include("mms_solution_generated.jl")
include("mms_model.jl")

# initial condition

function run(mpicomm, dim, topl, warpfun, N, timeend, FT, dt)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = warpfun,
                                         )
  dg = DGModel(MMSModel{dim}(),
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

  Q = init_ode_state(dg, FT(0))

  lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, ))
  # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  # Print some end of the simulation information
  engf = norm(Q)
  Qe = init_ode_state(dg, FT(timeend))

  engfe = norm(Qe)
  errf = euclidean_distance(Q, Qe)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  norm(Q - Qe)            = %.16e
  norm(Q - Qe) / norm(Qe) = %.16e
  """ engf engf/eng0 engf-eng0 errf errf / engfe
  errf
end

using Test
let
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 4
  base_num_elem = 4

  expected_result = [1.5606126271800713e-01 5.3315235468489059e-03 2.2572701271271977e-04;
                     2.5754410198969766e-02 1.1781217145186925e-03 6.1752962472852435e-05]
  lvls = integration_testing ? size(expected_result, 2) : 1

    for FT in (Float64,) #Float32)
      result = zeros(FT, lvls)
      for dim = 2:3
        for l = 1:lvls
          if dim == 2
            Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
            brickrange = (range(FT(0); length=Ne[1]+1, stop=1),
                          range(FT(0); length=Ne[2]+1, stop=1))
            topl = BrickTopology(mpicomm, brickrange,
                                 periodicity = (false, false))
            dt = 1e-2 / Ne[1]
            warpfun = (x1, x2, _) -> begin
              (x1 + sin(x1*x2), x2 + sin(2*x1*x2), 0)
            end

          elseif dim == 3
            Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
            brickrange = (range(FT(0); length=Ne[1]+1, stop=1),
                          range(FT(0); length=Ne[2]+1, stop=1),
            range(FT(0); length=Ne[2]+1, stop=1))
            topl = BrickTopology(mpicomm, brickrange,
                                 periodicity = (false, false, false))
            dt = 5e-3 / Ne[1]
            warpfun = (x1, x2, x3) -> begin
              (x1 + (x1-1/2)*cos(2*π*x2*x3)/4,
               x2 + exp(sin(2π*(x1*x2+x3)))/20,
              x3 + x1/4 + x2^2/2 + sin(x1*x2*x3))
            end
          end
          timeend = 1
          nsteps = ceil(Int64, timeend / dt)
          dt = timeend / nsteps

          @info (ArrayType, FT, dim)
          result[l] = run(mpicomm, dim, topl, warpfun,
                          polynomialorder, timeend, FT, dt)
          @test result[l] ≈ FT(expected_result[dim-1, l])
        end
        if integration_testing
          @info begin
            msg = ""
            for l = 1:lvls-1
              rate = log2(result[l]) - log2(result[l+1])
              msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
            end
            msg
          end
        end
      end
    end
end

nothing

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text: println(Chord(:triad))


text:
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end
using Test
import Thermodynamics
const TD = Thermodynamics

import UnPack
import BenchmarkTools

import RootSolvers
const RS = RootSolvers

import CLIMAParameters
const CP = CLIMAParameters

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const param_set = EarthParameterSet()

ArrayType = Array{Float64}
profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
UnPack.@unpack e_int, T, ρ, p, θ_liq_ice, q_tot = profiles

function thermo_state_ρeq_newton()
    ts =
        TD.PhaseEquil_dev_only.(
            param_set,
            ρ,
            e_int,
            q_tot;
            sat_adjust_method = RS.NewtonsMethod,
        )
end

function thermo_state_ρeq_regula_falsi()
    ts =
        TD.PhaseEquil_dev_only.(
            param_set,
            ρ,
            e_int,
            q_tot;
            sat_adjust_method = RS.RegulaFalsiMethod, maxiter = 20,
        )
end

function thermo_state_ρeq_secant()
    ts =
        TD.PhaseEquil_dev_only.(
            param_set,
            ρ,
            e_int,
            q_tot;
            sat_adjust_method = RS.SecantMethod, maxiter = 50,
        )
end

function thermo_state_pθq()
    ts = TD.PhaseEquil_pθq.(param_set, p, θ_liq_ice, q_tot)
end

function thermo_state_pTq()
    ts = TD.PhaseEquil_pTq.(param_set, p, T, q_tot)
end


text:
using LinearAlgebra
using JLD2

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

function timesteptest( cfl = 0.1, testcase = 2, order = 2, verbose = false)

    N = 2 # vector dimension
    IN = Matrix{Float64}(I, N, N)

    if testcase == 1 || testcase == 2
	K0 = [0 1; 1 0]
	S0 = [0 0; 0 0]
    elseif testcase == 0 || testcase == 3
	K0 = [0 0; 0 0]
	S0 = [0 1; -1 0]
    end

    #Final time
    period = 1
    T = 5pi
    omega = 2*pi/period

    lamb = eigvals(K0 .+ S0)
    maxeig = maximum(abs.(lamb))

    #time step
    dt = cfl/maxeig
    nsteps = ceil(Int64,T/dt)
    dt = T/nsteps

    if verbose
        println("Testcase: ", testcase, "  Cfl: ", cfl, " Final time: ", T, " nsteps: ", nsteps )
    end

    # Initial conditions
    u = [1.0; 0.0]
    v = [0.0; 0.0]

    t = 0.0

    # timefunctions and forcing
    if testcase == 1
	# timefunc1(t) = 0.5*(sin(0.5*omega*(t)))^2;
	timefunc1(t) = 0.25*(1.0 - cos(omega*t));
	uforce1(t::Float64) = [0.0; 0.0]
	vforce1(t::Float64) = [0.0; 0.0]

	timefunc =timefunc1
	uforce = uforce1
	vforce = vforce1

    elseif testcase == 0
	timefunc0(t::Float64) = 0.25*(1-sin(omega*t))
	uforce0(t::Float64) = [0.0; 0.0]
	vforce0(t::Float64) = [0.0; 0.0]

	timefunc =timefunc0
	uforce = uforce0
	vforce = vforce0
    elseif testcase == 2
	timefunc2(t::Float64) = 4/T^2 *t*(T-t)
	phi12(t::Float64) = 0.25*(t - sin(omega*t)/omega)
	phidot2(t::Float64) = 0.5*(sin(0.5*omega*(t)))^2
	uforce2(t::Float64) = [(timefunc2(t) - phidot2(t))*sin(phi12(t)); 0.0]
	vforce2(t::Float64) = [0.0; -(timefunc2(t) - phidot2(t)) * cos(phi12(t))]

	timefunc =timefunc2
	uforce = uforce2
	vforce = vforce2
    elseif testcase == 3
	timefunc3(t::Float64) =  4/T^2 *t*(T-t)
	phi13(t::Float64) = 0.25*(t - sin(omega*t)/omega)
	phidot3(t::Float64) = 0.5*(sin(0.5*omega*(t)))^2
	uforce3(t::Float64) = [-phidot3(t)*sin(phi13(t)); timefunc3(t)*cos(phi13(t))]
	vforce3(t::Float64) = [-timefunc3(t)*sin(phi13(t)); phidot3(t)*cos(phi13(t))]

	timefunc =timefunc3
	uforce = uforce3
	vforce = vforce3
    end

    K(t::Float64) = timefunc(t)*K0
    S(t::Float64) = timefunc(t)*S0

    #Create time stepper
    gamma, stages = Juqbox.getgamma(order)
    timestepper = Juqbox.svparams(K,S,IN)

    #Time integration
    usave = u
    vsave = -v
    tsave = t

    # Stormer-Verlet
    start = time()

    for ii in 1:nsteps
	for jj in 1:stages
	    t, u, v, v05 =  Juqbox.step(timestepper,t,u,v,dt*gamma[jj],uforce,vforce)
	    # t, u, v =  Juqbox.magnus(timestepper,t,u,v,dt,uforce,vforce)
	end
	usave = [usave u]
	vsave = [vsave -v]
	tsave = [tsave t]
    end
    elapsed = time() - start

    if verbose
        @show(elapsed)
    end

    # Evaluate analytical solution
    if testcase == 1 || testcase == 2 || testcase==3
#        phi = 0.25.*(tsave - 1/omega.*sin.(omega.*tsave))
        phi = 0.25*(t - 1.0/omega*sin(omega*t))
        cg = cos(phi)
        ce = -im*sin(phi)
    elseif testcase == 0
#        phi = 0.25.*( tsave + 1/omega.*(cos.(omega.*tsave) .- 1) );
        phi = 0.25*( t + 1/omega*(cos(omega*t) - 1.0) );
        cg = cos(phi);
        ce = -sin(phi);
    end

    #compute errors @ final time
    # cg_err = sqrt( (usave[1,end]-real(cg) )^2 + (vsave[1,end]-imag(cg) )^2 );
    # ce_err = sqrt( (usave[2,end]-real(ce) )^2 + (vsave[2,end]-imag(ce) )^2 );

    cg_err = sqrt( (u[1]-real(cg) )^2 + (v[1]+imag(cg) )^2 ); # negative imaginary part in v
    ce_err = sqrt( (u[2]-real(ce) )^2 + (v[2]+imag(ce) )^2 );

    if verbose
        println("cg-err = " , cg_err, " ce-err = ", ce_err)
    end

    return t,u,v,dt,cg_err,ce_err
end

function timestep_convergence( errFileName:: String, writeFile:: Bool=false)
    CFL_vec = 10.0.^(-1.0:-0.5:-2.0)

    ntests = 4
    err_mat  = zeros(length(CFL_vec),2,ntests)
    order    = 2

    verbose = false

    for j = 1:ntests
	for i = 1:length(CFL_vec)
	    t,u,v,dt,cg_err,ce_err = timesteptest(CFL_vec[i], j-1, order, verbose)
	    err_mat[i,1,j] = cg_err
	    err_mat[i,2,j] = ce_err
	end
    end

    if writeFile
        @save errFileName err_mat
        println("Saved final errors on file: ", errFileName)
        return true
    else
        # compare results with those on the reference file
        err_mat0 = copy(err_mat)
        @load errFileName err_mat
        # check sizes
        if length(err_mat0) != length(err_mat)
            printf("timestep_convergence: size missmatch with reference solution!")
            return false
        end
        # compare
        max_diff = maximum(abs.(err_mat - err_mat0))
        if verbose
            println("Max abs diff = ", max_diff)
        end
        success = (max_diff <= 1e-13)
        return success

        # for j = 1:ntests
	#     for i = 1:length(CFL_vec)
        #         println("test #", j, " cfl = ", CFL_vec[i], " diff (cg) = ", abs(err_mat[i,1,j] - err_mat0[i,1,j]), " diff (ce) = ", abs(err_mat[i,2,j] - err_mat0[i,2,j]) )
        #     end
        # end
    end

end