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

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

# ╔═╡ 0507f8ed-8b64-48af-a6c4-ff7c4211b9e9
begin
	using Pkg
	Pkg.activate(".")
end

# ╔═╡ 3310d2f7-450a-4cd1-9c3a-46d02d23a7c6
using Revise

# ╔═╡ d27b2d72-d264-11eb-0be5-13dcacfd2adc
using DeconvOptim, TestImages, ImageShow, Plots, LinearAlgebra, IndexFunArrays, Noise, FourierTools, SpecialFunctions, FFTW, LaTeXStrings, PlutoUI, Images, Tullio

# ╔═╡ 952b251c-207b-4412-b6e1-268fce1647d9
begin
	img = Float32.(testimage("fabio_gray"));
	img_1D = img[:, 200]
end;

# ╔═╡ fa8cd9c9-a2fd-495c-8d22-ada7bb9c39f6
otf(x, Δx=1) = begin
	x = abs(x)
	if x <= Δx
		SpecialFunctions.jinc(x * Δx *(1-x/Δx)) .* 2 / π * (acos(x/Δx) - x/Δx * sqrt(1-(x/Δx)^2))
	else
		zero(x)
	end
end

# ╔═╡ 9a07bc88-be76-4531-bc91-df0d20c3221c
begin
	x = range(-1.5, 1.5, length=size(img, 1))
	freqs = fftshift(fftfreq(size(img_1D, 1), 1))
	psf = Float32.(DeconvOptim.generate_psf(size(img), 20))
	psf_1D = psf[1, :]
	psf_1D ./= sum(psf_1D)
	otf_1D = abs.(ffts(psf_1D))
end;

# ╔═╡ 24fc86da-0b7c-493c-8977-2a19ef6dc133
img_n = Float32.(poisson(DeconvOptim.conv(img, psf), 1000));

# ╔═╡ 3bda922d-552b-42ec-9055-33a141b5841a
blur(x, otf=otf_1D) = iffts(ffts(x) .* otf)

# ╔═╡ 308e2562-5a20-4082-8a6d-9fb135738c49
md"
Mathematically:

$(S * \text{PSF})(\mathbf r) = \int_{-\infty}^{\infty} S(\mathbf r - \mathbf x) \cdot \text{PSF}(\mathbf x) \, \mathrm d \mathbf x$
"

# ╔═╡ b542187a-3ae7-4430-a81d-f968d9000427
img_blurry = DeconvOptim.conv(img, psf);

# ╔═╡ 88d83629-f1b2-4d69-928f-d4acfbc76b70
reg_1D = TV(num_dims=1);

# ╔═╡ d6d0f436-48fd-4e8a-863c-3e9d3850cc91
begin
	reg_tik = Tikhonov()
	reg_TV = TV()
	reg_GR = DeconvOptim.GR()
end

# ╔═╡ 061faf49-662c-4508-83b4-ddcf0970ed0d
md"### DeconvOptim.jl: Microscopy Image Deconvolution
"

# ╔═╡ 23829201-9756-4ef8-90c9-3917b761fe4b
load("../docs/src/assets/logo.png")

# ╔═╡ 30f21bb8-6d09-4fce-9d2a-568bfaf3ff7a
md"
* **Felix Wechsler:** Master Student at the Leibniz Institute of Photonic Technology in  Jena, Germany
* https://github.com/roflmaostc/DeconvOptim.jl
* `]add DeconvOptim`
"

# ╔═╡ 7a0a44c9-fb07-44e6-9a8b-8f720b84e6f6
md"""### Image Convolution

* Typical description of isotropic blur of an image
* The blurring kernel describes blur
    * In optics/microscopy a finite sized dot called Point Spread Function (**PSF**)
    * often a Gaussian function used in image processing
    * cigarre shaped object for motion blur

Discrete version:

$(S * \text{PSF})[i] = \sum_{m} S[i-m] \cdot \text{PSF}[m]$
"""

# ╔═╡ 49686d9a-1683-428f-83ba-a9131c2ad432
[Gray.(img) Gray.(DeconvOptim.conv(img, psf))]

# ╔═╡ 491ec7cb-1663-4acf-b81f-9acafba2b63d
md"## Convolution Theorem

$(S * \text{PSF})(\mathbf r) = \mathcal{F}^{-1}\bigg[ \mathcal{F}[S] \cdot \mathcal{F}[\text{PSF}]  \bigg]$

* we can express the convolution with a Fast Fourier Transform (FFT) which only takes $\mathcal O(N \log(N))$ operations

* For large kernels (especially in 3D), sliding kernels are slower

*  $\mathcal{F}[\text{PSF}]$ is called the $\text{OTF}$

"

# ╔═╡ b6f9e42e-5a23-40ad-9e73-8f02da48f69c
md"### Optical System act as low pass filter

*  $\text{OTF}$ shows the frequency throughput

"

# ╔═╡ e31f4704-5bce-4c8d-b3a7-873a3460d9f8
plot(x, otf.(x), xlabel="frequency / maximum frequency", ylabel="contrast")

# ╔═╡ dc8cf4e2-a87c-47ea-8247-3ae772852241
md"## Frequency spectrum of blurred sample $Y(\mathbf r)$

Blurred sample:
$Y(\mathbf r) = (S * \text{PSF})(\mathbf r)$
"

# ╔═╡ f8033d13-faee-41f5-bdd6-ae8721e8b8a8
begin
	plot(freqs, abs.(ffts(img_blurry)[:, 128]), yaxis=:log, ylabel="real part of FFT output in AU", xlabel="frequency in 1/px", ylims=(1e-4, 1e2), label="blurred")
	plot!(freqs, abs.(ffts(img)[:, 128]), ylabel="abs of FFT output in AU",
				xlabel="frequency in 1/px", yaxis=:log,
				ylims=(1e-4, 1e4), label="ground truth")
	#plot!(freqs, abs.(ffts(DeconvOptim.conv(img_1D, psf_1D))))
end

# ╔═╡ fe7e0292-31f6-43b5-83a5-a38698a87563
md"## Deconvolution Pipeline

* based on:
    * Zygote.jl
    * Optim.jl
    * Tullio.jl
    * CUDA.jl
"

# ╔═╡ e9ef5ba4-56c0-4595-bc28-e04882f44a9a
load("../docs/src/assets/tex/pipeline.png")

# ╔═╡ 90e5708c-05a4-46e4-b1e4-9a61c96dae32
TV_by_hand(x) = @tullio r = sqrt(1f-8 + abs2(x[i, j] - x[i+1, j]) +
						 	     abs2(x[i, j] - x[i, j+1]))

# ╔═╡ 03139ac5-3525-4fad-abf1-84421492b763
DeconvOptim.generate_TV(4, [1,2, 3], [1,1, 1], 1, 0)[1]

# ╔═╡ 4e84e739-9c59-4939-8b04-aec7dc069d67
md"
### Deconvolve with DeconvOptim.jl

"

# ╔═╡ 9d9a5da2-14df-46e6-b7e6-5a33aade1754
@bind reg_list2 Select(["1" => ("Tikhonov"), "2" => ("Total Variation TV"), "3" => ("Good's Roughness GR")])

# ╔═╡ 5ae1a6a3-4505-4123-9f1e-8a1d4ac0b4e1
reg = [reg_tik, reg_TV, reg_GR][parse(Int, reg_list2)]

# ╔═╡ 15d54e5e-64f4-4a1d-8cc8-9334b2e3784f
md"
iterations =
$(@bind iter Slider(0:50, show_value=true))

λ = $(@bind λ Slider(0:0.001:0.3, show_value=true))

regularizer = $(@bind reg_bool CheckBox())"

# ╔═╡ 803368e6-53fd-4413-b3f5-ffe46ee8983e
img_deconv, res_img = deconvolution(img_blurry, psf, regularizer=reg_bool ? reg : nothing, iterations=iter, λ=λ);

# ╔═╡ e58e1f63-2c81-48fb-866a-4bb70bd428a6
Gray.(img_deconv)

# ╔═╡ 438a6639-bd35-464f-a81d-d98eb65e006e
res_1D, o = deconvolution(real(blur(img_1D)), psf_1D, iterations=iter, regularizer=reg_1D, λ=0.01);

# ╔═╡ 7a16df73-95ad-47f5-907c-6fd23c6000cf
[Gray.(img) Gray.(img_blurry) Gray.(img_deconv)]

# ╔═╡ 5c8030b7-8805-4771-9329-23abb2744544
begin
	plot(freqs, abs.(ffts(img_blurry)[:, 128]), yaxis=:log, ylabel="abs of FFT output in AU", xlabel="frequency in 1/px", ylims=(1e-4, 1e2), label="blurred")
	plot!(freqs, abs.(ffts(img)[:, 128]),  yaxis=:log,
				ylims=(1e-4, 1e4), label="ground truth")
	plot!(freqs, abs.(ffts(img_deconv)[:, 128]), label="deconvolved image")
end

# ╔═╡ ff9af79c-f06c-42ac-9c8d-6f09d2ff4056
[Gray.(img_1D); Gray.(res_1D)];

# ╔═╡ d8aee845-e922-41da-a17e-37ffa3e692f0
md"### Real Microscopy Data"

# ╔═╡ fc93fa3d-8599-463b-b9e6-8043b90e9d63
load("figures/real_data_large.png")

# ╔═╡ 97c45bfd-d1ef-49ad-908d-7360c03b0170
md"Image taken from [DeconvolutionLab2](http://bigwww.epfl.ch/deconvolution/deconvolutionlab2/)."

# ╔═╡ bf37664a-f726-46b6-9592-da419165af91
md"## Conclusion - DeconvOptim.jl
"

# ╔═╡ f8117250-bac8-43a1-aae4-5c8bab3a522d
[Gray.(ones(130, 012)) load("../docs/src/assets/logo.png")]

# ╔═╡ b5a70276-e9b9-46e8-8c67-0cad2cfa19da
md"* Flexible Image Deconvolution Software
* N-dimensional signal deconvolution
* Works both on CPU and GPUs
    * GPUs usually 5-15x speed improvement
"

# ╔═╡ Cell order:
# ╠═3310d2f7-450a-4cd1-9c3a-46d02d23a7c6
# ╠═0507f8ed-8b64-48af-a6c4-ff7c4211b9e9
# ╠═d27b2d72-d264-11eb-0be5-13dcacfd2adc
# ╠═952b251c-207b-4412-b6e1-268fce1647d9
# ╠═24fc86da-0b7c-493c-8977-2a19ef6dc133
# ╠═3bda922d-552b-42ec-9055-33a141b5841a
# ╟─fa8cd9c9-a2fd-495c-8d22-ada7bb9c39f6
# ╠═9a07bc88-be76-4531-bc91-df0d20c3221c
# ╟─308e2562-5a20-4082-8a6d-9fb135738c49
# ╠═b542187a-3ae7-4430-a81d-f968d9000427
# ╠═88d83629-f1b2-4d69-928f-d4acfbc76b70
# ╠═e58e1f63-2c81-48fb-866a-4bb70bd428a6
# ╠═d6d0f436-48fd-4e8a-863c-3e9d3850cc91
# ╠═803368e6-53fd-4413-b3f5-ffe46ee8983e
# ╠═5ae1a6a3-4505-4123-9f1e-8a1d4ac0b4e1
# ╠═438a6639-bd35-464f-a81d-d98eb65e006e
# ╟─061faf49-662c-4508-83b4-ddcf0970ed0d
# ╟─23829201-9756-4ef8-90c9-3917b761fe4b
# ╟─30f21bb8-6d09-4fce-9d2a-568bfaf3ff7a
# ╟─7a0a44c9-fb07-44e6-9a8b-8f720b84e6f6
# ╟─49686d9a-1683-428f-83ba-a9131c2ad432
# ╟─491ec7cb-1663-4acf-b81f-9acafba2b63d
# ╟─b6f9e42e-5a23-40ad-9e73-8f02da48f69c
# ╟─e31f4704-5bce-4c8d-b3a7-873a3460d9f8
# ╟─dc8cf4e2-a87c-47ea-8247-3ae772852241
# ╟─f8033d13-faee-41f5-bdd6-ae8721e8b8a8
# ╟─fe7e0292-31f6-43b5-83a5-a38698a87563
# ╟─e9ef5ba4-56c0-4595-bc28-e04882f44a9a
# ╠═90e5708c-05a4-46e4-b1e4-9a61c96dae32
# ╠═03139ac5-3525-4fad-abf1-84421492b763
# ╟─4e84e739-9c59-4939-8b04-aec7dc069d67
# ╟─9d9a5da2-14df-46e6-b7e6-5a33aade1754
# ╟─15d54e5e-64f4-4a1d-8cc8-9334b2e3784f
# ╟─7a16df73-95ad-47f5-907c-6fd23c6000cf
# ╟─5c8030b7-8805-4771-9329-23abb2744544
# ╟─ff9af79c-f06c-42ac-9c8d-6f09d2ff4056
# ╟─d8aee845-e922-41da-a17e-37ffa3e692f0
# ╟─fc93fa3d-8599-463b-b9e6-8043b90e9d63
# ╟─97c45bfd-d1ef-49ad-908d-7360c03b0170
# ╟─bf37664a-f726-46b6-9592-da419165af91
# ╟─f8117250-bac8-43a1-aae4-5c8bab3a522d
# ╟─b5a70276-e9b9-46e8-8c67-0cad2cfa19da


text:
# Note: The "EHR Launch" workflow is also known as the "Embedded Launch" workflow.

function _get_endpoints_from_extension(extension)
    endpoints = Dict{String, String}()
    for e ∈ extension
        value_uri = e.valueUri
        url = e.url
        endpoints[url] = value_uri
    end
    return endpoints
end

"""
    provider_ehr_launch(config::ProviderEHRLaunchConfig, uri_string::AbstractString; kwargs...)

## Optional Keyword Arguments:
- `scope::AbstractString`. Default value: `"launch"`.
"""
function provider_ehr_launch(config::ProviderEHRLaunchConfig,
                             uri_string::AbstractString;
                             kwargs...)
    uri = URIs.URI(uri_string)
    return provider_ehr_launch(config, uri; kwargs...)
end

