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

using Markdown
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.peek, el) ? Base.peek(el) : missing
        el
    end
end

# ╔═╡ 77489ff6-bd82-11ea-2d09-8737634cc579
using Random

# ╔═╡ e9b93d52-bd82-11ea-15c9-699eeebaa243
using Plots

# ╔═╡ 7eab739e-bdba-11ea-3ff3-fd88b08946d3
using LaTeXStrings

# ╔═╡ 508fed7e-bdbe-11ea-1d75-8935a47d06d8
using Statistics

# ╔═╡ fc08751c-bd79-11ea-072b-33a239fa0605
md"# Stochastic Thinking"

# ╔═╡ f895f5da-bd7c-11ea-2250-efdcd9d49a38
md"## Modelling the recovery process using probability
Recall code: $I_{n+1} = (1 + c) I_n = \lambda * I_n$"

# ╔═╡ 65f53264-bd7d-11ea-03cc-9513654c222c
function run_infection(I₀, λ, T=20)
	I = zeros(T)
	I[1] = I₀

	for n ∈ 1:T-1
		I[n+1] = λ * I[n]
	end

	return I
end

# ╔═╡ 2dce78b0-bd7e-11ea-0381-d16883397ebd
run_infection(1.0, 0.1)

# ╔═╡ 38ad903e-bd7e-11ea-08b6-f9a9a086260f
md"## Adding new data to extend a vector"

# ╔═╡ c7c0a170-bd7f-11ea-3b68-a9ae0c8ccfde
v = [1.0]

# ╔═╡ e40b10e0-bd7f-11ea-2883-bf485b88d6ad
push!(v, 7.0) # ! means the function modifies its first argument

# ╔═╡ b9db886e-bd80-11ea-1347-5d5415a22dd5
v

# ╔═╡ c427a4fc-bd80-11ea-365a-59fd43802b47
[v; 10]

# ╔═╡ 92666d14-bd81-11ea-2ba7-e181ffa9f476
pushfirst!(v, 8) # deque - double endedn queue

# ╔═╡ 0f3c4b1c-bd80-11ea-0fd0-df1c82c3ed1c
function run_infection_vectorized(I₀, λ, T=20)
	Is = [I₀]
	I = I₀

	for n ∈ 1:T-1
		I_next = λ * I

		push!(Is, I_next)

		I = I_next
	end

	return Is
end

# ╔═╡ 71759792-bd81-11ea-1cbe-35ca8456b6a1
run_infection_vectorized(1.0, 1.1)

# ╔═╡ 7a610d46-bd81-11ea-0ff9-5b97563bce1e
md"## Randomness"

# ╔═╡ 29e900e6-bd82-11ea-0344-2b2276b99422
rand()

# ╔═╡ 8cec500a-bd82-11ea-2fbb-03dec9457d7e
Random.seed!(3);

# ╔═╡ a1c76a96-bd82-11ea-22a0-d9110840ecf3
rand()

# ╔═╡ d39e0c64-bd82-11ea-188d-6f762d4dd688
rand(10)

# ╔═╡ ed82a2b8-bd82-11ea-1d5c-8759b29ce1cc
scatter(rand(50))

# ╔═╡ 0682ca0c-bd83-11ea-0afa-8794eb7bd9b8
begin
	local r = rand(50)

	scatter(r, 0.5 * one.(r), ylim=(0, 1))
end

# ╔═╡ be2d6b30-bd83-11ea-2683-49e64cbccb3c
@bind num_points html"<input type='range' min=1 max=100>"

# ╔═╡ 811488a0-bdab-11ea-0f9d-b1bc9f854e07
#TODO coloring exercise
begin
	local r = rand(num_points)

	scatter(r[1:num_points],0.5 .* one.(r[1:num_points]),ylim=(0, 1),xlim=(0, 1))
end

# ╔═╡ 89984d06-bdab-11ea-3ff4-55786243da18
md"## Event with probability $p$"

# ╔═╡ a8196c14-bdac-11ea-3b3d-f744310c289d
begin
	local p = 0.25
	local r = rand()

	if r < p
		true
	else
		false
	end
end

# ╔═╡ cfac1f86-bdac-11ea-316c-a9abea7ed509
md"Bernoulli Trial"

# ╔═╡ 46de61a4-bdad-11ea-3aaf-a5f9b96e533f
function bernoulli(p)
	r = rand()

	return r < p
end

# ╔═╡ 8a20f6b6-bdad-11ea-0ef5-af788315bd82
2 < 3 ? true : false # ternary operator

# ╔═╡ bc82b248-bdad-11ea-3ad5-019463cf2ef1
2 < 3

# ╔═╡ c7552f0e-bdad-11ea-3fe8-fd3cf143ff30
begin
	local p = 0.25

	trials = [bernoulli(p) for i ∈ 1:100]
end

# ╔═╡ fe9569d2-bdad-11ea-0f72-09084f78f692
scatter(trials)

# ╔═╡ 38c6f580-bdae-11ea-2cbe-935e598ee04c
trials

# ╔═╡ 3a406aea-bdae-11ea-0c69-b13a4340a01a
count(trials)

# ╔═╡ 1b778d18-bdaf-11ea-3ee4-0f3a16536b9f
function bernoulli_experiment(p, N=100)
	trials = [bernoulli(p) for i ∈ 1:N]
	return count(trials)
end

# ╔═╡ 45a5c6a4-bdaf-11ea-1812-359502135397
count(trials .== false) + count(trials)

# ╔═╡ 65bf3b8a-bdaf-11ea-1878-bb69724e1a7b
count(.!(trials)) # ! is "not" in Julia

# ╔═╡ 84062240-bdaf-11ea-3494-8d49391f642a
bernoulli_experiment(0.25)

# ╔═╡ a9b7c836-bdaf-11ea-1691-9d38af43f16a
bernoulli_experiment(0.25)

# ╔═╡ be407118-bdaf-11ea-2973-853fc1a266e1
md"A quantity that has different **outcomes** in different runs is called a **random variable**"

# ╔═╡ de5398ea-bdaf-11ea-1d81-5df4eb22e68b
md"## Monte Carlo simulation"

# ╔═╡ 73e79708-bdb0-11ea-123f-9fbb2c18dc12
md"Run the same random process a lot of times and look at the results"

# ╔═╡ 85f366ac-bdb0-11ea-16b7-2bec6b2f5356
begin
	local p = 0.25
	local N = 20 # number of trials
	local num_expts = 1000

	results = [bernoulli_experiment(p, N) for i in 1:num_expts]
end

# ╔═╡ af9e25e6-bdb0-11ea-23aa-7db6e35a4447
scatter(results, zcolor=results)

# ╔═╡ 87eaac94-bdb1-11ea-117a-ffa693e0b098
md" $X$ = number of heads

**Probability distribution** of a random variable $X$

$\text{Prob}(X = x)$
"

# ╔═╡ 976c1838-bdb1-11ea-30b1-f5c5eb8e2360
md"Probability that $X = x$ = proportion of time that the result was $x$"

# ╔═╡ f411c49a-bdb2-11ea-00b9-dfcc68181d37
md"## Counting"

# ╔═╡ 0922211a-bdb3-11ea-1e49-3518de406ace
md"Need to count how many of the trials come out as 5"

# ╔═╡ 286acb30-bdb3-11ea-1be2-c51955c53ab8
md"Need data structure where we store the counts:

- `Dict` (dictionary)
- `Vector`"

# ╔═╡ 50d0c656-bdb3-11ea-0a32-f7d63d08f64a
results

# ╔═╡ 99aab29c-bdb3-11ea-3dd0-57d1e619f762
maximum(results) #not 'max': max(1, 2, 3)

# ╔═╡ 2beb3cb0-bdb6-11ea-1451-05c1e3087165
minimum(results)

# ╔═╡ 6e8a4d8a-bdb7-11ea-0fd5-6d942ceb5283
l = maximum(results) + 1

# ╔═╡ 8e80e7fc-bdb7-11ea-1967-bbff2030d913
counts = zeros(Int, l)

# ╔═╡ 9b648744-bdb7-11ea-3e0e-7307a4084f9a
for score ∈ results
	counts[score + 1] += 1
end

# ╔═╡ e2635850-bdb7-11ea-21bc-5beadadfa675
counts

# ╔═╡ 7aecb5b0-bdba-11ea-026f-09ebc21a36a6
import Pkg; Pkg.add("LaTeXStrings")

# ╔═╡ 96b1bbb2-bdb8-11ea-3e81-1bca33878a99
begin
	plot(0:maximum(results), counts, m=:o)
	ylabel!("frequency of having n heads")
	xlabel!("n")
	#xlabel!("n")
end

# ╔═╡ 7090316e-bdba-11ea-27a6-dbc5ec194d02
function count_them(results)
	counts = zeros(Int, maximum(results) + 1)

	for score ∈ results
		counts[score + 1] += 1
	end

	return counts
end

# ╔═╡ c56f3320-bdbb-11ea-2904-6197a97efbbf
function run_experiments(p, N, num_expts=1000)
	results = [bernoulli_experiment(p, N) for i ∈ 1:num_expts]

	return results
end

# ╔═╡ 4c85f43c-bdbc-11ea-1ba1-810cadb75f6e
begin
	local data = count_them(run_experiments(0.25, 20, 10000))
	plot(0:length(data)-1, data, m=:o)
end

# ╔═╡ 8c2374a2-bdbc-11ea-12a9-eda80a7e4e21
md"## Probability distribution"

# ╔═╡ a00fa58a-bdbc-11ea-14a2-fd43cfedf7bb
md"Probability = relative frequency

Divide by the number of experiments"

# ╔═╡ bc8cf4ec-bdbc-11ea-33f3-d9aa8f26cd60
begin
	local data = run_experiments(0.25, 20, 10^5) |> count_them
	plot(0:length(data)-1, data ./ sum(data), m=:o)
end

# ╔═╡ 0839ca58-bdbd-11ea-2a3e-c99ad45c7121
begin
	num_expts = 10^5
	data = run_experiments(0.25, 20, num_expts) |> count_them
	probs = data .// num_expts
end

# ╔═╡ cb4d967a-bdbd-11ea-35c8-ad831d395195
sum(probs)

# ╔═╡ cf802b2c-bdbd-11ea-0d10-390f5ea9fe19
x = 3 // 4

# ╔═╡ fb862028-bdbd-11ea-3c26-e5e8ebd9630d
typeof(x)

# ╔═╡ 04180cb0-bdbe-11ea-2a64-a33a4f778d68
x + x

# ╔═╡ 08a5bf84-bdbe-11ea-18d4-dd7255a35b56
data

# ╔═╡ 741484a8-bdbe-11ea-3fd8-e55fe46723f8
begin
	local results = run_experiments(0.25, 20, 10^5)
	mean(results)
end

# ╔═╡ 5d7c00b8-bdbe-11ea-2b46-6d1117a42368
md"Expected mean value = N * p"

