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

using Markdown
using InteractiveUtils

# ╔═╡ 4daac2cf-dd9c-4bff-861a-fdcb828df03c
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ b673d52b-0d39-4dba-b2f6-6d3c2c207e3a
begin
	using BenchmarkTools
	using QuantumStateBase
	using QuantumStatePlots
	using Plots
	wf = WignerFunction(-10:0.1:10, -10:0.1:10)
end;

# ╔═╡ 79ca1104-0218-4426-8722-0024b92a0eef
md"
# Quantum State Base

JingYu Ning
"

# ╔═╡ 4d8b01d4-dcbe-4d1a-9c22-da4e581d8143
md"
## Fock state

Fock state is thhe common basis for quantum mechanics. In this prokect, we construct most of the states in Fock basis. For example, `CoherentState`, `SqueezedState` and `SqueezedThermalState` etc.
"

# ╔═╡ 16771043-30ec-4d89-8ff7-46161f33b8b3
FockState(0)

# ╔═╡ c31bba4f-e23e-4f26-aae7-c02c79f4f6b5
plot_wigner(wf(FockState(0)), Contour)

# ╔═╡ 4974bb09-0901-4721-ac9e-37269c0242b1
md"
Also, the `QuantumStateBase` accept various alias like `VacuumState`, `SinglePhotonState`, and `NumberState` etc.
"

# ╔═╡ b0c1c790-8dcc-4d76-803c-bf6a3bd61d95
vec(NumberState(5)) == vec(FockState(5))

# ╔═╡ 4a26a842-9907-49df-a954-ced09352f39d
vec(VacuumState()) == vec(FockState(0))

# ╔═╡ ac4e9d5c-ae31-4b60-a1d7-21017db9da2d
vec(SinglePhotonState()) == vec(FockState(1))

# ╔═╡ 4dce2059-882d-40b1-b2ec-31c9e44307aa
md"
## Coherent state

Coherent state is defined as:

$| \alpha \rangle = exp(-\frac{1}{2}|\alpha|^2)
\sum^\infty_{n=0} \frac{\alpha^n}{\sqrt{n!}} | n \rangle$

We can therefore define displacement operator and implies that

$| \alpha \rangle = \hat{D}(\alpha) | 0 \rangle$

One can use the pre-defined constructor to declare a coherent state.
"

# ╔═╡ ab919e88-14b8-4ff8-896c-dcfe706a0ef7
CoherentState(α(5., π/4))

# ╔═╡ d6ad648e-be5e-408e-a93b-3b24e82001d2
md"
It is also recommend to construct a coherent state by apply the displacement operator onto a vacuum state
"

# ╔═╡ e61381a9-1291-4bd9-adb9-584128607060
displace!(VacuumState(), α(5., π/4))

# ╔═╡ f20ea730-270d-4a61-a8ba-11fbc6e38358
plot_wigner(wf(CoherentState(α(5., π/4))), Contour)

# ╔═╡ 3fbb03ea-5a92-4589-9748-63e6a0486240
md"
## Squeezed state

Squeezed state is defined as:

$exp(\frac{1}{2} (\xi^* \hat{a}^2 - \xi \hat{a}^{\dagger 2})) | 0 \rangle$

We implies Squeezing operator

$\hat{S}(\xi) = exp(\frac{1}{2} (\xi^* \hat{a}^2 - \xi \hat{a}^{\dagger 2}))$

One can use the pre-defined constructor to declare a squeezed state.
"

# ╔═╡ d5cf7b3f-c9a2-4450-8bee-cf510100968f
SqueezedState(ξ(0.8, π/4))

# ╔═╡ 360282a0-47e1-4501-ae88-ce08c5958c7a
md"
It is also recommend to construct a squeezed state by apply the squeezing operator onto a vacuum state
"

# ╔═╡ ed4b3452-e969-4e37-a9d0-aaec89e59120
squeeze!(VacuumState(), ξ(0.8, π/4))

# ╔═╡ 1c765e8a-c47e-4133-87ac-810a55c257c9
plot_wigner(wf(SqueezedState(ξ(0.8, π/4))), Contour)

# ╔═╡ 5e5052e9-9745-40f2-99fa-7eed9e3d5197
md"
## Other State

There are numerous of pre-defined constructor to construct state. For example, `ThermalState`, `SqueezedThermalState` etc.

The most powerful featuure is that one can construct their own state by operators such as $a$ and $a^{\dagger}$, $\hat{S}$ and $\hat{D}$ and so on.

The following sectioon is about constructing a cat state.
"

# ╔═╡ 85eee830-aa57-45b1-b412-83faa69c0508
md"
## Cat State in single mode

Cat state is the superposition of two coherent state with opposite phase.

The following formulas are coherent states defined in the Fock state.

$\begin{align*}
| \alpha \rangle &= exp(-\frac{1}{2}|\alpha|^2)
\sum^\infty_{n=0} \frac{\alpha^n}{\sqrt{n!}} | n \rangle \\

&= \hat{D}(\alpha) | 0 \rangle \\

| -\alpha \rangle &= exp(-\frac{1}{2}|-\alpha|^2)
\sum^\infty_{n=0} \frac{(-\alpha)^n}{\sqrt{n!}} | n \rangle \\

&= \hat{D}(-\alpha) | 0 \rangle
\end{align*}$

The even cat state is therefore defined as:

$| cat_e \rangle = | \alpha \rangle + | -\alpha \rangle$

And odd cat state is therefore defined as:

$| cat_o \rangle = | \alpha \rangle - | -\alpha \rangle$
"

# ╔═╡ 92f0daf2-ee46-11eb-27f0-5f90c0bcb61d
begin
	αₚ = displace!(VacuumState(), α(2., π/4))
	αₙ = displace!(VacuumState(), α(-2., π/4))

	catₑ = StateVector(vec(αₚ) + vec(αₙ), QuantumStateBase.DIM)
	catₒ = StateVector(vec(αₚ) - vec(αₙ), QuantumStateBase.DIM)
end;

# ╔═╡ 7c5728ce-5d9e-4f9f-b036-baa09b34b310
plot_wigner(wf(catₑ), Contour)

# ╔═╡ 8d56e643-aef3-426e-b098-bdb45a6f20b3
plot_wigner(wf(catₒ), Contour)

# ╔═╡ 06cc70ba-3de8-40b2-9ace-00ab2c8656dd
begin
	cat_points = rand(catₑ, 4096)

	scatter(
		cat_points[1, :], cat_points[2, :],
		title="Even cat state", xlabel="θ", ylabel="x", legend=false,
		ylim=(-5, 5), size=(800, 400),
	)
end

# ╔═╡ b39f1784-95ec-4529-b922-fe6bed27eae2
md"
## The blazing fast non-Gaussian data sampler

It is always a difficult topic when it comes to sampling data from arbitrary prabability density function.

In this project, we introduse a non-Gaussian prabability density function sampler. The sampler is implemented by rejection method. And can sample about 4000 points from a pdf in few seconds.
"

# ╔═╡ 1b8d62a2-3bdc-4a22-9f21-bbb5198f8b7b
@benchmark rand(catₑ, 4096)

# ╔═╡ Cell order:
# ╟─79ca1104-0218-4426-8722-0024b92a0eef
# ╟─4daac2cf-dd9c-4bff-861a-fdcb828df03c
# ╠═b673d52b-0d39-4dba-b2f6-6d3c2c207e3a
# ╟─4d8b01d4-dcbe-4d1a-9c22-da4e581d8143
# ╠═16771043-30ec-4d89-8ff7-46161f33b8b3
# ╠═c31bba4f-e23e-4f26-aae7-c02c79f4f6b5
# ╟─4974bb09-0901-4721-ac9e-37269c0242b1
# ╠═b0c1c790-8dcc-4d76-803c-bf6a3bd61d95
# ╠═4a26a842-9907-49df-a954-ced09352f39d
# ╠═ac4e9d5c-ae31-4b60-a1d7-21017db9da2d
# ╟─4dce2059-882d-40b1-b2ec-31c9e44307aa
# ╠═ab919e88-14b8-4ff8-896c-dcfe706a0ef7
# ╟─d6ad648e-be5e-408e-a93b-3b24e82001d2
# ╠═e61381a9-1291-4bd9-adb9-584128607060
# ╠═f20ea730-270d-4a61-a8ba-11fbc6e38358
# ╟─3fbb03ea-5a92-4589-9748-63e6a0486240
# ╠═d5cf7b3f-c9a2-4450-8bee-cf510100968f
# ╟─360282a0-47e1-4501-ae88-ce08c5958c7a
# ╠═ed4b3452-e969-4e37-a9d0-aaec89e59120
# ╠═1c765e8a-c47e-4133-87ac-810a55c257c9
# ╟─5e5052e9-9745-40f2-99fa-7eed9e3d5197
# ╟─85eee830-aa57-45b1-b412-83faa69c0508
# ╠═92f0daf2-ee46-11eb-27f0-5f90c0bcb61d
# ╠═7c5728ce-5d9e-4f9f-b036-baa09b34b310
# ╠═8d56e643-aef3-426e-b098-bdb45a6f20b3
# ╠═06cc70ba-3de8-40b2-9ace-00ab2c8656dd
# ╟─b39f1784-95ec-4529-b922-fe6bed27eae2
# ╠═1b8d62a2-3bdc-4a22-9f21-bbb5198f8b7b


text:
module Lecture2

using Test

include("../../lib/KTBC.jl")
import .KTBC.CreateKTBC

include("../../lib/SpecialVectors.jl")
using .SpecialVectors: constants, linear, squares, rampAtK, sines, cosines, exponentials

include("../../lib/DifferencesMatrix.jl")
using .DifferencesMatrix: backwardDiff, forwardDiff

dims = 4

K, T, B, C = CreateKTBC(dims)

Kₐ = convert(Array{Int}, K)

cV = constants(dims)

@test (K * cV)[2:3] == [
  0.0
  0.0
]

lV = linear(dims)

@test (K * lV)[2:3] == [
  0.0
  0.0
]

sV = squares(dims)

@test (-K * sV)[2:3] == [
 2.0
 2.0
]

rV = rampAtK(2, dims)

@test (-K * rV)[1:3] == [
  0.0
  1.0
  0.0
]

t = π / 4
multiplier = round(2 * cos(t) - 2, digits=2)

sinV = round.(sines(t, dims), digits=2)
@test round.((-K * sinV)[1:3], digits=1) == round.((multiplier * sinV)[1:3], digits=1)

cosV = round.(cosines(t, dims), digits=2)
@test round.((-K * cosV)[2:3], digits=1) == round.((multiplier * cosV)[2:3], digits=1)

expV = round.(exponentials(t, dims), digits=2)
@test round.((-K * expV)[2:3], digits=1) == round.((multiplier * expV)[2:3], digits=1)

# Matlab experiment
function mathlabExperiment(n::Int)
    h = 1 / (n + 1)
    u = cos.(π * [1:n...] * h / 2)
    c = (π / 2)^2
    f = c * u

    K, T, B, C = CreateKTBC(n)

    U = (h^2) * T \ f
    e = 1 - U[1]
    g = [c / 2; f]

    T = backwardDiff(n + 1) * forwardDiff(n + 1)

    V = (h^2) * T \ g
    E = 1 - V[1]

    return f, e, T, V, E
end

f, e, T, V, E = mathlabExperiment(3)
f, e, T, V, E = mathlabExperiment(7)
f, e, T, V, E = mathlabExperiment(15)

end

ext: jl


text:
"""
    fluence_DA_inf_TD(t, β::Array{Float64,1}, ρ::Float64, nmed::Float64)

Compute the time-domain fluence in an infinite medium with Eqn. 3 of Patterson. et al. 1989.

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `nmed::Float64`: the sample medium's index of refraction

# Examples
julia> fluence_DA_inf_TD(0:1:5, [0.1,10.0], 1.0, 1.0)
"""
function fluence_DA_inf_TD(t, β::Array{Float64,1}, ρ::Float64, nmed::Float64 = 1.0)
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
    ν::Float64 = 29.9792458/nmed

    ϕ = Array{Float64}(undef, length(t))

    Threads.@threads for n in eachindex(t)

        ϕ[n] = -(ρ^2/(4D*ν*t[n]))
        ϕ[n] = ϕ[n] - μa*ν*t[n]
        ϕ[n] = ν*exp(ϕ[n])/((4π*D*ν*t[n])^(3/2))

        if isnan(ϕ[n])
            ϕ[n] = 0
        end

    end

    return ϕ
end

"""
    fluence_DA_inf_CW(ρ::Float64, β::Array{Float64,1})

Compute the fluence for a steady-state source in an infinite medium.

# Arguments
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)

# Examples
julia> fluence_DA_inf_CW(1.0, [0.1,10.0])
"""
function fluence_DA_inf_CW(ρ::Float64, β::Array{Float64,1})
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp

    ϕ = exp(-sqrt(3*μsp*μa)*ρ)/(4*π*ρ*D)

    if isnan(ϕ)
        ϕ = 0
    end

    return ϕ
end

"""
    fluence_DA_inf_FD(ρ::Float64, β::Array{Float64,1})

Compute the fluence for a steady-state source in an infinite medium.

# Arguments
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)

# Examples
julia> fluence_DA_inf_FD(1.0, [0.1,10.0])
"""
function fluence_DA_inf_FD(ρ, β, ω, nmed)
    μa = β[1]
    μsp = β[2]
    D = 1/3μsp
    ν = 29.9792458/nmed

    ϕ = exp(-ρ*sqrt(μa/D + im*ω/(ν*D)))/(4*π*ρ*D)

    if isnan(ϕ)
        ϕ = 0
    end

    return ϕ
end

N = 512
T = 10e-12
ω = 2*π*Vector(0:N-1)/(N*T)

ω = 0.0977*Vector(1:512)

map((ω) -> fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0), ω)
quadgk((ω) -> exp(im*ω*1)*fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0), -Inf, Inf)

function TD_test(t, N, ub)

    Rt = Array{Float64,1}(undef, length(t))

    x, w = gauss(N, 0, ub)

    map!(t -> sum(FD.(x, t) .* w), Rt, t)

return Rt./pi
end