"""
    provider_ehr_launch(config::ProviderEHRLaunchConfig, uri::URIs.URI; kwargs...)

## Optional Keyword Arguments:
- `scope::AbstractString`. Default value: `"launch"`.
"""
function provider_ehr_launch(config::ProviderEHRLaunchConfig,
                             uri::URIs.URI;
                             scope::AbstractString = "launch")
    queryparams = URIs.queryparams(uri)
    iss          = queryparams["iss"]
    launch_token = queryparams["launch"]
    iss_metadata_endpoint = "$(iss)/Metadata"
    metadata_response = HTTP.request(
        "GET",
        iss_metadata_endpoint;
        headers = Dict(
            "Accept" => "application/fhir+json",
            "Epic-Client-ID" => config.client_id,
        ),
    )
    metadata_response_json = JSON3.read(String(metadata_response.body))
    extension = metadata_response_json.rest[1].security.extension[1].extension
    endpoints = _get_endpoints_from_extension(extension)
    authorize_endpoint = endpoints["authorize"]
    token_endpoint = endpoints["token"]
    authorize_uri_with_querystring_params = URIs.URI(
        URIs.URI(authorize_endpoint);
        query = Dict(
            "client_id" => config.client_id,
            "launch" => launch_token,
            "redirect_uri" => config.redirect_uri,
            "response_type" => "code",
            "scope" => scope,
        )
    )
    authorize_response = HTTP.request(
        "GET",
        authorize_uri_with_querystring_params;
        redirect = false,
    )
    headers = Dict(authorize_response.headers)
    location = headers["Location"]
    location_uri = URIs.URI(location)
    location_queryparams = URIs.queryparams(location_uri)
    authorization_code = location_queryparams["code"]

    authz_code_info = _AuthorizationCodeInformation(;
        authorization_code = authorization_code,
        client_id          = config.client_id,
        redirect_uri       = config.redirect_uri,
        token_endpoint     = token_endpoint,
    )

    access_token_info = authorization_code_to_access_token(authz_code_info)

    launch_token_is_jwt, launch_token_jwt_decoded = try_decode_jwt(launch_token)

    ehr_launch_result = ProviderEHRLaunchResult(;
        authorization_code             = access_token_info.authorization_code,
        authorization_code_is_jwt      = access_token_info.authorization_code_is_jwt,
        authorization_code_jwt_decoded = access_token_info.authorization_code_jwt_decoded,
        access_token                   = access_token_info.access_token,
        access_token_is_jwt            = access_token_info.access_token_is_jwt,
        access_token_jwt_decoded       = access_token_info.access_token_jwt_decoded,
        access_token_response          = access_token_info.access_token_response,
        launch_token                   = launch_token,
        launch_token_is_jwt            = launch_token_is_jwt,
        launch_token_jwt_decoded       = launch_token_jwt_decoded,
    )

    return ehr_launch_result
end

ext: jl


text:
using StatisticalNonlocality, Distributions, Random, LinearAlgebra
import StatisticalNonlocality: ou_transition_matrix
import StatisticalNonlocality: uniform_phase
import Distributions: Uniform

# Set random seed for reproducibility
Random.seed!(10001)

# simulate a continuous time markov process
T = ou_transition_matrix(4)
T = uniform_phase(4)

γ = 0.01
eT = exp(γ * T) # get the transition probabilities
# each column gives the transition probability
# column i, means, given that I am in state i, each row j gives the probability to transition to state j

ceT = cumsum(eT, dims = 1)

# Define the jump map
function next_state(current_state_index::Int, cT)
    vcT = view(cT, :, current_state_index)
    u = rand(Uniform(0, 1))
    # choose a random uniform variable and decide next state
    # depending on where one lies on the line with respect to the probability
    # of being found between given probabilities
    for i in eachindex(vcT)
        if u < vcT[i]
            return i
        end
    end

    return i
end

M = 1000
markov_chain = zeros(Int, M)
markov_chain[1] = 4
for i = 2:M
    markov_chain[i] = next_state(markov_chain[i-1], ceT)
end

using GLMakie
fig2 = scatter(markov_chain)

fig = hist(markov_chain)

# E[u(x⃗, t + τ) * u(x⃗, t) ] - E[u(x⃗,t)]E[u(x⃗,t+τ)] =  ∑ₘ( Pₘ(t) uₘ(x⃗) * ∑ₙ( exptT[n, m] uₙ(x⃗) ) ) - (∑ₘ Pₘ(t) uₘ) (∑ₘ Pₘ(t + τ) uₘ)
# under more assumptions
# ∑ₘₙ ( Pₘ(t) ( exptT[n, m]uₘ(x⃗) uₙ(x⃗) ) ) = u⃗' exp(T) * diag(P) * u⃗

# anti-deriv  T⁻¹ ( I - exp(t T) ) *

exp(T) * ones(4) ./ 4

exp(3 * T) * [1 -1 1 -1]'
##
basic(α) = [α 1-α; 1-α α]

# {-1, 1}
α = 2 / 3
tmp = basic(α)

# E[X[n]X[n+d]] = (2α - 1)ᵈ
# α = 2/3 => 3⁻ᵈ
sum(tmp * [1, 0] * Diagonal([-1 1])) + sum(tmp * [0, 1] * Diagonal([-1 1]))

# 1/2 * {-1} * (-1 * α + 1 * (1-α)) + 1/2 * {1} * (-1 * (1-α) + 1 * α)
1 / 2 * (-1) * (-1 * α + 1 * (1 - α)) + 1 / 2 * (1) * (-1 * (1 - α) + 1 * α)
[-1 1] * tmp^100 * [-1 / 2, 1 / 2]

[-1 , 1 ]' * tmp * [-1 / 2, 1 / 2]
λ, V = eigen(tmp)
Λ = Diagonal(λ)
V⁻¹ = inv(V)
tmpmat = zeros(size(tmp))
for i = 1:length(λ)
    tmpmat += λ[i] * V[:, i] * transpose(V⁻¹[i, :])
end

##
V⁻¹ = inv(V)
V * Λ * V⁻¹ - tmp
Diagonal([0, 1])
V * Λ * V⁻¹
V * Diagonal([0, 1]) * V⁻¹
V[:, 2] * V⁻¹[2, :]'

##
λ, V = eigen(T)
Λ = Diagonal(λ)
V⁻¹ = inv(V)
tmp = zeros(length(λ))
tmp[end] = 1
Λ∞ = Diagonal(tmp)
V * Λ * V⁻¹ - T
real.(V * Λ∞ * V⁻¹)
real.(V[:, end] * transpose(V⁻¹[end, :]))

tmpmat = zeros(ComplexF64, size(T))

for i = 1:length(λ)
    tmpmat += λ[i] * V[:, i] * transpose(V⁻¹[i, :])
end

T̃ = real.(tmpmat)
norm(T̃ - T)

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
using LinearAlgebra, Statistics
using Distributions, QuasiMonteCarlo, Random

n = 10000
d = 10

# synthetic d-dimensional least squares problem
wgen = rand(d);
w0 = zeros(d);

# Quasi-random distributions
qrngs = Dict([("IID Uniform"     , UniformSample()),
              ("Sobol"           , SobolSample()),])

qrng_names = collect(keys(qrngs))

# data distribution
# x ~ N(0,I_d)
# y ~ N(x'wopt, 1)
# goal is to learn w, which is to minimize E[0.5*(x'w - y)^2]

# online mode:  draw samples from this integral at every SGD iteration
# offline mode: draw n samples iid uniformly from this integral
#               to form a training set,
#               then draw one sample (possibly following a low-discrepancy sequence)
#               from the training set at every SGD iteration