# ╔═╡ Cell order:
# ╟─fc08751c-bd79-11ea-072b-33a239fa0605
# ╟─f895f5da-bd7c-11ea-2250-efdcd9d49a38
# ╠═65f53264-bd7d-11ea-03cc-9513654c222c
# ╠═2dce78b0-bd7e-11ea-0381-d16883397ebd
# ╟─38ad903e-bd7e-11ea-08b6-f9a9a086260f
# ╠═c7c0a170-bd7f-11ea-3b68-a9ae0c8ccfde
# ╠═e40b10e0-bd7f-11ea-2883-bf485b88d6ad
# ╠═b9db886e-bd80-11ea-1347-5d5415a22dd5
# ╠═c427a4fc-bd80-11ea-365a-59fd43802b47
# ╠═92666d14-bd81-11ea-2ba7-e181ffa9f476
# ╠═0f3c4b1c-bd80-11ea-0fd0-df1c82c3ed1c
# ╠═71759792-bd81-11ea-1cbe-35ca8456b6a1
# ╟─7a610d46-bd81-11ea-0ff9-5b97563bce1e
# ╠═29e900e6-bd82-11ea-0344-2b2276b99422
# ╠═77489ff6-bd82-11ea-2d09-8737634cc579
# ╠═8cec500a-bd82-11ea-2fbb-03dec9457d7e
# ╠═a1c76a96-bd82-11ea-22a0-d9110840ecf3
# ╠═d39e0c64-bd82-11ea-188d-6f762d4dd688
# ╠═e9b93d52-bd82-11ea-15c9-699eeebaa243
# ╠═ed82a2b8-bd82-11ea-1d5c-8759b29ce1cc
# ╠═0682ca0c-bd83-11ea-0afa-8794eb7bd9b8
# ╟─be2d6b30-bd83-11ea-2683-49e64cbccb3c
# ╠═811488a0-bdab-11ea-0f9d-b1bc9f854e07
# ╟─89984d06-bdab-11ea-3ff4-55786243da18
# ╠═a8196c14-bdac-11ea-3b3d-f744310c289d
# ╟─cfac1f86-bdac-11ea-316c-a9abea7ed509
# ╠═46de61a4-bdad-11ea-3aaf-a5f9b96e533f
# ╠═8a20f6b6-bdad-11ea-0ef5-af788315bd82
# ╠═bc82b248-bdad-11ea-3ad5-019463cf2ef1
# ╠═c7552f0e-bdad-11ea-3fe8-fd3cf143ff30
# ╠═fe9569d2-bdad-11ea-0f72-09084f78f692
# ╠═38c6f580-bdae-11ea-2cbe-935e598ee04c
# ╠═3a406aea-bdae-11ea-0c69-b13a4340a01a
# ╠═1b778d18-bdaf-11ea-3ee4-0f3a16536b9f
# ╠═45a5c6a4-bdaf-11ea-1812-359502135397
# ╠═65bf3b8a-bdaf-11ea-1878-bb69724e1a7b
# ╠═84062240-bdaf-11ea-3494-8d49391f642a
# ╠═a9b7c836-bdaf-11ea-1691-9d38af43f16a
# ╟─be407118-bdaf-11ea-2973-853fc1a266e1
# ╟─de5398ea-bdaf-11ea-1d81-5df4eb22e68b
# ╟─73e79708-bdb0-11ea-123f-9fbb2c18dc12
# ╠═85f366ac-bdb0-11ea-16b7-2bec6b2f5356
# ╠═af9e25e6-bdb0-11ea-23aa-7db6e35a4447
# ╟─87eaac94-bdb1-11ea-117a-ffa693e0b098
# ╟─976c1838-bdb1-11ea-30b1-f5c5eb8e2360
# ╟─f411c49a-bdb2-11ea-00b9-dfcc68181d37
# ╟─0922211a-bdb3-11ea-1e49-3518de406ace
# ╟─286acb30-bdb3-11ea-1be2-c51955c53ab8
# ╠═50d0c656-bdb3-11ea-0a32-f7d63d08f64a
# ╠═99aab29c-bdb3-11ea-3dd0-57d1e619f762
# ╠═2beb3cb0-bdb6-11ea-1451-05c1e3087165
# ╠═6e8a4d8a-bdb7-11ea-0fd5-6d942ceb5283
# ╠═8e80e7fc-bdb7-11ea-1967-bbff2030d913
# ╠═9b648744-bdb7-11ea-3e0e-7307a4084f9a
# ╠═e2635850-bdb7-11ea-21bc-5beadadfa675
# ╠═7aecb5b0-bdba-11ea-026f-09ebc21a36a6
# ╠═7eab739e-bdba-11ea-3ff3-fd88b08946d3
# ╠═96b1bbb2-bdb8-11ea-3e81-1bca33878a99
# ╠═7090316e-bdba-11ea-27a6-dbc5ec194d02
# ╠═c56f3320-bdbb-11ea-2904-6197a97efbbf
# ╠═4c85f43c-bdbc-11ea-1ba1-810cadb75f6e
# ╟─8c2374a2-bdbc-11ea-12a9-eda80a7e4e21
# ╟─a00fa58a-bdbc-11ea-14a2-fd43cfedf7bb
# ╠═bc8cf4ec-bdbc-11ea-33f3-d9aa8f26cd60
# ╠═0839ca58-bdbd-11ea-2a3e-c99ad45c7121
# ╠═cb4d967a-bdbd-11ea-35c8-ad831d395195
# ╠═cf802b2c-bdbd-11ea-0d10-390f5ea9fe19
# ╠═fb862028-bdbd-11ea-3c26-e5e8ebd9630d
# ╠═04180cb0-bdbe-11ea-2a64-a33a4f778d68
# ╠═08a5bf84-bdbe-11ea-18d4-dd7255a35b56
# ╠═508fed7e-bdbe-11ea-1d75-8935a47d06d8
# ╠═741484a8-bdbe-11ea-3fd8-e55fe46723f8
# ╟─5d7c00b8-bdbe-11ea-2b46-6d1117a42368

max_stars_repo_name: elihugarret/my-notebooks


text:
abstract type AggOperation <: Integer end
Base.Int32(d::AggOperation) = Int32(Integer(d))
struct Mean <: AggOperation end
struct Median <: AggOperation end
struct Each <: AggOperation end
Base.Integer(::Mean) = 0
Base.Integer(::Median) = 1
Base.Integer(::Each) = 2
const MEAN = Mean()
const MEDIAN = Median()
const EACH = Each()

abstract type FilterType <: Integer end
Base.Int32(d::FilterType) = Int32(Integer(d))
struct Butterworth <: FilterType end
struct Chebyshev <: FilterType end
struct Bessel <: FilterType end
Base.Integer(::Butterworth) = 0
Base.Integer(::Chebyshev) = 1
Base.Integer(::Bessel) = 2
const BUTTERWORTH = Butterworth()
const CHEBYSHEV_TYPE_1 = Chebyshev()
const BESSEL = Bessel()

abstract type WindowFunction <: Integer end
Base.Int32(d::WindowFunction) = Int32(Integer(d))
struct NoWindow <: WindowFunction end
struct Hanning <: WindowFunction end
struct Hamming <: WindowFunction end
struct BlackmanHarris <: WindowFunction end
Base.Integer(::NoWindow) = 0
Base.Integer(::Hanning) = 1
Base.Integer(::Hamming) = 2
Base.Integer(::BlackmanHarris) = 3
const NO_WINDOW = NoWindow()
const HANNING = Hanning()
const HAMMING = Hamming()
const BLACKMAN_HARRIS = BlackmanHarris()

abstract type DetrendOperation <: Integer end
Base.Int32(d::DetrendOperation) = Int32(Integer(d))
struct Constant <: DetrendOperation end
struct Linear <: DetrendOperation end
Base.Integer(::Constant) = 1
Base.Integer(::Linear) = 2
const CONSTANT = Constant()
const LINEAR = Linear()

abstract type NoiseType <: Integer end
Base.Int32(d::NoiseType) = Int32(Integer(d))
struct Fifty <: NoiseType end
struct Sixty <: NoiseType end
Base.Integer(::Fifty) = 1
Base.Integer(::Sixty) = 2
const FIFTY = Fifty()
const SIXTY = Sixty()