FD(ω, t) = real.(exp(im*ω*t)*fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0))
x, w = gauss(1000, 0, 500)
sum(FD.(x, 1.0) .* w)

a = map(ω -> fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0)*exp(im*ω*t), [0,1,2, 3, 4, -3, -2, -1])

r = map(t -> (quadgk((ω) -> exp(im*ω*t)*fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0), -Inf, Inf))[1], 0.01:0.05:3)

map(t -> (quadgk((ω) -> 1/2pi*exp(im*ω*t)*fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0), 0, Inf))[1], 0.01:0.05:3)

map(t -> (quadgk((ω) -> exp(im*ω*t)*fluence_FD(1.0, 0., [0.1,0.1], [10.,10.], 1.0, ω, 1.0), 0, Inf))[1], 0.01:0.05:3)
r1 = map(t -> (quadgk((ω) -> exp(im*ω*t)*fluence_FD(1.0, 0., [0.1,0.1], [10.,10.], 1.0, ω, 1.0), -1000, 1000))[1], 1:1:3)

map(t -> (quadgk((ω) -> exp(im*ω*t)*fluence_DA_inf_FD(1.0, [0.1,10.0], ω, 1.0), 0, 1024))[1], 0.01:0.05:3)

function FD(ω)
    D = 1/30
    ν = 29.9792458 #speed of light in cm/ns

    ϕ = exp(-sqrt(0.1/D + im*ω/(ν*D)))
end

using QuadGK

tind = 0.1:0.1:3
r = map(t -> (quadgk((ω) -> exp(im*ω*t)*FD(ω), -Inf, Inf))[1], tind)

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
ENV["GKS_DOUBLE_BUF"] = "1"

using GR

const x = [ 0, 0.850651, 0.850651, -0.850651, -0.850651, -0.525731, 0.525731, 0.525731, -0.525731, 0, 0, 0 ]
const y = [ -0.525731, 0, 0, 0, 0, 0.850651, 0.850651, -0.850651, -0.850651, -0.525731, 0.525731, 0.525731 ]
const z = [ 0.850651, 0.525731, -0.525731, -0.525731, 0.525731, 0, 0, 0, 0, -0.850651, -0.850651, 0.850651 ]

const connections = [
  [ 3, 3, 7, 2 ],
  [ 3, 8, 3, 2 ],
  [ 3, 5, 6, 4 ],
  [ 3, 4, 9, 5 ],
  [ 3, 6, 12, 7 ],
  [ 3, 7, 11, 6 ],
  [ 3, 11, 3, 10 ],
  [ 3, 10, 4, 11 ],
  [ 3, 9, 10, 8 ],
  [ 3, 8, 1, 9 ],
  [ 3, 1, 2, 12 ],
  [ 3, 12, 5, 1 ],
  [ 3, 3, 11, 7 ],
  [ 3, 7, 12, 2 ],
  [ 3, 6, 11, 4 ],
  [ 3, 5, 12, 6 ],
  [ 3, 8, 10, 3 ],
  [ 3, 2, 1, 8 ],
  [ 3, 10, 9, 4 ],
  [ 3, 9, 1, 5 ]
]

# use GR's distinct colors
const colors = signed.([inqcolor(979+i) | 0xdf000000 for i in 1:size(connections)[1]])

setviewport(0, 1, 0, 1)
setwindow(-1, 1, -1, 1)
setwindow3d(-1, 1, -1, 1, -1, 1)

for angle in 0:360
  clearws()
  setspace3d(30 + angle, 80, 0, 0)

  polygonmesh3d(x, y, z, vcat(connections...), colors)

  updatews()
  sleep(0.01)
end


text:
######################################################## model
F(z::BorderedArray,θ::AbstractVector) = F(z.u,(θ=θ,p=z.p))
function F(u::AbstractVector,parameters::NamedTuple)

	@unpack θ,p = parameters
	μ₁,μ₂, a₁,a₂, k = θ

	f = first(u)*first(p)*first(θ)
	F = similar(u,typeof(f))

	F[1] = ( 10^a₁ + (p*u[2])^2 ) / ( 1 + (p*u[2])^2 ) - u[1]*10^μ₁
	F[2] = ( 10^a₂ + (k*u[1])^2 ) / ( 1 + (k*u[1])^2 ) - u[2]*10^μ₂

	return F
end

######################################################### targets and initial guess
X = StateSpace( 2, 3:0.001:7, [4,5] )
θ = SizedVector{5}(0.0,0.0,-1.0,-1.0,2.0)


text:
## arquivo modelo

using IL_simulations, PDBTools

top_dir     =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/VSIL/ITP"
pdb_dir     =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/VSIL/PDB"
input_dir   =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/ubiquitin_files"
top_ion_dir =  "/home/viniciusp/Documents/doutorado/ANALYSE/Repository/IL_simulations.jl/ion_files"

atoms   = readPDB("$input_dir/ubq.pdb")
MMP     = mass(atoms)

cation1 = "BMIM"
anion1  = "NO3"
anion2  = "DCA"

# function to calculate the mean MDDF giving a set of MD trajectories
#function molarmass(anion, cation, pdb_dir)
#  atcation = readPDB("$pdb_dir/$(cation)_VSIL.pdb")
#  atanion  = readPDB("$pdb_dir/$(anion)_VSIL.pdb")
#  MMIL = mass(atcation) + mass(atanion)
#  println("molar mass of the ânion is = ", mass(atanion))
#  println("molar mass of the cation is = ", mass(atcation))
#  return MMIL
#end

MMIL  = molarmass(anion1, cation1, pdb_dir)
#MMIL2 = molarmass(anion2, cation1, pdb_dir)
c       =  0.10

## test with one IL/salt and a protein/polymer
#data = IL_simulations.Data_il(protein="$input_dir/ubq.pdb", MMP = MMP, cation = cation1, anion = anion1,MM = MMIL,c=c)
#nions, nwater, sides = prot_elec(data)
#IL_simulations.top(dict, top_dir, input_dir, data, nions, nwater)
#pack_input(data,pdb_dir,nions,nwater,sides)

função para escrever o arquivo

# test with one IL/salt
data = IL_simulations.Data_elec(cation = cation1, anion = anion1, MM = MMIL,c=c)
nions, nwater, sides = IL_simulations.sol_elec(data)
IL_simulations.topelec(dict, top_dir, top_ion_dir, data, nions, nwater)
IL_simulations.pack_input_sol(data,pdb_dir,nions,nwater,sides)
IL_simulations.analyzeIN(pdb_dir,data)

## test with two ILs/salts and a protein/polymer
#data = IL_simulations.Data_il(protein="$input_dir/ubq.pdb", MMP = MMP, cation = cation1, anion = anion1,MM = MM1,c=c)
#nions, nwater, sides = prot_melec(data)
#IL_simulations.top(dict, top_dir, input_dir, data, nions, nwater)
#pack_input(data,pdb_dir,nions,nwater,sides)
#
## test with two ILs/salts
#data = IL_simulations.Data_il(protein="$input_dir/ubq.pdb", MMP = MMP, cation = cation1, anion = anion1,MM = MM1,c=c)
#nions, nwater, sides = sol_melec(data)
#IL_simulations.top(dict, top_dir, input_dir, data, nions, nwater)
#pack_input(data,pdb_dir,nions,nwater,sides)


text:
##
using QuasiMonteCarlo: LowDiscrepancySample, sample
using Plots

function kern_sqexp(x)
    exp(-x^2/2)
end

function cov!(K, kernel, x1, x2, l)
    for (k1, xv1) in enumerate(x1)
        for (k2, xv2) in enumerate(x2)
        end
    end
end

lb = [0.0, -0.5]  # Lower bounds
ub = [4.0, 3.5]   # Upper bounds
n = 100           # Number of samples

halton = LowDiscrepancySample([3, 5])  # Halton with prime number sequence
xtrain = sample(n, lb, ub, halton)

scatter(xtrain[1,:], xtrain[2,:])
##

# f(x) = cos.(x[1,:]).*sin.(x[2,:])
# ytrain = f(xtrain)

# #scatter3d(xtrain[1,:], xtrain[2,:], ytrain, marker_z=ytrain)
# #savefig("test.png")

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
using Match
using Base.Test

include("testtypes.jl")

import Base: show

# Type matching
test1(item) = @match item begin
    n::Int                       => "Integers are awesome!"
    str::AbstractString          => "Strings are the best"
    m::Dict{Int, AbstractString} => "Ints for Strings?"
    d::Dict                      => "A Dict! Looking up a word?"
    _                            => "Something unexpected"
end

d = Dict{Int,AbstractString}(1=>"a",2=>"b")

@test test1(66)     == "Integers are awesome!"
@test test1("abc")  == "Strings are the best"
@test test1(d)      == "Ints for Strings?"
@test test1(Dict()) == "A Dict! Looking up a word?"
@test test1(2.0)    == "Something unexpected"

# Pattern extraction
# inspired by http://thecodegeneral.wordpress.com/2012/03/25/switch-statements-on-steroids-scala-pattern-matching/

# type Address
#     street::AbstractString
#     city::AbstractString
#     zip::AbstractString
# end

# type Person
#     firstname::AbstractString
#     lastname::AbstractString
#     address::Address
# end

test2(person) = @match person begin
    Person("Julia", lastname,  _) => "Found Julia $lastname"
    Person(firstname, "Julia", _) => "$firstname Julia was here!"
    Person(firstname, lastname ,Address(_, "Cambridge", zip)) => "$firstname $lastname lives in zip $zip"
    _::Person  => "Unknown person!"
end

@test test2(Person("Julia", "Robinson", Address("450 Serra Mall", "Stanford", "94305")))         == "Found Julia Robinson"
@test test2(Person("Gaston", "Julia",   Address("1 rue Victor Cousin", "Paris", "75005")))       == "Gaston Julia was here!"
@test test2(Person("Edwin", "Aldrin",   Address("350 Memorial Dr", "Cambridge", "02139")))       == "Edwin Aldrin lives in zip 02139"
@test test2(Person("Linus", "Pauling",  Address("1200 E California Blvd", "Pasadena", "91125"))) == "Unknown person!"  # Really?

# Guards, pattern extraction
# translated from Scala Case-classes http://docs.scala-lang.org/tutorials/tour/case-classes.html

##
## Untyped lambda calculus definitions
##

# abstract Term

# immutable Var <: Term
#     name::AbstractString
# end

# immutable Fun <: Term
#     arg::AbstractString
#     body::Term
# end

# immutable App <: Term
#     f::Term
#     v::Term
# end

# scala defines these automatically...
import Base.==

==(x::Var, y::Var) = x.name == y.name
==(x::Fun, y::Fun) = x.arg == y.arg && x.body == y.body
==(x::App, y::App) = x.f == y.f && x.v == y.v

# Not really the Julian way
function show(io::IO, term::Term)
    @match term begin
       Var(n)    => print(io, n)
       Fun(x, b) => begin
                        print(io, "^$x.")
                        show(io, b)
                    end
       App(f, v) => begin
                        print(io, "(")
                        show(io, f)
                        print(io, " ")
                        show(io, v)
                        print(io, ")")
                    end
    end
end

# Guard test is here
function is_identity_fun(term::Term)
   @match term begin
     Fun(x, Var(y)), if x == y end => true
     _ => false
   end
end

id = Fun("x", Var("x"))
t = Fun("x", Fun("y", App(Var("x"), Var("y"))))

let io = IOBuffer()
    show(io, id)
    @assert takebuf_string(io) == "^x.x"
    show(io, t)
    @assert takebuf_string(io) == "^x.^y.(x y)"
    @assert is_identity_fun(id)
    @assert !is_identity_fun(t)
end

# Test single terms

myisodd(x::Int) = @match(x, i => i%2==1)
@assert filter(myisodd, 1:10) == filter(isodd, 1:10) == [1,3,5,7,9]

# Alternatives, Guards

function parse_arg(arg::AbstractString, value::Any=nothing)
   @match (arg, value) begin
      ("-l",              lang),   if lang != nothing end => "Language set to $lang"
      ("-o" || "--optim", n::Int),      if 0 < n <= 5 end => "Optimization level set to $n"
      ("-o" || "--optim", n::Int)                         => "Illegal optimization level $(n)!"
      ("-h" || "--help",  nothing)                        => "Help!"
      bad                                                 => "Unknown argument: $bad"
   end
end

@assert parse_arg("-l", "eng")  == "Language set to eng"
@assert parse_arg("-l")         == "Unknown argument: (\"-l\",nothing)"
@assert parse_arg("-o", 4)      == "Optimization level set to 4"
@assert parse_arg("--optim", 5) == "Optimization level set to 5"
@assert parse_arg("-o", 0)      == "Illegal optimization level 0!"
@assert parse_arg("-o", 1.0)    == "Unknown argument: (\"-o\",1.0)"

@assert parse_arg("-h") == parse_arg("--help") == "Help!"

# Regular Expressions

# TODO: Fix me!

# Note: the following test only works because Ipv4Addr and EmailAddr
# are (module-level) globals!

# const Ipv4Addr = r"(\d{1,3})\.(\d{1,3})\.(\d{1,3})\.(\d{1,3})"
# const EmailAddr = r"\b([A-Z0-9._%+-]+)@([A-Z0-9.-]+\.[A-Z]{2,4})\b"i

# function regex_test(str)
#     ## Defining these in the function doesn't work, because the macro
#     ## (and related functions) don't have access to the local
#     ## variables.

#     # Ipv4Addr = r"(\d{1,3})\.(\d{1,3})\.(\d{1,3})\.(\d{1,3})"
#     # EmailAddr = r"\b([A-Z0-9._%+-]+)@([A-Z0-9.-]+\.[A-Z]{2,4})\b"i

#     @match str begin
#         Ipv4Addr(_, _, octet3, _),       if int(octet3) > 30 end => "IPv4 address with octet 3 > 30"
#         Ipv4Addr()                                               => "IPv4 address"

#         EmailAddr(_,domain), if endswith(domain, "ucla.edu") end => "UCLA email address"
#         EmailAddr                                                => "Some email address"

#         r"MCM.*"                                                 => "In the twentieth century..."