function gen_xy(q::Vector{Float64})
    xi = quantile(Normal(), q[1:d])
    yi = quantile(Normal(xi'wgen, 1), q[end])
    return (xi, yi)
end

function gen_syn_data(num_examples::Int64, qrng_name::String, dim::Int64, seed::Int64)
    Random.seed!(seed)

    lb = zeros(dim)
    ub = ones(dim)

    distrib = qrngs[qrng_name]
    qmc_samples = QuasiMonteCarlo.sample(num_examples, lb, ub, distrib)

    syn_data = []
    for i in 1:num_examples
        q = mod.(qmc_samples[:,i], 1)
        push!(syn_data, gen_xy(q))

    end
    return syn_data
end

function sgd(T::Int64, w0::Vector{Float64}, wopt::Vector{Float64},
	samples::Vector{Any}, ordering::Vector{Int64})
	subopt = zeros(T)
	wk = w0
	# diminishing step size optimized for IID Uniform
	alpha = 1/(d + 2 + d/norm(wgen)^2);
	for t in 1:T
		(xi, yi) = samples[ordering[t]]
		alpha = 1/(1/alpha + (1 - alpha*(d+2))/(1 - alpha))
		wk -= alpha * (xi'wk - yi) * xi
		subopt[t] = norm(wk - wopt).^2
	end
	return subopt
end

max_stars_repo_name: EugeneLYC/qmc-ordering

max_issues_repo_name: EugeneLYC/qmc-ordering

max_forks_repo_name: EugeneLYC/qmc-ordering


text:
export colorBlack1, colorBlack2, colorBlack3, colorBlack4, colorBlack5,
    lineStyle1, lineStyle2, lineStyle3, lineStyle4,
    lineWidth1, lineWidth2, lineWidth3, lineWidth4,
    marker1, marker2, marker3, marker4,
    markerSize1, markerSize2, markerSize3, markerSize4,
    legendFontSize, arrowaxes, removeaxes, lengthdimension

const colorBlack1 = "0"
const colorBlack2 = "0.25"
const colorBlack3 = "0.5"
const colorBlack4 = "0.6"
const colorBlack5 = "0.85"

const lineStyle1 = "-"
const lineStyle2 = linestyle = (0, (4, 1))
const lineStyle3 = linestyle = (0, (2, 0.81))
const lineStyle4 = linestyle = (0, (8, 1))

const lineWidth1 = 1
const lineWidth2 = 1
const lineWidth3 = 1.5
const lineWidth4 = 1

const marker1 = "D"
const marker2 = "o"
const marker3 = "s"
const marker4 = "*"

const markerSize1 = 4
const markerSize2 = 4
const markerSize3 = 4
const markerSize4 = 5

const legendFontSize  =  10

@doc raw"""
# Function call

```
arrowaxes(fig=gcf(), ax=gca(); xmin, xmax, ymin, ymax,
    xa=0, ya=0, xlabel="", ylabel="",
    xneg = false, yneg = false,
    color="black", backgroundcolor="none", axisoverhang = 0.18,
    linewidth = 0.75, headwidth = 0.06, headlength = 0.09, overhang = 0.1,
    labelsep = 0.06, left=0.2, right=0.85, bottom=0.20, top=0.85)
```

# Description

Creates a plot with a horizontal and a vertical axis instead of a frame.

# Variables

`fig` Figure handle; by default the figure handle of a actual figure is used

`xmin` Minimum of horizontal axis; default value = actual scaling

`xmax` Maximum of horizontal axis; default value = actual scaling

`ymin` Minimum of vertical axis; default value = actual scaling

`ymax` Maximum of vertical axis; default value = actual scaling

`ax` Axes handle; by default the axes handle of a actual figure is used

`xa` Horizontal location of the vertical axis; default value = 0

`ya` Vertical location of the horizontal axis; default value = 0

`xlabel` Label of x-axis; default value = ""

`ylabel` Label of y-axis; default value = ""

`xneg` If true, the horizontal arrow is drawn from right to left;
default value = `false`

`yneg` If true, the vertical arrow is drawn from top to bottom;
default value = `false`

`color` Color of the axes; default value = "black"

`backgroundcolor` Background color of the label; if labelrsep is equal to 0, the
background color "white" can be used; default value = "none"

`axisoverhang` Overhang of the axis, relative to plot area;
default value = 0.18 (18% of plot width)

`linewidth` Line width of the axes; default value = 0.75

`headwidth` Width of head, relative to plot area;
default value = 0.045 (4.5% of horizontal plot width)

`headlength` Length of head, relative to plot area;
default value = 0.07 (7% of horizontal plot width)

`overhang` Overhang of arrow head; default value = 0.0

`labelsep` Location of labels from axis, relative to plot area;
default value = 0.05 (5% of plot width)

`left` Left side of the figure; default value = 0.2

`right` Right side of the figure; default value = 0.85

`bottom` Bottom side of the figure; default value = 0.2

`top` Top side of the figure; default value = 0.85

# Examples

Copy and paste the following code:

```julia
using Unitful,Unitful.DefaultSymbols,ElectricalEngineering,PyPlot
figure(figsize=(3.3, 2.5))
rc("text", usetex=true); rc("font", family="serif")
x=collect(0.0:0.1:5.0); y=exp.(sin.(x));
plot(x,y,color=colorBlack1, linewidth=lineWidth1, linestyle=lineStyle1)
xlim(0,5); ylim(0,3); arrowaxes(xlabel=L"$x$",ylabel=L"$y$")
# save3fig("arrowaxes", crop=true)
```
"""
function arrowaxes(fig=gcf(), ax=gca();
    xmin=ax.get_xlim()[1], xmax=ax.get_xlim()[2],
    ymin=ax.get_ylim()[1], ymax=ax.get_ylim()[2],
    xa=0, ya=0, xlabel="", ylabel="",
    xneg = false, yneg = false,
    color="black", backgroundcolor="none", axisoverhang = 0.18,
    linewidth = 0.6, headwidth = 5, headlength = 10, overhang = 0,
    labelsep = 0.05,
    left=0.2, right=0.85, bottom=0.20, top=0.85)
    # The basic idea of this implementation is taken from:
    # http://www.yueshen.me/2015/1011.html
    # https://stackoverflow.com/questions/33737736/matplotlib-axis-arrow-tip/44138298#44138298

    # Create enough space around the plot area of fit in arrows
    fig.subplots_adjust(left=left, right=right, top=top, bottom=bottom)
    # fig[:subplots_adjust](left=left, right=right, top=top, bottom=bottom)
    ax.spines["right"].set_visible(false)
    ax.spines["top"].set_visible(false)
    ax.yaxis.set_ticks_position("left") # Only show vertical ticks on left
    ax.xaxis.set_ticks_position("bottom") # Only show horizontal ticks on bottom
    # Set position of horizontal axis:
    ax.spines["bottom"].set_position(("data",ya))
    # Remove horizontal axis, since it will be overwritten by arrow
    # ax[:spines]["bottom"][:set_color]("none")
    ax.spines["bottom"].set_color("none")
    # Set position of vertical axis
    ax.spines["left"].set_position(("data",xa))
    # Remove vertical axis since it will be overwritten by arrow
    ax.spines["left"].set_color("none")

    # Axis limits
    (xMIN,xMAX) = ax.get_xlim() # Horizontal limits
    dx = xMAX-xMIN # Difference of horizontal limits
    (yMIN,yMAX) = ax.get_ylim() # Vertical limits
    dy = yMAX-yMIN # Difference of vertical limits
    # Get width and height of axes object to compute
    # matching arrowhead length and width
    dps = fig.dpi_scale_trans.inverted()
    bbox = ax.get_window_extent().transformed(dps)
    width = bbox.width # width of bbox
    height = bbox.height # height of bbox
    # Physical arrow dimensions
    hwx = headwidth*dy # Width of x-axis arrow head
    hlx = headlength*dx # Length of x-axis arrow length
    hwy = headwidth*dx * height/width # Width of y-axis arrow head
    hly = headlength*dy * width/height # Length of y-axis arrow length

    xbeg = xneg ? xmax : xmin
    xend = xneg ? xmin-axisoverhang*dx : xmax+axisoverhang*dx
    ybeg = yneg ? ymax : ymin
    yend = yneg ? ymin-axisoverhang*dy : ymax+axisoverhang*dy

    # Set new limits
    xlim(xbeg,xend)
    ylim(ybeg,yend)
    # Horizontal arrow created by annotate instead of axis.arrow
    annotate("", xy=(xend, ya),
        xytext=(xbeg, ya), xycoords="data",
        arrowprops=Dict("edgecolor"=>color,
            "facecolor"=>color,
            "width"=>0.2,
            "linewidth"=>linewidth,
            "linestyle"=>"-",
            "headlength"=>headlength,
            "headwidth"=>headwidth,
            "joinstyle"=>"round"),
        annotation_clip=false)
    # Vertical arrow created by annotate instead of axis.arrow
    annotate("", xy=(xa,yend),
        xytext=(xa,ybeg), xycoords="data",
        arrowprops=Dict("edgecolor"=>color,
            "facecolor"=>color,
            "width"=>0.2,
            "linewidth"=>linewidth,
            "linestyle"=>"-",
            "headlength"=>headlength,
            "headwidth"=>headwidth,
            "joinstyle"=>"round"),
        annotation_clip=false)
    # https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.text.html#matplotlib.axes.Axes.text
    # Horizontal label
    ha = xneg ? "left" : "right"
    ax.text(xend, ya+dy*labelsep, xlabel,
        ha=ha, va="bottom", backgroundcolor=backgroundcolor)
    # Vertical label
    va = yneg ? "bottom" : "top"
    ax.text(xa+dx*labelsep, yend, ylabel,
        ha="left",va=va, backgroundcolor=backgroundcolor)
end

"""
# Function call

`removeaxes(ax=gca())`

# Description

Removes the axis of the actual plot

# Variables

`ax` Axes handle; by default the axes handle of a actual figure is used
"""
function removeaxes(ax=gca())
    #
    setproperty!(ax,"set_axis_off",true)
    xticks([])
    yticks([])
    ax.spines["right"].set_visible(false)
    ax.spines["top"].set_visible(false)
    ax.spines["bottom"].set_visible(false)
    ax.spines["left"].set_visible(false)
end

@doc raw"""
# Function call

```
lengthdimension(x1=0, y1=0, x2=1, y2=1;
    label = "", labeltsep = 0.1, labelrsep=0.5, labelrelrot=false,
    labelrelangle=0, ha = "center", va = "center",
    color="black", backgroundcolor = "none",
    arrowstyle1 = "<|-", arrowstyle2 = "-|>", linewidth=0.6, linestyle = "-",
    width=0.2, headlength=5, headwidth=2.5,
    par=0, paroverhang = 0.02, parcolor = "black",
    parlinewidth=0.6, parlinestyle = "-")

```

# Description

This function draws length dimension with label based on rectangular x- and
y-coordinates. The length dimension arrow can be shifted parallel to the
specified  coordinates by means the single parameter `par`. Additionally,
auxiliary lines  from the specified coordinates to the length dimension are
drawn.

# Variables

`x1` x-component of the begin of the length dimension; default value = 0

`y1` y-component of the begin of the length dimension; default value = 0

`x2` x-component of the end of the length dimension; default value = 1

`y2` y-component of the end of the length dimension; default value = 1

`label` Label of the angle; default value =""

`labelphisep` Angular separation of the label with respect to the arc; if
`labelphisep == 0`, the label is located at angle `phi1` and if `labelphisep ==
1`, the label is located at angle `phi2`; default value = 0.5, right in the
middle between `phi1` and `phi2`

labelrsep` Radial per unit location of label (in direction of the phasor):
`labelrsep = 0` locates the label right on the arc. A positive value locates the
label outside the arc, a negative value locates the label inside the arc;
default value = 0.1

`labelrelrot` Relative rotation of label; if `labelrelrot == false` (default
value) then the label is not rotated relative to the center of the arc;
otherwise the label is rotated relative to the  angle `labelrelangle`

`labelrelangle` Relative angle of label with respect to center of the arc; this
angle is only applied, it `labelrelrot == true`; this angle indicates the
relative orientation of the label with respect to the center of the arc;
default value = 0

`ha` Horizontal alignment of label; actually represents the radial alignment of
the label; default value = "center"

`va` Vertical alignment of label; actually represents the tangential alignment
of the label; default value = "center"

`color` Color of the arc; default value = "black"

`backgroundcolor` Background color of the label; if labelrsep is equal to 0, the
background color "white" can be used; default value = "none"

`arrowstyle1` Arrow style of the begin of the arc; default value = "."; valid
strings are:

- `.` dot marker
- `<|-` arrow

`arrowstyle2` Arrow style of the end of the arc; default value = "-|>"; valid
strings are:

- `.` dot marker
- `-|>` arrow

`headlength` Length of arrow head; default value = 5

`headwidth` Width of arrow head; default value = 2.5

`par` In order to be able to draw the length dimension parallel to the specified
coordiantes `x1`, `y1`, `x2`, `y2`, `par` is used to specify the per unit
tangential shift (offset) of the length dimension; default value = 0 (no shift)

`paroverhang` The auxiliary lines, drawn in case of `par != 0`, show the absolute
overhang `paroverhang`; default value = 0.02

`parcolor` Color of the auxiliary lines; default value = "black"

`parlinewidth` Line width of the auxiliary line; default value = 0.6

`parlinestyle` Line style of the auxiliary line; default value = "-"

# Example

Copy and paste code:

```julia
using Unitful, Unitful.DefaultSymbols, PyPlot, ElectricalEngineering
figure(figsize=(3.3, 2.5))
rc("text", usetex=true); rc("font", family="serif")
t1=0.2;t2=0.3;ymax=1
t=[0,t1,t1+t2,2*t1+t2,2*(t1+t2),3*t1+2*t2]
y=[0,ymax,0,ymax,0,ymax]
step(t,y,linewidth=1,color="black")
grid(true); xlim(0,1), ylim(-0.1,1.3);
xlabel(L"$t$\,(s)"); ylabel(L"$y$")
lengthdimension(0,0.2,t1,0.2,label=L"$t_1$",labeltsep=0,backgroundcolor="white")
lengthdimension(t1,0.2,t1+t2,0.2,label=L"$t_2$",labeltsep=0,backgroundcolor="white")
lengthdimension(t1,ymax,2*t1+t2,ymax,label=L"$T$",par=0.1,labeltsep=0.05,labelrsep=0.7)
# save3fig("lengthdimension",crop=true)
```
"""
function lengthdimension(
    x1 = 0,
    y1 = 0,
    x2 = 1,
    y2 = 1;
    label = "",
    labeltsep = 0.1,
    labelrsep = 0.5,
    labelrelrot = false,
    labelrelangle = 0,
    ha = "center",
    va = "center",
    color="black",
    backgroundcolor = "none",
    arrowstyle1 = "<|-",
    arrowstyle2 = "-|>",
    linewidth = 0.6,
    linestyle = "-",
    width = 0.2,
    headlength = 5,
    headwidth = 2.5,
    par = 0,
    paroverhang = 0.02,
    parcolor = "black",
    parlinewidth = 0.6,
    parlinestyle = "-")

    # Length of line
    dr = sqrt((x2-x1)^2 + (y2-y1)^2)
    # x-difference of line
    drx = (x2 - x1)/dr
    # y-difference of line
    dry = (y2 - y1)/dr
    # Angle of line
    absangle = atan(dry, drx)
    # Orientation tangential to phasor (lagging by 90°)
    # x-component of tangential component with repsect to length
    dtx = +dry
    # y-component of tangential component with repsect to length
    dty = -drx
    # x-component of parallel shift of line
    dpx = -par*dtx
    # y-component of parallel shift of line
    dpy = -par*dty
    # x-component of auxiliary line indicating parallel shift of line
    dox = -paroverhang*dtx*sign(par)
    # y-component of auxiliary line indicating parallel shift of line
    doy = -paroverhang*dty*sign(par)

    # Draw line
    plot([x1+dpx,x2+dpx], [y1+dpy,y2+dpy],
        color=color, linewidth=linewidth,linestyle=linestyle)
    # Draw line indicating parallel shift
    if par!=0
        plot([x1,x1+dpx+dox], [y1,y1+dpy+doy],
            color=parcolor, linewidth=parlinewidth, linestyle=parlinestyle)
        plot([x2,x2+dpx+dox], [y2,y2+dpy+doy],
            color=parcolor, linewidth=parlinewidth, linestyle=parlinestyle)
    end

    # Draw arrows at begin and end of line, depending on the specified arrow
    # styles

    # Begin of line
    if arrowstyle1 == "<|-"
        # Arrow head
        annotate("", xy=(x1+dpx, y1+dpy),
            xytext=(x1+dpx+0.001*drx, y1+dpy+0.001*dry),
            xycoords="data",
            arrowprops=Dict("edgecolor"=>color, "facecolor"=>color,
                "width"=>width, "linestyle"=>"-",
                "headlength"=>headlength,
                "headwidth"=>headwidth),
            annotation_clip=false)
    elseif arrowstyle1 == "."
        # Dot marker
        plot(x1+dpx, y1+dpy, marker=".", color=color, clip_on=false)
    end

    # End of line
    if arrowstyle2 == "-|>"
        # Arrow head
        annotate("", xy=(x2+dpx, y2+dpy),
            xytext=(x2+dpx-0.001*drx, y2+dpy-0.001*dry),
            xycoords="data",
            arrowprops=Dict("edgecolor"=>color, "facecolor"=>color,
                "width"=>width, "linestyle"=>"-",
                "headlength"=>headlength,
                "headwidth"=>headwidth),
            annotation_clip=false)
    elseif arrowstyle2 == "."
        # Marker dot
        plot(x2+dpx, y2+dpy ,marker=".", color=color, clip_on=false)
    end

    # Plot label
    if labelrelrot == false
        # Without relative rotation of label
        text(x1 + dr*drx*labelrsep - dtx*labeltsep + dpx,
            y1 + dr*dry*labelrsep - dty*labeltsep + dpy,
            label, ha=ha, va=va, rotation=labelrelangle*180/pi,
            backgroundcolor=backgroundcolor)
    else
        # Applying relative rotation of label
        text(x1 + dr*drx*labelrsep - dtx*labeltsep + dpx,
            y1 + dr*dry*labelrsep - dty*labeltsep + dpy,
            label, ha=ha, va=va, rotation=(absangle+labelrelangle)*180/pi,
            backgroundcolor=backgroundcolor)
    end
end


text:
# A simple reference counting util

# import this function in Base so that we do not have conflicting definition
import Base.dec

export RefCounter, inc, dec, ref

type RefCounter
  count :: Int
end

RefCounter() = RefCounter(0)

function inc(r :: RefCounter)
  r.count += 1
end
function dec(r :: RefCounter)
  r.count -= 1
  return r.count
end

function ref(r :: RefCounter)
  inc(r)
  r
end


text:
using StaticArrays
using CLIMA.VariableTemplates
import CLIMA.DGmethods: BalanceLaw,
                        vars_aux, vars_state, vars_gradient, vars_diffusive,
                        flux_nondiffusive!, flux_diffusive!, source!,
                        gradvariables!, diffusive!,
                        init_aux!, init_state!,
                        boundary_state!, wavespeed, LocalGeometry,
                        num_state, num_gradient
using CLIMA.DGmethods.NumericalFluxes: NumericalFluxNonDiffusive,
                                       NumericalFluxDiffusive,
                                       NumericalFluxGradient
import CLIMA.DGmethods.NumericalFluxes: boundary_flux_diffusive!

abstract type AdvectionDiffusionProblem end
struct AdvectionDiffusion{dim, P, fluxBC} <: BalanceLaw
  problem::P
  function AdvectionDiffusion{dim}(problem::P
                                  ) where {dim, P <: AdvectionDiffusionProblem}
    new{dim, P, false}(problem)
  end
  function AdvectionDiffusion{dim, fluxBC}(problem::P
                                  ) where {dim, P <: AdvectionDiffusionProblem,
                                           fluxBC}
    new{dim, P, fluxBC}(problem)
  end
end

# Stored in the aux state are:
#   `coord` coordinate points (needed for BCs)
#   `u` advection velocity
#   `D` Diffusion tensor
vars_aux(::AdvectionDiffusion, FT) = @vars(coord::SVector{3, FT},
                                           u::SVector{3, FT},
                                           D::SMatrix{3, 3, FT, 9})
#
# Density is only state
vars_state(::AdvectionDiffusion, FT) = @vars(ρ::FT)

# Take the gradient of density
vars_gradient(::AdvectionDiffusion, FT) = @vars(ρ::FT)

# The DG auxiliary variable: D ∇ρ
vars_diffusive(::AdvectionDiffusion, FT) = @vars(σ::SVector{3,FT})

"""
    flux_nondiffusive!(m::AdvectionDiffusion, flux::Grad, state::Vars,
                       aux::Vars, t::Real)

Computes non-diffusive flux `F` in:

```
∂ρ
-- = - ∇ • (u ρ - σ) = - ∇ • F
∂t
```
Where

 - `u` is the advection velocity
 - `ρ` is the advected quantity
 - `σ` is DG auxiliary variable (`σ = D ∇ ρ` with D being the diffusion tensor)
"""
function flux_nondiffusive!(m::AdvectionDiffusion, flux::Grad, state::Vars,
                            aux::Vars, t::Real)
  ρ = state.ρ
  u = aux.u
  flux.ρ += u * ρ
end

"""
flux_diffusive!(m::AdvectionDiffusion, flux::Grad, auxDG::Vars)

Computes diffusive flux `F` in:

```
∂ρ
-- = - ∇ • (u ρ - σ) = - ∇ • F
∂t
```
Where

 - `u` is the advection velocity
 - `ρ` is the advected quantity
 - `σ` is DG auxiliary variable (`σ = D ∇ ρ` with D being the diffusion tensor)
"""
function flux_diffusive!(m::AdvectionDiffusion, flux::Grad, auxDG::Vars)
  σ = auxDG.σ
  flux.ρ += -σ
end
flux_diffusive!(m::AdvectionDiffusion, flux::Grad, state::Vars, auxDG::Vars,
                aux::Vars, t::Real) = flux_diffusive!(m, flux, auxDG)

"""
    gradvariables!(m::AdvectionDiffusion, transform::Vars, state::Vars,
                   aux::Vars, t::Real)

Set the variable to take the gradient of (`ρ` in this case)
"""
function gradvariables!(m::AdvectionDiffusion, transform::Vars, state::Vars,
                        aux::Vars, t::Real)
  transform.ρ = state.ρ
end

"""
    diffusive!(m::AdvectionDiffusion, transform::Vars, gradvars::Vars,
               aux::Vars)

Set the variable to take the gradient of (`ρ` in this case)
"""
function diffusive!(m::AdvectionDiffusion, auxDG::Vars, gradvars::Grad,
                    aux::Vars)
  ∇ρ = gradvars.ρ
  D = aux.D
  auxDG.σ = D * ∇ρ
end
diffusive!(m::AdvectionDiffusion, auxDG::Vars, gradvars::Grad,
           state::Vars, aux::Vars, t::Real) = diffusive!(m, auxDG, gradvars,
                                                         aux)

"""
    source!(m::AdvectionDiffusion, _...)

There is no source in the advection-diffusion model
"""
source!(m::AdvectionDiffusion, _...) = nothing

"""
    wavespeed(m::AdvectionDiffusion, nM, state::Vars, aux::Vars, t::Real)

Wavespeed with respect to vector `nM`
"""
function wavespeed(m::AdvectionDiffusion, nM, state::Vars, aux::Vars, t::Real)
  u = aux.u
  abs(dot(nM, u))
end

"""
    init_aux!(m::AdvectionDiffusion, aux::Vars, geom::LocalGeometry)

initialize the auxiliary state
"""
function init_aux!(m::AdvectionDiffusion, aux::Vars, geom::LocalGeometry)
  aux.coord = geom.coord
  init_velocity_diffusion!(m.problem, aux, geom)
end

function init_state!(m::AdvectionDiffusion, state::Vars, aux::Vars,
                     coords, t::Real)
  initial_condition!(m.problem, state, aux, coords, t)
end

Neumann_data!(problem, ∇state, aux, x, t) = nothing
Dirichlet_data!(problem, state, aux, x, t) = nothing

function boundary_state!(nf, m::AdvectionDiffusion, stateP::Vars, auxP::Vars,
                         nM, stateM::Vars, auxM::Vars, bctype, t, _...)
  if bctype == 1 # Dirichlet
    Dirichlet_data!(m.problem, stateP, auxP, auxP.coord, t)
  elseif bctype ∈ (2, 4) # Neumann
    stateP.ρ = stateM.ρ
  elseif bctype == 3 # zero Dirichlet
    stateP.ρ = 0
  end
end

function boundary_state!(nf::CentralNumericalFluxDiffusive,
    m::AdvectionDiffusion,
    state⁺::Vars, diff⁺::Vars, aux⁺::Vars,
    n⁻::SVector,
    state⁻::Vars, diff⁻::Vars, aux⁻::Vars,
    bctype, t,
    _...)

  if bctype ∈ (1,3) # Dirchlet
    # Just use the minus side values since Dirchlet
    diff⁺.σ = diff⁻.σ
  elseif bctype == 2 # Neumann with data
    FT = eltype(diff⁺)
    ngrad = num_gradient(m, FT)
    ∇state = Grad{vars_gradient(m, FT)}(similar(parent(diff⁺), Size(3, ngrad)))
    # Get analytic gradient
    Neumann_data!(m.problem, ∇state, aux⁻, aux⁻.coord, t)
    diffusive!(m, diff⁺, ∇state, aux⁻)
    # compute the diffusive flux using the boundary state
  elseif bctype == 4 # zero Neumann
    FT = eltype(diff⁺)
    ngrad = num_gradient(m, FT)
    ∇state = Grad{vars_gradient(m, FT)}(similar(parent(diff⁺), Size(3, ngrad)))
    # Get analytic gradient
    ∇state.ρ = SVector{3, FT}(0, 0, 0)
    # convert to auxDG variables
    diffusive!(m, diff⁺, ∇state, aux⁻)
  end
  nothing
end

function boundary_flux_diffusive!(nf::CentralNumericalFluxDiffusive,
                                  m::AdvectionDiffusion{dim, P, true},
                                  F,
                                  state⁺, diff⁺, aux⁺, n⁻,
                                  state⁻, diff⁻, aux⁻,
                                  bctype, t,
                                  _...) where {dim, P}

  # Default initialize flux to minus side
  if bctype ∈ (1,3) # Dirchlet
    # Just use the minus side values since Dirchlet
    flux_diffusive!(m, F, state⁻, diff⁻, aux⁻, t)
  elseif bctype == 2 # Neumann data
    FT = eltype(diff⁺)
    ngrad = num_gradient(m, FT)
    ∇state = Grad{vars_gradient(m, FT)}(similar(parent(diff⁺), Size(3, ngrad)))
    # Get analytic gradient
    Neumann_data!(m.problem, ∇state, aux⁻, aux⁻.coord, t)
    # get the diffusion coefficient
    D = aux⁻.D
    # exact the exact data
    ∇ρ = ∇state.ρ
    # set the flux
    F.ρ = - D * ∇ρ
  elseif bctype == 4 # Zero Neumann
    FT = eltype(diff⁺)
    F.ρ = SVector{3, FT}(0, 0, 0)
  end
  nothing
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
__precompile__(true)
module IntervalMatrices

using LinearAlgebra

import Random
using Random: AbstractRNG, GLOBAL_RNG, seed!

using Reexport
@reexport using IntervalArithmetic

# ================================
# Type defining an interval matrix
# ================================
include("matrix.jl")

# =================================
# Arithmetic for interval matrices
# =================================
include("arithmetic.jl")

# =======================================================
# Methods to handle the exponential of an interval matrix
# =======================================================
include("exponential.jl")

# ============================================================
# Method to compute the correction hull of an interval matrix
# ============================================================
include("correction_hull.jl")

# ========
# Exports
# ========
export AbstractIntervalMatrix,
       IntervalMatrix

export inf, sup,
       rand, sample,
       opnorm,
       quadratic_expansion,
       expm_overapproximation,
       expm_underapproximation

end # module


text:
using Test, FastChebInterp, StaticArrays

# similar to ≈, but acts elementwise on tuples
≈′(a::Tuple, b::Tuple; kws...) where {N} = length(a) == length(b) && all(xy -> isapprox(xy[1],xy[2]; kws...), zip(a,b))

@testset "1d test" begin
    lb,ub = -0.3, 0.9
    f(x) = exp(x) / (1 + 2x^2)
    f′(x) = f(x) * (1 - 4x/(1 + 2x^2))
    @test_throws ArgumentError chebpoints(-1, lb, ub)
    x = chebpoints(48, lb, ub)
    interp = chebinterp(f.(x), lb, ub)
    @test ndims(interp) == 1
    x1 = 0.2
    @test interp(x1) ≈ f(x1)
    @test chebgradient(interp, x1) ≈′ (f(x1), f′(x1))
end

@testset "2d test" begin
    lb, ub = [-0.3,0.1], [0.9,1.2]
    f(x) = exp(x[1]+2*x[2]) / (1 + 2x[1]^2 + x[2]^2)
    ∇f(x) = f(x) * SVector(1 - 4x[1]/(1 + 2x[1]^2 + x[2]^2), 2 - 2x[2]/(1 + 2x[1]^2 + x[2]^2))
    x = chebpoints((48,39), lb, ub)
    interp = chebinterp(f.(x), lb, ub)
    interp0 = chebinterp(f.(x), lb, ub, tol=0)
    @test ndims(interp) == 2
    x1 = [0.2, 0.3]
    @test interp(x1) ≈ f(x1)
    @test interp(x1) ≈ interp0(x1) rtol=1e-15
    @test all(n -> n[1] < n[2], zip(size(interp.coefs), size(interp0.coefs)))
    @test chebgradient(interp, x1) ≈′ (f(x1), ∇f(x1))

    # univariate function in 2d should automatically drop down to univariate polynomial
    f1(x) = exp(x[1]) / (1 + 2x[1]^2)
    interp1 = chebinterp(f1.(x), lb, ub)
    @test interp1(x1) ≈ f1(x1)
    @test size(interp1.coefs, 2) == 1 # second dimension should have been dropped

    # complex and vector-valued interpolants:
    f2(x) = [f(x), cis(x[1]*x[2] + 2x[2])]
    ∇f2(x) = vcat(transpose(∇f(x)), transpose(SVector(im*x[2], im*(x[1] + 2)) * cis(x[1]*x[2] + 2x[2])))
    interp2 = chebinterp(f2.(x), lb, ub)
    @test interp2(x1) ≈ f2(x1)
    @test chebjacobian(interp2, x1) ≈′ (f2(x1), ∇f2(x1))

    # chebinterp_v1
    av1 = Array{ComplexF64}(undef, 2, size(x)...)
    av1[1,:,:] .= f.(x)
    av1[2,:,:] .= (x -> f2(x)[2]).(x)
    interp2v1 = chebinterp_v1(av1, lb, ub)
    @test interp2v1(x1) ≈ f2(x1)
    @test chebjacobian(interp2v1, x1) ≈′ (f2(x1), ∇f2(x1))
end

@testset "1d regression" begin
    x = range(1,2, length=200)
    y = @. sin.(x + 0.2) * x
    c = chebregression(x, y, 2)
    c2 = x.^(0:2)' \ y # Vandermonde-style fit
    @test c(1.1) ≈ c2[1] + 1.1 * c2[2] + 1.1^2 * c2[3] rtol=1e-13
end

@testset "2d regression" begin
    x = vec(SVector.(range(-3, 2, length=50), range(0,1, length=40)'))
    fr(x::SVector{2}) = sin(x[1]/3 + 1)*exp(x[2])
    c = chebregression(x, fr.(x), (2,3))

    # Vandermonde-style fit:
    powsvec(x) = vec([x[1]^i * x[2]^j for i=0:2, j=0:3])'
    c2 = mapreduce(powsvec, vcat, x) \ fr.(x)

    p = SVector(1.12345, sqrt(0.5))
    @test c(p) ≈ powsvec(p)*c2 rtol=1e-13

    cM = chebregression(mapreduce(transpose∘Vector, vcat, x), fr.(x), (2,3))
    @test cM.coefs == c.coefs

    fr4(x::SVector{2}) = cos(x[1]/3 + 1)*sin(x[2])
    c4 = chebregression(x, fr4.(x), (2,3))
    F14 = [SVector(fr(x),fr4(x)) for x in x]
    c14 = chebregression(x, F14, (2,3))
    @test first.(c14.coefs) ≈ c.coefs rtol=1e-13
    @test last.(c14.coefs) ≈ c4.coefs rtol=1e-13
    c14M = chebregression(x, mapreduce(transpose∘Vector, vcat, F14), (2,3))
    @test c14M.coefs == c14.coefs
end

@testset "degree 1" begin
    lb,ub = -0.3, 0.9
    x = chebpoints(1, lb, ub)
    interp = chebinterp((x -> 2x+3).(x), lb, ub)
    @test interp(0) ≈ 3
    @test interp(0.1) ≈ 3.2
    @test chebgradient(interp, 0)[2] ≈ 2 ≈ chebgradient(interp, 0.1)[2]

    lb, ub = [-0.3,-0.1], [0.9,1.2]
    x = chebpoints((1,1), lb, ub)
    interp = chebinterp((x -> (2x[1]+3)*(4x[2]+5)).(x), lb, ub)
    @test interp([0,0]) ≈ 3*5
    @test interp([0.1,0]) ≈ 3.2*5
    @test interp([0,0.1]) ≈ 3*5.4
    @test chebgradient(interp, [0,0])[2] ≈ [2*5, 3*4]
end

@testset "degree 0" begin
    lb,ub = -0.3, 0.9
    x = chebpoints(0, lb, ub)
    interp = chebinterp(sin.(x), lb, ub)
    @test interp(0) ≈ interp(0.1) ≈ sin((lb+ub)/2)
    @test chebgradient(interp, 0)[2] == 0 == chebgradient(interp, 0.1)[2]

    lb, ub = [-0.3,-0.1], [0.9,1.2]
    x = chebpoints((0,0), lb, ub)
    interp = chebinterp((x -> sin(x[1])*cos(x[2])).(x), lb, ub)
    @test interp([0,0]) ≈ interp([0.1,0.2]) ≈ sin(0.3)*cos(0.55)
    @test chebgradient(interp, [0.1,0.2])[2] == [0, 0]

    lb, ub = [-0.3,-0.1], [0.9,1.2]
    x = chebpoints((1,0), lb, ub)
    interp = chebinterp((x -> (2x[1]+3)*cos(x[2])).(x), lb, ub)
    @test interp([0,0]) ≈ 3*cos(0.55)
    @test interp([0.1,0.22]) ≈ 3.2*cos(0.55)
    @test chebgradient(interp, [0.1,0.2])[2] == [2*cos(0.55), 0]
end


text:
@testset "Vocabulary reduction" begin
    d = 5
    vocab = string.('a':'z')
    vectors = rand(d, length(vocab))
    vectors./=permutedims(map(norm, (vectors[:, i] for i in 1:size(vectors,2))))
    wv = WordVectors(vocab, vectors)

    rvocab = ["a", "z"]
    rwv = vocab_reduction(wv, rvocab, 1)
    @test rwv.vocab == rvocab
    @test size(rwv) == (d, length(rvocab))

    rvocab = ["a", "z"]
    rwv = vocab_reduction(wv, rvocab, Dict(r=>1 for r in rvocab))
    @test rwv.vocab == rvocab
    @test size(rwv) == (d, length(rvocab))

    rvocab = ["a", "z"]
    rwv = vocab_reduction(wv, [rvocab..., "missing"], Dict(r=>1 for r in rvocab))
    @test rwv.vocab == rvocab
    @test size(rwv) == (d, length(rvocab))
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
import PathWeightSampling
using Test

system = PathWeightSampling.cooperative_chemotaxis_system(dtimes=0:0.1:10)
conf = PathWeightSampling.generate_full_configuration(system)
@test length(conf.s_traj[:, 1]) == 1
@test length(conf.r_traj[:, 1]) == 40
@test length(conf.x_traj[:, 1]) == 2

mconf = PathWeightSampling.generate_configuration(system)
@test length(mconf.s_traj[:, 1]) == 1
@test length(mconf.x_traj[:, 1]) == 2

for u in conf.r_traj.u
    @test sum(u) == sum(conf.r_traj.u[1]) # this is true because the total number of receptors must be a constant
end

for u in conf.x_traj.u
    @test sum(u) == sum(conf.x_traj.u[1]) # this is true because the total number of CheY proteins must be a constant
end

cond_ens = PathWeightSampling.ConditionalEnsemble(system)
for t = 0:0.1:10
    local u0 = PathWeightSampling.sample_initial_condition(cond_ens)
    u_new, weight = PathWeightSampling.propagate(mconf, cond_ens, u0, (0.0, t))
    @test size(u0) == size(u_new)
    @test sum(u0[2:end]) == sum(u_new[2:end]) # this is true because the total number of receptors must be a constant
    @test !isinf(weight)
end

algorithm = PathWeightSampling.SMCEstimate(16)
result = PathWeightSampling.mutual_information(system, algorithm, num_samples=1)
@test !all(isinf.(result.MutualInformation[1]))

system_direct_agg = PathWeightSampling.cooperative_chemotaxis_system(dtimes=0:0.1:10, dist_aggregator=PathWeightSampling.GillespieDirect())
@test typeof(system.dist.aggregator) <: PathWeightSampling.DepGraphAggregator
@test typeof(system_direct_agg.dist.aggregator) <: PathWeightSampling.DirectAggregator

cens1 = PathWeightSampling.ConditionalEnsemble(system)
cens2 = PathWeightSampling.ConditionalEnsemble(system_direct_agg)

a = PathWeightSampling.energy_difference(conf, cens1)
b = PathWeightSampling.energy_difference(conf, cens2)
@test a == b

mens1 = PathWeightSampling.MarginalEnsemble(system)
mens2 = PathWeightSampling.MarginalEnsemble(system_direct_agg)

mconf = PathWeightSampling.marginal_configuration(conf)
a = PathWeightSampling.energy_difference(mconf, mens1)
b = PathWeightSampling.energy_difference(mconf, mens2)
@test a == b


text:
module MultiComponentFlash
    using LinearAlgebra, ForwardDiff, StaticArrays, Roots
    # Types that define specific cubic equations of state
    export SoaveRedlichKwong, RedlichKwong, PengRobinson, ZudkevitchJoffe
    # The generic cubic form that supports the above
    export GenericCubicEOS, number_of_components
    # Flash interfaces
    export flash_2ph, flash_2ph!, flash_storage
    export stability_2ph, stability_2ph!
    # Algorithms for flash
    export SSIFlash, NewtonFlash, SSINewtonFlash
    # Mixtures and their molecular makeup
    export MolecularProperty, MultiComponentMixture
    # K-values
    export wilson_estimate, wilson_estimate!, initial_guess_K, initial_guess_K!
    # Vapor-liquid equilibrium
    export solve_rachford_rice

    export liquid_mole_fraction, vapor_mole_fraction
    export component_fugacity, mixture_fugacities, mixture_compressibility_factor

    export force_scalars, force_coefficients, force_coefficients!
    export critical_pressure, critical_temperature, critical_volume, acentric_factor, molar_weight
    export cubic_benchmark
    export single_phase_label

    # Types and utilities for coupling flash to other codes.
    export SinglePhaseLiquid, SinglePhaseVapor, TwoPhaseLiquidVapor, FlashedPhase, UnknownPhaseState, FlashedMixture2Phase
    export phase_saturations, mass_density, mass_densities, molar_volume, phase_is_present
    export inverse_flash_update!
    export lbc_viscosity, lbc_viscosities
    export set_partials, set_partials_phase_mole_fractions!, set_partials_vapor_fraction

    const MINIMUM_COMPOSITION = 1e-10
    const IDEAL_GAS_CONSTANT = 8.3144598
    include("mixtures.jl")
    include("kvalues.jl")
    include("rachford_rice.jl")
    include("eos.jl")
    include("flash.jl")
    include("derivatives.jl")
    include("stability.jl")
    include("tables.jl")

    include("flow_coupler.jl")
    include("utils.jl")
end


text:
# Constants
const whichRounding = Dict{Int, Function}(
    0 => simplerounding!,
    1 => recursiverounding!,
)

const whichProjection = Dict{Int, Function}(
    0 => minimumdistance,
    1 => minimumdeltaalpha,
    2 => minimumdistance_reduce,
)

const whichPerturb = Dict{Int, Function}(
    -1 => noperturb,
    0 => perturb!,
    1 => randomperturb!,
    2 => frequencyperturb!,

    100 => perturbRestart!,
)

const whichRestart = Dict{Int, Function}(
    -1 => norestart,
    0 => perturbRestart!,

    100 => perturb!,
    101 => randomperturb!,
    102 => frequencyperturb!,
)

const __paramToInt = Dict{String, Int}([
    "logLevel" => 1,
    "timelim" => 2,
    "roundingMethod" => 10,
    "recursiveRounding" => 11,
    "projectionMethod" => 20,
    "perturbationMethod" => 30,
    "perturbationTT" => 31,
    "perturbationTTmin" => 32,
    "perturbationTTmax" => 33,
    "checkCycle" => 34,
    "restart" => 35,
    "presolve" => 41,
    "improve" => 42,
    "feasibilityCheck" => 100
])

const __paramToString = Dict{Int, String}([
    1 => "logLevel",
    2 => "timelim",
    10 => "roundingMethod",
    11 => "recursiveRounding",
    20 => "projectionMethod",
    30 => "perturbationMethod",
    31 => "perturbationTT",
    32 => "perturbationTTmin",
    33 => "perturbationTTmax",
    34 => "checkCycle",
    35 => "restart",
    41 => "presolve",
    42 => "improve",
    100 => "feasibilityCheck"
])

max_stars_repo_name: remi-garcia/FeasibilityPump

max_issues_repo_name: remi-garcia/FeasibilityPump

max_forks_repo_name: remi-garcia/FeasibilityPump


text:
##########
##########
# Functions to perform perturbations on the food webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

# set up =======================================================================
using DrWatson
@quickactivate "trait-based-rewiring"
# logger set up
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())
using BenchmarkTools

include(srcdir("01_eqn_component_functions.jl"))
include(srcdir("00_species_list_creation.jl"))
include(srcdir("03_initial_sims_functions.jl"))
include(srcdir("04_perturbation_sims_functions.jl"))

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# pulse perturbations ==========================================================

# little line to set the network id (this happens from the shell line)
network_id = parse(Int64, ARGS[1])

species_list = readdlm(datadir("./intermediate-files/species-lists-pre-perturbation/species_list_0001.csv"), ',', Any)
network = create_interaction_matrix(species_list)

# set pulse amounts
c = 210
x = 100
f = 1
z = c/(x*(1/f))
# set vector callback
extinction_callback =
    VectorContinuousCallback(extinction_threshold,
                            extinction_affect!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))
steady_callback = TerminateSteadyState(1e-8,
                                        1e-6,
                                        DiffEqCallbacks.allDerivPass,
                                        min_t = 1205)

# set callback set
"""
    affect_inc!(integrator)

Affect for the condition to work on. Increase the mortality by making the amount
(the fourth element in p) equal to z.
"""
function affect_inc!(integrator, z)
  integrator.p[4] = z
end

"""
    affect_dec!(integrator)

Affect for the condition to work on. Decrease the fourth element in p,
returning it to 1.0.
"""
function affect_dec!(integrator)
  integrator.p[4] = nothing
end

save_positions1 = (true, true)
cb_inc = DiscreteCallback(inc_mort_cond, affect_inc!,
                            save_positions=save_positions1)
save_positions2 = (false, true)
cb_dec = DiscreteCallback(dec_mort_cond, affect_dec!,
                            save_positions=save_positions2)
callback_set = CallbackSet(extinction_callback,
                            #steady_callback, #keep this out until fixed
                            cb_inc,
                            cb_dec)

#  allocate x0
u0 = Float64[]
for i in 1:size(species_list, 1)
    append!(u0, species_list[i, biomass_pos])
end

# set perturbation tspan
tspan = (1000.0, 1100.0)

# set constant stop values
const stop1 = [1010.0]
const stop2 = [1080.0]
const tstop = [1010.;1080.]

# loop through random selections of species to perturb
#@info{"starting loop..."}

# set times to save from the simulation
save_times = [collect(range(1000.0, 1300.0, step = 0.1));
                collect(range(1250, 4800, step = sample(200:400, 1)[1]));
                collect(range(4800.0, 5000.0, step = 0.1))]

for j in 1:100

    # pick set of random species to perturb
    all = Array(1:size(species_list, 1))
    perturb_sp = rand(Array(1:size(species_list, 1)),
                        Int(round(size(species_list, 1)*(1/2))))
    non_perturb = setdiff(Set(all), Set(perturb_sp))

    # set p
    amount = nothing
    p = [species_list, network, 0.3, amount, perturb_sp, non_perturb]

    prob_pulse = ODEProblem(biomass_perturb_random,
                            u0,
                            tspan,
                            p,
                            tstops=tstop)
    @benchmark sol_prob_pulse = solve(prob_pulse,
                                        OwrenZen5(),
                                        callback = callback_set,
                                        reltol = 1e-4,
                                        # save every part for the perturbation and
                                        # the first 150 time steps after the
                                        # perturbation, but then just grab a few
                                        # from there to the end
                                        saveat = save_times,
                                        progress = true,
                                        progress_steps = 1)
    # save output of the solver
    sol_name = string("sol-pulse-net-", lpad(network_id,4,"0"), "-sim-", lpad(j, 4, "0"), ".csv")
    writedlm(datadir("./simulation-outpus/pulse-sims/output/$name"), sol_prob_pulse, ",")
    #plot(sol_prob_pulse)
    u_matrix = zeros(length(sol_prob_pulse.u), length(sol_prob_pulse.u[1]))
    for i in 1:size(u_matrix, 1)
        u_matrix[i,:] = sol_prob_pulse.u[i]
    end
    output = hcat(sol_prob_pulse.t, u_matrix)
    # write out the biomass outputs
    output_name = string("pulse-net-", lpad(network_id,4,"0"), "-sim-", lpad(j, 4, "0"), ".csv")
    writedlm(datadir("./simulation-outpus/pulse-sims/output/$name"), output_name, ",")
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
n = 10000
x = .25 * randn(Float64, (10000, 3))
ccall((:dvz_demo_scatter, "../../build/libdatoviz.so"), Int32, (Cint, Ref{Cdouble}), 10000, x)

max_stars_repo_name: djhoese/datoviz

max_issues_repo_name: djhoese/datoviz

max_forks_repo_name: djhoese/datoviz


text:
module Highlight

module SQL

using Crayons

const SYNTAX = Dict(
  :ACTIONS   =>   ["ALTER", "CREATE", "DELETE", "DROP", "INSERT", "EXPLAIN", "SELECT", "TRUNCATE", "UPDATE"],
  :KEYWORDS  =>   ["BEGIN", "CONSTRAINT", "FROM", "GROUP", "HAVING", "ORDER",
                  "ROLLBACK", "SET", "TABLE", "TRANSACTION", "UNION", "VIEW"],
  :OPERATORS =>   ["AND", "BETWEEN", "IN", "XOR", "OR ", "!=", "<>", ">=", "<=", ">", "<", "/"],
  :OTHER =>       ["ASC", "AS ", "BY", "DESC", "OFFSET", "ON"],
  :CONDITIONAL => ["CASE", "ELSE", "ELSEIF", "EXISTS", "IF", "NOT", "THEN", "WHERE"],
  :VALUES =>      ["FALSE", "KEY", "NULL", "PRIMARY", "TRUE", "VALUE"],
  :JOINS =>       ["INNER", "JOIN", "LEFT", "OUTER", "RIGHT"],
  :TYPES =>       ["BOOLEAN", "DATE", "SERIAL", "TEXT", "VARCHAR"],
)

const CRAYONS = Dict(
  :ACTIONS      => crayon"light_magenta",
  :KEYWORDS     => crayon"cyan",
  :OPERATORS    => crayon"light_red",
  :OTHER        => crayon"blue",
  :CONDITIONAL  => crayon"magenta",
  :VALUES       => crayon"green",
  :JOINS        => crayon"yellow",
  :SPECIAL      => crayon"red",
  :TYPES        => crayon"light_green",
)

end

using Crayons
using .SQL

function highlight(s::String, m::Module = SQL)
  for (key, word_list) in m.SYNTAX
    for word in word_list
      s = replace(s, Regex(word, "i") => string(m.CRAYONS[key], word, crayon"default"))
    end
  end

  s = replace(replace(replace(s, r"\"[a-zA-Z]*\"" => s"---->\0<----"), "---->\"" => string("\"", crayon"green")), "\"<----" => string("\"", crayon"default"))
  s = replace(replace(replace(s, r"\"[a-zA-Z_]*\"" => s"---->\0<----"), "---->\"" => string("\"", crayon"magenta")), "\"<----" => string("\"", crayon"default"))

  s = replace(s, "\"" => string(m.CRAYONS[:SPECIAL], "\"", crayon"default"))
  s = replace(s, "'" => string(m.CRAYONS[:SPECIAL], "'", crayon"default"))
  s = replace(s, "`" => string(m.CRAYONS[:SPECIAL], "`", crayon"default"))
  s = replace(s, "*" => string(m.CRAYONS[:SPECIAL], "*", crayon"default"))
  s = replace(s, "." => string(m.CRAYONS[:SPECIAL], ".", crayon"default"))

  s
end

end


text:
using LinearAlgebra: Matrix
using CompressingSolvers
using LinearAlgebra
using Test
using Plots

@testset "CompressingSolvers.jl" begin
    # Write your tests here.
    # Testing domain.jl
    @testset "domain.jl" begin
        include("test_domains.jl")
    end
    # Testing creating_problems.jl
    @testset "create_problems.jl" begin
        include("./test_create_problems.jl")
    end
    @testset "reconstruction.jl" begin
        include("./test_reconstruction.jl")
    end
end


text:
include("functions/GenerateCareerPaths.jl")
using Compose, GraphPlot

# Deinition on constants which are the vertices of our graph
const in_world  = 1
const _3D_BDL   = 2
const _3D_B     = 3
const _2C_BDL   = 4
const _2C_B     = 5
const _2B_BDL   = 6
const _2B_B     = 7
const _1B_BDL   = 8
const _1B_B     = 9
const _1A_BDL   = 10
const _1A_B     = 11
const out_world = 12

FirstLevel = Vector{AbstractLevel}()
push!(FirstLevel, AcademicLevel("3D-BDL"))
push!(FirstLevel, AcademicLevel("3D-B"))
push!(FirstLevel, AcademicLevel("2C-BDL"))
push!(FirstLevel, AcademicLevel("2C-B"))
push!(FirstLevel, AcademicLevel("2B-BDL"))
push!(FirstLevel, AcademicLevel("2B-B"))
push!(FirstLevel, AcademicLevel("1B-BDL"))
push!(FirstLevel, AcademicLevel("1B-B"))
push!(FirstLevel, AcademicLevel("1A-BDL"))
push!(FirstLevel, AcademicLevel("1A-B"))

g = DiGraph(12)

# Recruitment edges i.e. statrting from in_world
add_edge!(g, in_world, _3D_BDL)
add_edge!(g, in_world, _2C_BDL)
add_edge!(g, in_world, _2B_BDL)
add_edge!(g, in_world, _1B_BDL)
add_edge!(g, in_world, _1A_BDL)

# Leaving the organization (retirement and attrition) i.e. connect all edges to the leaving vertex
for i in 2:11
    add_edge!(g, i, out_world)
end

# Connecting nodes with even index to all nodes with a higher index
for i in 1:5
    for j in (2*i + 1):11
        add_edge!(g, 2*i, j)
    end
end

# Connecting nodes with odd index to nodes with a higher odd index
for i in 1:5
    for j in (i+1):5
        add_edge!(g, 2*i + 1, 2*j + 1)
    end
end
levelnames = [FirstLevel[i].Name for i in 1:length(FirstLevel)]
draw(PNG("InitGraph.png", 50cm, 50cm), gplot(g, nodelabel=["In_world"; levelnames ;"Out_world"]))

Paths = ConstructCareerPaths(g, [FirstLevel])

nbVertices = 0
for i in 1:length(Paths)
    nbVertices += length(Paths[i].Path)
end
PlotG = DiGraph(length(Paths) + nbVertices + length(Paths))

NodeLabels = Array{String}(length(Paths) + nbVertices + length(Paths))
i = 1
for path in Paths
    NodeLabels[i] = "World"
    i += 1
    for status in path.Path
        NodeLabels[i] = status.AcademicLevel.Name
        add_edge!(PlotG, i - 1, i)
        i += 1
    end
    add_edge!(PlotG, i - 1, i)
    NodeLabels[i] = "World"
    i += 1
end
draw(SVG("CareerPaths.svg", 50cm, 50cm), gplot(PlotG, nodelabel=NodeLabels))
draw(PNG("CareerPaths.png", 50cm, 50cm), gplot(PlotG, nodelabel=NodeLabels))

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
export LaurentDirichlet

##
# LaurentDirichlet represents functions in the basis
#       1, z^(-1)+1,z+1,z^(-2)+z^(-1),z^2+z,…
# which incorporates decay at z = -1
# TODO: restructure as SumSpace made out of HardyDirichlet
##

immutable LaurentDirichlet{DD} <: UnivariateSpace{ComplexBasis,DD}
    domain::DD
    LaurentDirichlet(d::DD)=new(d)
end
LaurentDirichlet{T<:Number}(d::Vector{T})=LaurentDirichlet(PeriodicDomain(d))
LaurentDirichlet(d::Domain)=LaurentDirichlet{typeof(d)}(d)
LaurentDirichlet()=LaurentDirichlet(PeriodicInterval())

spacescompatible(a::LaurentDirichlet,b::LaurentDirichlet)=domainscompatible(a,b)

canonicalspace(S::LaurentDirichlet)=Laurent(domain(S))

Conversion{DD}(a::LaurentDirichlet{DD},b::Laurent{DD})=ConcreteConversion(a,b)
bandinds{DD}(::ConcreteConversion{LaurentDirichlet{DD},Laurent{DD}})=0,2

function getindex{DD}(C::ConcreteConversion{LaurentDirichlet{DD},Laurent{DD}},k::Integer,j::Integer)
    if k==1 && j==2
        one(eltype(C))
    elseif k==j || j==k+2
        one(eltype(C))
    else
        zero(eltype(C))
    end
end

conversion_rule{DD}(b::LaurentDirichlet,a::Laurent{DD})=b

differentiate{DD}(f::Fun{LaurentDirichlet{DD}}) = differentiate(Fun(f,Laurent))

for op in (:+,:-,:(.+),:(.-),:(.*))
    @eval begin
        $op{DD}(f::Fun{Laurent{DD}},g::Fun{LaurentDirichlet{DD}}) = $op(f,Fun(g,Laurent))
        $op{DD}(f::Fun{LaurentDirichlet{DD}},g::Fun{Laurent{DD}}) = $op(Fun(f,Laurent),g)
        $op{DD}(f::Fun{LaurentDirichlet{DD}},g::Fun{LaurentDirichlet{DD}}) = $op(Fun(f,Laurent),Fun(g,Laurent))
    end
end

Base.real{DD}(f::Fun{LaurentDirichlet{DD}}) = real(Fun(f,Laurent))
Base.imag{DD}(f::Fun{LaurentDirichlet{DD}}) = imag(Fun(f,Laurent))

coefficients(v::Vector,::Laurent,::LaurentDirichlet)=laurentdirichlettransform!(copy(v))
coefficients(v::Vector,::LaurentDirichlet,::Laurent)=laurentidirichlettransform!(copy(v))

function laurentdirichlettransform!(w::Vector)
    if length(w) > 1
        for k=length(w)-2:-1:1
            @inbounds w[k] -= w[k+2]
        end
        @inbounds w[1] -= w[2]
    end

    w
end

function laurentidirichlettransform!(w::Vector)
    if length(w) == 1
    elseif length(w) == 2
        @inbounds w[1] += w[2]
    else
        @inbounds w[1] += w[2]+w[3]
        for k=4:length(w)
            @inbounds w[k-2] += w[k]
        end
    end

    w
end

##
# CosDirichlet represents
# 1,cos(θ)+1,cos(2θ)+cos(θ),…
#

immutable CosDirichlet{DD} <: RealUnivariateSpace{DD}
    domain::DD
end

CosDirichlet()=CosDirichlet(PeriodicInterval())

spacescompatible(a::CosDirichlet,b::CosDirichlet) = domainscompatible(a,b)

canonicalspace(S::CosDirichlet)=CosSpace(domain(S))

Conversion(a::CosSpace,b::CosDirichlet) = ConcreteConversion(a,b)
bandinds{CS<:CosSpace,CD<:CosDirichlet}(::ConcreteConversion{CD,CS})=0,1
getindex{CS<:CosSpace,CD<:CosDirichlet}(C::ConcreteConversion{CD,CS},k::Integer,j::Integer) =
    (k==j||j==k+1)?one(eltype(C)):zero(eltype(C))

conversion_rule(b::CosDirichlet,a::CosSpace)=b


text:
module LNR

using Compat
using Lazy
import Base: seek, position, eof, peek

export LineNumberingReader, line, column, Cursor, cursor, seekline, seekcol

# Cursor Type

immutable Cursor
  line::Int
  column::Int
end

Base.isless(x::Cursor, y::Cursor) =
  x.line < y.line || (x.line == y.line && x.column < y.column)

# LNR

immutable LineNumberingReader{T<:IO} <: IO
  io::T
  lines::Vector{Int} # The byte index of the first char in each line
end

LineNumberingReader(io::IO) = LineNumberingReader(io, [1])
LineNumberingReader(s::AbstractString) = LineNumberingReader(IOBuffer(s))

eof(r::LineNumberingReader) = eof(r.io)
position(r::LineNumberingReader) = position(r.io)
peek(r::LineNumberingReader) = Base.peek(r.io)

scannedindex(r::LineNumberingReader, i) = i < r.lines[end]

function Base.read(r::LineNumberingReader, ::Type{UInt8})
  c = read(r.io, UInt8)
  c == '\n' && !scannedindex(r, position(r)) && !eof(r) &&
     push!(r.lines, position(r)+1)
  return c
end

function Base.skip(r::LineNumberingReader, n::Integer)
  if n > 0 && !scannedindex(r, position(r) + n)
    @dotimes n read(r, UInt8)
  else
    skip(r.io, n)
  end
end

# Seeking

function seek(r::LineNumberingReader, pos)
  scannedindex(r, pos) ? seek(r.io, pos) : skip(r, pos-position(r))
  return r
end

function seekline(r::LineNumberingReader, line::Integer)
  line ≤ length(r.lines) && return seek(r, r.lines[line]-1)
  while length(r.lines) < line && !eof(r)
    readline(r)
  end
  return r
end

seekline(r::LineNumberingReader) = seekline(r, line(r))

function seekcol(r::LineNumberingReader, col::Integer)
  seekline(r)
  c = 1
  while c < col && !eof(r)
    read(r, Char) == '\n' && return skip(r,-1)
    c += 1
  end
  return r
end

seek(r::LineNumberingReader, c::Cursor) =
  @> r seekline(c.line) seekcol(c.column)

# Cursor finding

function line(s::LineNumberingReader)
  p = position(s)+1 # after the cursor
  for i = 1:length(s.lines)
    s.lines[i] > p && return i - 1
  end
  return length(s.lines)
end

function column(r::LineNumberingReader)
  pos = position(r)
  col = 1
  seekline(r)
  while position(r) < pos
    read(r, Char)
    col += 1
  end
  return col
end

cursor(s::LineNumberingReader) = Cursor(line(s), column(s))
cursor(l, c) = Cursor(l, c)

function position(r::LineNumberingReader, p::Cursor)
  withstream(r) do
    seek(r, p)
    position(r)
  end
end

# Position is 0-based
index(r::LineNumberingReader, p::Cursor) = position(r, p) + 1

function peekchar(r::LineNumberingReader)
  withstream(r) do
    read(r, Char)
  end
end

# Util

function withstream(f, io)
  pos = position(io)
  try
    f()
  finally
    seek(io, pos)
  end
end

end


text:
# this script is used for testing hover functions

include(pwd()*"/src/const.jl")
include(pwd()*"/src/mathfunctions.jl")
include(pwd()*"/src/aoaget.jl")
include(pwd()*"/src/clcdget.jl")
include(pwd()*"/src/rotorforce.jl")
include(pwd()*"/hovertest/hoverfunc.jl")

hovertmp = uitest2hover(2.*π/180,0.,2.*π/180,rb)
print("the result is : $(hovertmp)")

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
import Libdl

"""

This is a modified version of the Clang.jl generated code for the tree-sitter
API. The generated code is not entirely correct, so some modifications were
needed. Those edits are clarified with a comment to make future updates easy.
"""
const DEPS = "$(@__DIR__)/../deps"
const ts_lib = "$(DEPS)/libmylang-parser.$(Libdl.dlext)"
const C99bool_t = Cuchar

import CEnum: @cenum

######################################################### Common

const TREE_SITTER_LANGUAGE_VERSION = 9
const TSSymbol = UInt16
const TSLanguage = Cvoid
const TSParser = Cvoid
const TSTree = Cvoid

@cenum(TSInputEncoding,
    TSInputEncodingUTF8 = 0,
    TSInputEncodingUTF16 = 1,
)
@cenum(TSSymbolType,
    TSSymbolTypeRegular = 0,
    TSSymbolTypeAnonymous = 1,
    TSSymbolTypeAuxiliary = 2,
)

struct TSPoint
    row::UInt32
    column::UInt32
end

struct TSRange
    start_point::TSPoint
    end_point::TSPoint
    start_byte::UInt32
    end_byte::UInt32
end

struct TSInput
    payload::Ptr{Cvoid}
    read::Ptr{Cvoid}
    encoding::TSInputEncoding
end

@cenum(TSLogType,
    TSLogTypearse = 0,
    TSLogTypeLex = 1,
)

struct TSLogger
    payload::Ptr{Cvoid}
    log::Ptr{Cvoid}
end

struct TSInputEdit
    start_byte::UInt32
    old_end_byte::UInt32
    new_end_byte::UInt32
    start_point::TSPoint
    old_end_point::TSPoint
    new_end_point::TSPoint
end

struct TSNode
    context::NTuple{4, UInt32}
    id::Ptr{Cvoid}
    tree::Ptr{TSTree}
end

struct TSTreeCursor
    tree::Ptr{Cvoid}
    id::Ptr{Cvoid}
    context::NTuple{2, UInt32}
end

######################################################### API

function ts_parser_new()
    ccall((:ts_parser_new, ts_lib), Ptr{TSParser}, ())
end

function ts_parser_delete(arg1)
    ccall((:ts_parser_delete, ts_lib), Cvoid, (Ptr{TSParser},), arg1)
end

function ts_parser_language(arg1)
    ccall((:ts_parser_language, ts_lib), Ptr{TSLanguage}, (Ptr{TSParser},), arg1)
end

# EDIT somehow the arguments of set_language were missing
function ts_parser_set_language(p, l)
    ccall((:ts_parser_set_language, ts_lib), C99bool_t, (Ptr{TSParser}, Ptr{TSLanguage}), p, l)
end

function ts_parser_logger(arg1)
    ccall((:ts_parser_logger, ts_lib), TSLogger, (Ptr{TSParser},), arg1)
end

function ts_parser_set_logger(arg1, arg2)
    ccall((:ts_parser_set_logger, ts_lib), Cvoid, (Ptr{TSParser}, TSLogger), arg1, arg2)
end

function ts_parser_print_dot_graphs(arg1, arg2)
    ccall((:ts_parser_print_dot_graphs, ts_lib), Cvoid, (Ptr{TSParser}, Cint), arg1, arg2)
end

function ts_parser_halt_on_error(arg1, bool)
    ccall((:ts_parser_halt_on_error, ts_lib), Cvoid, (Ptr{TSParser}, Cint), arg1, bool)
end

function ts_parser_parse(arg1, arg2, arg3)
    ccall((:ts_parser_parse, ts_lib), Ptr{TSTree}, (Ptr{TSParser}, Ptr{TSTree}, TSInput), arg1, arg2, arg3)
end

function ts_parser_parse_string(arg1, arg2, arg3, arg4)
    ccall((:ts_parser_parse_string, ts_lib), Ptr{TSTree}, (Ptr{TSParser}, Ptr{TSTree}, Cstring, UInt32), arg1, arg2, arg3, arg4)
end

function ts_parser_parse_string_encoding(arg1, arg2, arg3, arg4, arg5)
    ccall((:ts_parser_parse_string_encoding, ts_lib), Ptr{TSTree}, (Ptr{TSParser}, Ptr{TSTree}, Cstring, UInt32, TSInputEncoding), arg1, arg2, arg3, arg4, arg5)
end

function ts_parser_enabled()
    ccall((:ts_parser_enabled, ts_lib), Cint, ())
end

function ts_parser_set_enabled(arg1, bool)
    ccall((:ts_parser_set_enabled, ts_lib), Cvoid, (Ptr{TSParser}, Cint), arg1, bool)
end

function ts_parser_operation_limit()
    ccall((:ts_parser_operation_limit, ts_lib), Cint, ())
end

function ts_parser_set_operation_limit(arg1, size_t)
    ccall((:ts_parser_set_operation_limit, ts_lib), Cvoid, (Ptr{TSParser}, Cint), arg1, size_t)
end

function ts_parser_reset(arg1)
    ccall((:ts_parser_reset, ts_lib), Cvoid, (Ptr{TSParser},), arg1)
end

function ts_parser_set_included_ranges(arg1, arg2, arg3)
    ccall((:ts_parser_set_included_ranges, ts_lib), Cvoid, (Ptr{TSParser}, Ptr{TSRange}, UInt32), arg1, arg2, arg3)
end

function ts_parser_included_ranges(arg1, arg2)
    ccall((:ts_parser_included_ranges, ts_lib), Ptr{TSRange}, (Ptr{TSParser}, Ptr{UInt32}), arg1, arg2)
end

function ts_tree_copy(arg1)
    ccall((:ts_tree_copy, ts_lib), Ptr{TSTree}, (Ptr{TSTree},), arg1)
end

function ts_tree_delete(arg1)
    ccall((:ts_tree_delete, ts_lib), Cvoid, (Ptr{TSTree},), arg1)
end

function ts_tree_root_node(arg1)
    ccall((:ts_tree_root_node, ts_lib), TSNode, (Ptr{TSTree},), arg1)
end

function ts_tree_edit(arg1, arg2)
    ccall((:ts_tree_edit, ts_lib), Cvoid, (Ptr{TSTree}, Ptr{TSInputEdit}), arg1, arg2)
end

function ts_tree_get_changed_ranges(arg1, arg2, arg3)
    ccall((:ts_tree_get_changed_ranges, ts_lib), Ptr{TSRange}, (Ptr{TSTree}, Ptr{TSTree}, Ptr{UInt32}), arg1, arg2, arg3)
end

function ts_tree_language(arg1)
    ccall((:ts_tree_language, ts_lib), Ptr{TSLanguage}, (Ptr{TSTree},), arg1)
end

function ts_node_start_byte(arg1)
    ccall((:ts_node_start_byte, ts_lib), UInt32, (TSNode,), arg1)
end

function ts_node_start_point(arg1)
    ccall((:ts_node_start_point, ts_lib), TSPoint, (TSNode,), arg1)
end

function ts_node_end_byte(arg1)
    ccall((:ts_node_end_byte, ts_lib), UInt32, (TSNode,), arg1)
end

function ts_node_end_point(arg1)
    ccall((:ts_node_end_point, ts_lib), TSPoint, (TSNode,), arg1)
end

function ts_node_symbol(arg1)
    ccall((:ts_node_symbol, ts_lib), TSSymbol, (TSNode,), arg1)
end

function ts_node_type(arg1)
    ccall((:ts_node_type, ts_lib), Cstring, (TSNode,), arg1)
end

function ts_node_string(arg1)
    ccall((:ts_node_string, ts_lib), Cstring, (TSNode,), arg1)
end

function ts_node_eq()
    ccall((:ts_node_eq, ts_lib), Cint, ())
end

function ts_node_is_null()
    ccall((:ts_node_is_null, ts_lib), Cint, ())
end

function ts_node_is_named()
    ccall((:ts_node_is_named, ts_lib), Cint, ())
end

# EDIT arguments were missing
function ts_node_is_missing(arg::TSNode)
    ccall((:ts_node_is_missing, ts_lib), C99bool_t, (TSNode,), arg)
end

# EDIT arguments were missing
function ts_node_is_extra(arg::TSNode)
    ccall((:ts_node_is_extra, ts_lib), C99bool_t, (TSNode,), arg)
end

# EDIT arguments were missing
function ts_node_is_ambiguous(arg)
    ccall((:ts_node_is_ambiguous, ts_lib), C99bool_t, (TSNode,), arg)
end

function ts_node_amb_count(arg1)
    ccall((:ts_node_amb_count, ts_lib), UInt32, (TSNode,), arg1)
end

function ts_node_amb_child_count(arg1, arg2)
    ccall((:ts_node_amb_child_count, ts_lib), UInt32, (TSNode, UInt32), arg1, arg2)
end

function ts_node_amb_child(arg1, arg2, arg3)
    ccall((:ts_node_amb_child, ts_lib), TSNode, (TSNode, UInt32, UInt32), arg1, arg2, arg3)
end

function ts_node_has_changes()
    ccall((:ts_node_has_changes, ts_lib), Cint, ())
end

# EDIT arguments were missing
function ts_node_has_error(arg)
    ccall((:ts_node_has_error, ts_lib), C99bool_t, (TSNode,), arg)
end

function ts_node_parent(arg1)
    ccall((:ts_node_parent, ts_lib), TSNode, (TSNode,), arg1)
end

function ts_node_child(arg1, arg2)
    ccall((:ts_node_child, ts_lib), TSNode, (TSNode, UInt32), arg1, arg2)
end

function ts_node_named_child(arg1, arg2)
    ccall((:ts_node_named_child, ts_lib), TSNode, (TSNode, UInt32), arg1, arg2)
end

function ts_node_child_count(arg1)
    ccall((:ts_node_child_count, ts_lib), UInt32, (TSNode,), arg1)
end

function ts_node_named_child_count(arg1)
    ccall((:ts_node_named_child_count, ts_lib), UInt32, (TSNode,), arg1)
end

function ts_node_next_sibling(arg1)
    ccall((:ts_node_next_sibling, ts_lib), TSNode, (TSNode,), arg1)
end

function ts_node_next_named_sibling(arg1)
    ccall((:ts_node_next_named_sibling, ts_lib), TSNode, (TSNode,), arg1)
end

function ts_node_prev_sibling(arg1)
    ccall((:ts_node_prev_sibling, ts_lib), TSNode, (TSNode,), arg1)
end

function ts_node_prev_named_sibling(arg1)
    ccall((:ts_node_prev_named_sibling, ts_lib), TSNode, (TSNode,), arg1)
end

function ts_node_first_child_for_byte(arg1, arg2)
    ccall((:ts_node_first_child_for_byte, ts_lib), TSNode, (TSNode, UInt32), arg1, arg2)
end

function ts_node_first_named_child_for_byte(arg1, arg2)
    ccall((:ts_node_first_named_child_for_byte, ts_lib), TSNode, (TSNode, UInt32), arg1, arg2)
end

function ts_node_descendant_for_byte_range(arg1, arg2, arg3)
    ccall((:ts_node_descendant_for_byte_range, ts_lib), TSNode, (TSNode, UInt32, UInt32), arg1, arg2, arg3)
end

function ts_node_named_descendant_for_byte_range(arg1, arg2, arg3)
    ccall((:ts_node_named_descendant_for_byte_range, ts_lib), TSNode, (TSNode, UInt32, UInt32), arg1, arg2, arg3)
end

function ts_node_descendant_for_point_range(arg1, arg2, arg3)
    ccall((:ts_node_descendant_for_point_range, ts_lib), TSNode, (TSNode, TSPoint, TSPoint), arg1, arg2, arg3)
end

function ts_node_named_descendant_for_point_range(arg1, arg2, arg3)
    ccall((:ts_node_named_descendant_for_point_range, ts_lib), TSNode, (TSNode, TSPoint, TSPoint), arg1, arg2, arg3)
end

function ts_node_edit(arg1, arg2)
    ccall((:ts_node_edit, ts_lib), Cvoid, (Ptr{TSNode}, Ptr{TSInputEdit}), arg1, arg2)
end

function ts_tree_cursor_new(arg1)
    ccall((:ts_tree_cursor_new, ts_lib), TSTreeCursor, (TSNode,), arg1)
end

function ts_tree_cursor_delete(arg1)
    ccall((:ts_tree_cursor_delete, ts_lib), Cvoid, (Ptr{TSTreeCursor},), arg1)
end

function ts_tree_cursor_reset(arg1, arg2)
    ccall((:ts_tree_cursor_reset, ts_lib), Cvoid, (Ptr{TSTreeCursor}, TSNode), arg1, arg2)
end

function ts_tree_cursor_current_node(arg1)
    ccall((:ts_tree_cursor_current_node, ts_lib), TSNode, (Ptr{TSTreeCursor},), arg1)
end

function ts_tree_cursor_goto_parent()
    ccall((:ts_tree_cursor_goto_parent, ts_lib), Cint, ())
end

function ts_tree_cursor_goto_next_sibling()
    ccall((:ts_tree_cursor_goto_next_sibling, ts_lib), Cint, ())
end

function ts_tree_cursor_goto_first_child()
    ccall((:ts_tree_cursor_goto_first_child, ts_lib), Cint, ())
end

function ts_tree_cursor_goto_first_child_for_byte(arg1, arg2)
    ccall((:ts_tree_cursor_goto_first_child_for_byte, ts_lib), Int64, (Ptr{TSTreeCursor}, UInt32), arg1, arg2)
end

function ts_language_symbol_count(arg1)
    ccall((:ts_language_symbol_count, ts_lib), UInt32, (Ptr{TSLanguage},), arg1)
end

function ts_language_symbol_name(arg1, arg2)
    ccall((:ts_language_symbol_name, ts_lib), Cstring, (Ptr{TSLanguage}, TSSymbol), arg1, arg2)
end

function ts_language_symbol_for_name(arg1, arg2)
    ccall((:ts_language_symbol_for_name, ts_lib), TSSymbol, (Ptr{TSLanguage}, Cstring), arg1, arg2)
end

function ts_language_symbol_type(arg1, arg2)
    ccall((:ts_language_symbol_type, ts_lib), TSSymbolType, (Ptr{TSLanguage}, TSSymbol), arg1, arg2)
end

function ts_language_version(arg1)
    ccall((:ts_language_version, ts_lib), UInt32, (Ptr{TSLanguage},), arg1)
end


text:
@noinline function _throw_conflicting_sgnum_and_bravais(bt, sgnum)
    throw(DomainError((bt, sgnum),
            "provided bravais type and space group number are mutually inconsistent"))
end
@noinline function _throw_triclinic_angles(Rs)
    throw(DomainError(Rs, "Triclinic Bravais lattice must be specificied with a "*
                          "lattice that is either all-acute or all-obtuse; provided "*
                          "lattice system is neither"))
end
@noinline function _throw_basis_required(Rs)
    throw(DomainError(Rs, "`Rs` must be supplied for the considered Bravais type "))
end

# --- 3D ---
# this is a translation and simplification (omitting several warnings / edge case handling)
# of SeeK's `get_path` extended bravais type branch table. It should be pretty much the
# same as Table 94 of the HPKOT paper (i.e. paper "behind" SeeK). It simply returns the
# extended Bravais symbol as a `Symbol`. Note the input direct basis is expected to
# be a conventional cell in ITA settings.
function extended_bravais(sgnum::Integer,
                          bt::String,
                          Rs::Union{Nothing, StaticVector{3, <:StaticVector{3, <:Real}}},
                          Dᵛ::Val{3})

    @boundscheck boundscheck_sgnum(sgnum, 3)

    if bt == "cP"
        if 195 ≤ sgnum ≤ 206
            return :cP1
        elseif 207 ≤ sgnum ≤ 230
            return :cP2
        else
            _throw_conflicting_sgnum_and_bravais(bt, sgnum)
        end

    elseif bt == "cF"
        if 195 ≤ sgnum ≤ 206
            return :cF1
        elseif 207 ≤ sgnum ≤ 230
            return :cF2
        else
            _throw_conflicting_sgnum_and_bravais(bt, sgnum)
        end

    elseif bt == "cI"
        return :cI1

    elseif bt == "tP"
        return :tP1

    elseif bt == "tI"
        a, _, c = basisnorms(Rs)
        if c ≤ a
            return :tI1
        else
            return :tI2
        end

    elseif bt == "oP"
        return :oP1

    elseif bt == "oF"
        a, b, c = basisnorms(Rs)
        if inv(a^2) > inv(b^2) + inv(c^2)
            return :oF1
        elseif inv(c^2) > inv(a^2) + inv(b^2)
            return :oF2
        else
            return :oF3
        end

    elseif bt == "oI"
        a, b, c = basisnorms(Rs)
        if c>a && c>b     # (c largest)
            return :oI1
        elseif a>b && a>c # (a largest)
            return :oI2
        else # b≥a && b≥c   (b largest, or equal)
            return :oI3
        end

    elseif bt == "oC"
        a, b, _ = basisnorms(Rs)
        if a ≤ b
            return :oC1
        else
            return :oC2
        end

    elseif bt == "oA"
        _, b, c = basisnorms(Rs)
        if b ≤ c
            return :oA1
        else
            return :oA2
        end

    elseif bt == "hP"
        if sgnum in (143:149..., 151, 153, 157, 159:163...)
            return :hP1
        else
            return :hP2
        end

    elseif bt == "hR"
        a, _, c = basisnorms(Rs)
        if sqrt(3)a ≤ sqrt(2)c
            return :hR1
        else
            return :hR2
        end

    elseif bt == "mP"
        return :mP1

    elseif bt == "mC"
        a, b, c = basisnorms(Rs)
        cosβ = dot(Rs[3], Rs[1])/(c*a)
        if b < a * sqrt(1 - cosβ^2)
            return :mC1
        else
            if -a*cosβ/c + a^2*(1 - cosβ^2) / b^2 ≤ 1
                return :mC2
            else
                return :mC3
            end
        end

    elseif bt == "aP"
        # SeeK-path only provides paths for all-acute or all-obtuse reciprocal aP cells
        # (more precisely, SeeK-path will transform to this setting automatically) so check:

        # get reciprocal basis inter-angles angles
        Gs = (Rs[2]×Rs[3], Rs[3]×Rs[1], Rs[1]×Rs[2]) # NB: "scale" prefactor omitted here
        aᴳ, bᴳ, cᴳ = basisnorms(Gs)
        cosαᴳ = dot(Gs[2], Gs[3]) / (bᴳ * cᴳ)
        cosβᴳ = dot(Gs[1], Gs[3]) / (aᴳ * cᴳ)
        cosγᴳ = dot(Gs[1], Gs[2]) / (aᴳ * bᴳ)

        # detect whether all-obtuse or all-acute reciprocal basis
        if cosαᴳ ≤ 0 && cosβᴳ ≤ 0 && cosγᴳ ≤ 0     # all-obtuse
            return :aP2
        elseif cosαᴳ ≥ 0 && cosβᴳ ≥ 0 && cosγᴳ ≥ 0 # all-acute
            return :aP3
        else
            _throw_triclinic_angles(Rs)
        end

    else
        throw(DomainError(bt, "undefined bravais type"))

    end
end # function

# --- 2D ---
function extended_bravais(sgnum::Integer,
                          bt::String,
                          Rs::Union{Nothing, StaticVector{2, <:StaticVector{2, <:Real}}},
                          Dᵛ::Val{2})
    @boundscheck boundscheck_sgnum(sgnum, 2)

    if bt == "hp"
        return :hp1
    elseif bt == "tp"
        return :tp1
    elseif bt == "op"
        if sgnum ∈ (3,4)
            return :op1
        elseif sgnum ∈ (6,7,8)
            return :op2
        else
            _throw_conflicting_sgnum_and_bravais(bt, sgnum)
        end
    elseif bt == "oc"
        a, b = basisnorms(Rs)
        if a ≤ b
            return :oc1
        else # a>b
            return :oc2
        end
    elseif bt == "mp"
        a, b = basisnorms(Rs)
        cosα = dot(Rs[1], Rs[2]) / (a * b)
        if cosα ≥ 0 # ∠(R₁, R₂) ≤ 90
            return :mp1
        else        # ∠(R₁, R₂) > 90
            return :mp2
        end
    else
        throw(DomainError(bt, "undefined bravais type"))
    end
end

# --- 1D ---
function extended_bravais(sgnum::Integer,
                          bt::String,
                          Rs::Union{Nothing, StaticVector{1, <:StaticVector{1, <:Real}}},
                          Dᵛ::Val{1})

    @boundscheck boundscheck_sgnum(sgnum, 1)
    return :lp # trivial case; no "extended" bravais types
end

function basisnorms(Rs::Union{AVec{V}, NTuple{3, V}}) where V <: StaticVector{3, <:Real}
    a = norm(Rs[1])
    b = norm(Rs[2])
    c = norm(Rs[3])
    return a, b, c
end
function basisnorms(Rs::Union{AVec{V}, NTuple{2, V}}) where V <: StaticVector{2, <:Real}
    a = norm(Rs[1])
    b = norm(Rs[2])
    return a, b
end
basisnorms(Rs::Nothing) = _throw_basis_required(Rs)


text:
using DifferenceLists
using Test

@test collect(dl(1, 2, 3)) == [1,2,3]
@test collect(dl()) == []
@test collect(dl(1)) == [1]
@test collect(dl(1, dl(2, 3), 4)) == [1, dl(2, 3), 4]
@test collect(todl([1, 2, 3])) == [1, 2, 3]
@test collect(push(push(dl(7, 8, 9), 1), 2)) == [7,8,9,1,2]
@test collect(pushfirst(pushfirst(dl(7, 8, 9), 2), 1)) == [1, 2, 7, 8, 9]
@test collect(concat(dl(1, 2), dl(3, 4))) == [1,2,3,4]
@test collect(dl(1, 2)(dl(3, 4), dl(5, 6, 7))) == [1, 2, 3, 4, 5, 6, 7]


text:
# ------------------------------------------------------------------------------------------
# # IMDB Training Notebook
# ------------------------------------------------------------------------------------------

using Knet
# Hyperparams LSTM
EPOCHS=3
BATCHSIZE=64
EMBEDSIZE=125
NUMHIDDEN=100
LR=0.0001
BETA_1=0.9
BETA_2=0.999
EPS=1e-08
MAXLEN=150 #maximum size of the word sequence
MAXFEATURES=30000 #vocabulary size
DROPOUT=0.35
SEED=1311194
gpu(0)
atype = gpu()<0 ? Array{Float32}:KnetArray{Float32}

#define model"
function initmodel()
    rnnSpec,rnnWeights = rnninit(EMBEDSIZE,NUMHIDDEN; rnnType=:lstm)
    inputMatrix = atype(xavier(Float32,EMBEDSIZE,MAXFEATURES))
    outputMatrix = atype(xavier(Float32,2,NUMHIDDEN))
    return rnnSpec,(rnnWeights,inputMatrix,outputMatrix)
end

function savemodel(weights,rnnSpec;localfile="model_imdb.jld")
    save(localfile,"weights",weights,"rnnSpec",rnnSpec)
end

# define loss and its gradient
function predict(weights, inputs, rnnSpec;train=false)
    rnnWeights, inputMatrix, outputMatrix = weights # (1,1,W), (X,V), (2,H)
    indices = hcat(inputs...)' # (B,T)
    rnnInput = inputMatrix[:,indices] # (X,B,T)
    if train
        rnnInput = dropout(rnnInput,DROPOUT)
    end
    rnnOutput = rnnforw(rnnSpec, rnnWeights, rnnInput)[1] # (H,B,T)
    if train
        rnnOutput = dropout(rnnOutput,DROPOUT)
    end
    return outputMatrix * rnnOutput[:,:,end] # (2,H) * (H,B) = (2,B)
end

loss(w,x,y,r;train=false)=nll(predict(w,x,r;train=train),y)
lossgradient = grad(loss);

# load data
include("imdb.jl")
@time (xtrn,ytrn,xtst,ytst,imdbdict)=imdb(maxlen=MAXLEN,maxval=MAXFEATURES,seed=SEED)
for d in (xtrn,ytrn,xtst,ytst); println(summary(d)); end
imdbarray = Array{String}(88584)
for (k,v) in imdbdict; imdbarray[v]=k; end

rnd = rand(1:length(xtrn))
println("Sample review:\n",join(imdbarray[xtrn[rnd]]," "),"\n")
println("Classification: ",join(ytrn[rnd]))

# prepare for training
weights = nothing; knetgc(); # Reclaim memory from previous run
rnnSpec,weights = initmodel()
optim = optimizers(weights, Adam; lr=LR, beta1=BETA_1, beta2=BETA_2, eps=EPS);

# 29s
info("Training...")
@time for epoch in 1:EPOCHS
    @time for (x,y) in minibatch(xtrn,ytrn,BATCHSIZE;shuffle=true)
        grads = lossgradient(weights,x,y,rnnSpec;train=true)
        update!(weights, grads, optim)
    end
end

info("Testing...")
@time accuracy(weights, minibatch(xtst,ytst,BATCHSIZE), (w,x)->predict(w,x,rnnSpec))

savemodel(weights,rnnSpec)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
#!/usr/bin/env julia

# Generate a set of model and residual visibilities and then write them to UVHDF5 file format
# More about this format, including scripts to convert to UVFITS and CASA Measurement set,
# can be found at https://github.com/Astrochem/UVHDF5

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
    "--out-model"
    help = "The file name to write out the model visibilities."
    default = "model.hdf5"
    "--out-resid"
    help = "The file name to write out the residual visibilities."
    default = "resid.hdf5"
end

parsed_args = parse_args(ARGS, s)

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.gridding

using HDF5

import YAML
config = YAML.load(open(parsed_args["config"]))

# The model writing is taken care of by `DJ_initialize.jl`, the RADMC-3D synthesis is taken care of by `DJ_plot_model.jl` and then this simply reads in the image, does FFT, downsamples, etc.

# read the wavelengths for all channels
fid = h5open(config["data_file"], "r")
nchan = length(read(fid["freqs"]))
close(fid)

# Read the parameters from the config file
pars = convert_dict(config["parameters"], config["model"])
grid = Grid(config["grid"])

# Mention the contribution of the prior to the lnprob
ln_prior = lnprior(pars, config["dpc_prior"]["mu"], config["dpc_prior"]["sig"], grid)

im = imread()
skim = imToSky(im, pars.dpc)
corrfun!(skim) # alpha = 1.0

# Determine dRA and dDEC from the image and distance
dRA = abs(skim.ra[2] - skim.ra[1])/2. # [arcsec] the half-size of a pixel
println("dRA is ", dRA, " arcsec")

# For *this purpose only*, read in the flagged data in addition to the unflagged data
# so that we can export a model for these otherwise flagged visibilities
dvarr = DataVis(config["data_file"], true)
# Do this as we do in `mach_three.jl`
for dset in dvarr
    # Conjugation is necessary for the SMA and ALMA
    visibilities.conj!(dset) # Swap UV convention
end

mvarr = Array(DataVis, nchan)
chi2s = Array(Float64, nchan)
lnprobs = Array(Float64, nchan)

for i=1:nchan
    dv = dvarr[i]

    # FFT the appropriate image channel
    vis_fft = transform(skim, i)

    # Interpolate the `vis_fft` to the same locations as the DataSet
    mvis = ModelVis(dv, vis_fft)

    # Apply the phase correction here, since there are fewer data points
    phase_shift!(mvis, pars.mu_RA + dRA, pars.mu_DEC - dRA)

    dvis = visibilities.ModelVis2DataVis(mvis)

    mvarr[i] = dvis

    chi2s[i] = visibilities.chi2(dv, mvis)
    lnprobs[i] = visibilities.lnprob(dv, mvis)
end

# Now generate the residual visibilities
rvarr = ResidVis(dvarr, mvarr)

# Now swap the model and residual visibilities back to ALMA/SMA convetion
for i=1:nchan
    visibilities.conj!(mvarr[i])
    visibilities.conj!(rvarr[i])
end

N = nchan * 2 * length(dvarr[1].VV)

# Only use the unmasked channels in the chi2 calculation
if haskey(config, "exclude")
    exclude = config["exclude"]
    # which channels of the dset to fit
    keylist = filter(x->(!in(x, exclude)), Int[i for i=1:nchan])
else
    keylist = Int[i for i=1:nchan]
end

println("Note: may include flagged visibilities!")
chi2s = chi2s[keylist]
# println("Chi^2 :", sum(chi2s))
println("Reduced Chi^2 ", sum(chi2s)/N)

# Calculate the lnprob between these two
lnprobs = lnprobs[keylist]

println("lnprior ", ln_prior)
println("lnlikelihood ", sum(lnprobs))
println("lnprob ", ln_prior + sum(lnprobs))

visibilities.write(mvarr, parsed_args["out-model"])
visibilities.write(rvarr, parsed_args["out-resid"])

# copy visibility flags from the original dataset to this new one
visibilities.copy_flags(config["data_file"], parsed_args["out-model"])
visibilities.copy_flags(config["data_file"], parsed_args["out-resid"])

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
module TestGlobalVariables
using AutoImports: global_variables
using Test

simple_code = """
using Somewhere: E
G = 1

\"\"\"
docstring of f
\"\"\"
function f(x)
    x + ((z) -> Y + z)(x)
end

function g(x)
    x
end

M.h(x) = nothing
"""

@testset "$label" for (label, code) in [
    ("simple", simple_code),
    ("module", """
               module MyModule
               $simple_code
               end
               """),
]
    vars = global_variables(code)

    if occursin("MyModule", code)
        @test vars.left == Set([:MyModule, :G, :E, :f, :g])
    else
        @test vars.left == Set([:G, :E, :f, :g])
    end

    @test Set([:Y, :M]) <= vars.right

    @testset for x in [:h, :x, :z]
        @test x ∉ vars.right
        @test x ∉ vars.left
    end
end

end  # module


text:
__precompile__()
module LittleScienceTools
using Statistics
using LinearAlgebra
using Printf
export Roots, Ising, Measuring, Random
include("Measuring/Measuring.jl")
# include("Roots/Roots.jl")
# include("Random/Random.jl")
# include("Ising/Ising.jl")

export interwine
include("utils.jl")

end # module LittleScienceTools

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
@testset "cluster.jl" begin
    @testset "constructor" begin
        # assertion error
        @test_throws AssertionError Cluster(10,11,[51],[85])
        @test_throws AssertionError Cluster(10,0,[51],[85; 85])

        # Test with known values
        c = Cluster(10,0,[51],[85])
        @test c.u₁ ≈ 10
        @test c.u₂ ≈ 0
        @test c.position ≈ [51]
        @test c.value ≈ [85]
    end

    @testset "getcluster(y, u₁, u₂)" begin

        # Test with known values
        y = zeros(Int64,10)
        y[1:3] .= 1
        y[5] = 1
        y[9:10] .= 1

        # assertion error
        @test_throws AssertionError getcluster(y,.5,1.0)

        cluster = getcluster(y,.5,0.0)

        @test cluster[1].position == collect(1:3)
        @test cluster[1].value == ones(Int64,3)

        @test cluster[2].position == collect(5:5)
        @test cluster[2].value == ones(Int64,1)

        @test cluster[3].position == collect(9:10)
        @test cluster[3].value == ones(Int64,2)
    end

    @testset "getcluster(y, u)" begin

        # Test with known values
        y = zeros(Int64,10)
        y[1:3] .= 1
        y[5] = 1
        y[9:10] .= 1

        cluster = getcluster(y,.5)

        @test cluster[1].position == collect(1:3)
        @test cluster[1].value == ones(Int64,3)

        @test cluster[2].position == collect(5:5)
        @test cluster[2].value == ones(Int64,1)

        @test cluster[3].position == collect(9:10)
        @test cluster[3].value == ones(Int64,2)
    end

    @testset "getcluster(y, u, runlength=r)" begin

        # Test with known values
        y = zeros(Int64,10)
        y[1:3] .= 1
        y[5] = 1
        y[9:10] .= 1

        cluster = getcluster(y,.5, runlength=1)

        @test cluster[1].position == [1;2;3]
        @test cluster[2].position == [5]
        @test cluster[3].position == [9;10]

        cluster = getcluster(y,.5, runlength=2)

        @test cluster[1].position == [1;2;3;5]
        @test cluster[2].position == [9;10]

        cluster = getcluster(y,.5, runlength=4)

        @test cluster[1].position == [1;2;3;5;9;10]

        cluster = getcluster(y[1:4],.5, runlength=4)
        @test cluster[1].position == [1;2;3]

        cluster = getcluster(y[6:7],.5, runlength=4)
        @test isempty(cluster)

    end

    @testset "length(c::Cluster)" begin

        # Test with known values
        y = zeros(Int64,10)
        y[1:3] .= 1
        y[5] = 1
        y[9:10] .= 1

        cluster = getcluster(y,.5)

        @test length.(cluster) == [3; 1; 2]
    end

    @testset "max(c::Cluster)" begin

        # Test with known values
        y = zeros(Int64,10)
        y[1] = 5
        y[2] = 2
        y[3] = 1
        y[5] = 3
        y[9:10] .= 1

        cluster = getcluster(y,.5)

        @test maximum.(cluster) == [5; 3; 1]
    end

    @testset "sum(c::Cluster)" begin

        # Test with known values
        y = zeros(Int64,10)
        y[1:3] .= 1
        y[5] = 1
        y[9:10] .= 1

        cluster = getcluster(y,.5)

        @test sum.(cluster) == [3; 1; 2]
    end

    @testset "merge(c₁::Cluster, c₂::Cluster)" begin

        # Test assertion
        c₁ = Cluster(1,0,[10],[100])
        c₂ = Cluster(0,0,[10],[100])
        c₃ = Cluster(1,1,[10],[100])

        @test_throws AssertionError Extremes.merge(c₁, c₂)
        @test_throws AssertionError Extremes.merge(c₁, c₃)

        # Test with known values
        y = zeros(Int64,10)
        y[1:3] .= 1
        y[5] = 1
        y[9:10] .= 1

        cluster = getcluster(y,.5)

        c = Extremes.merge(cluster[1], cluster[2])

        @test c.u₁ ≈ cluster[1].u₁
        @test c.u₂ ≈ cluster[1].u₁
        @test c.position == [1;2;3;5]
        @test c.value == [1;1;1;1]
    end

end


text:
function  poisson(theta, y, x)
	lambda = exp.(x*theta)
	logdensity = y .* (x*theta) .- lambda .- lgamma.(y.+1)
end