@brainflow_rethrow function perform_lowpass(data, sampling_rate::Integer, cutoff::Float64, order::Integer,
    filter_type::Integer, ripple::Float64)
    ccall((:perform_lowpass, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    return
end

@brainflow_rethrow function perform_highpass(data, sampling_rate::Integer, cutoff::Float64, order::Integer,
    filter_type::Integer, ripple::Float64
)
    ccall((:perform_highpass, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    return
end

@brainflow_rethrow function perform_bandpass(data, sampling_rate::Integer, center_freq::Float64,
    band_width::Float64, order::Integer, filter_type::Integer, ripple::Float64)
    ccall((:perform_bandpass, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    return
end

@brainflow_rethrow function perform_bandstop(data, sampling_rate::Integer, center_freq::Float64,
    band_width::Float64, order::Integer, filter_type::Integer, ripple::Float64)
    ccall((:perform_bandstop, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    return
end

@brainflow_rethrow function remove_environmental_noise(data, sampling_rate::Integer, noise_type::Integer)
    ccall((:remove_environmental_noise, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Cint),
            data, length(data), Int32(sampling_rate), Int32(noise_type))
    return
end

@brainflow_rethrow function perform_rolling_filter(data, period::Integer, operation::Integer)
    ccall((:perform_rolling_filter, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Cint),
            data, length(data), Int32(period), Int32(operation))
    return
end

@brainflow_rethrow function detrend(data, operation)
    ccall((:detrend, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint),
            data, length(data), Int32(operation))
    return
end

@brainflow_rethrow function perform_wavelet_denoising(data, wavelet::String, decomposition_level::Integer)
    ccall((:perform_wavelet_denoising, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
            data, length(data), wavelet, Int32(decomposition_level))
    return
end

@brainflow_rethrow function perform_downsampling(data, period::Integer, operation::Integer)
    len = Integer(floor(length(data) / period))
    downsampled_data = Vector{Float64}(undef, len)
    ccall((:perform_downsampling, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
            data, length(data), Int32(period), Int32(operation), downsampled_data)
    return downsampled_data
end

@brainflow_rethrow function write_file(data, file_name::String, file_mode::String)
    shape = size(data)
    flatten = transpose(vcat(data))
    flaten = reshape(flatten, (1, shape[1] * shape[2]))
    flatten = copy(flatten)
    ccall((:write_file, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
            flatten, shape[1], shape[2], file_name, file_mode)
    return
end

@brainflow_rethrow function get_nearest_power_of_two(value::Integer)
    power_of_two = Vector{Cint}(undef, 1)
    ccall((:get_nearest_power_of_two, DATA_HANDLER_INTERFACE), Cint, (Cint, Ptr{Cint}),
                Int32(value), power_of_two)
    return power_of_two[1]
end

@brainflow_rethrow function get_num_elements_in_file(file_name::String)
    num_elements = Vector{Cint}(undef, 1)
    ccall((:get_num_elements_in_file, DATA_HANDLER_INTERFACE), Cint, (Ptr{UInt8}, Ptr{Cint}),
                file_name, num_elements)
    return num_elements
end

@brainflow_rethrow function read_file(file_name::String)
    num_elements = get_num_elements_in_file(file_name)
    data = Vector{Float64}(undef, num_elements[1])
    num_cols = Vector{Cint}(undef, 1)
    num_rows = Vector{Cint}(undef, 1)
    ccall((:read_file, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
                data, num_rows, num_cols, file_name, num_elements[1])
    value = transpose(reshape(data, (num_cols[1], num_rows[1])))
    return value
end

@brainflow_rethrow function perform_wavelet_transform(data, wavelet::String, decomposition_level::Integer)
    wavelet_coeffs = Vector{Float64}(undef, length(data) + 2 * (40 + 1))
    lengths = Vector{Cint}(undef, decomposition_level + 1)
    ccall((:perform_wavelet_transform, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
            data, length(data), wavelet, Int32(decomposition_level), wavelet_coeffs, lengths)
    return wavelet_coeffs[1:sum(lengths)], lengths
end

@brainflow_rethrow function perform_inverse_wavelet_transform(wavelet_output, original_data_len::Integer, wavelet::String, decomposition_level::Integer)
    original_data = Vector{Float64}(undef, original_data_len)
    ccall((:perform_inverse_wavelet_transform, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
            wavelet_output[1], Int32(original_data_len), wavelet, Int32(decomposition_level), wavelet_output[2], original_data)
    return original_data
end

@brainflow_rethrow function get_csp(data, labels)
    n_epochs = size(data, 3)
    n_channels = size(data, 1)
    n_times = size(data, 2)

    temp_data1d = Vector{Float64}(undef, Integer(n_epochs * n_channels * n_times))
    for e=1:n_epochs
        for c=1:n_channels
            for t=1:n_times
                temp_data1d[(e-1) * n_channels * n_times + (c-1) * n_times + t] = data[c, t, e]
            end
        end
    end

    temp_filters = Vector{Float64}(undef, Integer(n_channels * n_channels))
    output_eigenvalues = Vector{Float64}(undef, Integer(n_channels))

    ccall((:get_csp, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}), temp_data1d, labels, Int32(n_epochs), Int32(n_channels), Int32(n_times), temp_filters, output_eigenvalues)

    output_filters = Array{Float64,2}(undef, n_channels, n_channels)
    for i=1:n_channels
        for j=1:n_channels
            output_filters[i, j] = temp_filters[(i-1) * n_channels + j]
        end
    end
    return output_filters, output_eigenvalues
end

@brainflow_rethrow function get_window(window_function::Integer, window_len::Integer)
    window_data = Vector{Float64}(undef, Integer(window_len))
    ccall((:get_window, DATA_HANDLER_INTERFACE), Cint, (Cint, Cint, Ptr{Float64}),
    Int32(window_function), Int32(window_len), window_data)
    return window_data
end

@brainflow_rethrow function perform_fft(data, window::Integer)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(length(data))
        throw(BrainFlowError(string("Data Len must be power of two ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
    end

    temp_re = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    temp_im = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    res = Vector{Complex}(undef, Integer(length(data) / 2) + 1)

    ccall((:perform_fft, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(window), temp_re, temp_im)
    for i in 1:Integer(length(data) / 2) + 1
        res[i] = Complex(temp_re[i], temp_im[i])
    end
    return res
end

@brainflow_rethrow function perform_ifft(data)

    temp_re = Vector{Float64}(undef, length(data))
    temp_im = Vector{Float64}(undef, length(data))
    res = Vector{Float64}(undef, 2 * (length(data) - 1))

    for i in 1:length(data)
        temp_re[i] = data[i].re
        temp_im[i] = data[i].im
    end

    ccall((:perform_ifft, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
            temp_re, temp_im, length(res), res)
    return res
end

@brainflow_rethrow function get_avg_band_powers(data, channels, sampling_rate::Integer, apply_filter::Bool)

    shape = size(data)
    data_1d = reshape(transpose(data[channels,:]), (1, size(channels)[1] * shape[2]))
    data_1d = copy(data_1d)

    temp_avgs = Vector{Float64}(undef, 5)
    temp_stddevs = Vector{Float64}(undef, 5)

    ccall((:get_avg_band_powers, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data_1d, size(channels)[1], shape[2], Int32(sampling_rate), Int32(apply_filter), temp_avgs, temp_stddevs)
    return temp_avgs, temp_stddevs
end

@brainflow_rethrow function get_psd(data, sampling_rate::Integer, window::Integer)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(length(data))
        throw(BrainFlowError(string("Data Len must be power of two ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
    end

    temp_ampls = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    temp_freqs = Vector{Float64}(undef, Integer(length(data) / 2) + 1)

    ccall((:get_psd, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    return temp_ampls, temp_freqs
end

@brainflow_rethrow function get_psd_welch(data, nfft::Integer, overlap::Integer, sampling_rate::Integer, window::Integer)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(nfft)
        throw(BrainFlowError(string("nfft must be power of two ", INVALID_ARGUMENTS_ERROR), Integer(INVALID_ARGUMENTS_ERROR)))
    end

    temp_ampls = Vector{Float64}(undef, Integer(nfft / 2) + 1)
    temp_freqs = Vector{Float64}(undef, Integer(nfft / 2) + 1)

    ccall((:get_psd_welch, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    return temp_ampls, temp_freqs
end

@brainflow_rethrow function get_band_power(psd, freq_start::Float64, freq_end::Float64)
    band_power = Vector{Float64}(undef, 1)
    ccall((:get_band_power, DATA_HANDLER_INTERFACE), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Float64, Float64, Ptr{Float64}),
            psd[1], psd[2], length(psd[1]), Float64(freq_start), Float64(freq_end), band_power)
    return band_power[1]
end

ext: jl


text:
module CumulantsUpdates
  using Cumulants
  using SymmetricTensors
  using StatsBase
  using FileIO
  using JLD2
  using Distributed
  import Cumulants: outerprodcum
  import SymmetricTensors: pyramidindices, getblockunsafe
  import LinearAlgebra: norm

  include("updates.jl")
  include("operations.jl")

  export dataupdat, momentupdat, momentarray
  export norm, moms2cums!, cums2moms, cumulantsupdate!, DataMoments, savedm, loaddm
end


text:
module Circuit

using DataStructures
using QXZoo.GateOps
using QXZoo.GateMap
using QXZoo.DefaultGates
using FileIO

include("utils/DLList.jl")
using .DS

export append!, push!
export export_circuit, load_circuit

# =========================================================================== #
#                  Circuit generation functions and structures
# =========================================================================== #

"""
    Circ

Structure for ordered quantum circuit gate operations.
Maintains mutable Doubly Linked List (DLList) of gate-calls and a set of the gates used.
"""
struct Circ
    circ_ops::DS.DLList{<:GateOps.AGateCall}
    gate_set::Set{<:GateOps.AGateSymbol}
    num_qubits::Int
end

"""
    Circ()

Default constructor for empty circuit. Registers Pauli and Hadamard gates during initialisation.
"""
function Circ()
    circ_ops = DS.DLList{Union{GateOps.GateCall1, GateOps.GateCall2}}()
    gate_set = Set{Union{GateOps.GateSymbol, GateOps.GateSymbolP}}()
    return Circ(circ_ops, gate_set, 0)
end
"""
    Circ(num_qubits::Int)

Circ constructor for given number of qubits circuit. Registers Pauli and Hadamard gates during initialisation.
"""
function Circ(num_qubits::Int)
    circ_ops = DS.DLList{Union{GateOps.GateCall1, GateOps.GateCall2}}()
    gate_set = Set{Union{GateOps.GateSymbol, GateOps.GateSymbolP}}()
    return Circ(circ_ops, gate_set, num_qubits)
end
"""
    Circ(num_qubits::Int, gate_set::Set{<:GateOps.AGateSymbol})

Circ constructor for given number of qubits circuit. Registers Pauli and Hadamard gates during initialisation, and the included gate_set.
"""
function Circ(num_qubits::Int, gate_set::Set{<:GateOps.AGateSymbol})
    circ_ops = DS.DLList{Union{GateOps.GateCall1, GateOps.GateCall2}}()
    return Circ(circ_ops, gate_set, num_qubits)
end

"""
    num_qubits(circuit::Circ)::Int

Return the number of qubits for the given circuit.
"""
function num_qubits(circuit::Circ)::Int
    return circuit.num_qubits
end

function Base.:append!(circ::Circ, gc)
    add_gatecall!(circ, gc);
end

function Base.:append!(circ1::Circ, circ2::Circ)
    if circ1.num_qubits == circ2.num_qubits
        append!(circ1.circ_ops, circ2.circ_ops);
        union!(circ1.gate_set, circ2.gate_set);
        empty!(circ2.gate_set);
    else
        error("Number of qubits incompatible between circuits: $(circ1.num_qubits), $(circ2.num_qubits)")
    end
end

"""
    add_gatecall!(circ::Circ, gc::GateOps.AGateCall)

Adds the given gate call `gc` to the circuit at the end position.

# Example 1
```julia-repl
julia> Circuit.add_gatecall!(circ, DefaultGates.x(4)) #Apply Paulix to qubit index 4
```
This functionality is also added by overloading `<<` and used as:
# Example 2
```julia-repl
julia> circ << GateOps.x(4) #Apply Paulix to qubit index 4
```
"""
function add_gatecall!(circ::Circ, gc::GateOps.AGateCall)
    if circ.num_qubits < gc.target
        error("Target is outside qubit register range")
    end
    if gc isa GateOps.GateCall2 && circ.num_qubits < gc.ctrl
        error("Ctrl is outside qubit register range")
    end
    if ~haskey(GateMap.gates, gc.gate_symbol)

        if (typeof(gc.gate_symbol) == GateOps.GateSymbolP) && (isdefined(GateMap, gc.gate_symbol.label) )
            GateMap.cache_gate!(gc.gate_symbol, ()->getproperty(GateMap, gc.gate_symbol.label)(gc.gate_symbol.param) )

        elseif haskey(GateMap.gates, GateOps.GateSymbol(gc.gate_symbol.label, 0, false) ) && (gc.gate_symbol.rt_depth != 0)
            GateMap.cache_gate!(gc.gate_symbol, ()->GateMap.gates[GateSymbol(gc.gate_symbol.label, 0, false)]()^(1/2^(gc.gate_symbol.rt_depth)) )

        elseif haskey(GateMap.gates, GateOps.GateSymbol(gc.gate_symbol.label, gc.gate_symbol.rt_depth, false) ) && (gc.gate_symbol.is_adj == true) #handle non-param adjoint gates
            GateMap.cache_gate!(gc.gate_symbol, ()->adjoint(GateMap.gates[GateOps.GateSymbol(gc.gate_symbol.label, gc.gate_symbol.rt_depth, false)]()) )

        else
            error("Gate $(gc.gate_symbol.label) not registered")
        end
    end
    if ~(gc.gate_symbol in circ.gate_set)
        push!(circ.gate_set, gc.gate_symbol)
    end
    push!(circ.circ_ops, gc);
end

function Base.push!(circ::Circ, gc::GateOps.AGateCall)
    if ~haskey(GateMap.gates, gc.gate_symbol)
        throw("Gate $(gc.gate_symbol) not registered")
    end
    if ~(gc.gate_symbol in circ.gate_set)
        push!(circ.gate_set, gc.gate_symbol)
    end
    push!(circ.circ_ops, gc);
end

Base.show(io::IO, circ::Circ) = print(io, "$(circ.num_qubits) qubits with $(circ.circ_ops.len) gate-calls using $(length(circ.gate_set)) unique gates.")
Base.show(io::IO, m::MIME"text/plain", circ::Circ) = print(io, "$(circ.num_qubits) qubits with $(circ.circ_ops.len) gate-calls using $(length(circ.gate_set)) unique gates.")

Base.:sqrt(U::Symbol) = Symbol("s"*string(U));
Base.:adjoint(U::Symbol) = Symbol("a"*string(U));
Base.:copy(circ::Circ) = Circ(circ.circ_ops, circ.gate_set, circ.num_qubits);
Base.:<<(circ::Circ, gc::GateOps.GateCall1) = add_gatecall!(circ, gc::GateOps.GateCall1)
Base.:<<(circ::Circ, gc::GateOps.GateCall2) = add_gatecall!(circ, gc::GateOps.GateCall2)
Base.:<<(circ1::Circ, circ2::Circ) = append!(circ1::Circ, circ2::Circ)

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

"""
    export_circuit(cct::Circ, file_name::String="circuit")

Exports the circuit to disk as a JLD2 file.
"""
function export_circuit(cct::Circ, file_name::String="circuit")
    save("$(file_name).jld2", Dict{String, Circ}("circuit"=>cct))
end

"""
    load_circuit(file_name::String="circuit")

Loads the circuit at `file_name.jld2` from disk.
"""
function load_circuit(file_name::String="circuit")
    load("$(file_name).jld2")
end

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

end


text:
# SPDX-License-Identifier: X11
# 2020-11-15 (Live)
# Billiards (200pt)

function main()
  Sx, Sy, Gx, Gy = parse.(Int, (split ∘ readline)())
  println((Gx - Sx) / (Gy/Sy + 1) + Sx)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module westbrick

using GLMakie
using FileIO
using LinearAlgebra
using StaticArrays

include("types.jl")
include("utils.jl")
include("visualize.jl")
include("pseudocontrol.jl")
include("action.jl")

#methods
export load_robot,
       load_object!,
       init_environment,
       init_visualization,
       visualize_trajectory!,
       visualize_environment!,
       translate!,
       turn!,
       go_to!,
       move!,
       pick!,
       place!

#types
export Object,
       Robot

end


text:
using Distributions
using JuliaWebAPI
using Logging
Logging.configure(level=DEBUG)

f(u) = exp(-u^2/2)/√(2pi)

const u1=Uniform(-5.0,+5.0)
const u2=Uniform(0.0, 0.5)

function simulate_pt()
    x=rand(u1)
    y=rand(u2)
    y<f(x)
end

function simulate(num::Int)
     hits=0
     for (i in 1:num)
        hits = hits + simulate_pt()
    end
    return hits/num*(0.5*10)
end

function simulateTime(num)
    tic()
    r = simulate(num)
    t = toc()
    return (r, t)
end

simulate(num::String) = simulate(parseint(num))

#process([(simulate,false), (simulateTime,false)], "tcp://127.0.0.1:9999"; bind=true)

max_stars_repo_name: JuliaComputing/JuliaInXL

max_issues_repo_name: JuliaComputing/JuliaInXL

max_forks_repo_name: JuliaComputing/JuliaInXL


text:
# Some useful wrappers
function get_function(model::MOI.ModelLike, ci::CI)
    return MOI.get(model, MOI.ConstraintFunction(), ci)
end

function get_set(model::MOI.ModelLike, ci::CI)
    return MOI.get(model, MOI.ConstraintSet(), ci)
end

function get_ci_row_dimension(model::MOI.ModelLike, ci::CI)
    return MOI.output_dimension(get_function(model, ci))
end

function get_scalar_term(
    model::MOI.ModelLike,
    ci::CI{VI,S},
) where {S<:MOI.AbstractScalarSet}
    return [-MOI.constant(get_set(model, ci))]
end

function get_scalar_term(
    model::MOI.ModelLike,
    ci::CI{F,S},
) where {F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}
    return [
        MOI.constant(get_function(model, ci)) -
        MOI.constant(get_set(model, ci)),
    ]
end

# This is used to fill the dual objective dictionary
function get_scalar_term(
    func::MOI.AbstractVectorFunction,
    ::MOI.AbstractVectorSet,
    i::Int,
)
    return MOI.constant(func)[i]
end

# This is used to fill the dual objective dictionary
function get_scalar_term(
    ::MOI.VariableIndex,
    set::MOI.AbstractScalarSet,
    i::Int,
)
    return -MOI.constant(set)
end

# This is used to fill the dual objective dictionary
function get_scalar_term(
    func::MOI.AbstractScalarFunction,
    set::MOI.AbstractScalarSet,
    i::Int,
)
    # In this case there i only one constant in the function and one in the set
    return MOI.constant(func) - MOI.constant(set)
end


text:
export randorth, randrealsv, randsparse

"""
    Construct random orthogonal matrix.
"""
function randorth(T::Type{<:Real}, m::Integer, n::Integer=m)
    tr = m < n
    n, m = minmax(m, n)
    A = randn(T, m, n)
    F = qrfact(A, Val(true))
    Q = copy(F.Q)[:,1:n]
    tr ? copy(Q') : Q
end
randorth(m::Integer, n::Integer=m) = randorth(Float64, m, n)

"""
    Construct random Matrix with given singular values
"""
function randrealsv(sv::AbstractVector{T}, m::Integer=0, n::Integer=0) where T<:Real
    TT = typeof(sqrt(one(eltype(sv))))
    ls = length(sv)
    m = max(m, ls)
    n = max(n, ls)
    Matrix((( randorth(TT, n, ls) * Diagonal(sv)) * randorth(TT, m, ls)')')
end

"""
    randsparse(m, n, r, p)

Construct sparse matrix with given dimensions and rank and p
"""
function randsparse(m::Integer, n::Integer, k::Integer, p::AbstractFloat)
    k <= min(m, n) || throw(ArgumentError("rank $k cannot be > min(m,n) = $(min(m, n))"))
    0 <= p <= 1 || throw(ArgumentError("probability $p in not in unit interval"))

    if m < n
        return copy(adjoint(randsparse(n, m, k, p)))
    end

    A = sprandn(m, k, p)

    if n > k
        e1(x) = x == zero(x) ? one(x) : expm1(x) / x # `(exp(x) - 1) / x` exactly
        α = log(one(p) - p)
        pp = e1(α / k) / e1(α) / k
        B = sprandn(k, n-k, pp)
        A = [A A * B]
        A[:,randperm(n)]
    else
        A
    end
end


text:
# test tree move NNI in a simple quartet
# Claudia February 2015

include("../src/types.jl")
include("../src/functions.jl")

using Base.Collections # for updateInCycle with priority queue

tree = "((1:0.1,2:0.2):0.5,3:0.3,4:0.4);"
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)

net = readTopology("prueba_tree.txt");
printEdges(net)
printNodes(net)

edge = chooseEdgeNNI(net);
println("$(edge.number)")
NNI!(net,edge)

for e in net.edge
    e.inCycle = -1
end

# Case 1
net.edge[1].inCycle = 2
net.edge[3].inCycle = 2
net.edge[4].inCycle = 2

flag,edge = chooseEdgeNNI(net,10);
NNI!(net,edge)

# Case 2
net.edge[1].inCycle = 2
net.edge[3].inCycle = 2
net.edge[4].inCycle = 2
net.edge[2].inCycle = 2

flag,edge = chooseEdgeNNI(net,10);
NNI!(net,edge)

# Case 3
net.edge[1].inCycle = 2
net.edge[2].inCycle = 3

flag,edge = chooseEdgeNNI(net,10);
NNI!(net,edge)

# Case 4
net.edge[1].inCycle = 2
net.edge[2].inCycle = 2

flag,edge = chooseEdgeNNI(net,10);
NNI!(net,edge)

# Case 5
net.edge[1].inCycle = 2
net.edge[2].inCycle = 2
net.edge[4].inCycle = 3
net.edge[5].inCycle = 3

flag,edge = chooseEdgeNNI(net,10);
NNI!(net,edge)

# Case 1
net.edge[2].inCycle = 2
net.edge[3].inCycle = 2
net.edge[4].inCycle = 2

net.node[3].inCycle = 2
net.node[6].inCycle = 2

flag,edge = chooseEdgeNNI(net,10)
NNI!(net,edge)
printEdges(net)
printNodes(net)


text:
# Implements a generic m, n, k game.
# https://en.wikipedia.org/wiki/M,n,k-game
# m: width
# n: height
# k: stones in a row required to win a game.

mutable struct MNKGame{M, N, K} <: Game
  board :: Vector{Int} # We may want https://github.com/JuliaArrays/StaticArrays.jl
  current_player :: Int
  status :: Status
  move_count :: Int
end

const TicTacToe = MNKGame{3, 3, 3}

function MNKGame{M, N, K}() where {M, N, K}
  MNKGame{M, N, K}(zeros(Int, M * N), 1, Status(), 0)
end

function Base.copy(s :: MNKGame{M, N, K}) :: MNKGame{M, N, K} where {M, N, K}
  MNKGame{M, N, K}(
    copy(s.board),
    s.current_player,
    s.status,
    s.move_count,
  )
end

current_player(game :: MNKGame) :: Int = game.current_player

# Returns a list of the Indices of all legal actions
function legal_actions(game :: MNKGame{M, N}) :: Vector{ActionIndex} where {M, N}
  if is_over(game)
    ActionIndex[]
  else
    result = Vector{ActionIndex}(undef, M*N - game.move_count)
    i = 1
    for j in 1:(M*N)
      if game.board[j] == 0
        result[i] = ActionIndex(j)
        i += 1
      end
    end
    @assert size(result)[1] + 1 == i "There should be " *
            "$(size(result)[1] + 1) legal actions but we found $(i)."

    result
  end
end

# A action is legal, if the board position is still empty
function is_action_legal( game :: MNKGame{M, N}
                        , index :: ActionIndex
                        ) :: Bool where {M, N}
  game.board[index] == 0 && !is_over(game)
end

function apply_action!( game :: MNKGame{M, N, K}
                      , index :: ActionIndex
                      ) :: MNKGame{M, N, K} where {M, N, K}

  @assert is_action_legal(game, index) "Action $index is not allowed."
  # Update the board state
  game.board[index] = game.current_player
  game.current_player = -game.current_player

  # Update the status cache
  game.move_count += 1
  game.status = tic_tac_mnk_status(game, index)
  game

end

# All possible directions for winning rows, in flat index notation.
const search_directions = [(1, 0), (0, 1), (1, 1), (1, -1)]

bound_check(m, n, fx, fy) = (fx > 0 && fx <= m && fy > 0 && fy <= n)

function tic_tac_mnk_status( game :: MNKGame{M, N, K}
                           , changed_index :: ActionIndex
                           ) where {M, N, K}

  ci = CartesianIndices((M, N))
  x, y = Tuple(ci[changed_index])
  moving = game.board[changed_index]

  for (dx, dy) in search_directions
    count = 1 # There is always the freshly placed token.
    # forward search along the direction
    # Seek using a focus point at (fx, fy).
    (fx, fy) = (x + dx, y + dy)
    while bound_check(M, N, fx, fy) && game.board[fx + (fy - 1) * M] == moving
      count += 1
      fx += dx
      fy += dy
    end
    # backward search along the direction
    (fx, fy) = (x - dx, y - dy)
    while bound_check(M, N, fx, fy) && game.board[fx + (fy - 1) * M] == moving
      count += 1
      fx -= dx
      fy -= dy
    end

    if count >= K
      return Status(moving)
    end
  end

  # Check if there are empty spaces left.
  if game.move_count == M*N
    Status(0)
  else
    Status()
  end

end

status(game :: MNKGame) :: Status = game.status

function policy_length(:: Type{<:MNKGame{M, N}}) :: Int where {M, N}
  M * N
end

# Size of the data representation of the game
function Base.size(:: Type{<:MNKGame{M, N}}) :: Tuple{Int, Int, Int} where {M, N}
  (M, N, 1)
end

# Data representation of the game as layered 2d image
function representation(game :: MNKGame{M, N}) :: Array{Float32, 3} where {M, N}
  reshape(game.current_player .* game.board, (M, N, 1))
end

# Augment single games
function augment(game :: MNKGame{M, N, K}) where {M, N, K}

  if M == N
    boards = apply_dihedral_group(reshape(game.board, (M, N)))
    map(boards) do b
      MNKGame{M, N, K}( reshape(b, (M * N,))
                      , game.current_player
                      , game.status
                      , game.move_count )
    end

  else
    boards = apply_klein_four_group(reshape(game.board, (M, N)))
    map(boards) do b
      MNKGame{M, N, K}( reshape(b, (M * N,))
                      , game.current_player
                      , game.status
                      , game.move_count )
    end

  end

end

# Augment game/label pairs
function augment( game :: MNKGame{M, N, K}
                , label :: Vector{Float32}
                ) where {M, N, K}
  if M == N

    matpol = reshape(label[2:end], (M, N))
    matpols = apply_dihedral_group(matpol)
    labels = [ vcat(label[1], reshape(mp, (M * N,))) for mp in matpols ]

  else

    matpol = reshape(label[2:end], (M, N))
    matpols = apply_klein_four_group(matpol)
    labels = [ vcat(label[1], reshape(mp, (M * N,))) for mp in matpols ]

  end

  augment(game), labels

end

function draw(game :: MNKGame{M, N}) :: Nothing where {M, N}
  board = reshape(game.board, (M,N))
  symbols = Dict(1 => "X", -1 => "O", 0 => "⋅")

  for i in 1:N
    for j in 1:M
      print(" $(symbols[board[j,i]])")
      if j != M
        print(" │")
      end
    end
    println()
    if i != N
      println("───$(repeat("┼───", M - 1))")
    end
  end
end

max_stars_repo_name: roSievers/jtac


text:
using .Zarr: ZArray, ZGroup, zgroup, zcreate,
to_zarrtype, zopen

struct ZarrDataset
  g::ZGroup
end
ZarrDataset(g::String) = ZarrDataset(zopen(g))

get_var_dims(ds::ZarrDataset,name) = reverse(ds[name].attrs["_ARRAY_DIMENSIONS"])
get_varnames(ds::ZarrDataset) = collect(keys(ds.g.arrays))
get_var_attrs(ds::ZarrDataset, name) = ds[name].attrs
get_global_attrs(ds::ZarrDataset) = ds.g.attrs
Base.getindex(ds::ZarrDataset, i) = ds.g[i]
Base.haskey(ds::ZarrDataset,k) = haskey(ds.g,k)

function add_var(p::ZarrDataset, T::Type{>:Missing}, varname, s, dimnames, attr; kwargs...)
  S = Base.nonmissingtype(T)
  add_var(p,S, varname, s, dimnames, attr; fill_value = defaultfillval(S), kwargs...)
end

function add_var(p::ZarrDataset, T::Type, varname, s, dimnames, attr;
  chunksize=s, kwargs...)
  attr["_ARRAY_DIMENSIONS"]=reverse(collect(dimnames))
  za = zcreate(T, p.g, varname, s...;attrs=attr,chunks=chunksize,kwargs...)
  za
end

#Special case for init with Arrays
function add_var(p::ZarrDataset, a::AbstractArray, varname, dimnames, attr;
  kwargs...)
  T = to_zarrtype(a)
  b = add_var(p,T,varname,size(a),dimnames,attr;kwargs...)
  b .= a
  a
end

create_empty(::Type{ZarrDataset}, path, gatts=Dict()) = ZarrDataset(zgroup(path, attrs=gatts))

backendlist[:zarr] = ZarrDataset
push!(backendregex, r"(.zarr$)|(.zarr/$)"=>ZarrDataset)

allow_parallel_write(::ZarrDataset) = true
allow_missings(::ZarrDataset) = true
to_dataset(g::ZGroup; kwargs...) = ZarrDataset(g)


text:
__precompile__()
module Delaunay

using StaticArrays
using LinearAlgebra
using PyCall

# strides(::Transpose) is not implemented. This workaround from
# https://github.com/JuliaPy/PyCall.jl/issues/555 should fix it.
using PyCall: PyObject
PyObject(x::Adjoint) = PyObject(copy(x))
PyObject(x::Transpose) = PyObject(copy(x))

const scipyspatial = PyNULL()

function __init__()
    copy!(scipyspatial, pyimport_conda("scipy.spatial", "scipy"))
end

function simplexindices_static(delaunaytriang::PyCall.PyObject)
    # The indices of simplices from the pycall object
    py_simplexinds = delaunaytriang."simplices"

    # Get the indices of the first simplex to determine
    # the number of dimensions.
    n_vertices = length(get(py_simplexinds, PyVector{Int32}, 0))

    # Convert to Julian, accounding for index differences
    # between Python and Julia
    n_simplices = length(py_simplexinds)
    juliainds = Vector{SVector{n_vertices, Int32}}(undef, n_simplices)
    for i = 1:n_simplices
        v = get(py_simplexinds, PyVector{Int32}, i - 1) .+ 1
        juliainds[i] = SVector{n_vertices}(v)
    end
    juliainds
end

function simplexindices(delaunaytriang::PyCall.PyObject)
    # The indices of simplices from the pycall object
    py_simplexinds = delaunaytriang."simplices"

    # Get the indices of the first simplex to determine
    # the number of dimensions.
    n_vertices = length(get(py_simplexinds, PyVector{Int32}, 0))

    # Convert to Julian, accounding for index differences
    # between Python and Julia
    n_simplices = length(py_simplexinds)
    juliainds = Vector{Vector{Int32}}(undef, n_simplices)
    for i = 1:n_simplices
        v = get(py_simplexinds, PyVector{Int32}, i - 1) .+ 1
        juliainds[i] = v
    end
    juliainds
end

function delaunay_static(points)
    triang = scipyspatial.Delaunay(points)
    simplexindices_static(triang)
end

function delaunay(points)
    triang = scipyspatial.Delaunay(points)
    simplexindices(triang)
end

function delaunayn(points)
    py = scipyspatial.Delaunay(points)
    indices = zeros(Int, length(py."simplices"), size(points, 2) + 1)
    pyarray_to_array!(py."simplices", indices, Int)
    return indices .+ 1 # Add 1 to account for base difference in indices
end

function pyarray_to_array!(pyobject, arr, T)
    for i = 1:length(pyobject)
        arr[i, :] = get(pyobject, PyVector{T}, i-1) # i-1 because of Python 0 indexing
    end
end

export delaunay, delaunay_static, delaunayn

end


text:
# C

function csetup()
    quote
        const c_filename = Ref{AbstractString}("code.c")
        function set_c_filename!(filename::AbstractString)
            return c_filename[] = filename
        end

        const c_chunks = Dict{AbstractString,Vector{String}}()
        function c_add_code!(c_code::AbstractString)
            return push!(get!(c_chunks, c_filename[], String[]), c_code)
        end
        function c_get_code()
            iobuffer = IOBuffer()
            allcode = Dict{AbstractString,String}()
            for (filename, chunks) in c_chunks
                for chunk in chunks
                    println(iobuffer, chunk)
                end
                allcode[filename] = String(take!(iobuffer))
            end
            return allcode
        end
        function c_write_code!()
            println("Generating C code:")
            allcode = c_get_code()
            for (filename, content) in allcode
                println("Generating $filename...")
                write(filename, content)
            end
            println("Done.")
            return nothing
        end
    end
end
export csetup

function cnewfile(filename::AbstractString, c_stmts::AbstractString)
    c_code = """
               #include <complex.h>
               #include <stdint.h>

               $c_stmts
               """
    quote
        set_c_filename!($filename)
        c_add_code!($c_code)
    end
end
export cnewfile

function ccode(c_stmts::AbstractString)
    c_code = """
               $c_stmts
               """
    quote
        c_add_code!($c_code)
    end
end
export ccode

function cvariable(name::VarName, type::VarType, c_expr::AbstractString)
    julia_code = quote
        const $(name.julia_name) = begin
            lib = CxxInterface.dlopen($(name.cxx_library))
            sym = CxxInterface.dlsym(lib, $(name.cxx_name); throw_error=false)
            if sym ≡ nothing
                libname = $(name.cxx_library)
                symname = $(name.cxx_name)
                println("Load error (library \"$libname\", symbol \"$symname\"). You can ignore this error when generating C++ code.")
                nothing
            else
                # ptr = cglobal(($(name.cxx_name), $(name.cxx_library)), $(type.julia_type))
                ptr = cglobal(sym, $(type.julia_type))
                res = unsafe_load(ptr)
                $(type.convert_to_final(:res))::$(type.final_julia_type)
            end
        end
    end
    julia_code = clean_code(julia_code)
    simple_julia_code = simplify_code(julia_code)

    c_code = """
        /*
        $(string(simple_julia_code))
        */
        const $(type.cxx_type) $(name.cxx_name) = $c_expr;
        """
    quote
        c_add_code!($c_code)
        $julia_code
    end
end
export cvariable

function cfunction(name::FnName, result::FnResult, arguments::AbstractVector{FnArg}, c_stmts::AbstractString)
    julia_code = quote
        function $(name.julia_name)($([:($(arg.julia_name)::$(arg.initial_julia_type)) for arg in arguments]...))
            res = ccall(($(name.cxx_name), $(name.cxx_library)), $(result.julia_type),
                        ($([arg.julia_type for arg in arguments if !arg.skip]...),),
                        $([arg.convert_from_initial(arg.julia_name) for arg in arguments if !arg.skip]...))
            return $(result.convert_to_final(:res))::$(result.final_julia_type)
        end
    end
    julia_code = clean_code(julia_code)
    simple_julia_code = simplify_code(julia_code)

    c_code = """
        /*
        $(string(simple_julia_code))
        */
        $(result.cxx_type) $(name.cxx_name)(
            $(join(["$(arg.cxx_type) $(arg.cxx_name)" for arg in arguments if !arg.skip], ",\n    "))
        ) {
            $c_stmts
        }
        """
    quote
        c_add_code!($c_code)
        $julia_code
    end
end
export cfunction


text:
using MicroLogging
using Compat
using Compat.Test

import MicroLogging: BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
    shouldlog, handle_message, min_enabled_level, catch_exceptions,
    configure_logging

if MicroLogging.core_in_base
    import Compat.Test: collect_test_logs, TestLogger, LogTestFailure

    function configure_logging(logger::TestLogger; min_level=Info)
        logger.min_level = min_level
        logger
    end
else
    using MicroLogging.LogTest
    import MicroLogging.LogTest: @test_logs, @test_deprecated
    import MicroLogging.LogTest: collect_test_logs, TestLogger, LogTestFailure

    # Copied from stdlib/Test/test/runtests.jl
    mutable struct NoThrowTestSet <: Test.AbstractTestSet
        results::Vector
        NoThrowTestSet(desc) = new([])
    end
    Test.record(ts::NoThrowTestSet, t::Test.Result) = (push!(ts.results, t); t)
    Test.finish(ts::NoThrowTestSet) = ts.results
end

#-------------------------------------------------------------------------------
@testset "Logging" begin

if !MicroLogging.core_in_base
    include("test_logs.jl")
    include("core.jl")
end

include("config.jl")

end


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 20, 2, 4, :y, Quadrilateral(8, 2), false),
  :properties => [1.0e5 0.0 1.0e-5 1.0e-5;],
  :x_coords => collect(range(0.0, stop=1.0, length=21)),
  :y_coords => collect(range(0.05, stop=-0.05, length=3)),
  :temp => [],
  :support => [
    (1, [0 0]),
    (2, [0 0]),
    (3, [0 0]),
    (4, [0 0]),
    (5, [0 0])
  ],
  :dtemp => [
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50
  ],
  :nspr => [
    (163, 2, 50.0)
  ]
)

@time m, dis_df, fm_df = p55(data)

@test m.loads[0:19] ≈ [0.0,-9.088858462889305e-10,7.925850199690076e-8,
  -5.429946419597694e-21,-4.382299263408736e-8,9.088858462780911e-10,
  7.925850199690076e-8,1.296940403436091e-9,9.125889568564588e-8,
  6.4039863087481865e-9,-4.512138001498236e-9,-1.0854382987425028e-20,
  -2.1036347161742548e-8,-6.403986308769904e-9,-4.512138001498197e-9,
  -1.2969404034577835e-9,9.125889568564593e-8,1.1885079699791057e-8,
  9.585609222810678e-8,-1.6280098183907838e-20]


text:
"""
    Exo1SplitBinion{T <: Real} = Parabolic{SplitBinion{T}}

An exo-perplex number is a parabolic Cayley-Dickson construct with perplex numbers.
"""
const Exo1SplitBinion{T <: Real} = Parabolic{SplitBinion{T}}

function Exo1SplitBinion(z::SplitBinion{T}) where T <: Real
    Exo1SplitBinion{T}(z, zero(SplitBinion{T}))
end

function Exo1SplitBinion(a::T, b::T, c::T, d::T) where T <: Real
    Exo1SplitBinion{T}(SplitBinion{T}(a, b), SplitBinion{T}(c, d))
end

function Exo1SplitBinion(a::T, b::T, c::T) where T <: Real
    Exo1SplitBinion{T}(SplitBinion{T}(a, b), SplitBinion{T}(c, zero(T)))
end

function Exo1SplitBinion(a::T, b::T) where T <: Real
    Exo1SplitBinion{T}(SplitBinion{T}(a, b), zero(SplitBinion{T}))
end

function Exo1SplitBinion(a::T) where T <: Real
    Exo1SplitBinion{T}(SplitBinion{T}(a, zero(T)), zero(SplitBinion{T}))
end

function Exo1SplitBinion(a::Real, b::Real)
    Exo1SplitBinion(promote(a, b)...)
end

function Exo1SplitBinion(a::Real, b::Real, c::Real)
    Exo1SplitBinion(promote(a, b, c)...)
end

function Exo1SplitBinion(a::Real, b::Real, c::Real, d::Real)
    Exo1SplitBinion(promote(a, b, c, d)...)
end

function Exo1SplitBinion(z::Exo1Real{T}) where T <: Real
    Exo1SplitBinion(z.l, zero(T), z.r, zero(T))
end

"""
    unreal(z::Exo1SplitBinion)

The unreal part of an exo-perplex number is a 3-vector.
"""
function unreal(z::Exo1SplitBinion)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Exo1SplitBinion)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", s: ")
    print(io, z.l.r)
    print(io, ", a: ")
    print(io, z.r.l)
    print(io, ", sa: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{Exo1SplitBinion{T}}) where T <: Real
    Exo1SplitBinion{T}(random(SplitBinion{T}), random(SplitBinion{T}))
end


text:
module TestOnlineLearning

using OnlineLearning, Base.Test
using Calculus

my_tests = ["test_glm.jl"]

for my_test in my_tests
  include(my_test)
end

end


text:
using Base.Iterators: product

"""
    pprint_time(s; meaningful_digits=0)

Returns a string that represents a given time (in seconds)
as a value scaled to the appropriate unit (minutes, hours, milliseconds etc)
and rounded to a given number of meaningful digits (if it is smaller than a minute).
If the latter is `0`, the result is rounded to an integer at all times.
"""
function pprint_time(s::Number; meaningful_digits::Int=0)
    if s >= 60
        pprint_large_time(s)
    else
        pprint_small_time(s, meaningful_digits)
    end
end

function pprint_large_time(s::Number)
    int_s = convert(Int, round(s))

    limits = [
        (24 * 3600, "d"),
        (3600, "h"),
        (60, "m"),
        (1, "s"),
    ]

    res = []

    for (limit, unit) in limits
        if int_s >= limit
            num_units = convert(Int, floor(int_s / limit))
            int_s -= num_units * limit
            push!(res, "$num_units$unit")
        end
    end

    join(res, " ")
end

"Round a given number to a certain number of meaningful digits."
function round_to_meaningful(s::Number, meaningful_digits::Int)
    multiplier = 10.0^(meaningful_digits - 1 - convert(Integer, floor(log10(s))))
    round(s * multiplier) / multiplier
end

function build_str(s, unit, meaningful_digits)
    if meaningful_digits == 0
        rounded_s = convert(Int, round(s))
    else
        rounded_s = round_to_meaningful(s, meaningful_digits)
    end
    "$rounded_s$unit"
end

function pprint_small_time(s::Number, meaningful_digits::Int)
    limits = [
        (1.0, "s"),
        (1e-3, "ms"),
        (1e-6, "us"),
        (1e-9, "ns")
    ]

    for (limit, unit) in limits[1:end-1]
        if s >= limit
            return build_str(s / limit, unit, meaningful_digits)
        end
    end

    limit, unit = limits[end]
    build_str(s / limit, unit, meaningful_digits)
end

# A simple lazy map function for use with rowmajor_product()

struct LazyMap
    func
    iter
end

lazymap(func, iter) = LazyMap(func, iter)

function Base.iterate(lm::LazyMap, state=nothing)
    if state === nothing
        pair = iterate(lm.iter)
    else
        pair = iterate(lm.iter, state)
    end

    if pair === nothing
        nothing
    else
        val, state = pair
        lm.func(val), state
    end
end

Base.length(lm::LazyMap) = length(lm.iter)

rowmajor_product(xss...) = lazymap(reverse, product(reverse(xss)...))

"""
    with_output_capture(func, pass_through=false)

Execute the callable `func` and capture its output (both `STDOUT` and `STDERR`) in a string.
Returns a tuple of the `func`'s return value and its output.
If `pass_through` is `true`, does not capture anything
and returns an empty string instead of the output.
"""
function with_output_capture(func, pass_through::Bool=false)

    if pass_through
        return func(), ""
    end

    stdout_old = stdout
    stderr_old = stderr

    rd, wr = redirect_stdout()
    redirect_stderr(wr)
    reader = @async read(rd, String)

    ret = nothing
    output = ""
    try
        ret = func()
    finally
        redirect_stdout(stdout_old)
        redirect_stderr(stderr_old)
        close(wr)
        output = fetch(reader)
        close(rd)
    end

    ret, output
end


text:
using Base.Test
using TestSetExtensions

@testset ExtendedTestSet "Rifraf" begin
    @includetests ARGS
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ONNXRuntime"
version = v"1.10.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/microsoft/onnxruntime.git", "0d9030e79888d1d5828730b254fedc53c7b640c1"),
    ArchiveSource("https://github.com/microsoft/onnxruntime/releases/download/v$version/onnxruntime-win-x64-$version.zip", "a0c6db3cff65bd282f6ba4a57789e619c27e55203321aa08c023019fe9da50d7"; unpack_target="onnxruntime-x86_64-w64-mingw32"),
    ArchiveSource("https://github.com/microsoft/onnxruntime/releases/download/v$version/onnxruntime-win-x86-$version.zip", "fd1680fa7248ec334efc2564086e9c5e0d6db78337b55ec32e7b666164bdb88c"; unpack_target="onnxruntime-i686-w64-mingw32")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir

if [[ $target == *-w64-mingw32* ]]; then
    chmod 755 onnxruntime-$target/onnxruntime-*/lib/*
    mkdir -p $includedir $libdir
    cp -av onnxruntime-$target/onnxruntime-*/include/* $includedir
    cp -av onnxruntime-$target/onnxruntime-*/lib/* $libdir
    install_license onnxruntime-$target/onnxruntime-*/LICENSE
else
    # Cross-compiling for aarch64-apple-darwin on x86_64 requires setting arch.: https://github.com/microsoft/onnxruntime/blob/v1.10.0/cmake/CMakeLists.txt#L186
    if [[ $target == aarch64-apple-darwin* ]]; then
        cmake_extra_args="-DCMAKE_OSX_ARCHITECTURES='arm64'"
    fi

    cd onnxruntime
    git submodule update --init --recursive
    mkdir -p build
    cd build
    cmake $WORKSPACE/srcdir/onnxruntime/cmake \
        -DCMAKE_INSTALL_PREFIX=$prefix \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_BUILD_TYPE=Release \
        -DONNX_CUSTOM_PROTOC_EXECUTABLE=$host_bindir/protoc \
        -Donnxruntime_BUILD_SHARED_LIB=ON \
        -Donnxruntime_BUILD_UNIT_TESTS=OFF \
        -Donnxruntime_DISABLE_RTTI=OFF \
        $cmake_extra_args
    make -j $nproc
    make install
    install_license $WORKSPACE/srcdir/onnxruntime/LICENSE
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
function platform_exclude_filter(p::Platform)
    libc(p) == "musl" ||
    p == Platform("i686", "Linux") || # No binary - and source build fails linking CXX shared library libonnxruntime.so
    Sys.isfreebsd(p)
end
platforms = supported_platforms(; exclude=platform_exclude_filter)
platforms = expand_cxxstring_abis(platforms; skip=!Sys.islinux)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libonnxruntime", "onnxruntime"], :libonnxruntime)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    HostBuildDependency(PackageSpec("protoc_jll", Base.UUID("c7845625-083e-5bbe-8504-b32d602b7110"), v"3.16.1"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version = v"8",
    julia_compat = "1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
using ITensors

using ITensors:
  blockdim,
  blockview,
  findfirstblock,
  tensor

# Get the range of the Index `i` for block `b`
# For example, `blockrange(Index([QN(0) => 2, QN(1) => 3]), 2) == 3:5`
function blockrange(i::Index, b::Int)
  start = 1
  for n in 1:b-1
    start += blockdim(i, n)
  end
  stop = start + blockdim(i, b) - 1
  return start:stop
end

"""
   scale_bases(M::MPO;
               normalize_op = "Id",
               combinedind_tags = ITensorLocalMPO.default_combinedind_tags)

Return a change of basis of the MPO that normalizes a specified operator basis element.

# Examples

```julia
using ITensors, ITensorLocalMPO

N = 20
s = siteinds("Electron", N; conserve_qns = true)

U = 1.0
a = AutoMPO()
for n in 1:N-1
  a .+= -1.0, "Cdagup", n, "Cup", n+1
  a .+= -1.0, "Cdagup", n+1, "Cup", n
  a .+= -1.0, "Cdagdn", n, "Cdn", n+1
  a .+= -1.0, "Cdagdn", n+1, "Cdn", n
end
for n in 1:N
  a .+= U, "Nupdn", n
end

C, X, invX = scale_bases(H; normalize_op = "Id")

# Combine the site indices
HX .*= C

# Convert to an MPS
HX = convert(MPS, copy(H))

# Apply the change of basis that scales
# the identity subspace
HX = apply(X, HX)

# Test inverting
Horig = apply(invX, HX)
Horig .*= dag.(C)
Horig = convert(MPO, Horig)

@show norm(H), norm(HX), norm(Horig)
```
"""
function scale_bases(H::MPO; kwargs...)
  s = dag.(firstsiteinds(H; plev = 0))
  return scale_bases(s; kwargs...)
end

default_combinedind_tags(n) = "Comb,Site,n=$n"

function scale_bases(s::Vector{<: Index}; normalize_op = "Id",
                     combinedind_tags = default_combinedind_tags)
  N = length(s)

  ops_to_scale = [op(normalize_op, sn) for sn in s]
  op_norms = norm.(ops_to_scale) .^ 2

  Cs = [combiner(s[n]', dag(s[n]); tags = default_combinedind_tags(n)) for n in 1:N]
  cs = combinedind.(Cs)

  ops_to_scale = ops_to_scale .* Cs

  Xs = [ITensor(cn', dag(cn)) for cn in cs]
  for X in Xs
    dX = minimum(dims(X))
    for n in 1:dX
      X[n, n] = 1
    end
  end
  invXs = deepcopy(Xs)

  # Determine the block of the QN() sector on each site
  qn0blocks = findfirstblock.(indblock -> qn(indblock) == QN(), cs)

  op_vectors = [parent(array(blockview(tensor(ops_to_scale[n]), qn0blocks[n]))) for n in 1:N]
  op_vectors ./= op_norms
  for n in 1:N
    op_vector = op_vectors[n]
    Q_qn0, _ = qr(op_vector)
    U_qn0 = Q_qn0 * I
    U_qn0[:, 1] .= op_vector
    U_qn0 = permutedims(U_qn0)
    range_qn0 = blockrange(cs[n], Int(qn0blocks[n]))
    Xs[n][range_qn0, range_qn0] = U_qn0
    invXs[n][range_qn0, range_qn0] = inv(U_qn0)
  end
  return Cs, Xs, invXs
end

"""
    combine_and_transform(H::MPO, C::Vector{<: ITensor}, X::Vector{<: ITensor}, invX::Vector{<: ITensor})

Combine and transform the sites of the MPO, convert to an MPS.
"""
function combine_and_transform(H::MPO, C::Vector{<: ITensor}, X::Vector{<: ITensor}, invX::Vector{<: ITensor})
  return apply(X, MPS(H .* C))
end

function inv_transform_and_uncombine(HX::MPS, C::Vector{<: ITensor}, X::Vector{<: ITensor}, invX::Vector{<: ITensor})
  return MPO(apply(invX, HX) .* dag.(C))
end

"""
    combine_and_transform(U::Vector{<: ITensor}, H::MPO, C::Vector{<: ITensor}, X::Vector{<: ITensor}, invX::Vector{<: ITensor})

Combine and transform the sites of the MPO as well as the gates. For gates `U`, gates are returned that are the tensor product of `U ⊗ Uᴴ`, with the indices combined and transformed.
"""
function combine_and_transform(U::Vector{<: ITensor}, H::MPO,
                               C::Vector{<: ITensor}, X::Vector{<: ITensor}, invX::Vector{<: ITensor})
  # Combine and transform the sites of the MPO, convert to an MPS
  HX = combine_and_transform(H, C, X, invX)

  UX = Vector{ITensor}(undef, length(U))
  for n in 1:length(U)
    Un = U[n]
    n1, n2 = findsites(H, Un)
    Undag = dag(swapprime(Un, 0 => 1))
    UXn = prime(Un) * prime(Undag, -1)

    C1 = C[n1]
    C2 = C[n2]
    c1 = combinedind(C1)
    c2 = combinedind(C2)

    UXn *= replaceprime(C1, 1 => 2, 0 => -1; inds = !c1)
    UXn *= replaceprime(C2, 1 => 2, 0 => -1; inds = !c2)
    UXn = prime(UXn; inds = (c1, c2))
    UXn *= dag(C1)
    UXn *= dag(C2)

    X1 = X[n1]
    X2 = X[n2]
    invX1 = invX[n1]
    invX2 = invX[n2]

    # TODO: use something like:
    #
    # apply([X1, X2], UXn, [invX1, invX2])
    #
    # for this
    UXn = replaceprime(UXn * prime(X1) * prime(X2), 0 => 1)
    UXn = UXn * invX1 * invX2
    UXn = replaceprime(UXn, 2 => 1)
    UX[n] = UXn
  end
  return UX, HX
end

max_stars_repo_name: mtfishman/ITensorLocalMPO

max_issues_repo_name: mtfishman/ITensorLocalMPO

max_forks_repo_name: mtfishman/ITensorLocalMPO


text:
#trotter expansion
function trotter(A::AbstractMatrix,B::AbstractMatrix,P::Int64)
    return (expm(full(A/(2*P)))*expm(full(B/P))*expm(full(A/(2*P))))^P
end

function (⊗)(A::AbstractMatrix,B::AbstractMatrix)
    @assert size(A)[1]==size(A)[2]
    @assert size(B)[1]==size(B)[2]

    return kron(A,B)
end

function (⊕)(A::AbstractMatrix,B::AbstractMatrix)
    return full(blkdiag(sparse(A),sparse(B)))
end

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
import ..update!, ..select_action, ..reset!

using StaticArrays
using Random
using LinearAlgebra: Diagonal

using Ramnet.Models.AltDiscriminators: FunctionalDiscriminator
using Ramnet.Optimizers
import Ramnet.Loss

mutable struct ContinousActionPolicyLoss <: Loss.AbstractLoss
    start_cov::Float64
    end_cov::Float64
    cov_decay::Int
    update_count::Int
end

function ContinousActionPolicyLoss(start_cov::Float64, end_cov::Float64, cov_decay::Int)
    ContinousActionPolicyLoss(start_cov, end_cov, cov_decay, 0)
end

# Policy approximator with only two admissible actions: {-1, 1}.
struct ContinousActionPolicy{OS,AS,T<:Real,O<:StaticArray{Tuple{OS},T,1},A<:StaticArray{Tuple{AS},T,1},E<:AbstractEncoder{T}}
    opt::FunctionalOptimizer{ContinousActionPolicyLoss}
    f::FunctionalDiscriminator{AS,T,E}
    rng::MersenneTwister
end

function ContinousActionPolicy(::Type{O}, ::Type{A}, n::Int, encoder::E; start_cov::Float64, end_cov::Float64, cov_decay::Int, start_learning_rate::Float64 = 0.1, end_learning_rate::Float64 = 1e5, learning_rate_decay::Int = 1000, epochs::Int = 1, partitioner::Symbol = :uniform_random, seed::Union{Nothing,UInt} = nothing) where {AS,OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},A<:StaticArray{Tuple{AS},T,1},E<:AbstractEncoder{T}}
    ContinousActionPolicy{OS,AS,T,O,A,E}(
        FunctionalOptimizer(ContinousActionPolicyLoss(start_cov, end_cov, cov_decay); start_learning_rate, end_learning_rate, learning_rate_decay, epochs),
        FunctionalDiscriminator{AS}(OS, n, encoder, partitioner; seed),
        isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)
    )
end

function cov(l::ContinousActionPolicyLoss)
    α = l.update_count / l.cov_decay
    UniformScaling((1 - α) * l.start_cov + α * l.end_cov)
end

function Loss.grad(l::ContinousActionPolicyLoss, action, pred_action)
    inv(cov(l)) * (pred_action - action)
end

function mode(p::ContinousActionPolicy{OS,AS,T,O,A,E}, observation::O) where {OS,AS,T,O,A,E}
    # println(typeof(predict(p.f, observation)))
    return SVector{AS,T}(predict(p.f, observation))
end

function update!(p::ContinousActionPolicy{OS,AS,T,O,A,E}, observation::O, action::A, G) where {OS,AS,T,O,A,E}
    # gradient = Optimizers.grad(p.opt, p, observation, action, G)

    gradient = Optimizers.grad(p.opt, action, predict(p.f, observation))
    train!(p.f, observation, learning_rate(p.opt) * gradient * G)

    p.opt.loss.update_count += 1

    return nothing
end

# TODO: The clamping of the action is specific to the hopper agent!
function select_action(p::ContinousActionPolicy{OS,AS,T,O,A,E}, observation::O; deterministic::Bool = False) where {OS,AS,T,O,A,E}
    if deterministic
        return mode(p, observation)
    end

    return A(cov(p.opt.loss) * randn(p.rng, AS) + predict(p.f, observation))
end

function reset!(p::ContinousActionPolicy; seed::Union{Nothing,UInt} = nothing)
    Ramnet.reset!(p.opt)
    p.opt.loss.update_count = 0
    Ramnet.reset!(p.f; seed)
    Random.seed!(p.rng, seed)

    return nothing
end


text:
using GPMLj, Test, Random

Random.seed!(123456)

include("gpflow/gpflow.jl")


text:
"""
$(TYPEDEF)

Rigid transform between two Pose2's, assuming (x,y,theta).

Related

Pose3Pose3, Point2Point2, MutablePose2Pose2Gaussian, DynPose2, InertialPose3
"""
struct Pose2Pose2{T} <: IncrementalInference.FunctorPairwise where {T <: IIF.SamplableBelief}
  z::T
  Pose2Pose2{T}() where {T <: IIF.SamplableBelief} = new{T}()
  Pose2Pose2{T}(z1::T) where {T <: IIF.SamplableBelief} = new{T}(z1)
end
Pose2Pose2(z::T) where {T <: IIF.SamplableBelief} = Pose2Pose2{T}(z)

getSample(s::Pose2Pose2{<:IIF.SamplableBelief}, N::Int=1) = (rand(s.z,N), )
function (s::Pose2Pose2{<:IIF.SamplableBelief})(
            res::Vector{Float64},
            userdata,
            idx::Int,
            meas::Tuple,
            wxi::Array{Float64,2},
            wxj::Array{Float64,2}  )
  #
  wXjhat = SE2(wxi[1:3,idx])*SE2(meas[1][1:3,idx])
  jXjhat = SE2(wxj[1:3,idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

# NOTE, serialization support -- will be reduced to macro in future
# ------------------------------------

"""
$(TYPEDEF)
"""
mutable struct PackedPose2Pose2  <: IncrementalInference.PackedInferenceType
  datastr::String
  PackedPose2Pose2() = new()
  PackedPose2Pose2(x::String) = new(x)
end
function convert(::Type{Pose2Pose2}, d::PackedPose2Pose2)
  return Pose2Pose2(extractdistribution(d.datastr))
end
function convert(::Type{PackedPose2Pose2}, d::Pose2Pose2)
  return PackedPose2Pose2(string(d.z))
end

function compare(a::Pose2Pose2,b::Pose2Pose2; tol::Float64=1e-10)
    TP = true
    TP = TP && norm(a.z.μ-b.z.μ) < (tol + 1e-5)
    TP = TP && norm(a.z.Σ.mat-b.z.Σ.mat) < tol
    return TP
end

## Deprecated
# function Pose2Pose2(mean::Array{Float64,1}, cov::Array{Float64,2}, w::Vector{Float64})
#   @warn "Pose2Pose2(mu,cov,w) is deprecated in favor of Pose2Pose2(T(...)) -- use for example Pose2Pose2(MvNormal(mu, cov))"
#   Pose2Pose2(MvNormal(mean, cov))
# end


text:
################################################################################
################################################################################
"""
    struct AngleAxis{T} <: Rotation{3,T}
    AngleAxis(Θ, x, y, z)

A 3×3 rotation matrix parameterized by a 3D rotation by angle θ about an
arbitrary axis `[x, y, z]`.

Note that the axis is not unique for θ = 0, and that this parameterization does
not continuously map the neighbourhood of the identity rotation (and therefore
might not be suitable for autodifferentation and optimization purposes).

Note: by default, the constructor will renormalize the input so that the axis
has length 1 (x² + y² + z² = 1).

Renormalization can be skipped by passing `false` as an additional constructor
argument, in which case the user provides the guarantee that the input arguments
represent a normalized rotation axis. Operations on an `AngleAxis` with a rotation
axis that does not have unit norm, created by skipping renormalization in this fashion,
are not guaranteed to do anything sensible.
"""
struct AngleAxis{T} <: Rotation{3,T}
    theta::T
    axis_x::T
    axis_y::T
    axis_z::T

    @inline function AngleAxis{T}(θ, x, y, z, normalize::Bool = true) where {T}
        if normalize
            # Not sure what to do with theta?? Should it become theta * norm ?
            norm = sqrt(x*x + y*y + z*z)
            new(θ, x/norm, y/norm, z/norm)
        else
            new(θ, x, y, z)
        end
    end
end

params(aa::AngleAxis) = SVector{4}(aa.theta, aa.axis_x, aa.axis_y, aa.axis_z)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 4 inputs (not 9).
@inline function AngleAxis(θ::Θ, x::X, y::Y, z::Z, normalize::Bool = true) where {Θ,X,Y,Z}
    AngleAxis{promote_type(promote_type(promote_type(Θ, X), Y), Z)}(θ, x, y, z, normalize)
end

# These functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{AA})(t::NTuple{9}) where {AA <: AngleAxis} = AA(QuatRotation(t)) # TODO: consider going directly from tuple (RotMatrix) to AngleAxis
@inline Base.getindex(aa::AngleAxis, i::Int) = QuatRotation(aa)[i]

@inline function Base.Tuple(aa::AngleAxis{T}) where T
    # Rodrigues' rotation formula.
    s, c = sincos(aa.theta)
    c1 = one(T) - c

    c1x2 = c1 * aa.axis_x^2
    c1y2 = c1 * aa.axis_y^2
    c1z2 = c1 * aa.axis_z^2

    c1xy = c1 * aa.axis_x * aa.axis_y
    c1xz = c1 * aa.axis_x * aa.axis_z
    c1yz = c1 * aa.axis_y * aa.axis_z

    sx = s * aa.axis_x
    sy = s * aa.axis_y
    sz = s * aa.axis_z

    # Note that the RotMatrix constructor argument order makes this look transposed:
    (one(T) - c1y2 - c1z2, c1xy + sz, c1xz - sy,
        c1xy - sz, one(T) - c1x2 - c1z2, c1yz + sx,
        c1xz + sy, c1yz - sx, one(T) - c1x2 - c1y2)
end

@inline function (::Type{Q})(aa::AngleAxis) where Q <: QuatRotation
    s, c = sincos(aa.theta / 2)
    return Q(c, s * aa.axis_x, s * aa.axis_y, s * aa.axis_z, false)
end

@inline function (::Type{AA})(q::QuatRotation) where AA <: AngleAxis
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3
    s2 = x * x + y * y + z * z
    sin_t2 = sqrt(s2)
    theta = 2 * atan(sin_t2, w)
    num_pert = eps(typeof(theta))^4
    inv_sin_t2 = 1 / (sin_t2 + num_pert)
    return principal_value(AA(theta, inv_sin_t2 * (x + num_pert), inv_sin_t2 * y, inv_sin_t2 * z, false))
end

# Trivial type conversions for RotX, RotY and RotZ
@inline function (::Type{AA})(r::RotX) where AA <: AngleAxis
    return AA(r.theta, 1, 0, 0)
end
@inline function (::Type{AA})(r::RotY) where AA <: AngleAxis
    return AA(r.theta, 0, 1, 0)
end
@inline function (::Type{AA})(r::RotZ) where AA <: AngleAxis
    return AA(r.theta, 0, 0, 1)
end

# Using Rodrigues formula on an AngleAxis parameterization (assume unit axis length) to do the rotation
# (implementation from: https://ceres-solver.googlesource.com/ceres-solver/+/1.10.0/include/ceres/rotation.h)
function Base.:*(aa::AngleAxis, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    w = rotation_axis(aa)
    st, ct = sincos(aa.theta)
    w_cross_pt = cross(w, v)
    m = dot(v, w) * (one(w_cross_pt[1]) - ct)
    T = promote_type(eltype(aa), eltype(v))
    return similar_type(v,T)(v[1] * ct + w_cross_pt[1] * st + w[1] * m,
                             v[2] * ct + w_cross_pt[2] * st + w[2] * m,
                             v[3] * ct + w_cross_pt[3] * st + w[3] * m)
end

@inline Base.:*(aa::AngleAxis, r::Rotation) = QuatRotation(aa) * r
@inline Base.:*(aa::AngleAxis, r::RotMatrix) = QuatRotation(aa) * r
@inline Base.:*(aa::AngleAxis, r::MRP) = QuatRotation(aa) * r
@inline Base.:*(r::Rotation, aa::AngleAxis) = r * QuatRotation(aa)
@inline Base.:*(r::RotMatrix, aa::AngleAxis) = r * QuatRotation(aa)
@inline Base.:*(r::MRP, aa::AngleAxis) = r * QuatRotation(aa)
@inline Base.:*(aa1::AngleAxis, aa2::AngleAxis) = QuatRotation(aa1) * QuatRotation(aa2)

@inline Base.inv(aa::AngleAxis) = AngleAxis(-aa.theta, aa.axis_x, aa.axis_y, aa.axis_z)
@inline Base.:^(aa::AngleAxis, t::Real) = AngleAxis(aa.theta*t, aa.axis_x, aa.axis_y, aa.axis_z)
@inline Base.:^(aa::AngleAxis, t::Integer) = AngleAxis(aa.theta*t, aa.axis_x, aa.axis_y, aa.axis_z) # to avoid ambiguity

# define identity rotations for convenience
@inline Base.one(::Type{AngleAxis}) = AngleAxis(0.0, 1.0, 0.0, 0.0)
@inline Base.one(::Type{AngleAxis{T}}) where {T} = AngleAxis{T}(zero(T), one(T), zero(T), zero(T))

# accessors
@inline rotation_angle(aa::AngleAxis) = aa.theta #  - floor((aa.theta+pi) / (2*pi)) * 2*pi
@inline rotation_axis(aa::AngleAxis) = SVector(aa.axis_x, aa.axis_y, aa.axis_z)

################################################################################
################################################################################
"""
    struct RotationVec{T} <: Rotation{3,T}
    RotationVec(sx, sy, sz)

Rodrigues vector parameterization of a 3×3 rotation matrix. The direction of the
vector [sx, sy, sz] defines the axis of rotation, and the rotation angle is
given by its norm.
"""
struct RotationVec{T} <: Rotation{3,T}
    sx::T
    sy::T
    sz::T
end

params(aa::RotationVec) = SVector{3}(aa.sx, aa.sy, aa.sz)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 4 inputs (not 9).
@inline RotationVec(x::X, y::Y, z::Z) where {X,Y,Z} = RotationVec{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{RV})(t::NTuple{9}) where {RV <: RotationVec} = RV(QuatRotation(t)) # TODO: go through AngleAxis once it's faster
@inline Base.getindex(aa::RotationVec, i::Int) = QuatRotation(aa)[i]
@inline Base.Tuple(rv::RotationVec) = Tuple(QuatRotation(rv))

function (::Type{AA})(rv::RotationVec) where AA <: AngleAxis
    # TODO: consider how to deal with derivative near theta = 0. There should be a first-order expansion here.
    theta = rotation_angle(rv)
    return theta > 0 ? AA(theta, rv.sx / theta, rv.sy / theta, rv.sz / theta, false) : AA(zero(theta), one(theta), zero(theta), zero(theta), false)
end

function (::Type{RV})(aa::AngleAxis) where RV <: RotationVec
    return RV(aa.theta * aa.axis_x, aa.theta * aa.axis_y, aa.theta * aa.axis_z)
end

function (::Type{Q})(rv::RotationVec) where Q <: QuatRotation
    return QuatRotation(AngleAxis(rv))
end

(::Type{RV})(q::QuatRotation) where {RV <: RotationVec} = RV(AngleAxis(q))

function Base.:*(rv::RotationVec{T1}, v::StaticVector{3, T2}) where {T1,T2}
    theta = rotation_angle(rv)
    if (theta > eps(T1)) # use eps here because we have the 1st order series expansion defined
        return AngleAxis(rv) * v
    else
        return similar_type(typeof(v), promote_type(T1,T2))(
                    v[1] + rv.sy * v[3] - rv.sz * v[2],
                    v[2] + rv.sz * v[1] - rv.sx * v[3],
                    v[3] + rv.sx * v[2] - rv.sy * v[1])
    end
end

@inline Base.:*(rv::RotationVec, r::Rotation) = QuatRotation(rv) * r
@inline Base.:*(rv::RotationVec, r::RotMatrix) = QuatRotation(rv) * r
@inline Base.:*(rv::RotationVec, r::MRP) = QuatRotation(rv) * r
@inline Base.:*(rv::RotationVec, r::AngleAxis) = QuatRotation(rv) * r
@inline Base.:*(r::Rotation, rv::RotationVec) = r * QuatRotation(rv)
@inline Base.:*(r::RotMatrix, rv::RotationVec) = r * QuatRotation(rv)
@inline Base.:*(r::MRP, rv::RotationVec) = r * QuatRotation(rv)
@inline Base.:*(r::AngleAxis, rv::RotationVec) = r * QuatRotation(rv)
@inline Base.:*(rv1::RotationVec, rv2::RotationVec) = QuatRotation(rv1) * QuatRotation(rv2)

@inline Base.inv(rv::RotationVec) = RotationVec(-rv.sx, -rv.sy, -rv.sz)
@inline Base.:^(rv::RotationVec, t::Real) = RotationVec(rv.sx*t, rv.sy*t, rv.sz*t)
@inline Base.:^(rv::RotationVec, t::Integer) = RotationVec(rv.sx*t, rv.sy*t, rv.sz*t) # to avoid ambiguity

# rotation properties
@inline rotation_angle(rv::RotationVec) = sqrt(rv.sx * rv.sx + rv.sy * rv.sy + rv.sz * rv.sz)
function rotation_axis(rv::RotationVec)     # what should this return for theta = 0?
    theta = rotation_angle(rv)
    return (theta > 0 ? SVector(rv.sx / theta, rv.sy / theta, rv.sz / theta) : SVector(one(theta), zero(theta), zero(theta)))
end

# define identity rotations for convenience
@inline Base.one(::Type{RotationVec}) = RotationVec(0.0, 0.0, 0.0)
@inline Base.one(::Type{RotationVec{T}}) where {T} = RotationVec{T}(zero(T), zero(T), zero(T))


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Plots, LaTeXStrings
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.12"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	mu, sig = 2, 3
	eta = sqrt(3)*sig/pi
	n, N = 15, 10^4
	dNormal   = Normal(mu, sig)
	dLogistic = Logistic(mu, eta)
	alphaUsed = 0.001:0.001:0.1
end;

# ╔═╡ 874034ba-3b2b-11eb-2b77-c1d2d28a1116
function alphaSimulator(dist, n, alpha)
    popVar        = var(dist)
    coverageCount = 0
    for _ in 1:N
        sVar = var(rand(dist, n))
        L = (n - 1) * sVar / quantile(Chisq(n-1),1-alpha/2)
        U = (n - 1) * sVar / quantile(Chisq(n-1),alpha/2)
        coverageCount +=  L < popVar && popVar < U
    end
    1 - coverageCount/N
end

# ╔═╡ 87407402-3b2b-11eb-36de-2b731fe48190
begin
	scatter(alphaUsed, alphaSimulator.(dNormal,n,alphaUsed),
		c=:blue, msw=0, label="Normal")
	scatter!(alphaUsed, alphaSimulator.(dLogistic, n, alphaUsed),
		c=:red, msw=0, label="Logistic")
	plot!([0,0.1],[0,0.1],c=:black, label="1:1 slope",
		xlabel=L"\alpha"*" used", ylabel=L"\alpha"*" actual",
		legend=:topleft, xlim=(0,0.1), ylims=(0,0.2))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.12"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═874034ba-3b2b-11eb-2b77-c1d2d28a1116
# ╠═87407402-3b2b-11eb-36de-2b731fe48190
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
"""
    MinimalDegreeOfSingerAlgebraCache

a dictionary that stores at the key `e` (an integer), if bound,
a dictionary that stores at the key `q` (an integer), if bound,
the value `m(q,e)`;
the values get stored by calls of `MinimalDegree`.
"""
MinimalDegreeOfSingerAlgebraCache = Dict{Integer,Dict{Integer,Integer}}();

"""
    VectorIterator( m, n, s )

Iterate over all vectors of length `n`
with entries in the set { 0, 1, 2, ... `m` }
and coefficient sum `s`.
The vectors are enumerated in lexicographical order.
(The functions change the vectors in place,
just collecting the results makes no sense.)
"""
mutable struct VectorIterator
    m::Int
    n::Int
    s::Int
end

# Auxiliary function for the iterator:
# Distribute `s` to the first `n` positions in the array `v`
# such that as many initial entries as possible are `m`.
function VectorIterator_ResetPrefix(v::Vector{Int}, m::Int, n::Int, s::Int)
    local rest::Int, i::Int, j::Int

    rest = s
    i = 1
    while m < rest
      v[i] = m
      rest = rest - m
      i = i + 1
    end
    v[i] = rest
    for j in (i+1):n
      v[j] = 0
    end

    return v
end

# Initialize the iterator.
function Base.iterate(vi::VectorIterator)
    local next::Vector{Int}

    if vi.s > vi.m * vi.n
      # The iterator is empty.
      return
    end
    next = VectorIterator_ResetPrefix(zeros(Int, vi.n), vi.m, vi.n, vi.s)
    return (next, next)
end

# Define the iteration step.
# Note that `state` is changed in place.
function Base.iterate(vi::VectorIterator, state::Vector{Int})
    local sum::Int, i::Int

    # Find the first position with a nonzero value
    # such that the value on the right is smaller than m.
    sum = -1
    for i in 1:(vi.n-1)
      sum = sum + state[i]
      if state[i] != 0 && state[ i+1 ] < vi.m
        state[i] = state[i] - 1
        state[i+1] = state[i+1] + 1
        VectorIterator_ResetPrefix(state, vi.m, i, sum)
        return (state, state)
      end
    end

    # There is no such position, we are done.
    return
end

"""
    MinimalDegreeCheap(q::Int, n::Int, e::T) where {T<:Integer}

Return either 0 or the number ``m(q,e)``,
which is the smallest number of powers of `q`
such that `e` divides the sum of these powers.

The return value 0 means that the cheap criteria from the two papers
do not suffice to determine ``m(q,e)``.

The argument `n` must be `OrderMod(q, e)`.
"""
function MinimalDegreeCheap(q::Int, n::Int, e::T) where {T<:Integer}
    if e == 1
      return 1
    end

    @assert powermod( q, n, e ) == 1 && OrderMod( q, e, n ) == n "$n is not the mult. order of $q mod $e"

    if q == 1
      # Example I.6.1.
      return Int(e)
    elseif iseven(n) && powermod(q, div(n, 2), e) == e-1
      # Decide the case m = 2, see Lemma I.6.3:
      return 2
    elseif n == 2
      # Prop. II.2.7.
      e1 = Int(gcd(e, q-1))
      e2 = Int(gcd(e, q+1))
      if e2 <= e1 || (iseven(e) && iseven(div(q^2-1, e)))
        return e1
      else
        return 2 * e1
      end
    end

    # Deal with those cases where e is a prime power,
    # and where we know the value of m(q,e).
    if IsPrimePowerInt(e)
      # Let e = p^k for a prime p.
      if mod(e, 2) == 0
        # see Prop. II.2.9, note that here we have e > 4 and q mod e <> q-1.
        if mod(q, 4) == 1
          return Int(gcd(e, q-1))
        else
          @assert (mod(q+1, e) != 0) "the case m = 2 cannot occur here"
          return 4
        end
      else
        # see Prop. II.2.10, note that q \equiv 1 \pmod{p}
        # if and only if \gcd( e, q-1 ) \not= 1 holds.
        m = Int(gcd(e, q-1))
        if m != 1
          return m
        elseif mod(n, 3) == 0
          # Here we know that m(q,e) > 2 and thus p \not= 3.
          # Hence \gcd( e, q^{n/3}-1 ) = 1 holds,
          # because otherwise q^{n/3} would be a p-element in (Z/eZ)^*.
          # Thus Lemma I.6.4 yields m(q,e) \leq 3.
          # (This is a generalization of Cor. II.2.11.)
          return 3
        elseif 2 * n == Phi(e)
          # Cor. II.2.14.
          return 3
        elseif mod(e, 11) == 0
          # Prop. II.2.16.
          return 5
        end
      end
    end

    # Deal with those cases where e is twice an odd prime power,
    # and where we know the value of m(q,e).
    if iseven(e) && IsPrimePowerInt(div(e, 2))
      # see Prop. II.2.17, note that the order of q mod e is a power of p
      # if and only if \gcd( e/2, q-1 ) \not= 1 holds.
      m = gcd(e, q-1)
      if m > 2
        return Int(m)
      end
    end

    # We give up.
    return 0
end;

"""
    MinimalDegreeHard(q::Int, n::Int, e::T) where {T<:Integer}

Return the number `m(q,e)`, which is the smallest number of powers of `q`
such that `e` divides the sum of these powers.
The argument `n` must be `OrderMod( q, e )`.
"""
function MinimalDegreeHard(q::Int, n::Int, e::T) where {T<:Integer}
    @assert powermod( q, n, e ) == 1 && OrderMod( q, e, n ) == n "$n is not the mult. order of $q mod $e"

    # If A[q,n,z] has small dimension and e is not very small
    # then enumerate its basis.
    # (For small e, we known that the m values are small,
    # so it is cheaper not to enumerate.)
    if e > 100
      z = div(BigInt(q)^n - 1, e)
      if z <= 10^5
        return LoewyStructureInfo(q, n, Int(z))[:m]
      end
    end

    # m(q,e) is divisible by \gcd( e, q-1 ), and m(q,e) \geq 3.
    e1 = Int(gcd(e, q-1))
    if e1 == 1
      m = 3
    elseif e1 == 2
      m = 4
    else
      m = e1
    end

    powers = [powermod(q, i, e) for i in 1:(n-1)]

    while true
      for a in 1:m
        for v in VectorIterator(a, n-1, m-a)
          if mod(a + sum(v .* powers), e) == 0
            return m
          end
        end
      end
      m = m + e1
    end
end;

"""
    MinimalDegree(q::Int, e::T) where {T<:Integer}
    MinimalDegree(q::Int, n::Int, e::T) where {T<:Integer}
    MinimalDegree(q::Int, n::Int, e::T, A::GAP.GapObj) where {T<:Integer}

Return the number `m(q,e)`, which is the smallest number of powers of `q`
such that `e` divides the sum of these powers.

If `n` is given, it must be the multiplicative order of `q` modulo `e`.
"""
function MinimalDegree(q::Int, e::T) where {T<:Integer}
    return MinimalDegree(q, OrderMod(q, e), e)
end;

function MinimalDegree(q::Int, n::Int, e::T) where {T<:Integer}
    q = Int(mod(q, e))

    if ! haskey(MinimalDegreeOfSingerAlgebraCache, e)
      cache_e = Dict{Integer,Integer}()
      MinimalDegreeOfSingerAlgebraCache[e] = cache_e
    else
      cache_e = MinimalDegreeOfSingerAlgebraCache[e]
      if haskey(cache_e, q)
        return cache_e[q];
      end
    end

    m = MinimalDegreeCheap(q, n, e)
    if m == 0
      m = MinimalDegreeHard(q, n, e)
    end

    # Extend the cache.
    cache_e[q] = m

    return m
end;

function MinimalDegree(q::Int, n::Int, e::T, A::GAP.GapObj) where {T<:Integer}
    # Compute the value.
    m = MinimalDegree(q, n, e)

    # Set the attribute value in A.
    GAP.Globals.SetMinimalDegreeOfSingerAlgebraJulia(A, m)

    return m
end;

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
module JSONTypeProvider

import JSON3

include("build_type.jl")
include("write_type.jl")

end


text:
module QSimulator

export ⊗

const ⊗ = kron

include("fourier.jl")
include("perturbative_transmon.jl")
include("basis_utils.jl")
include("systems.jl")
include("operators.jl")
include("composite_systems.jl")
include("time_evolution.jl")
end


text:
using Test

ENV["JULIA_PKG_SERVER"] = ""
using JuliaZH
@test !isempty(ENV["JULIA_PKG_SERVER"])

@testset "check ENV" begin
    @test ENV["REPL_LOCALE"] == "zh_CN"
    JuliaZH.set_mirror("BFSU")
    @test ENV["JULIA_PKG_SERVER"] == "https://mirrors.bfsu.edu.cn/julia"
end


text:
# ---
# title: 860. Lemonade Change
# id: problem860
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Greedy
# link: <https://leetcode.com/problems/lemonade-change/description/>
# hidden: true
# ---
#
# At a lemonade stand, each lemonade costs `/$/5`.
#
# Customers are standing in a queue to buy from you, and order one at a time (in
# the order specified by `bills`).
#
# Each customer will only buy one lemonade and pay with either a `/$/5`, `/$/10`, or
# `/$/20` bill.  You must provide the correct change to each customer, so that the
# net transaction is that the customer pays /$/5.
#
# Note that you don't have any change in hand at first.
#
# Return `true` if and only if you can provide every customer with correct
# change.
#
#
#
# **Example 1:**
#
#
#
#     Input: [5,5,5,10,20]
#     Output: true
#     Explanation:
#     From the first 3 customers, we collect three /$/5 bills in order.
#     From the fourth customer, we collect a /$/10 bill and give back a /$/5.
#     From the fifth customer, we give a /$/10 bill and a /$/5 bill.
#     Since all customers got correct change, we output true.
#
#
# **Example 2:**
#
#
#
#     Input: [5,5,10]
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: [10,10]
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: [5,5,10,10,20]
#     Output: false
#     Explanation:
#     From the first two customers in order, we collect two /$/5 bills.
#     For the next two customers in order, we collect a /$/10 bill and give back a /$/5 bill.
#     For the last customer, we can't give change of /$/15 back because we only have two /$/10 bills.
#     Since not every customer received correct change, the answer is false.
#
#
#
#
# **Note:**
#
#   * `0 <= bills.length <= 10000`
#   * `bills[i]` will be either `5`, `10`, or `20`.
#
#
## @lc code=start
using LeetCode

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