#         _                                                        => "No match"
#     end
# end

# @assert regex_test("128.97.27.37")                 == "IPv4 address"
# @assert regex_test("96.17.70.24")                  == "IPv4 address with octet 3 > 30"

# @assert regex_test("beej@cs.ucla.edu")             == "UCLA email address"
# @assert regex_test("beej@uchicago.edu")            == "Some email address"

# @assert regex_test("MCMLXXII")                     == "In the twentieth century..."
# @assert regex_test("Open the pod bay doors, HAL.") == "No match"

# Pattern extraction from arrays

# extract first, rest from array
# (b is a subarray of the original array)
@assert @match([1:4;], [a,b...])                             == (1,[2,3,4])
@assert @match([1:4;], [a...,b])                             == ([1,2,3],4)
@assert @match([1:4;], [a,b...,c])                           == (1,[2,3],4)

# match particular values at the beginning of a vector
@assert @match([1:10;], [1,2,a...])                          == [3:10;]
@assert @match([1:10;], [1,a...,9,10])                       == [2:8;]

# match / collect columns
@assert @match([1 2 3; 4 5 6], [a b...])                    == ([1,4] , [2 3; 5 6])
@assert @match([1 2 3; 4 5 6], [a... b])                    == ([1 2; 4 5] , [3,6])
@assert @match([1 2 3; 4 5 6], [a b c])                     == ([1,4] , [2,5] , [3,6])
@assert @match([1 2 3; 4 5 6], [[1,4] a b])                 == ([2,5] , [3,6])

@assert @match([1 2 3 4; 5 6 7 8], [a b... c])              == ([1,5] , [2 3; 6 7] , [4,8])

if VERSION < v"0.5.0-dev"
    # match / collect rows
    @assert @match([1 2 3; 4 5 6], [a, b])                      == ([1 2 3], [4 5 6])
    @assert @match([1 2 3; 4 5 6], [[1 2 3], a])                ==  [4 5 6]             # TODO: don't match this
    @assert @match([1 2 3; 4 5 6], [1 2 3; a])                  ==  [4 5 6]

    @assert @match([1 2 3; 4 5 6; 7 8 9], [a, b...])            == ([1 2 3], [4 5 6; 7 8 9])
    @assert @match([1 2 3; 4 5 6; 7 8 9], [a..., b])            == ([1 2 3; 4 5 6], [7 8 9])
    @assert @match([1 2 3; 4 5 6; 7 8 9], [1 2 3; a...])        ==  [4 5 6; 7 8 9]

    @assert @match([1 2 3; 4 5 6; 7 8 9; 10 11 12], [a,b...,c]) == ([1 2 3], [4 5 6; 7 8 9], [10 11 12])

    # match invidual positions
    @assert @match([1 2; 3 4], [1 a; b c])                      == (2,3,4)
    @assert @match([1 2; 3 4], [1 a; b...])                     == (2,[3 4])

    @assert @match([ 1  2  3  4
                     5  6  7  8
                     9 10 11 12
                    13 14 15 16
                    17 18 19 20 ],

                    [1      a...
                     b...
                     c... 15 16
                     d 18 19 20])                               == ([2 3 4], [5 6 7 8; 9 10 11 12], [13 14], 17)

else
    # match / collect rows
    @assert @match([1 2 3; 4 5 6], [a, b])                      == ([1,2,3], [4,5,6])
    @assert @match([1 2 3; 4 5 6], [[1,2,3], a])                ==  [4,5,6]             # TODO: don't match this
    #@assert @match([1 2 3; 4 5 6], [1 2 3; a])                  ==  [4,5,6]

    @assert @match([1 2 3; 4 5 6; 7 8 9], [a, b...])            == ([1,2,3], [4 5 6; 7 8 9])
    @assert @match([1 2 3; 4 5 6; 7 8 9], [a..., b])            == ([1 2 3; 4 5 6], [7,8,9])
    #@assert @match([1 2 3; 4 5 6; 7 8 9], [1 2 3; a...])        ==  [4 5 6; 7 8 9]

    #@assert @match([1 2 3; 4 5 6; 7 8 9; 10 11 12], [a,b...,c]) == ([1,2,3], [4 5 6; 7 8 9], [10 11 12])

    # match invidual positions
    #@assert @match([1 2; 3 4], [1 a; b c])                      == (2,3,4)
    #@assert @match([1 2; 3 4], [1 a; b...])                     == (2,[3,4])

    # @assert @match([ 1  2  3  4
    #                  5  6  7  8
    #                  9 10 11 12
    #                 13 14 15 16
    #                 17 18 19 20 ],
    #
    #                 [1      a...
    #                  b...
    #                  c... 15 16
    #                  d 18 19 20])                               == ([2,3,4], [5 6 7 8; 9 10 11 12], [13,14], 17)

end

# match 3D arrays
m = reshape([1:8;], (2,2,2))
@assert @match(m, [a b])                                    == ([1 3; 2 4], [5 7; 6 8])
if VERSION < v"0.5.0-dev"
    @assert @match(m, [[1 a; b c] d])                           == (3,2,4,[5 7; 6 8])
end

# match against an expression
function get_args(ex::Expr)
    @match ex begin
        Expr(:call, [:+, args...], _) => args
        _ => "None"
    end
end

@assert get_args(Expr(:call, :+, :x, 1)) == [:x, 1]

# Zach Allaun's fizzbuzz (https://github.com/zachallaun/Match.jl#awesome-fizzbuzz)

function fizzbuzz(range::Range)
    io = IOBuffer()
    for n in range
        @match (n%3, n%5) begin
            (0,0) => print(io, "fizzbuzz ")
            (0,_) => print(io, "fizz ")
            (_,0) => print(io, "buzz ")
            (_,_) => print(io, n, ' ')
        end
    end
    takebuf_string(io)
end

@assert fizzbuzz(1:15) == "1 2 fizz 4 buzz fizz 7 8 fizz buzz 11 fizz 13 14 fizzbuzz "

# Zach Allaun's "Balancing Red-Black Trees" (https://github.com/zachallaun/Match.jl#balancing-red-black-trees)

abstract RBTree

immutable Leaf <: RBTree
end

immutable Red <: RBTree
    value
    left::RBTree
    right::RBTree
end

immutable Black <: RBTree
    value
    left::RBTree
    right::RBTree
end

function balance(tree::RBTree)
    @match tree begin
        ( Black(z, Red(y, Red(x, a, b), c), d)
         || Black(z, Red(x, a, Red(y, b, c)), d)
         || Black(x, a, Red(z, Red(y, b, c), d))
         || Black(x, a, Red(y, b, Red(z, c, d)))) => Red(y, Black(x, a, b),
                                                         Black(z, c, d))
        tree => tree
    end
end

@assert balance(Black(1, Red(2, Red(3, Leaf(), Leaf()), Leaf()), Leaf())) ==
            Red(2, Black(3, Leaf(), Leaf()),
                Black(1, Leaf(), Leaf()))

function num_match(n)
    @match n begin
        0      => "zero"
        1 || 2 => "one or two"
        3:10   => "three to ten"
        _      => "something else"
    end
end

@assert num_match(0) == "zero"
@assert num_match(2) == "one or two"
@assert num_match(4) == "three to ten"
@assert num_match(12) == "something else"
@assert num_match("hi") == "something else"
@assert num_match('c') == "something else"

# Interpolation of matches in quoted expressions
test_interp(item) = @match item begin
    [a, b] => :($a + $b)
end
@assert test_interp([1, 2]) == :(1 + 2)


text:
module NCCL

using Printf

import CUDAdrv

const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
isfile(ext) || error("NCCL.jl has not been built, please run Pkg.build(\"NCCL\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const libnccl = nothing
end

include("base.jl")
include("communicator.jl")

@show id = UniqueID()
@show comm = Communicator(length(CUDAdrv.devices()), id, 0)
@show comm = Communicator(CUDAdrv.devices())

end


text:
module Stan

using Mamba

"""The directory which contains the executable `bin/stanc`. Inferred
from `Main.CMDSTAN_HOME` or `ENV["CMDSTAN_HOME"]` when available. Use
`set_CMDSTAN_HOME!` to modify."""
CMDSTAN_HOME=""

function __init__()
    global CMDSTAN_HOME = if isdefined(Main, :CMDSTAN_HOME)
        eval(Main, :CMDSTAN_HOME)
    elseif haskey(ENV, "CMDSTAN_HOME")
        ENV["CMDSTAN_HOME"]
    else
        warn("Environment variable CMDSTAN_HOME not found. Use set_CMDSTAN_HOME!.")
        ""
    end
end

"""Set the path for `CMDSTAN`.

Example: `set_CMDSTAN_HOME!(homedir() * "/src/src/cmdstan-2.11.0/")`"""
set_CMDSTAN_HOME!(path) = global CMDSTAN_HOME=path

if !isdefined(Main, :JULIA_SVG_BROWSER)
    JULIA_SVG_BROWSER = ""
    try
        JULIA_SVG_BROWSER = ENV["JULIA_SVG_BROWSER"]
    catch e
        println("Environment variable JULIA_SVG_BROWSER not found.")
        JULIA_SVG_BROWSER = ""
    end
end

include("stanmodel.jl")
include("stancode.jl")
include("utilities.jl")

export
# from this file
set_CMDSTAN_HOME!,
# From stancode.jl
stan,
stan_summary,
read_stanfit,
read_stanfit_samples,
CMDSTAN_HOME,
JULIA_SVG_BROWSER,
# From stanmodel.jl
Stanmodel,
Data,
RNG,
Output,
# From sampletype.jl
Sample,
Hmc,
diag_e,
unit_e,
dense_e,
Engine,
Nuts,
Static,
Metrics,
SampleAlgorithm,
Fixed_param,
Adapt,
# From optimizetype.jl
Optimize,
Lbfgs,
Bfgs,
Newton,
# From diagnosetype.jl
Diagnose,
Diagnostics,
Gradient,
# From variationaltype.jl
Variational,
# From cmdline.jl
cmdline

end # module


text:
"""
String utilities and performance improvements

Copyright 2016 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module StringUtils
export @u_str, @sinterpolate
export s_unescape_string, s_escape_string, s_print_unescaped, s_print_escaped

include("literals.jl")
include("unicodenames.jl")
include("htmlnames.jl")

# From Formatting.jl
import Base.show

export FormatSpec, FormatExpr, printfmt, printfmtln, format, generate_formatter
export pyfmt, cfmt, fmt, fmt_default, fmt_default!, reset!, defaultSpec

include("cformat.jl" )
include("fmtspec.jl")
include("fmtcore.jl")
include("formatexpr.jl")
include("fmt.jl")

end


text:
## floating-point functions ##

copysign(x::Float64, y::Float64) = box(Float64,copysign_float(unbox(Float64,x),unbox(Float64,y)))
copysign(x::Float32, y::Float32) = box(Float32,copysign_float(unbox(Float32,x),unbox(Float32,y)))
copysign(x::Float32, y::Real) = copysign(x, float32(y))
copysign(x::Float64, y::Real) = copysign(x, float64(y))
@vectorize_2arg Real copysign

flipsign(x::Float64, y::Float64) = box(Float64,xor_int(unbox(Float64,x),and_int(unbox(Float64,y),0x8000000000000000)))
flipsign(x::Float32, y::Float32) = box(Float32,xor_int(unbox(Float32,x),and_int(unbox(Float32,y),0x80000000)))
flipsign(x::Float32, y::Real) = flipsign(x, float32(y))
flipsign(x::Float64, y::Real) = flipsign(x, float64(y))
@vectorize_2arg Real flipsign

signbit(x::Float64) = signbit(reinterpret(Int64,x))
signbit(x::Float32) = signbit(reinterpret(Int32,x))
signbit(x::Float16) = signbit(reinterpret(Int16,x))

maxintfloat(::Type{Float64}) = 9007199254740992.
maxintfloat(::Type{Float32}) = float32(16777216.)
maxintfloat(::Type{Float16}) = float16(2048f0)
maxintfloat{T<:FloatingPoint}(x::T)  = maxintfloat(T)
maxintfloat() = maxintfloat(Float64)

isinteger(x::FloatingPoint) = (trunc(x)==x)&isfinite(x)

num2hex(x::Float16) = hex(reinterpret(Uint16,x), 4)
num2hex(x::Float32) = hex(box(Uint32,unbox(Float32,x)),8)
num2hex(x::Float64) = hex(box(Uint64,unbox(Float64,x)),16)

function hex2num(s::String)
    if length(s) <= 8
        return box(Float32,unbox(Int32,parseint(Int32,s,16)))
    end
    return box(Float64,unbox(Int64,parseint(Int64,s,16)))
end

@vectorize_1arg Real iround
@vectorize_1arg Real itrunc
@vectorize_1arg Real ifloor
@vectorize_1arg Real iceil

@vectorize_1arg Number abs
@vectorize_1arg Number abs2
@vectorize_1arg Number angle

@vectorize_1arg Number isnan
@vectorize_1arg Number isinf
@vectorize_1arg Number isfinite

iround{T<:Integer,R<:Real}(::Type{T}, x::AbstractArray{R,1}) = [ iround(T, x[i]) for i = 1:length(x) ]
iround{T<:Integer,R<:Real}(::Type{T}, x::AbstractArray{R,2}) = [ iround(T, x[i,j]) for i = 1:size(x,1), j = 1:size(x,2) ]
iround{T<:Integer,R<:Real}(::Type{T}, x::AbstractArray{R}) = reshape([ iround(T, x[i]) for i = 1:length(x) ], size(x))

# adapted from Matlab File Exchange roundsd: http://www.mathworks.com/matlabcentral/fileexchange/26212
# for round, og is the power of 10 relative to the decimal point
# for signif, og is the absolute power of 10
# digits and base must be integers, x must be convertable to float

function _signif_og(x, digits, base)
    if base == 10
        oftype(x, 10. ^ floor(log10(abs(x)) - digits + 1.))
    elseif base == 2
        oftype(x, 2. ^ floor(log2(abs(x)) - digits + 1.))
    else
        oftype(x, float(base) ^ floor(log2(abs(x))/log2(base) - digits + 1.))
    end
end

function signif(x, digits::Integer, base::Integer=10)
    if digits < 0
        throw(DomainError())
    end
    x = float(x)
    if x==0 || !isfinite(x)
        return x
    end
    og = _signif_og(x, digits, base)
    round(x/og) * og
end

for f in (:round, :ceil, :floor, :trunc)
    @eval begin
        function ($f)(x, digits::Integer, base::Integer=10)
            x = float(x)
            og = oftype(eltype(x),base)^digits
            ($f)(x * og) / og
        end
    end
end

# isapprox: Tolerant comparison of floating point numbers
function isapprox(x::FloatingPoint, y::FloatingPoint; rtol::Real=rtoldefault(x,y), atol::Real=atoldefault(x,y))
    (isinf(x) || isinf(y)) ? x == y : abs(x-y) <= atol + rtol.*max(abs(x), abs(y))
end

# promotion of non-floats
isapprox(x::Real, y::FloatingPoint; rtol::Real=rtoldefault(x, y), atol::Real=atoldefault(x, y)) = isapprox(promote(x, y)...; rtol=rtol, atol=atol)
isapprox(x::FloatingPoint, y::Real; rtol::Real=rtoldefault(x, y), atol::Real=atoldefault(x, y)) = isapprox(promote(x, y)...; rtol=rtol, atol=atol)

# other real numbers
isapprox(x::Real, y::Real; rtol::Real=0, atol::Real=0) = abs(x-y) <= atol

# complex numbers
isapprox(z::Complex, w::Complex; rtol::Real=rtoldefault(abs(z), abs(w)), atol::Real=atoldefault(abs(z), abs(w))) = abs(z-w) <= atol + rtol*max(abs(z), abs(w))

# real-complex combinations
isapprox(x::Real, z::Complex; rtol::Real=rtoldefault(x, abs(z)), atol::Real=atoldefault(x, abs(z))) = isapprox(complex(x), z; rtol=rtol, atol=atol)
isapprox(z::Complex, x::Real; rtol::Real=rtoldefault(x, abs(z)), atol::Real=atoldefault(x, abs(z))) = isapprox(complex(x), z; rtol=rtol, atol=atol)

# default tolerance arguments
rtoldefault(x::FloatingPoint, y::FloatingPoint) = cbrt(max(eps(x), eps(y)))
atoldefault(x::FloatingPoint, y::FloatingPoint) = sqrt(max(eps(x), eps(y)))

# promotion of non-floats
for fun in (:rtoldefault, :atoldefault)
    @eval begin
        ($fun)(x::Real, y::FloatingPoint) = ($fun)(promote(x,y)...)
        ($fun)(x::FloatingPoint, y::Real) = ($fun)(promote(x,y)...)
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
@testset "c\"... (*)(...)\"" begin
	@eval c``

	asprintf = @eval c"extern int asprintf(char const **str, char *fmt, ...);"
	asprintf = asprintf[]  # get pointer to function
	@test eltype(asprintf) <: eltype(c"int (*)(char const **, char *, ...)")

	f(x) = asprintf(x, "%s i%c %d great test of CBinding.jl v%3.1f%c\n", "This", 's', 0x01, 0.1, '!')
	str = Ref(c"char const *"(C_NULL))
	expect = "This is 1 great test of CBinding.jl v0.1!\n"
	bytes = 0
	for i in 1:100
		b = @allocated len = f(str)
		bytes += i > 1 ? b : 0
		@test len == length(expect)
		@test unsafe_string(str[]) == expect
		Libc.free(str[])
		str[] = C_NULL
	end
	@test_broken bytes == 0
	@test_throws MethodError asprintf(1234)
	@test_throws MethodError asprintf(1234, "still wrong")

	f2a = @eval c"int jl_ver_major();"
	f2a = f2a[]
	@test eltype(f2a) <: eltype(c"int (*)()")
	@test typeof(f2a()) === Cint
	@test f2a() isa Cint
	@test f2a() == Base.VERSION.major
	@test_throws MethodError f2a("no arguments, please!")

	f2b = @eval c"int jl_ver_minor();"
	f2b = f2b[]
	@test eltype(f2b) <: eltype(c"int (*)()")
	@test typeof(f2b()) === Cint
	@test f2b() isa Cint
	@test f2b() == Base.VERSION.minor
	@test_throws MethodError f2b("no arguments, please!")

	function add(val1::Cint, val2::Cint)::Cint
		return val1+val2
	end
	Cadd = c"int (*)(int, int)"(add)
	@test eltype(Cadd) <: eltype(c"int (*)(int, int)")
	@test typeof(Cadd(Cint(10), Cint(3))) === typeof(add(Cint(10), Cint(3)))
	@test Cadd(Cint(10), Cint(3)) == add(Cint(10), Cint(3))

	c"CB" = c"int (*)(int, int)"
	Capply = c"int (*)(CB, int, int)"() do func::c"CB", val1::Cint, val2::Cint
		return func(val1, val2)::Cint
	end
	@test eltype(Capply) <: eltype(c"int (*)(CB, int, int)")
	@test typeof(Capply(Cadd, Cint(10), Cint(3))) === typeof(add(Cint(10), Cint(3)))
	@test Capply(Cadd, Cint(10), Cint(3)) == add(Cint(10), Cint(3))
end


text:
export StackFrames

using CircularArrayBuffers: CircularArrayBuffer

"""
    StackFrames(::Type{T}=Float32, d::Int...)

Use a pre-initialized [`CircularArrayBuffer`](@ref) to store the latest several states specified by `d`. Before processing any observation, the buffer is filled with `zero{T}
by default.
"""
struct StackFrames{T,N} <: AbstractArray{T,N}
    buffer::CircularArrayBuffer{T,N}
end

Base.size(x::StackFrames) = size(x.buffer)
Base.getindex(x::StackFrames, I...) = getindex(x.buffer, I...)
Base.IndexStyle(x::StackFrames) = IndexStyle(x.buffer)

StackFrames(d::Int...) = StackFrames(Float32, d...)

function StackFrames(::Type{T}, d::Vararg{Int,N}) where {T,N}
    b = CircularArrayBuffer{T}(d...)
    append!(b, zeros(T, d...))
    StackFrames(b)
end

function (p::StackFrames{T,N})(state::AbstractArray) where {T,N}
    push!(p.buffer, state)
    p
end

function RLBase.reset!(p::StackFrames{T,N}) where {T,N}
    fill!(p.buffer, zero(T))
    p
end

"""
When pushing a `StackFrames` into a `CircularArrayBuffer` of the same dimension,
only the latest frame is pushed. If the `StackFrames` is one dimension lower,
then it is treated as a general `AbstractArray` and is pushed in as a frame.
"""
function Base.push!(cb::CircularArrayBuffer{T,N}, p::StackFrames{T,N}) where {T,N}
    push!(cb, selectdim(p.buffer, N, size(p.buffer, N)))
end


text:
#####
##### Informal benchmarks
#####

using FunctionalTables

N = 10000
a = rand(1:100, N)
b = rand('a':'z', N)
c = rand(Float64, N)

ft = sort(FunctionalTable((a = a, b = b, c = c)), (:a, :b))

m(ft) = map((_, ft) -> ft, by(ft, (:a, :b)))

@code_warntype m(ft)
@time m(ft)
# f77a599: 0.215s
# 0.09 after by optimization
# 0.07 after collect_columns rewrite

f(ft) = iterate(by(ft, (:a, :b)))
@code_warntype f(ft)

using FunctionalTables
using FunctionalTables: collect_columns, SINKCONFIG, empty_sinks, collect_columns!, start_sinks

@code_warntype collect_columns(SINKCONFIG, [(a = 1, b = 1)], TrustOrdering(()), NamedTuple{(),Tuple{}})

@code_warntype start_sinks(SINKCONFIG, (a = 1, ), NamedTuple{(:a, ), Tuple{Int}})

A = fill((a = 1, ), 3)
elt, state = iterate(A)

@code_warntype collect_columns(SINKCONFIG, A, TrustOrdering(), typeof(NamedTuple()))
@code_warntype collect_columns!(sinks, 0, SINKCONFIG, A, VerifyOrdering(), elt, state)


text:
using CSV, DataFrames, Plots
using Statistics
# %% Compare ILC results
ilc = CSV.read("ilc_results_sosatccatp_300sims_seed5147.csv", DataFrame)
println("ILC:")
println(
    "- SAT only: r = ",
    mean(ilc.r_SOonly_wo_FGmarg),
    " ± ",
    std(ilc.r_SOonly_wo_FGmarg),
)
println("- SAT+FYST: r = ", mean(ilc.r_wo_FGmarg), " ± ", std(ilc.r_wo_FGmarg))
p = scatter(
    ilc.r_wo_FGmarg,
    ilc.r_SOonly_wo_FGmarg,
    xlab = "Tensor-to-scalar ratio, SO-SAT + FYST",
    ylab = "Tensor-to-scalar ratio, SO-SAT only",
    lab = "",
    legend = :topright,
    xlims = [-0.003, 0.01],
    ylims = [-0.003, 0.01],
    aspect_ratio = :equal,
    title = "Internal Linear Combination",
)
p = plot!(
    ilc.r_wo_FGmarg,
    ilc.r_wo_FGmarg,
    lw = 2,
    lab = "r(SAT only) = r(SAT+FYST)",
)
p = plot!(
    ilc.r_wo_FGmarg,
    ilc.r_wo_FGmarg .+ 1.14e-3,
    lw = 2,
    lab = "r(SAT only) = r(SAT+FYST) + 1.14e-3",
)
p = histogram!(
    ilc.r_SOonly_wo_FGmarg,
    inset = (1, bbox(0.5, 0.5, 0.5, 0.4)),
    yticks = nothing,
    subplot = 2,
    bg_inside = nothing,
    lab = "SAT only",
    α = 0.7,
    legend = :best,
    legendfontsize = 5,
    xlims = [-0.005, 0.01],
)
p = histogram!(p[2], ilc.r_wo_FGmarg, α = 0.7, lab = "SATFYST")
p = vline!([0], c = :grey, lw = 2, ls = :dot, lab = "")
p = hline!([0], c = :grey, lw = 2, ls = :dot, lab = "")
savefig("scatterplot_so_vs_soccatp_ilc.pdf")
display(p)

# %% Compare Parametric Maximum Likelihood results
milca = CSV.read("milca_results_sosatccatp_300sims_seed5147.csv", DataFrame)
println("MILCA:")
println(
    "- SAT only: r = ",
    mean(milca.r_SOonly_wo_FGmarg),
    " ± ",
    std(milca.r_SOonly_wo_FGmarg),
)
println(
    "- SAT+FYST: r = ",
    mean(milca.r_wo_FGmarg),
    " ± ",
    std(milca.r_wo_FGmarg),
)
p = scatter(
    milca.r_wo_FGmarg,
    milca.r_SOonly_wo_FGmarg,
    xlab = "Tensor-to-scalar ratio, SO-SAT + FYST",
    ylab = "Tensor-to-scalar ratio, SO-SAT only",
    lab = "",
    legend = :topright,
    xlims = [-0.003, 0.01],
    ylims = [-0.003, 0.01],
    aspect_ratio = :equal,
    title = "Parametric Maximum Likelihood",
)
p = plot!(
    milca.r_wo_FGmarg,
    milca.r_wo_FGmarg,
    lw = 2,
    lab = "r(SAT only) = r(SAT+FYST)",
)
p = plot!(
    milca.r_wo_FGmarg,
    milca.r_wo_FGmarg .+ 1.095e-3,
    lw = 2,
    lab = "r(SAT only) = r(SAT+FYST) + 1.10e-3",
)
p = histogram!(
    milca.r_SOonly_wo_FGmarg,
    inset = (1, bbox(0.5, 0.5, 0.5, 0.4)),
    yticks = nothing,
    subplot = 2,
    bg_inside = nothing,
    lab = "SAT only",
    α = 0.7,
    legend = :best,
    legendfontsize = 5,
    xlims = [-0.005, 0.01],
)
p = histogram!(p[2], milca.r_wo_FGmarg, α = 0.7, lab = "SATFYST")
p = vline!([0], c = :grey, lw = 2, ls = :dot, lab = "")
p = hline!([0], c = :grey, lw = 2, ls = :dot, lab = "")
savefig("scatterplot_so_vs_soccatp_milca.pdf")
display(p)

# %% Compare ILC and Parametric Maximum Likelihood results
p = scatter(
    ilc.r_wo_FGmarg,
    milca.r_wo_FGmarg,
    xlab = "Tensor-to-scalar ratio, ILC",
    ylab = "Tensor-to-scalar ratio, ML",
    lab = "",
    legend = :topright,
    xlims = [-0.003, 0.01],
    ylims = [-0.003, 0.01],
    aspect_ratio = :equal,
    title = "ILC vs Parametric ML",
)
p = plot!(ilc.r_wo_FGmarg, ilc.r_wo_FGmarg, lw = 2, lab = "r(ML) = r(ILC)")
p = histogram!(
    milca.r_wo_FGmarg,
    inset = (1, bbox(0.5, 0.5, 0.5, 0.4)),
    yticks = nothing,
    subplot = 2,
    bg_inside = nothing,
    lab = "ML",
    α = 0.7,
    legend = :best,
    legendfontsize = 5,
    xlims = [-0.005, 0.01],
)
p = histogram!(p[2], ilc.r_wo_FGmarg, α = 0.7, lab = "ILC")
p = vline!([0], c = :grey, lw = 2, ls = :dot, lab = "")
p = hline!([0], c = :grey, lw = 2, ls = :dot, lab = "")
savefig("scatterplot_ilc_vs_milca.pdf")
display(p)


text:
#SpiceData: Show functions
#-------------------------------------------------------------------------------

Base.show(io::IO, ::BigEndian) = print(io, "BigEndian")
Base.show(io::IO, ::LittleEndian) = print(io, "LittleEndian")

_showcompact(io::IO, ::SpiceFormat) = print(io, "Format:Unknown")
_showcompact(io::IO, ::Format_9601) = print(io, "SPICE:9601")
_showcompact(io::IO, ::Format_9602) = print(io, "CppSim:9602")
_showcompact(io::IO, ::Format_2001) = print(io, "SPICE:2001")
_showcompact(io::IO, ::Format_2013) = print(io, "SPICE:2013")

Base.show(io::IO, fmt::Format_Unknown) = _showcompact(io, fmt)

function Base.show(io::IO, fmt::SpiceFormat)
	_showcompact(io, fmt)
	print(io, " (x: $(xtype(fmt))[], y: $(ytype(fmt))[])")
end

function _show(io::IO, r::DataReader, compact::Bool = false)
	#Base (compact) information:
	print(io, DataReader, "(")
	print(io, basename(r.filepath))
	print(io, ", nsig=", length(r.signalnames))
	print(io, ", npts=", length(r.sweep))
	print(io, ", ")
	print(io, r.format)
	print(io, ")")
	if compact; return; end

	#Extended information:
	println(io)
	print(io, ">> (", r.endianness, ")")
	print(io, " sweep = '", r.sweepname, "'")
	println(io)
	tags = r.tags
	println(io, ">> ", tags.date, " (", tags.time, ")")
	println(io, ">> ", tags.id)
	println(io, ">> ", tags.comments)
end

Base.show(io::IO, r::DataReader) = _show(io, r)
Base.show(io::IOContext, r::DataReader) = _show(io, r, haskey(io.dict, :compact))

#End


text:
module Solution
  export oddcount

  oddcount(n) = n ÷ 2

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
export minConf_PQN, pqn_options

mutable struct PQN_params
	verbose
	optTol
	progTol
	maxIter
	maxProject
	suffDec
	corrections
	adjustStep
	bbInit
	SPGoptTol
	SPGprogTol
	SPGiters
	SPGtestOpt
end

function pqn_options(;verbose=3,optTol=1f-5,progTol=1f-7,
					 maxIter=20,maxProject=100000, suffDec=1f-4,
					 corrections=10, adjustStep=false,bbInit=false,
					 SPGoptTol=1f-6,SPGprogTol=1f-7,
					 SPGiters=10,SPGtestOpt=false)
	return PQN_params(verbose,optTol,progTol,maxIter,
	                  maxProject,suffDec,corrections,
					  adjustStep,bbInit,SPGoptTol,
					  SPGprogTol,SPGiters,SPGtestOpt)
end

function  minConf_PQN(funObj,x,funProj,options)
# function [x,f] = minConf_PQN(funObj,funProj,x,options)
#
# Function for using a limited-memory projected quasi-Newton to solve problems of the form
#   min funObj(x) s.t. x in C
#
# The projected quasi-Newton sub-problems are solved the spectral projected
# gradient algorithm
#
#   @funObj(x): function to minimize (returns gradient as second argument)
#   @funProj(x): function that returns projection of x onto C
#
#   options:
#       verbose: level of verbosity (0: no output, 1: final, 2: iter (default), 3:
#       debug)
#       optTol: tolerance used to check for optimality (default: 1e-5)
#       progTol: tolerance used to check for progress (default: 1e-9)
#       maxIter: maximum number of calls to funObj (default: 15)
#       maxProject: maximum number of calls to funProj (default: 100000)
#       numDiff: compute derivatives numerically (0: use user-supplied
#       derivatives (default), 1: use finite differences, 2: use complex
#       differentials)
#       suffDec: sufficient decrease parameter in Armijo condition (default: 1e-4)
#       corrections: number of lbfgs corrections to store (default: 10)
#       adjustStep: use quadratic initialization of line search (default: 0)
#       bbInit: initialize sub-problem with Barzilai-Borwein step (default: 1)
#       SPGoptTol: optimality tolerance for SPG direction finding (default: 1e-6)
#       SPGiters: maximum number of iterations for SPG direction finding (default:10)

	fsave = 0;
	nVars = length(x);
	#iter_save = zeros(nVars,1);

	# Output Parameter Settings
	if options.verbose >= 3
	   @printf("Running PQN...\n");
	   @printf("Number of L-BFGS Corrections to store: %d\n",options.corrections);
	   @printf("Spectral initialization of SPG: %d\n",options.bbInit);
	   @printf("Maximum number of SPG iterations: %d\n",options.SPGiters);
	   @printf("SPG optimality tolerance: %.2e\n",options.SPGoptTol);
	   @printf("SPG progress tolerance: %.2e\n",options.SPGprogTol);
	   @printf("PQN optimality tolerance: %.2e\n",options.optTol);
	   @printf("PQN progress tolerance: %.2e\n",options.progTol);
	   @printf("Quadratic initialization of line search: %d\n",options.adjustStep);
	   @printf("Maximum number of function evaluations: %d\n",options.maxIter);
	   @printf("Maximum number of projections: %d\n",options.maxProject);
	end

	# Output Log
	if options.verbose >= 2
	        @printf("%10s %10s %10s %15s %15s %15s\n","Iteration","FunEvals","Projections","Step Length","Function Val","Opt Cond");
	end

	# Make objective function (if using numerical derivatives)
	funEvalMultiplier = 1;

	# Project initial parameter vector
	x = funProj(x);
	projects = 1;

	# Evaluate initial parameters
	f, g = funObj(x);
	funEvals = 1;

	# Check Optimality of Initial Point
	projects = projects+1;
	if maximum(abs.(funProj(x-g)-x)) < options.optTol
	    if options.verbose >= 1
	        @printf("First-Order Optimality Conditions Below optTol at Initial Point\n");
	    end
	    fsave = f;
	    return;
	end

	i = 1;
	while funEvals <= options.maxIter

	    # Compute Step Direction
	    if i == 1

	        p = funProj(x-g);

	        projects = projects+1;
	        S = zeros(Float32, nVars,0);
	        Y = zeros(Float32, nVars,0);
	        Hdiag = 1;
	    else
	        y = g-g_old;
	        s = x-x_old;
	        S, Y, Hdiag = lbfgsUpdate(y,s,options.corrections,options.verbose,S,Y,Hdiag);

	        # Make Compact Representation
	        k = size(Y,2);
	        L = zeros(Float32, k, k);
	        for j = 1:k
	            L[j+1:k,j] = S[:,j+1:k]'*Y[:,j];
	        end
	        N = [S/Hdiag Y];
	        M = [S'*S/Hdiag L;L' -diagm(diag(S'*Y))];
	        HvFunc(v) = lbfgsHvFunc2(v,Hdiag,N,M);

	        if options.bbInit
	            # Use Barzilai-Borwein step to initialize sub-problem
	            alpha = dot(s,s)/dot(s,y);
	            if alpha <= 1e-10 || alpha > 1e10
	                alpha = min(1,1/sum(abs.(g)));
	            end

	            # Solve Sub-problem
	            xSubInit = x-alpha*g;
	            feasibleInit = false;
	        else
	            xSubInit = x;
	            feasibleInit = true;
	        end
	        # Solve Sub-problem
	        p, subProjects = solveSubProblem(x,g,HvFunc,funProj,options.SPGoptTol,options.SPGprogTol,options.SPGiters,options.SPGtestOpt,feasibleInit,xSubInit);
	        projects = projects+subProjects;
	    end
	    d = p-x;
	    g_old = g;
	    x_old = x;

	    # Check that Progress can be made along the direction
	    gtd = dot(g,d);
	    if gtd > -options.progTol && i>(options.corrections/2)
	        if options.verbose >= 1
	            @printf("Directional Derivative below progTol\n");
	        end
	        break;
	    end

	    # Select Initial Guess to step length
	    if i == 1 || options.adjustStep == 0
	       t = 1;
	    else
	        t = min(1,2*(f-f_old)/gtd);
	    end

	    # Bound Step length on first iteration
	    if i == 1
	        t = min(1,1/sum(abs.(d)));
	    end

	    # Evaluate the Objective and Gradient at the Initial Step
	    if t == 1
	        x_new = p;
	    else
	        x_new = x + Float32(t)*d;
	    end
	    f_new, g_new = funObj(x_new);
	    funEvals = funEvals+1;

	    # Backtracking Line Search
	    f_old = f;
	    while f_new > f + options.suffDec*dot(g,(x_new-x)) || ~isLegal(f_new) || ~isLegal(g_new)
	        temp = t;

	        # Backtrack to next trial value
	        if ~isLegal(f_new) || ~isLegal(g_new)
	            if options.verbose == 3
	                @printf("Halving Step Size\n");
	            end
	            t = t/2;
	        else
	            if options.verbose == 3
	                @printf("Cubic Backtracking\n");
	            end
	            t = polyinterp([0 f gtd; t f_new g_new'*d]);
	        end

	        # Adjust if change is too small/large
	        if t < temp*1e-3
	            if options.verbose == 3
	                @printf("Interpolated value too small, Adjusting\n");
	            end
	            t = temp*1e-3;
	        elseif t > temp*0.6
	            if options.verbose == 3
	                @printf("Interpolated value too large, Adjusting\n");
	            end
	            t = temp*0.6;
	        end

	        # Check whether step has become too small
	        if sum(abs.(t*d)) < options.progTol || t == 0
	            if options.verbose == 3
	                @printf("Line Search failed\n");
	            end
	            t = 0;
	            f_new = f;
	            g_new = g;
	            break;
	        end

	        # Evaluate New Point
	        f_prev = f_new;
	        t_prev = temp;
	        x_new = x + Float32(t)*d;
	        f_new, g_new = funObj(x_new);
	        funEvals = funEvals+1;

	        if funEvals > options.maxIter
	            break
	        end

	    end

	    # Take Step
	    x = x_new;
	    f = f_new;
	    fsave = [fsave ; f];
	    # iter_save = [iter_save x];
	    g = g_new;

	    optCond = maximum(abs.(funProj(x-g)-x));
	    projects = projects+1;

	    # Output Log
	    if options.verbose >= 2
	            @printf("%10d %10d %10d %15.5e %15.5e %15.5e\n",i,funEvals*funEvalMultiplier,projects,t,f,optCond);
	    end

	    # Check optimality
	        if optCond < options.optTol && (options.verbose >= 1)
	            @printf("First-Order Optimality Conditions Below optTol\n");
	            break;
	        end

	    if (maximum(abs.(t*d)) < options.progTol) && (funEvals*funEvalMultiplier > options.maxIter/2)
	        if options.verbose >= 1
	            @printf("Step size below progTol\n");
	        end
	        break;
	    end

	    if abs.(f-f_old) < options.progTol
	        if options.verbose >= 1
	            @printf("Function value changing by less than progTol\n");
	        end
	        break;
	    end

	    if funEvals*funEvalMultiplier > options.maxIter
	        if options.verbose >= 1
	            @printf("Function Evaluations exceeds maxIter\n");
	        end
	        break;
	    end

	    if projects > options.maxProject
	        if options.verbose >= 1
	            @printf("Number of projections exceeds maxProject\n");
	        end
	        break;
	    end

	    i = i + 1;
	#    pause
	end

	return x, fsave, funEvals
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
using Lynx
using Test
using Observables
using Colors

var"@asset_str"(::LineNumberNode, ::Module, path) = joinpath(@__DIR__, "..", "assets", path)

include("input.jl")
include("widgets.jl")

@testset "Examples" begin
    include("../examples/basic.jl")
    include("../examples/events.jl")
    include("../examples/images.jl")
    include("../examples/random.jl")
    include("../examples/widgets.jl")
    include("../examples/spritesheet.jl")
end


text:
using CRF

tests = [
    "base.jl",
    "features.jl",
]

for test in tests
    println(">>> $(test)")
    include(test)
end


text:
# Based on http://web.mit.edu/puzzle/www/2015/puzzle/ukacd/

@testset "ukacd" begin
    missing_words = String["hip-flask","infamise","facade","furnaced","ill fame","loofah","odoriferously","focalisation","deceitful","modify","defoliate","definabilities","exemplification","unfeminine","icefield","boastfulness","cuttlefishes","zinckification","folkland","unfeeling","flanges","have a short fuse","pelletifying","fable","efts","state of play","fangled","fitchews","fasciated","soulful","flaunchings","rust-proofing","Mustafa","carfuffle","far-out","goff","quinquefoliate","usufruct","denazified","thoughtfulness","pacificists","calves'-foot jelly","field of view","sweet flag","fandangos","sail-fish","Cadfael","fag","The Return of the Native","gold-foil","feldspathoid","caftan","road-fund licence","purfling","febrifugal","fish eagle","old-fogey","heffalump","confutation","flimping","sell-off","defuzed","battle fatigue","effluences","chaudfroid","flavorous","floury","elf","fusses","dog-faced","felt","law of excluded middle","Schottky effect","faddy","superficializes","tuner amplifiers","off-plan","counter-flory","unselfishness","sound effects","Fawkes","Mafikeng","Laffer curve","felines","foundering","fixed assets","foins","flatlets","dog-fights","lift-off","Fluellen","pouf","suffused","Aesop's Fables","manful","falsifying","kifs","custard coffin","officials","vote of no confidence"]

    fs = Collective.common_features(corpus, [replace(word, r"[^a-zA-Z]", "") for word in missing_words])
    @test any(f.description == "contains at least 1 'f's" for f in fs)
end


text:
module TestAbstractArray
	using Base.Test
	using DataArrays

	unsorted_dv = @data [2, 1, NA]
	sorted_dv = @data [NA, 1, 2]

	@assert isequal(sort(unsorted_dv), sorted_dv)
	@assert isequal(sortperm(unsorted_dv), [3, 2, 1])
	# TODO: Make this work
	# tiedrank(dv)

	@assert first(unsorted_dv) == 2
	@assert isna(last(unsorted_dv))
end


text:
using Pkg, Test
Pkg.activate("../../.")
include("../../src/Toolips.jl")
using Main.Toolips

function test()

    new_webapp("TESTAPP")
end

test()


text:
using PowerSystems

sys = System(0, enable_compression=true)
add_component!(sys, Bus(nothing))
for i in 1:3
    plant = RenewableDispatch(nothing)
    set_name!(plant, "plant$i")
    add_component!(sys, plant)
end

# required_metadata

metadata_file = "data/siip_lookahead_metadata.json"
add_time_series!(sys, metadata_file)
to_json(sys, "lookahead_system/sys.json", force=true)

using Plots
pyplot()
td = get_time_series(
    Deterministic,
    get_component(RenewableDispatch, sys, "plant1"),
    "max_active_power"
)

for window in iterate_windows(td)
    plot!(window, legend=false, color=:black)
end

savefig("lookahead_system/plant1_max_active_power.png")

max_stars_repo_name: NREL-SIIP/reV-PowerSystems

max_issues_repo_name: NREL-SIIP/reV-PowerSystems

max_forks_repo_name: NREL-SIIP/reV-PowerSystems


text:
using Test
using BackwardsLinalg
using Random, LinearAlgebra

@testset "svd grad U" begin
    function loss(A)
        M, N = size(A)
        U, S, V = svd(A)
        psi = U[:,1]
        Random.seed!(2)
        H = randn(ComplexF64, M, M)
        H+=H'
        real(psi'*H*psi)[]
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        a = randn(ComplexF64, M, N)
        @test gradient_check(loss, a)
    end
end

@testset "svd grad V" begin
    function loss_v(A)
        M, N = size(A)
        U, S, V = svd(A)
        Random.seed!(2)
        H = randn(ComplexF64, N, N)
        H+=H'
        psi = V[:,1]
        real(psi'*H*psi)[]
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        a = randn(ComplexF64, M,N)
        @show loss_v(a)
        @test gradient_check(loss_v, a)
    end
end

@testset "svd grad U,V" begin
    function loss_uv(A)
        M, N = size(A)
        U, S, V = svd(A)
        psi = V[1,1]
        psi_l = U[1,1]
        real(conj(psi_l)*psi)[]
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        a = randn(ComplexF64, M,N)
        @show loss_uv(a)
        @test gradient_check(loss_uv, a)
    end
end

@testset "svd grad U,V imag diag" begin
    function loss_uv(A)
        M, N = size(A)
        U, S, V = svd(A)
        psi = V[1,1]
        psi_l = U[1,1]
        real(conj(psi_l)*psi)[]
    end

    A = [-1+1im 2+1im;1-2im 3+0.8im]
    @show loss_uv(A)
    da = [0 0; 1 0im]
    ndiff = (loss_uv(A .+ 1e-4*da) - loss_uv(A .- 1e-4*da)) ./ 2e-4 + im*(loss_uv(A .+ 1e-4im*da) - loss_uv(A .- 1e-4im*da)) ./ 2e-4
    grad = loss_uv'(A)
    @show grad[2,1], ndiff
    @test gradient_check(loss_uv, A)
    @test isapprox(grad[2,1], ndiff, atol=1e-3)
end

@testset "svd grad S" begin
    function loss(A)
        U, S, V = svd(A)
        S |> sum
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        H1 = randn(ComplexF64, M, M)
        H1 += H1'
        a = randn(ComplexF64, M, N)
        @test gradient_check(loss, a)
    end
end

@testset "rsvd" begin
    for shape in [(100, 30), (30, 30), (30, 100)]
        A = randn(ComplexF64, shape...)
        U, S, V = rsvd(A, 30)
        @test U*Diagonal(S)*V' ≈ A
    end

    A = randn(100, 30) * randn(30, 70)
    U, S, V = rsvd(A, 30)
    @test U*Diagonal(S)*V' ≈ A
end

@testset "rsvd grad U" begin
    function loss(A)
        M, N = size(A)
        U, S, V = rsvd(A)
        psi = U[:,1]
        Random.seed!(2)
        H = randn(ComplexF64, M, M)
        H+=H'
        real(psi'*H*psi)[]
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        a = randn(ComplexF64, M, N)
        @test gradient_check(loss, a)
    end
end

@testset "rsvd grad V" begin
    function loss_v(A)
        M, N = size(A)
        U, S, V = rsvd(A)
        Random.seed!(2)
        H = randn(ComplexF64, N, N)
        H+=H'
        psi = V[:,1]
        real(psi'*H*psi)[]
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        a = randn(ComplexF64, M,N)
        @test gradient_check(loss_v, a)
    end
end

@testset "rsvd grad S" begin
    function loss(A)
        U, S, V = rsvd(A)
        S |> sum
    end

    for (M, N) in [(6, 3), (3, 6), (3,3)]
        K = min(M, N)
        H1 = randn(ComplexF64, M, M)
        H1 += H1'
        a = randn(ComplexF64, M, N)
        @test gradient_check(loss, a)
    end
end


text:
################################################################################
# Augmented Random Search ARS
################################################################################
using LinearAlgebra
using Statistics
using Folds

import LinearAlgebra.normalize
import Dojo.GeometryBasics.update

# ARS options: hyper parameters
Base.@kwdef struct HyperParameters{T}
    main_loop_size::Int = 100
    horizon::Int = 200
    step_size::T = 0.02
    n_directions::Int = 16
    b::Int = 16
    noise::T = 0.03
    seed::Int = 1
    env_name::String = "halfcheetah"
end
# assert self.b<=self.n_directions, "b must be <= n_directions"

# observation filter
mutable struct Normalizer{T}
    n::Vector{T}
    mean::Vector{T}
    mean_diff::Vector{T}
    var::Vector{T}
end

function Normalizer(num_inputs::Int)
    n = zeros(num_inputs)
    mean = zeros(num_inputs)
    mean_diff = zeros(num_inputs)
    var = zeros(num_inputs)
    return Normalizer{eltype(n)}(n, mean, mean_diff, var)
end

function observe(normalizer::Normalizer{T}, x::AbstractVector{T}) where T
    normalizer.n .+= 1
    last_mean = deepcopy(normalizer.mean)
    normalizer.mean .+= (x .- normalizer.mean) ./ normalizer.n
    normalizer.mean_diff .+= (x .- last_mean) .* (x .- normalizer.mean)
    normalizer.var .= max.(1e-2, normalizer.mean_diff ./ normalizer.n)
end

function normalize(normalizer, inputs)
    obs_mean = normalizer.mean
    obs_std = sqrt.(normalizer.var)
    return (inputs .- obs_mean) ./ obs_std
end

# linear policy
mutable struct Policy{T}
    hp::HyperParameters{T}
    θ::Matrix{T}
end

function Policy(input_size::Int, output_size::Int, hp::HyperParameters{T}; scale=1.0e-1) where T
    return Policy{T}(hp, scale * randn(output_size, input_size))
end

function evaluate(policy::Policy{T}, input::AbstractVector{T}) where T
    return policy.θ * input
end

function positive_perturbation(policy::Policy{T}, input::AbstractVector{T}, δ::AbstractMatrix{T}) where T
    return (policy.θ + policy.hp.noise .* δ) * input
end

function negative_perturbation(policy::Policy{T}, input::AbstractVector{T}, δ::AbstractMatrix{T}) where T
    return (policy.θ - policy.hp.noise .* δ) * input
end

function sample_δs(policy::Policy{T}) where T
    return [randn(size(policy.θ)) for i = 1:policy.hp.n_directions]
end

function update(policy::Policy, rollouts, σ_r)
    stp = zeros(size(policy.θ))
    for (r_pos, r_neg, d) in rollouts
        stp += (r_pos - r_neg) * d
    end
    policy.θ += policy.hp.step_size * stp ./ (σ_r * policy.hp.b)
    return nothing
end

function sample_policy(policy::Policy{T}) where T
    δ = [randn(size(policy.θ)) for i = 1:policy.hp.n_directions]
    θp = [policy.θ + policy.hp.noise .* δ[i] for i = 1:policy.hp.n_directions]
    θn = [policy.θ - policy.hp.noise .* δ[i] for i = 1:policy.hp.n_directions]
    return [θp..., θn...], δ
end

function eval_sample_policy(θ::Matrix{T}, input::AbstractVector{T}) where T
    return θ * input
end

function rollout_policy(θ::Matrix, env::Environment, normalizer::Normalizer, hp::HyperParameters;
        reset=reset)
    state = reset(env)
    rewards = 0.0
    done = false
    num_plays = 0
    while !done && num_plays < hp.horizon
        observe(normalizer, state)
        state = normalize(normalizer, state)
        action = eval_sample_policy(θ, state)
        state, reward, done, _ = step(env, action)
        reward = max(min(reward, 1), -1)
        rewards += reward
        num_plays += 1
    end
    return rewards
end

function train(env::Environment, policy::Policy{T}, normalizer::Normalizer{T},
        hp::HyperParameters{T}; distributed=false, usefolds=false, foldsexec=Folds.ThreadedEx(;basesize=1)) where T
    println("Training linear policy with Augmented Random Search (ARS)\n ")
    if distributed
        envs = [deepcopy(env) for i = 1:(2 * hp.n_directions)]
        normalizers = [deepcopy(normalizer) for i = 1:(2 * hp.n_directions)]
        hps = [deepcopy(hp) for i = 1:(2 * hp.n_directions)]
        print("  $(nprocs()) processors")
    elseif usefolds
        envs = [deepcopy(env) for i = 1:(2*hp.n_directions)]
        print("  $(Threads.nthreads()) threads with Folds")
    else
        envs = [deepcopy(env) for i = 1:(Threads.nthreads())]
        print(" $(Threads.nthreads()) ")
    end

    # pre-allocate for rewards
    rewards = zeros(2 * hp.n_directions)

    for episode = 1:hp.main_loop_size
        # initialize deltas and rewards
        θs, δs = sample_policy(policy)

        # evaluate policies
        roll_time = @elapsed begin
            if distributed
                rewards .= pmap(rollout_policy, θs, envs, normalizers, hps)
            elseif usefolds
                @assert length(envs) == size(θs, 1) "$(length(envs))"
                Folds.map!(rewards, θs, envs, foldsexec) do θ, env
                    rollout_policy(θ, env, normalizer, hp)
                end
            else
                Threads.@threads for k = 1:(2 * hp.n_directions)
                    rewards[k] = rollout_policy(θs[k], envs[Threads.threadid()], normalizer, hp)
                end
            end
        end
        # reward evaluation
        r_max = [max(rewards[k], rewards[hp.n_directions + k]) for k = 1:hp.n_directions]
        σ_r = std(rewards)
        order = sortperm(r_max, rev = true)[1:hp.b]
        rollouts = [(rewards[k], rewards[hp.n_directions + k], δs[k]) for k = order]

        # policy update
        update(policy, rollouts, σ_r)

        # finish, print:
        println("episode $episode reward_evaluation $(mean(rewards)). Took $(roll_time) seconds")
    end

    return nothing
end

# display learned policy
function display_policy(env::Environment, policy::Policy, normalizer::Normalizer, hp::HyperParameters; rendering = false, reset=reset)
    obs = reset(env)
    traj = [copy(env.state)]
    done = false
    num_plays = 1.
    reward_evaluation = 0
    while !done && num_plays < hp.horizon
        rendering && render(env)
        sleep(env.mechanism.timestep)
        observe(normalizer, obs)
        push!(traj, copy(env.state))
        obs = normalize(normalizer, obs)
        action = evaluate(policy, obs)
        obs, reward, done, _ = step(env, action)
        reward_evaluation += reward
        num_plays += 1
    end
    Dojo.close(env)

    return traj
end

# display learned policy
function display_random_policy(env::Environment, hp::HyperParameters; rendering = false)
    obs = reset(env)
    traj = [copy(env.state)]
    done = false
    num_plays = 1.
    reward_evaluation = 0
    while !done && num_plays < hp.horizon
        rendering && render(env)
        sleep(env.mechanism.timestep)
        push!(traj, copy(env.state))
        action = sample(env.input_space)
        obs, reward, done, _ = step(env, action)
        reward_evaluation += reward
        num_plays += 1
    end
    Dojo.close(env)

    return traj
end


text:
module InsuranceContractsController
using Genie.Renderer.Html

import Base: @kwdef
import TimeZones
using TimeZones
import ToStruct
using ToStruct
import JSON
using JSON
import SearchLight
using SearchLight
import BitemporalPostgres
using BitemporalPostgres
include("InsuranceContracts.jl")
using .InsuranceContracts
export Contract,
    ContractRevision,
    ContractPartnerRef,
    ContractPartnerRefRevision,
    ProductItem,
    ProductItemRevision,
    ProductItemTariffRef,
    ProductItemTariffRefRevision,
    ProductItemPartnerRef,
    ProductItemPartnerRefRevision
include("InsurancePartners.jl")
using .InsurancePartners
export Partner, PartnerRevision
include("InsuranceTariffs.jl")
using .InsuranceTariffs
export Tariff, TariffRevision
export ContractSection, PartnerSection, TariffSection, csection, tsection, psection
export insurancecontracts_view

@kwdef mutable struct PartnerSection
    tsdb_validfrom::TimeZones.ZonedDateTime = now(tz"UTC")
    tsw_validfrom::TimeZones.ZonedDateTime = now(tz"UTC")
    ref_history::SearchLight.DbId = DbId(InfinityKey)
    ref_version::SearchLight.DbId = MaxVersion
    partner_revision::PartnerRevision = PartnerRevision()
end

@kwdef mutable struct TariffSection
    tsdb_validfrom::TimeZones.ZonedDateTime = now(tz"UTC")
    tsw_validfrom::TimeZones.ZonedDateTime = now(tz"UTC")
    ref_history::SearchLight.DbId = DbId(InfinityKey)
    ref_version::SearchLight.DbId = MaxVersion
    tariff_revision::TariffRevision = TariffRevision()
end

@kwdef mutable struct ContractSection
    tsdb_validfrom::TimeZones.ZonedDateTime = now(tz"UTC")
    tsw_validfrom::TimeZones.ZonedDateTime = now(tz"UTC")
    ref_history::SearchLight.DbId = DbId(InfinityKey)
    ref_version::SearchLight.DbId = MaxVersion
    contract_revision::ContractRevision = ContractRevision()
    contract_partnerref_revision::ContractPartnerRefRevision = ContractPartnerRefRevision()
    productitem_revision::ProductItemRevision = ProductItemRevision(position = 0)
    productitem_tariffref_revision::ProductItemTariffRefRevision =
        ProductItemTariffRefRevision()
    productitem_partnerref_revision::ProductItemPartnerRefRevision =
        ProductItemPartnerRefRevision()
    ref_entities::Dict{DbId,Union{PartnerSection,ContractSection,TariffSection}} =
        Dict{DbId,Union{PartnerSection,ContractSection,TariffSection}}()
end

function csection()
    psection = PartnerSection()
    ContractSection(ref_entities = (Dict(DbId(1) => psection)))
end

function insurancecontracts_view()
    html(:insurancecontracts, :insurancecontracts, contracts = all(Contract))
end

function renderhistory(history_id::Int)
    renderhforest(
        mkforest(
            DbId(history_id),
            MaxDate,
            ZonedDateTime(1900, 1, 1, 0, 0, 0, 0, tz"UTC"),
            MaxDate,
        ),
        0,
    )
end

function renderhforest(f, i)
    Html.ul() do
        for_each(f) do node
            renderhnode(node, i + 1)
        end
    end
end

function renderhnode(node, i)
    Html.li() do
        """ Version $(string(node.interval.ref_version.value))
            as of $(node.interval.tsworld_validfrom) created $(node.interval.tsdb_validfrom)
        """ *
        if length(node.shadowed) > 0
            Html.span(class = "caret") * Html.ul(class = "nested") do
                renderhforest(node.shadowed, i + 1)
            end
        else
            ""
        end
    end

end

function historyforest_view(history_id::Int)
    html(
        :insurancecontracts,
        :historyforest,
        hforest = renderhistory(history_id),
        hid = history_id,
        entitytype = "Contract",
    )
end

function get_revision(
    ctype::Type{CT},
    rtype::Type{RT},
    hid::DbId,
    vid::DbId,
) where {CT<:Component,RT<:ComponentRevision}
    find(
        rtype,
        SQLWhereExpression(
            "ref_component=? and ref_valid  @> BIGINT ?",
            find(ctype, SQLWhereExpression("ref_history=?", hid))[1].id,
            vid,
        ),
    )[1]
end

function csection(history_id::Integer, version_id::Integer)::ContractSection
    ContractSection(
        ref_history = DbId(history_id),
        ref_version = DbId(version_id),
        contract_revision = get_revision(
            Contract,
            ContractRevision,
            DbId(history_id),
            DbId(version_id),
        ),
        contract_partnerref_revision = get_revision(
            ContractPartnerRef,
            ContractPartnerRefRevision,
            DbId(history_id),
            DbId(version_id),
        ),
        productitem_revision = get_revision(
            ProductItem,
            ProductItemRevision,
            DbId(history_id),
            DbId(version_id),
        ),
        productitem_tariffref_revision = get_revision(
            ProductItemTariffRef,
            ProductItemTariffRefRevision,
            DbId(history_id),
            DbId(version_id),
        ),
        productitem_partnerref_revision = get_revision(
            ProductItemPartnerRef,
            ProductItemPartnerRefRevision,
            DbId(history_id),
            DbId(version_id),
        ),
        ref_entities = Dict{DbId,Union{PartnerSection,ContractSection,TariffSection}}(),
    )
end

function psection(history_id::Integer, version_id::Integer)::PartnerSection
    PartnerSection(
        partner_revision = get_revision(
            Partner,
            PartnerRevision,
            DbId(history_id),
            DbId(version_id),
        ),
    )
end

function tsection(history_id::Integer, version_id::Integer)::Tariffsection
    TariffSection(
        tariff_revision = get_revision(
            Tariff,
            TariffRevision,
            DbId(history_id),
            DbId(version_id),
        ),
    )
end

function csection_view(history_id::Int,version_id::Int)
    html(
        :insurancecontracts,
        :insurancecontract,
        csect = csection(history_id,version_id)
    )
end

end #module

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
__precompile__()

module ODEInterfaceDiffEq

using Reexport
@reexport using DiffEqBase

using ODEInterface, Compat, DataStructures, FunctionWrappers
using LinearAlgebra

import DiffEqBase: solve

const warnkeywords =
    (:save_idxs, :d_discontinuities, :unstable_check, :tstops,
     :calck, :progress, :dense,:save_start)

function __init__()
    global warnlist = Set(warnkeywords)
end

const KW = Dict{Symbol,Any}

include("algorithms.jl")
include("integrator_types.jl")
include("integrator_utils.jl")
include("solve.jl")

export ODEInterfaceAlgorithm, dopri5, dop853, odex, seulex, radau, radau5, rodas,
       ddeabm, ddebdf

end # module


text:
using Sesam
using ModelingToolkit, DifferentialEquations
using Test
using DistributionFits
using MTKHelpers

@testset "sesam3_protect" begin
    include("test_sesam3_protect.jl")
end;

@testset "sesam3" begin
    #include("test/test_sesam3.jl")
    include("test_sesam3.jl")
end;

@testset "sesam3_revMM" begin
    include("test_sesam3_revMM.jl")
end;

@testset "sesam3_CN" begin
    #include("test/test_sesam3CN.jl")
    include("test_sesam3CN.jl")
end;

@testset "seam3" begin
    include("test_seam3.jl")
end;

@testset "plant_face" begin
    include("test_plant_face.jl")
end;

@testset "plant_face_fluct" begin
    include("test_plant_face_fluct.jl")
end;

@testset "util" begin
    include("test_util.jl")
end;

# @testset "parupdater" begin
#     include("test_parupdater.jl")
# end;


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export maximummatching, isperfectmatching

function backtrack(pred, n)
    path = [n]
    x = n
    while pred[x] != x
        x = pred[x]
        pushfirst!(path, x)
    end
    return path
end

function findaugpath(G::UndirectedGraph, matching)
    head = nodeset(G)  # Initial exposed node set
    for m in matching
        (u, v) = getedge(G, m)
        setdiff!(head, u, v)
    end
    @debug "head: $(head)"
    pred = Dict{Int,Int}()  # F tree
    for v in head
        pred[v] = v
    end
    while !isempty(head)
        v = pop!(head)
        @debug "  v: $(v)"
        unmarked = setdiff(edgeset(G), copy(matching))
        for e in intersect(incidences(G, v), unmarked)
            w = neighbors(G, v)[e]
            @debug "    w: $(w)"
            if !(w in keys(pred))
                e2s = intersect(incidences(G, w), matching)  # TODO use only() in Julia 1.4
                @assert length(e2s) == 1
                e2 = pop!(e2s)
                x = neighbors(G, w)[e2]
                @debug "      x: $(x)"
                pred[w] = v
                pred[x] = w
                union!(head, x)
            else
                p1 = backtrack(pred, v)
                p2 = backtrack(pred, w)
                length(p2) % 2 == 1 || continue
                if p1[1] != p2[1]
                    P = vcat(p1, reverse(p2))
                    @debug "      P: $(P) (extend tree)"
                    return P
                end
                # Contraction
                root = intersect(p1, p2)[end]
                bpath = union([root], setdiff(p1, p2), reverse(setdiff(p2, p1)))
                Gcont = edgecontraction(G, bpath)
                b = Gcont.vnode
                @debug "      Recursion: $(Gcont.collapsed) -> $(b)"
                Mcont = setdiff(matching, edgeset(nodesubgraph(G, bpath)))
                P = findaugpath(Gcont, Mcont)
                @assert length(P) % 2 == 0
                pos = findfirst(x -> x == b, P)
                if pos === nothing
                    @debug "      P: $(P) (no blossom in augpath)"
                    return P
                end
                # Lifting
                if pos == 1 || pos == length(P)
                    nbrpos = pos == 1 ? 2 : length(P) - 1
                    inc = findedgekey(Gcont, P[pos], P[nbrpos])
                    dst = neighbors(G, P[nbrpos])[inc]
                    if root == dst
                        P[pos] = root
                        @debug "      P: $(P) (blossom root is exposed)"
                        return P
                    end
                else
                    inedge = findedgekey(Gcont, P[pos], P[pos-1])
                    outedge = findedgekey(Gcont, P[pos], P[pos+1])
                    inlet = neighbors(G, P[pos-1])[inedge]
                    outlet = neighbors(G, P[pos+1])[outedge]
                    @assert inlet == root || outlet == root
                    dst = inlet == root ? outlet : inlet
                end
                # Build blossom graph from edges
                bgedges = Int[]
                push!(bgedges, findedgekey(G, bpath[1], bpath[end]))
                for i in 1:(length(bpath) - 1)
                    push!(bgedges, findedgekey(G, bpath[i], bpath[i + 1]))
                end
                bgraph = edgesubgraph(G, bgedges)
                # Find path in the blossom
                r = shortestpathnodes(bgraph, root, dst)
                if length(r) % 2 == 0
                    # find even length path
                    redges = shortestpathedges(bgraph, root, dst)
                    cobg = edgesubgraph(bgraph, setdiff(bgedges, redges))
                    r = shortestpathnodes(cobg, root, dst)
                end
                baug = pos % 2 == 0 ? reverse(r) : r
                P = union(P[1:pos-1], baug, P[pos+1:end])
                @debug "      P: $(P) (lifted)"
                return P
            end
        end
    end
    return Int[]
end

function findmaxmatch(G::UndirectedGraph, matching)
    P = findaugpath(G, matching)
    @assert length(P) % 2 == 0
    if !isempty(P)
        Pedges = Set{Int}()
        for i in 1:length(P)-1
            push!(Pedges, findedgekey(G, P[i], P[i+1]))
        end
        augment = symdiff(Pedges, matching)
        @assert length(matching) + 1 == length(augment)
        return findmaxmatch(G, augment)
    else
        return matching
    end
end

"""
    maximummatching(G::UndirectedGraph; method=:Blossom) -> Set{Int}

Compute maximum cardinality matching by Edmonds' blossom algorithm and return the set of matched edges.
"""
maximummatching(G::UndirectedGraph; method=:Blossom) = findmaxmatch(G, Set{Int}())

"""
    isperfectmatching(G::UndirectedGraph) -> Bool
    isperfectmatching(G::UndirectedGraph, matching::Set{Int}) -> Bool

Return if the given graph has a perfect matching.
"""
isperfectmatching(G::UndirectedGraph) = length(maximummatching(G)) * 2 == nodecount(G)
isperfectmatching(G::UndirectedGraph, matching::Set{Int}) = length(matching) * 2 == nodecount(G)


text:
cd(@__DIR__); include("setups/simpletree.jl");
pyplot(dpi = 200)

## (a)
plot(size = (200,500))
gplot!(W, X, width = 1)
scatter_gplot!(X[ib1, :]; c = :pink)
scatter_gplot!(X[ib2, :]; c = :orange)
scatter_gplot!(X[ib3, :]; c = :green)
scatter_gplot!(X[ib4, :]; c = :yellow)
scatter_gplot!(X[ijc, :]; c = :red)
scatter_gplot!(X[ir, :]; c = :grey)
plot!(framestyle = :none)
plot!([0, 0], [0.3, 2], line = :arrow, c = :black, lw = 1)
plt = annotate!(0, 3, text("root", 9))
savefig(plt, "../figs/SimpleTree.png")

## (b)
plot(size = (200, 500))
gplot!(Ws, Xs, width = 1)
scatter_gplot!(Xs[11:11, :]; c = :pink, ms = 4)
scatter_gplot!(Xs[10:10, :]; c = :orange, ms = 4)
scatter_gplot!(Xs[13:13, :]; c = :green, ms = 4)
scatter_gplot!(Xs[12:12, :]; c = :yellow, ms = 4)
scatter_gplot!(Xs[[2,4,6,8], :]; c = :red, ms = 4)
scatter_gplot!(Xs[[1,3,5,7,9], :]; c = :grey, ms = 4)
plt = plot!(framestyle = :none, xlims = [minimum(X[:, 1]), maximum(X[:, 1])],
            ylims = [minimum(X[:, 2]), maximum(X[:, 2])])
savefig(plt, "../figs/SimpleTree_simplified.png")


text:
# dependencies
using LFAToolkit
using LinearAlgebra
using CSV
using DataFrames

# setup
dimension = 1
numbercomponents = 1
mesh = []
if dimension == 1
    mesh = Mesh1D(1.0)
elseif dimension == 2
    mesh = Mesh2D(1.0, 1.0)
elseif dimension == 3
    mesh = Mesh3D(1.0, 1.0, 1.0)
end
convergencefactors = DataFrame()

# test range
for fineP = 1:4
    println("fine_p = ", 2^fineP)
    for coarseP = 0:fineP-1
        println("  coarse_p = ", 2^coarseP)
        # setup
        # -- bases
        coarsep = 2^coarseP
        finep = 2^fineP

        ctofbasis = TensorH1LagrangePProlongationBasis(
            coarsep + 1,
            finep + 1,
            numbercomponents,
            dimension,
        )

        # -- diffusion operators
        finediffusion = GalleryOperator("diffusion", finep + 1, finep + 1, mesh)
        coarsediffusion = GalleryOperator("diffusion", coarsep + 1, finep + 1, mesh)

        # -- smoothers
        identity = IdentityPC(finediffusion)
        jacobi = Jacobi(finediffusion)
        chebyshev = Chebyshev(finediffusion)

        # -- p-multigrid preconditioner
        multigrid = PMultigrid(finediffusion, coarsediffusion, identity, [ctofbasis])

        # compute smoothing factor
        # -- setup
        numbersteps = 8
        θ_min = -π / 2
        θ_max = 3π / 2
        θ_step = 2π / numbersteps
        θ_range = θ_min:θ_step:(θ_max-θ_step)
        if dimension == 3
            numbersteps = 8
            θ_step = 2π / numbersteps
            θ_range = θ_min:θ_step:(θ_max-θ_step)
        end

        # -- compute
        maxeigenvalue = zeros(4)
        θ_maxeigenvalue = -1 * ones(4, dimension)
        # -- 1D --
        if dimension == 1
            for i = 1:numbersteps
                θ = [θ_range[i]]
                if abs(θ[1]) > π / 128
                    M = computesymbols(multigrid, [0], [0, 0], θ)
                    for k = 1:4
                        S = computesymbols(chebyshev, [k], θ)
                        eigenvalues = [abs(val) for val in eigvals(S * M * S)]
                        currentmaxeigenvalue = max(eigenvalues...)
                        if (currentmaxeigenvalue > maxeigenvalue[k])
                            maxeigenvalue[k] = currentmaxeigenvalue
                            θ_maxeigenvalue[k, :] = θ / π
                        end
                    end
                end
            end
            # -- 2D --
        elseif dimension == 2
            for i = 1:numbersteps, j = 1:numbersteps
                θ = [θ_range[i], θ_range[j]]
                if sqrt(abs(θ[1])^2 + abs(θ[2])^2) > π / 128
                    M = computesymbols(multigrid, [0], [0, 0], θ)
                    for k = 1:4
                        S = computesymbols(chebyshev, [k], θ)
                        eigenvalues = [abs(val) for val in eigvals(S * M * S)]
                        currentmaxeigenvalue = max(eigenvalues...)
                        if (currentmaxeigenvalue > maxeigenvalue[k])
                            maxeigenvalue[k] = currentmaxeigenvalue
                            θ_maxeigenvalue[k, :] = θ / π
                        end
                    end
                end
            end
            # -- 3D --
        elseif dimension == 3
            for i = 1:numbersteps, j = 1:numbersteps, k = 1:numbersteps
                θ = [θ_range[i], θ_range[j], θ_range[k]]
                if sqrt(abs(θ[1])^2 + abs(θ[2])^2 + abs(θ[3])^2) > π / 128
                    M = computesymbols(multigrid, [0], [0, 0], θ)
                    for k = 1:4
                        eigenvalueestimates =
                            [0.0, 1.3393, 1.9893, 2.1993, 2.2932, 2.3463, 2.3823, 2.4102]
                        eigenvalue = eigenvalueestimates[finep]
                        S = []
                        if k == 1
                            S = computesymbols(jacobi, [1.0], θ)
                        else
                            S = computesymbols(
                                chebyshev,
                                [k, eigenvalue / 10.0, eigenvalue],
                                θ,
                            )
                        end
                        eigenvalues = [abs(val) for val in eigvals(S * M * S)]
                        currentmaxeigenvalue = max(eigenvalues...)
                        if (currentmaxeigenvalue > maxeigenvalue[k])
                            maxeigenvalue[k] = currentmaxeigenvalue
                            θ_maxeigenvalue[k, :] = θ / π
                        end
                    end
                end
            end
        end
        for k = 1:4
            append!(
                convergencefactors,
                DataFrame(
                    finep = finep,
                    coarsep = coarsep,
                    k = k,
                    θ = θ_maxeigenvalue[k, :],
                    ρ = maxeigenvalue[k],
                ),
            )
        end
    end
end

CSV.write("two-grid-chebyshev.csv", convergencefactors)


text:
# Tests for mutual information function
# (NB These tests will fail with the shrinkage estimator, because the 1D,
# 2D and 3D distributions on the right will shrink towards different
# targets by default, whereas the mutual information formula overcomes
# by estimating the probabilities of the 3D distribution, then projecting)

mi = entropy1 + entropy2 - entropy12

# Test mutual information
@test get_mutual_information(arr1, arr2) ≈ mi
println("Mutual information passed")

# Test mutual information formula
@test apply_mutual_information_formula(entropy1, entropy2, entropy12) ≈ mi
println("Mutual information formula passed.")


text:
"""
    Triweight(μ, σ)
"""
struct Triweight{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    Triweight{T}(µ::T, σ::T) where {T} = new{T}(µ, σ)
end

function Triweight(μ::T, σ::T; check_args::Bool=true) where {T <: Real}
    @check_args Triweight (σ, σ > zero(σ))
    return Triweight{T}(μ, σ)
end

Triweight(μ::Real, σ::Real; check_args::Bool=true) = Triweight(promote(μ, σ)...; check_args=check_args)
Triweight(μ::Integer, σ::Integer; check_args::Bool=true) = Triweight(float(μ), float(σ); check_args=check_args)
Triweight(μ::Real=0.0) = Triweight(μ, one(μ); check_args=false)

@distr_support Triweight d.μ - d.σ d.μ + d.σ

## Conversions

convert(::Type{Triweight{T}}, μ::Real, σ::Real) where {T<:Real} = Triweight(T(μ), T(σ))
Base.convert(::Type{Triweight{T}}, d::Triweight) where {T<:Real} = Triweight{T}(T(d.μ), T(d.σ))
Base.convert(::Type{Triweight{T}}, d::Triweight{T}) where {T<:Real} = d

## Parameters

location(d::Triweight) = d.μ
scale(d::Triweight) = d.σ
params(d::Triweight) = (d.μ, d.σ)
@inline partype(d::Triweight{T}) where {T<:Real} = T

## Properties
mean(d::Triweight) = d.μ
median(d::Triweight) = d.μ
mode(d::Triweight) = d.μ

var(d::Triweight) = d.σ^2 / 9
skewness(d::Triweight{T}) where {T<:Real} = zero(T)
kurtosis(d::Triweight{T}) where {T<:Real} = T(1)/33 - 3

## Functions
function pdf(d::Triweight{T}, x::Real) where T<:Real
    u = abs(x - d.μ)/d.σ
    u >= 1 ? zero(T) : 1.09375*(1 - u*u)^3/d.σ
end
logpdf(d::Triweight, x::Real) = log(pdf(d, x))

function cdf(d::Triweight{T}, x::Real) where T<:Real
    u = (x - d.μ)/d.σ
    u ≤ -1 ? zero(T) : u ≥ 1 ? one(T) : 0.03125*(1 + u)^4*@horner(u,16,-29,20,-5)
end

function ccdf(d::Triweight{T}, x::Real) where T<:Real
    u = (d.μ - x)/d.σ
    u ≤ -1 ? zero(T) : u ≥ 1 ? one(T) : 0.03125*(1 + u)^4*@horner(u,16,-29,20,-5)
end

@quantile_newton Triweight

function mgf(d::Triweight{T}, t::Float64) where T<:Real
    a = d.σ*t
    a2 = a*a
    a == 0 ? one(T) : 105*exp(d.μ*t)*((15/a2+1)*cosh(a)-(15/a2-6)/a*sinh(a))/(a2*a2)
end

function cf(d::Triweight{T}, t::Float64) where T<:Real
    a = d.σ*t
    a2 = a*a
    a == 0 ? complex(one(T)) : 105*cis(d.μ*t)*((1-15/a2)*cos(a)+(15/a2-6)/a*sin(a))/(a2*a2)
end


text:
# Chains for NNs:
#
# (c) A. Dominik, 2021

"""
    abstract type DNN end

Mother type for DNN hierarchy with implementation for a chain of layers.

### Signatures:
```Julia
(m::DNN)(x) = (for l in m.layers; x = l(x); end; x)
(m::DNN)(x,y) = m(x,y)
(m::DNN)(d::Knet.Data) = mean( m(x,y) for (x,y) in d)
(m::DNN)(d::Tuple) = mean( m(x,y) for (x,y) in d)
(m::DNN)(d::NNHelferlein.DataLoader) = mean( m(x,y) for (x,y) in d)
```
"""
abstract type DNN
end
(m::DNN)(x) = (for l in m.layers; x = l(x); end; x)
(m::DNN)(x,y) = m(x,y)
(m::DNN)(d::Knet.Data) = mean( m(x,y) for (x,y) in d)
(m::DNN)(d::Tuple) = mean( m(x,y) for (x,y) in d)
(m::DNN)(d::NNHelferlein.DataLoader) = mean( m(x,y) for (x,y) in d)

"""
    struct Classifier <: DNN

Classifyer with nll loss.

### Signatures:
    (m::Classifier)(x,y) = nll(m(x), y)
"""
struct Classifier <: DNN
    layers
    Classifier(layers...) = new(layers)
end
(m::Classifier)(x,y) = Knet.nll(m(x), y)

"""
    struct Regressor

Regression network with square loss.

### Signatures:
    (m::Regression)(x,y) = sum(abs2.( m(x) - y))
"""
struct Regressor <: DNN
    layers
    Regressor(layers...) = new(layers)
end
(m::Regressor)(x,y) = sum(abs2, m(x) .- y)

"""
    struct Chain

Simple wrapper to chain layers and execute them one after another.
"""
struct Chain <: DNN
    layers
    Chain(layers...) = new(layers)
end

function Base.summary(mdl::DNN; indent=0)
    n = get_n_params(mdl)
    ls = length(mdl.layers)
    s1 = "$(typeof(mdl)) with $ls layers,"
    return print_summary_line(indent, s1, n)
end

function print_network(mdl::DNN; n=0, indent=0)

    top = indent == 0
    if top
        println("Neural network summary:")
        println(summary(mdl))
        println("Details:")
    else
        println(summary(mdl, indent=indent))
    end

    indent += 4
    println(" ")
    for l in mdl.layers
        if l isa DNN
            n = print_network(l, n=n, indent=indent)
            println(" ")
        else
            println(summary(l, indent=indent))
            n += 1
        end
    end

    if top
        println(" ")
        println("Total number of layers: $n")
        println("Total number of parameters: $(get_n_params(mdl))")
    end
    return n
end

function print_summary_line(indent, line, params)

    LIN_LEN = 60
    s1 = " "^indent * line
    len = length(s1)
    gap = " "
    if len < LIN_LEN
        gap = " "^(LIN_LEN-len)
    end

    s2 = @sprintf("%8d params", params)

    return "$s1 $gap $s2"
end


text:
using Gen
using Gen2DAgentMotion
import Distributions
using Plots
import Random

scene = Gen2DAgentMotion.example_apartment_floorplan()
couch = Point(0.15, 0.85)
planner_params = PlannerParams(1000, 1.0, 0.05, 1000, 0.01, 0.02)#PlannerParams(400, 3.0, 200, 0.02)
obs_params = ObsModelParams(0.05, 0.2, 0.02)

@gen function model(T::Int)
    destination ~ uniform_coord()
    start = couch
    path = Point[start]
    (path_rest, failed, tree) = plan_and_optimize_path(scene, start, Point(destination), planner_params)
    append!(path, path_rest)
    points_along_path_and_alignment = ({:observations} ~ path_observation_model(path, obs_params, T))
    points_along_path = points_along_path_and_alignment[1]
    alignment = points_along_path_and_alignment[2]
    return (scene, start, path, failed, points_along_path, alignment)
end

get_path(trace) = get_retval(trace)[3]
get_T(trace) = get_args(trace)[1]
get_points_along_path(trace) = get_retval(trace)[5]
get_alignment(trace) = get_retval(trace)[6]

function generate_synthetic_data()
    T = 10
    trace = simulate(model, (T,))
    (_, _, path, failed) = get_retval(trace)
    observations = Vector{Point}(undef, T)
    for t in 1:T
        observations[t] = Point(trace[:observations => (:x, t)], trace[:observations => (:y, t)])
    end
    return observations
end

function infer(
        observations::Vector{Point}, num_particles::Int, num_samples::Int)
    T = length(observations)
    constraints = choicemap()
    for t in 1:T
        constraints[:observations => (:x, t)] = observations[t].x
        constraints[:observations => (:y, t)] = observations[t].y
    end
    (traces, log_weights, _) = importance_sampling(model, (T,), constraints, num_particles)
    weights = exp.(log_weights)
    idx = Vector{Int}(undef, num_samples)
    Distributions.rand!(Distributions.Categorical(weights / sum(weights)), idx)
    return traces[idx]
end

function draw_scene!()
    for wall in scene.walls
        plot!([wall.a.x, wall.b.x], [wall.a.y, wall.b.y], color="black", label=nothing, aspect_ratio=:equal)
    end
end

function draw_paths(traces)
    for trace in traces
        path = get_path(trace)
        for i in 1:(length(path)-1)
            plot!([path[i].x, path[i+1].x], [path[i].y, path[i+1].y], color="gray", label=nothing)
        end
    end
end

function draw_trace(trace)
    p = scatter([couch.x], [couch.y], color="blue", label=nothing, markerstrokewidth=0)
    scatter!([trace[:destination][1]], [trace[:destination][2]], color="red", label=nothing, markerstrokewidth=0)
    draw_scene!()
    path = get_path(trace)
    for i in 1:(length(path)-1)
        plot!([path[i].x, path[i+1].x], [path[i].y, path[i+1].y], color="gray", label=nothing)
    end
    xs = [trace[:observations => (:x, t)] for t in 1:get_T(trace)]
    ys = [trace[:observations => (:y, t)] for t in 1:get_T(trace)]
    scatter!(xs, ys, color="black", label=nothing, markerstrokewidth=0)
    return p
end

function draw_simulated_traces()
    Random.seed!(2)
    n = 10
    plots = []
    for i in 1:n
        trace = simulate(model, (10,))
        push!(plots, draw_trace(trace))
    end
    plot(plots...)
    savefig("simulated.png")
end

draw_simulated_traces()

function draw_inferences()
    Random.seed!(2)
    plots = []
    for i in 1:4
        observations = generate_synthetic_data()
        traces = infer(observations, 5000, 100)
        dests = [tr[:destination] for tr in traces]
        p = scatter([couch.x], [couch.y], color="blue", label=nothing, markerstrokewidth=0)
        draw_paths(traces)
        xs = [dest[1] for dest in dests]
        ys = [dest[2] for dest in dests]
        scatter!(xs, ys, color="red", label=nothing, markerstrokewidth=0)
        xs = [obs.x for obs in observations]
        ys = [obs.y for obs in observations]
        scatter!(xs, ys, color="black", label=nothing, markerstrokewidth=0)
        draw_scene!()
        push!(plots, p)
    end
    plot(plots...)
    savefig("inferences.png")
end

draw_inferences()


text:
# Domain
x_max = 0.005

# Models
sigma = 0.0004
filter = (x, y) -> exp(-(x^2 + y^2)/2/sigma^2)
filter_dx = (x, y) -> -x/sigma^2*filter(x, y)
filter_dy = (x, y) -> -y/sigma^2*filter(x, y)
K = 2
v_max = 0.05
tau = 1/200

# Noise
sigma_noise = 0.01

# Medium
n_im = 500
dx = 0.0002

# Inverse problem
lp_norm = 1
jump_threshold = 0.2

# Sequence
p = 0.05
initial_position_generator =
() -> x_max*(1.2 * rand() - 0.1)