text:
export svdvals_gkl

"""
Compute the largest singular values of a matrix A using the Golub-Kahan-Lanczos
bidiagonalization method `\cite{Golub1965}` implemented in `svdvals_gkl()`.

This implementation uses full one-sided reorthogonalization as described in
`\cite{Simon2000}`.

Inputs
------

- `A`    : The matrix or matrix-like object whose truncated SVD is desired.
           `A` must support `size()`, `A*v` and `A'*v` methods.

- `nvals`: Number of singular values requested.
           Default: 6

- `v0`   : Initial guess vector.
           Default: a randomized unit vector.

Keyword arguments
-----------------

- `maxiter`: Maximum number of iterations.
             Default: the smaller dimension of A.

- `βth`    : The threshold value of β below which an invariant subspace is
             deemed to be found.
             Default: `0.1*√eps(eltype(A))`

- `σth`    : The threshold value below which a Ritz value estimate of the
             singular value is considered to be converged.
             Default: `0.1*√eps(eltype(A))`

Outputs
-------

- `converged_values`   : The requested singular values

- `convergence_history`: A `Dict{Symbol,Any}` containing the following keys:
    - `:isconverged` : Did the calculation converge?
    - `:iters`       : Number of iterations run
    - `:reorth`      : Number of reorthogonalizations
    - `:mvps`        : Number of matrix-vector products computed
    - `:B`           : The `Bidiagonal` matrix computed during the
                       bidiagonalization process
    - `:β`           : Norm of Lanczos iterate
    - `:ω²`          : The Frobenius norm of the difference between A and the
                       low rank approximation to A computable from the
                       currently computed singular values [Simon2000]
    - `:vals`        : Ritz values computed at each iteration
    - `:valerrs`     : Error bounds on Ritz values at each iteration

Side effects
------------

If an invariant subspace is found which smaller than the either dimension of A,
an informational message is printed and only the singular values corresponding
to this subspace are returned.

References
----------

```bibtex
@article{Golub1965,
    author = {Golub, G. and Kahan, W.},
    doi = {10.1137/0702016},
    journal = {Journal of the Society for Industrial and Applied Mathematics
        Series B Numerical Analysis},
    volume = 2,
    number = 2,
    pages = {205--224},
    title = {Calculating the Singular Values and Pseudo-Inverse of a Matrix},
    year = 1965
}

@article{Simon2000,
    author = {Simon, Horst D. and Zha, Hongyuan},
    doi = {10.1137/S1064827597327309},
    journal = {SIAM Journal on Scientific Computing},
    number = 6,
    pages = {2257--2274},
    title = {Low-Rank Matrix Approximation Using the {L}anczos Bidiagonalization
        Process with Applications},
    volume = 21,
    year = 2000
}
```
"""
function svdvals_gkl(A, nvals::Int=6, v0=randn(size(A,2));
    maxiter::Int=minimum(size(A)),
    βth::Real = 0.1*√eps(eltype(A)),
    σth::Real = 0.1*√eps(eltype(A)),
    krestart::Int = 20,
    ktrunc::Int = nvals
    )

    m, n = size(A)
    T = eltype(A)
    Tσ= eltype(svdfact(A[1,1])[:S])
    @assert length(v0) == n
    p = v0

    αs = T[]
    βs = T[]

    α = Inf
    r = Array(T, m)
    u = Array(T, m) #Same as r
    v = Array(T, n) #Same as p

    V = Array(T, min(m, n), 0) #List of converged right vectors

    ω² = ω²₀ = vecnorm(A)^2

    convergence_history = Dict()
    convergence_history[:isconverged] = false
    convergence_history[:ω²] = [ω²]
    convergence_history[:vals] = Vector{Tσ}[]
    convergence_history[:valerrs] = Vector{Tσ}[]
    convergence_history[:mvps] = 0
    convergence_history[:reorth] = 0

    converged_vectors = Int[]
    local β, k, oldσ, S, Δσ, converged_values, converged_values_errors
    fasterror = true
    for k=1:maxiter
        #Reorthogonalize right vectors [Simon2000]
        #Here we use doubly reorthogonalized classical Gram-Schmidt
        if m ≥ n && k>1
            p -= V*(V'p)
            p -= V*(V'p)
            convergence_history[:reorth] += 2size(V,2)
        end

        β = norm(p)

        v[:] = p/β
        #r = A*v
        #k>1 && (r -= β*u)
        A_mul_B!(r, A, v)
        k>1 && (axpy!(-β, u, r))

        #Reorthogonalize left vectors [Simon2000]
        if m < n
            r -= V*(V'r)
            r -= V*(V'r)
            convergence_history[:reorth] += 2size(V,2)
        end

        α = norm(r)
        u[:] = r/α
        #p = A'u - α*v
        Ac_mul_B!(p, A, u)
        axpy!(-α, v, p)

        β = norm(p)
        push!(αs, α)
        push!(βs, β)
        convergence_history[:mvps] += 2

        #Update Simon-Zha approximation error
        ω² -= α^2
        length(βs) > 1 && (ω² -= βs[end-1]^2)
        push!(convergence_history[:ω²], ω²)

        #Compute error bars on singular values
        d = √(α*β)

        #if fasterror
        #    σ = svdvals(Bidiagonal(αs, βs[1:end-1], false))
        #    Δσ = k==1 ? [d] : svdvals_error!(Δσ, oldσ, σ, d)
        #    fasterror = fasterror && all(isfinite(Δσ)) && all(Δσ.>√eps(Tσ))
        #    oldσ = σ
        #end

        #if !fasterror
            S = svdfact(Bidiagonal(αs, βs[1:end-1], false))
            σ = svdvals(S)
            e1= abs(S[:U][end,:])
            e2= abs(S[:Vt][:,end])
            Δσ= Tσ[d*min(e1[i], e2[i]) for i in eachindex(σ)]
        #end

        #Debug svdvals_error!
        #println("Iteration $k: $(√eps())")
        #for i=1:min(k,10)
        #    println(σ[i], '\t', Δσ[i], '\t', Δσ2[i], '\t')
        #end
        #println()
        #end

        push!(convergence_history[:vals], σ)
        push!(convergence_history[:valerrs], Δσ)

        #Check number of converged values
        converged_values = Array(Tσ, 0)
        converged_values_errors = Array(Tσ, 0)
        converged_vectors = Array(Int, 0)
        for i in eachindex(σ)
            if Δσ[i] ≤ σth
                push!(converged_values, σ[i])
                push!(converged_values_errors, Δσ[i])
                push!(converged_vectors, i)
            end
        end

        #true: do complete reorthogonalization
        if true
            V = [V m≥n?v:u]
        end

        #If invariant subspace has been found, stop
        if β ≤ βth
            if k != minimum(size(A))
                #In exact arithmetic, Lanczos is guaranteed to find an
                #invariant subspace corresponding to the entire range of the
                #matrix. For small test matrices it is entirely possible to
                #attain this limit, so that the Krylov subspace is of rank k.
                #For other matrices this may indicate to continue with a
                #different choice of starting vector.
                info("Invariant subspace of dimension $k found")
            end
            convergence_history[:isconverged] = true
            converged_values = σ
            break
        end

        #If at least n converged values have been found, stop
        if length(converged_values) ≥ nvals
            convergence_history[:isconverged] = true
            break
        end
    end

    convergence_history[:iters] = k
    convergence_history[:B] = Bidiagonal(αs, βs[1:end-1], false)
    convergence_history[:β] = βs
    @assert issorted(converged_values, rev=true)
    converged_values, convergence_history
end

"""
Compute approximate errors in the Rayleigh-Ritz approximation to singular
values in the Golub-Kalan-Lanczos bidiagonalization method implemented in
`svdvals_gkl()`.

Warning
-------
The approximations here are only good to ~`√eps()`, since differences smaller
than that cannot be captured reliably in floating point arithmetic.

Inputs
------

- `Δθ`: Old error vector. Will get appended to.
- `τ` : Old Ritz values, sorted either in ascending or descending order.
- `θ` : Current Ritz values, sorted in the same order as `τ`
- `β` : The norm of the current residual vector

Outputs
-------

- `Δθ`: The new error vector, updated in place.

Implementation notes
--------------------

The notation used corresponds to `\cite[2/e, Corollary 7.9.2]{Parlett1980}`.

A given Ritz value ``θ`` has an associated error bound which is the product of
the norm of the residual and last entry of its associated eigenvector
`\cite[Section 13.2]{Parlett1980}`.

`\cite{Hill1992}` provides an approximation to the last component of each
eigenvector given the current Ritz values and the previous Ritz values, using
the Cauchy interlacing theorem.

This function implements the upper bound of Theorem 3 in `\cite{Hill1992}`,
adapted to the computation of error bounds on Ritz values corresponding to
singular values.

The simple modification is to square the Ritz values to get estimates of
(nonzero) eigenvalues of ``A'A`` and apply Theorem 3 of `\cite{Hill1992}` to
the Ritz values of the eigenvalues, thus deriving the final formulae used here:

```math
\Delta\theta_{j}<\beta\times
\begin{cases}
    \sqrt{\frac{\tau_{1}-\theta_{1}}{\theta_{2}-\theta_{1}}
          \frac{\tau_{1}+\theta_{1}}{\theta_{2}+\theta_{1}}}     & i=1\\
    \sqrt{\frac{\theta_{j}-\tau_{j-1}}{\theta_{j}-\theta_{j-1}}
          \frac{\tau_{j}-\theta_{j}}{\theta_{j+1}-\theta_{j}}
          \frac{\theta_{j}+\tau_{j-1}}{\theta_{j}+\theta_{j-1}}
          \frac{\tau_{j}+\theta_{j}}{\theta_{j+1}+\theta_{j}}}   & i=2,\dots,k-1\\
    \sqrt{\frac{\theta_{k}-\tau_{k-1}}{\theta_{k}-\theta_{k-1}}
          \frac{\theta_{k}+\tau_{k-1}}{\theta_{k}+\theta_{k-1}}} & i=k
\end{cases}
```

We also sneak in an `abs()` before taking the outermost square root to avoid
roundoff error when computing very small error bounds. It also has the
advantage of rendering the formulae agnostic to sort order. So long as `τ` and
`θ` are sorted the same way, the same formulae apply since reversing the sort
order merely reverses the indices in the Cauchy interlacing theorem.

References
----------

```bibtex
@book{Parlett1980,
    author = {Parlett, Beresford N},
    address = {Philadelphia, PA},
    doi = {10.1137/1.9781611971163},
    publisher = {SIAM},
    title = {The symmetric eigenvalue problem},
    year = 1980
}

@article{Hill1992,
	Author = {Hill, R O, Jr. and B N Parlett},
	Doi = {10.1137/0613019},
	Journal = {SIAM Journal on Matrix Analysis and Applications},
	Number = 1,
	Pages = {239-247},
	Title = {Refined Interlacing Properties},
	Volume = 13,
	Year = 1992
}
```
"""
function svdvals_error!{T}(
        Δθ::AbstractVector, τ::AbstractVector{T},
        θ::AbstractVector{T}, β::Real=1.0)

    n = length(θ)
    @assert (issorted(τ, rev=true) && issorted(θ, rev=true)) ||
            (issorted(τ) && issorted(θ))
    @assert length(Δθ)==length(τ)==n-1

    #Check Cauchy interlacing property
    #for i in eachindex(τ)
    #    if !(θ[i]+length(θ)^2*eps(T) ≥ τ[i] ≥ θ[i+1]-length(θ)^2*eps(T))
    #        warn("Interlacing violated $i: $(θ[i]) ≥ $(τ[i]) ≥ $(θ[i+1])")
    #        println(min(abs(τ[i] - θ[i])/eps(T), abs(θ[i+1] - τ[i])/eps(T)))
    #    end
    #end

    Δθ[1] = β*√abs((τ[1]-θ[1])*(τ[1]+θ[1])/
               ((θ[2]-θ[1])*(θ[2]+θ[1])))

    for j=2:n-1
        Δθ[j] = β*√abs(((θ[j]-τ[j-1])/(θ[j]-θ[j-1])) *
                    ((θ[j]+τ[j-1])/(θ[j]+θ[j-1])) *
                    ((τ[j]-θ[j])/(θ[j+1]-θ[j])) *
                    ((τ[j]+θ[j])/(θ[j+1]+θ[j])))
    end
    if n>1
        push!(Δθ, β*√abs((θ[n]-τ[n-1])*(θ[n]+τ[n-1])/
                        ((θ[n]-θ[n-1])*(θ[n]+θ[n-1]))))
    end
    Δθ
end


text:
include("11_two_site_state_preparation.jl")

println("
#######################################################
# Tutorial 12                                         #
#                                                     #
# n-site circuit optimization with MPS                #
#######################################################
")

# Beyond full state, use approximations
n = 30
i = [Index(2, "S=1/2") for j in 1:n]

# Product state |Z+Z+Z+…Z+⟩
Zp = MPS(i, "Z+")

# Product state |Z-Z-Z-…Z-⟩
Zm = MPS(i, "Z-")

inner(Zp, Zp)
inner(Zm, Zp)

# Do some arithmatic with them, like tensors:
# TODO: fix this when normalize PR is merged!
#Cat = (Zp + Zm) / √2
#
# Performs a truncated MPS addition,
# returns an MPS approximate up to some tolerance.
Cat = (Zp + Zm) * (1 / √2)
@show maxlinkdim(Zp)
@show maxlinkdim(Zm)
@show maxlinkdim(Cat)
@show norm(Cat)
@show inner(Cat, Zp)
@show inner(Cat, Zm)

# Apply operators:
j = n ÷ 2
Xⱼ = op("X", i[j])
XⱼZp = apply(Xⱼ, Zp)

XⱼZp = MPS(i, [k == j ? "Z-" : "Z+" for k in 1:n])
@show inner(XⱼZp, apply(Xⱼ, Zp))

function ising(n; h)
  # Symbolic representation
  # of the Hamiltonian ("operator sum")
  H = OpSum()
  for j in 1:(n - 1)
    H -= "Z", j, "Z", j + 1
  end
  for j in 1:n
    H += h, "X", j
  end
  return H
end

h = 0.5
H = MPO(ising(n; h=h), i)

# Performs a truncated MPO*MPS contraction,
# returns an MPS approximate up to some tolerance.
HZp = H * Zp
@show inner(Zp', HZp)
@show inner(Zp', H, Zp)

# Gradient energy minimization
# Already written above! New Hamiltonian.
function E(ψ)
  ψHψ = inner(ψ', H, ψ)
  ψψ = inner(ψ, ψ)
  return ψHψ / ψψ
end

function minimize(f, ∂f, x; nsteps, γ, kwargs...)
  for n in 1:nsteps
    println("n = ", n, ", f_x = ", f(x))
    x = -(x, γ * ∂f(x); kwargs...)
  end
  return x
end

ψ⁰ = MPS(i, "Z+")
∂E(x) = gradient(E, x)[1]
ψ = minimize(E, ∂E, ψ⁰; nsteps=50, γ=0.1, maxdim=10, cutoff=1e-5)

# TODO: change to:
# ψ /= norm(ψ)
ψ *= inv(norm(ψ))

Eᵈᵐʳᵍ, ψᵈᵐʳᵍ = dmrg(H, ψ⁰; nsweeps=10, maxdim=10, cutoff=1e-5)

@show maxlinkdim(ψ⁰)
@show maxlinkdim(ψ)
@show maxlinkdim(ψᵈᵐʳᵍ)
@show E(ψ⁰), norm(∂E(ψ⁰))
@show E(ψ), norm(∂E(ψ))
@show E(ψᵈᵐʳᵍ), norm(∂E(ψᵈᵐʳᵍ))

#
# Circuit optimization
#

# Form the circuit from parameters
Ry_layer(θ, i) = [op("Ry", i[j]; θ=θ[j]) for j in 1:n]
CX_layer(i) = [op("CX", i[j], i[j+1]) for j in 1:2:(n-1)]

function U(θ, i; nlayers)
  n = length(i)
  Uᶿ = Ry_layer(θ[1:n], i)
  for l in 1:(nlayers - 1)
    Uᶿ = [Uᶿ; CX_layer(i)]
    Uᶿ = [Uᶿ; Ry_layer(θ[(1:n) .+ l * n], i)]
  end
  return Uᶿ
end

nlayers = 6
maxdim = 10
cutoff = 1e-5

# Find the circuit `U(θ)` that minimizes:
# E(θ) = ⟨0|U(θ)† H U(θ)|0⟩ = ⟨θ|H|θ⟩
function E(θ)
  # Apply the circuit:
  # |θ⟩ = U(θ)|0⟩
  ψᶿ = apply(U(θ, i; nlayers=nlayers), ψ⁰; maxdim=maxdim, cutoff=cutoff)

  # Compute the expecation value: ⟨θ|H|θ⟩
  # No need to normalize!
  return inner(ψᶿ', H, ψᶿ)
end

println()
println("Circuit optimization")
θ⁰ = zeros(nlayers * n)
∂E(ψ) = gradient(E, ψ)[1]
θ = minimize(E, ∂E, θ⁰; nsteps=20, γ=0.1)

ψᶿ = apply(U(θ, i; nlayers=nlayers), ψ⁰; maxdim=maxdim, cutoff=cutoff)

@show maxlinkdim(ψ⁰)
@show maxlinkdim(ψᶿ)
@show E(θ⁰), norm(∂E(θ⁰))
@show E(θ), norm(∂E(θ))

# TODO: Finish this implementation
# Example of preparing the ground state
# Find `U(θ)` that minimizes
#
# F(θ) = -|⟨ψ|U(θ)|0⟩|²
function F(θ)
  # Apply the circuit:
  # |θ⟩ = U(θ)|0⟩
  ψᶿ = apply(U(θ, i; nlayers=nlayers), ψ⁰; maxdim=maxdim, cutoff=cutoff)

  # -|⟨ψ|θ⟩|²
  return -abs(inner(ψ, ψᶿ))^2
end

println()
println("State preparation")
θ⁰ = zeros(nlayers * n)
∂F(ψ) = gradient(F, ψ)[1]
θ = minimize(F, ∂F, θ⁰; nsteps=20, γ=0.1)

ψᶿ = apply(U(θ, i; nlayers=nlayers), ψ⁰; maxdim=maxdim, cutoff=cutoff)

@show maxlinkdim(ψ⁰)
@show maxlinkdim(ψᶿ)
@show F(θ⁰), norm(∂F(θ⁰))
@show F(θ), norm(∂F(θ))

ext: jl


text:
include("unit/test_core.jl")
include("unit/test_helpers.jl")
include("unit/test_locales.jl")
using .test_core
using .test_helpers
using .test_locales

# TEST CORE
test_core.test_equality()

test_core.test_ordering()

test_core.test_fromisoformat()

test_core.test_today()

test_core.test_year()

test_core.test_month()

test_core.test_day()

test_core.test_datetuple()

test_core.test_isoformat()

test_core.test_dmyformat()

test_core.test_month_length()

test_core.test_month_name()

test_core.test_weekday()

test_core.test_iso_weekday()

test_core.test_day_name()

test_core.test_notation()

test_core.test_to_julian()

test_core.test_to_gregorian()

test_core.test_month_index()

# TEST HELPERS
test_helpers.test_helper_functions()

# TEST LOCALES
test_locales.test_locale_data_structure()

test_locales.test_locale_map()

max_stars_repo_name: ceferisbarov/hijri-converter

max_issues_repo_name: ceferisbarov/hijri-converter


text:
##############################################################################
#
# build.jl
#
# Build for cvortex.jl. Downloads correct precompiled binary.
#
# Copyright 2019-2020 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
##############################################################################

using BinDeps
using Libdl

# version of cvortex package to use
cvortexver="0.3.7"

# binaries url prefix
url="https://github.com/hjabird/cvortex/releases/download/$cvortexver/"

tagfile = "installed_vers"
if !isfile(tagfile) || readchomp(tagfile) != "$cvortexver $WORD_SIZE"
    @info("Installing CVortex $cvortexver library.")
    working_dll = true
    opencl_dll = true
    if Sys.iswindows()
		@info("Downloading $url/cvortex_Win_x64_release.dll")
        run(download_cmd("$url/cvortex_Win_x64_release.dll", "libcvortex.dll"))
    elseif Sys.isapple()
        error("Sorry, no apple release of the binaries yet!")
        global working_dll = false
    elseif Sys.islinux()
		@info("Downloading $url/cvortex_Linux_x64_release.dll")
        run(download_cmd("$url/cvortex_Linux_x64_release.so", "libcvortex.so"))
    end

    try
        dlopen("./libcvortex")
		@info("Successfully linked to library using dlopen.")
        open(tagfile, "w") do f
            println(f, "$cvortexver")
        end
    catch
		@info("Failed to successfully link to library using dlopen. Downloading "*
			"CVortex without OpenMP.")
        global opencl_dll = false
        if Sys.iswindows()
		@info("Downloading $url/cvortex_Win_x64_release_NoOCL.dll")
            run(download_cmd("$url/cvortex_Win_x64_release_NoOCL.dll", "libcvortex.dll"))
        elseif Sys.islinux()
		@info("Downloading $url/cvortex_Linux_x64_release_NoOCL.so")
            run(download_cmd("$url/cvortex_Linux_x64_release_NoOCL.so", "libcvortex.so"))
        end
        try
            dlopen("./libcvortex")
			@info("Successfully linked to library using dlopen. (No OpenCL version)")
            open(tagfile, "w") do f
                println(f, "$cvortexver NoOpenCL")
            end
        catch
			@info("Failed to successfully link to library using dlopen.")
            global working_dll = false
        end
    end

    if !working_dll
        error("Could not install working binary.")
    elseif !opencl_dll
        @warn("The OpenCL accelerated CVortex binary did not work. "*
            "An OpenMP only (no GPU acceleration) version was installed"*
            " instead. If you're expecting GPU acceleration to work, "*
            "check that OpenCL is properly installed on your platform.")
    end
else
    @info("CVortex $cvortexver is already installed.")
end


text:
# %%
using Tmatrix

# inputs
rx, rz = 1e-6, 1.3e-6
n_max = 3
k1_r = 1e7;
k1_i = 0.0;
k2_r = 1.5e7;
k2_i = 1e3;

n_θ_points = 10
n_ϕ_points = 20

rotationally_symmetric = false

θ_array, ϕ_array = meshgrid_θ_ϕ(
    n_θ_points,
    n_ϕ_points;
    min_θ = 1e-16,
    min_ϕ = 1e-16,
    rotationally_symmetric = rotationally_symmetric,
)

r_array, n̂_array = ellipsoid(rx, rz, θ_array)


text:
struct FastSetIndex <: AbstractOptic
    path::NTuple{N, Int} where N
end
FastSetIndex(i::Int) = FastSetIndex((i,))

path_getindex(obj, ::Tuple{}) = obj
path_getindex(obj, path::Tuple) =
    path_getindex(getindex(obj, first(path)), tail(path))
left(f::FastSetIndex, obj) =
    (obj, path_getindex(obj, f.path))

path_setindex(obj, update, (idx,)::Tuple{Int}) =
    setindex(obj, update, idx)
path_setindex(obj, update, path) =
    setindex(obj,
        path_setindex(getindex(obj, first(path)),
        update, tail(path)),
        first(path))
function right(f::FastSetIndex, obj, update)
    path_setindex(obj, update, f.path)
end

(a::FastSetIndex ⨟ b::FastSetIndex) =
    FastSetIndex((a.path..., b.path...))


text:
using Pkg

cd(@__DIR__)
Pkg.activate(".")
Pkg.instantiate()
#Pkg.precompile()

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# misc functions
function epoch(; data, tbl, τ, sfreq, kwargs...)
    epoch(data, tbl, τ, sfreq; kwargs...)
end

function epoch(
    data::Array{T,1},
    tbl::DataFrame,
    τ::Tuple{Number,Number},
    sfreq;
    kwargs...,
) where {T<:Union{Missing,Number}}
    data_r = reshape(data, (1, :))
    epoch(data_r, tbl, τ, sfreq; kwargs...)
end

epoch(data::Matrix, tbl::DataFrame, τ::Vector, sfreq; eventtime) =
    epoch(data, tbl, Tuple(τ...), sfreq; eventtime = eventtime)
epoch(
    data::Matrix,
    tbl::DataFrame,
    τ::Tuple{Number,Number},
    sfreq;
    eventtime::String = "latency",
) = epoch(data, tbl, τ, sfreq; eventtime = Symbol(eventtime))

function epoch(
    data::Array{T,2},
    tbl::DataFrame,
    τ::Tuple{Number,Number},
    sfreq;
    eventtime::Symbol = :latency,
) where {T<:Union{Missing,Number}}
    # data: channels x times

    # partial taken from EEG.jl

    numEpochs = size(tbl, 1)
    τ = round_times(τ, sfreq)

    times = range(τ[1], stop = τ[2], step = 1 ./ sfreq)
    lenEpochs = length(times)
    numChans = size(data, 1)
    epochs = Array{Union{Missing,Float64}}(
        missing,
        Int(numChans),
        Int(lenEpochs),
        Int(numEpochs),
    )

    # User feedback
    @debug "Creating epochs: $numChans x $lenEpochs x $numEpochs"

    for si = 1:size(tbl, 1)
        #eventonset = tbl[si,eventtime] # in samples
        #d_start = eventonset
        d_start = Int(round(tbl[si, eventtime] - sum(times .< 0)))
        d_end = Int(round(tbl[si, eventtime] + lenEpochs - 1 - sum(times .< 0)))
        e_start = 1
        e_end = lenEpochs
        #println("d: $(size(data)),e: $(size(epochs)) | $d_start,$d_end,$e_start,$e_end | $(tbl[si,eventtime])")
        if d_start < 1
            e_start = e_start + (-d_start + 1)
            d_start = 1
        end
        if d_end > size(data, 2)
            e_end = e_end - (d_end - size(data, 2))
            d_end = size(data, 2)
        end
        #println("d: $(size(data)),e: $(size(epochs)) | $d_start,$d_end,$e_start,$e_end | $(tbl[si,eventtime])")
        epochs[:, e_start:e_end, si] = data[:, d_start:d_end]
    end
    return (epochs, times)
end

function round_times(τ, sfreq)
    # function to round τ to sfreq samples. This specifies the epoch length.
    # its a function to be the same for epoch & timeexpanded analyses
    return round.(τ .* sfreq) ./ sfreq

end
function dropMissingEpochs(X, y)
    missingIx = .!any(ismissing.(y), dims = (1, 2))
    print(size(missingIx))
    goodIx = dropdims(missingIx, dims = (1, 2))
    return X[goodIx, :], Array{Float64}(y[:, :, goodIx])
end

"""
$(SIGNATURES)
Flatten a 1D array from of a 2D/3D array. Also drops the empty dimension
"""
function linearize(x::AbstractArray)
    # used in condense_long to generate the long format
    return dropdims(reshape(x, :, 1), dims = 2)
end
function linearize(x::String)
    return x
end

"""
$(SIGNATURES)
Equates the length of data and designmatrix by cutting the shorter one

The reason we need this is because when generating the designmatrix, we do not know how long the data actually are. We only assume that event-latencies are synchronized with the data
"""
function zeropad(X, data::AbstractArray{T,2}) where {T<:Union{Missing,<:Number}}
    @debug("2d zeropad")
    if size(X, 1) > size(data, 2)
        X = X[1:size(data, 2), :]
    else
        data = data[:, 1:size(X, 1)]
    end
    return X, data
end
function zeropad(X, data::AbstractVector{T}) where {T<:Union{Missing,<:Number}}
    @debug("1d zeropad")
    if size(X, 1) > length(data)
        X = X[1:length(data)]
    else
        data = data[1:size(X, 1)]
    end
    return X, data
end

function zeropad(X, data::AbstractArray{T,3}) where {T<:Union{Missing,<:Number}}
    @debug("3d zeropad")
    if size(X, 1) > size(data, 3)
        X = X[1:size(data, 3), :]
    else
        data = data[:, :, 1:size(X, 1)]
    end
    return X, data
end

function clean_data(
    data::AbstractArray{T,2},
    winrej::AbstractArray{<:Number,2},
) where {T<:Union{Float64,Missing}}
    data = Array{Union{Float64,Missing}}(data)
    for row = 1:size(winrej, 1)
        data[:, Int.(winrej[row, 1]:winrej[row, 2])] .= missing
    end
    return data
end


text:
#!/usr/local/bin/julia -q

arg = ARGS[1]
targetPath = joinpath(pwd(), arg)
targetFiles = readdir(targetPath)

for file in targetFiles
  symlink(joinpath(targetPath, file), joinpath(homedir(), file))
end

max_stars_repo_name: moldray/juliaStow

max_issues_repo_name: moldray/juliaStow

max_forks_repo_name: moldray/juliaStow


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.2
#   kernelspec:
#     display_name: Julia 1.5.0
#     language: julia
#     name: julia-1.5
# ---

using LinearAlgebra

include("SOneTo.jl")

"""
    abstract type StaticArrays{S,T,N} <: AbstractArray{T, N} end
    StaticScalar{T} = StaticArray{Tuple{}, T, 0}
    StaticVector{N,T} = StaticArray{Tuple{N}, T, 1}
    StaticMatrix{N, T} = StaticArray{Tuple{N, M}, T, 2}
`StaticArray`s are Julia arrays with fixed, known size.
"""

abstract type StaticArray{S <: Tuple, T, N} <: AbstractArray{T, N} end
const StaticScalar{T} = StaticArray{Tuple{}, T, 0}
const StaticVector{N, T} = StaticArray{Tuple{N}, T, 1}
const StaticMatrix{N, M, T} = StaticArray{Tuple{N, M}, T, 2}
const StaticVecOrMat{T} = Union{StaticVector{<:Any, T}, StaticMatrix{<:Any, <:Any, T}}

# +
# Being a member of StaticMatrixLike, StaticVecOrMatLike, or StaticArrayLike implies that Size(A)
# returns a static Size instance (none of the dimensions are Dynamic). The converse may not be true.
# These are akin to aliases like StridedArray and in similarly bad taste, but the current approach
# in Base necessitates their existence.
const StaticMatrixLike{T} = Union{
    StaticMatrix{<:Any, <:Any, T},
    Transpose{T, <:StaticVecOrMat{T}},
    Adjoint{T, <:StaticVecOrMat{T}},
    Symmetric{T, <:StaticMatrix{<:Any, <:Any, T}},
    Hermitian{T, <:StaticMatrix{<:Any, <:Any, T}},
    Diagonal{T, <:StaticVector{<:Any, T}},
    # We specifically list *Triangular here rather than using
    # AbstractTriangular to avoid ambiguities in size() etc.
    UpperTriangular{T, <:StaticMatrix{<:Any, <:Any, T}},
    LowerTriangular{T, <:StaticMatrix{<:Any, <:Any, T}},
    UnitUpperTriangular{T, <:StaticMatrix{<:Any, <:Any, T}},
    UnitLowerTriangular{T, <:StaticMatrix{<:Any, <:Any, T}}
}
const StaticVecOrMatLike{T} = Union{StaticVector{<:Any, T}, StaticMatrixLike{T}}
const StaticArrayLike{T} = Union{StaticVecOrMatLike{T}, StaticArray{<:Tuple, T}}

const AbstractScalar{T} = AbstractArray{T, 0} # not exported, but useful none-the-less
const StaticArrayNoEltype{S, N, T} = StaticArray{S, T, N}

# +
# utils.jl
# -

# For convenience
TupleN{T,N} = NTuple{N,T}

@inline convert_ntuple(::Type{T}, d::T) where {T} = T # For zero-dimensional arrays
@inline convert_ntuple(::Type{T}, d::NTuple{N,T}) where {N,T} = d
@generated function convert_ntuple(::Type{T}, d::NTuple{N,Any}) where {N,T}
    exprs = ntuple(i->:(convert(T, d[$i])), Val(N))
    return quote
        @_inline_meta
        $(Expr(:tuple, exprs...))
    end
end

# Base gives up on tuples for promote_eltype... (TODO can we improve Base?)
@generated function promote_tuple_eltype(::Union{T,Type{T}}) where T <: Tuple
    t = Union{}
    for i = 1:length(T.parameters)
        tmp = T.parameters[i]
        if tmp <: Vararg
            tmp = tmp.parameters[1]
        end
        t = :(promote_type($t, $tmp))
    end
    return quote
        @_inline_meta
        $t
    end
end

# The ::Tuple variants exist to make sure anything that calls with a tuple
# instead of a Tuple gets through to the constructor, so the user gets a nice error message
Base.@pure tuple_length(T::Type{<:Tuple}) = length(T.parameters)
Base.@pure tuple_length(T::Tuple) = length(T)
Base.@pure tuple_prod(T::Type{<:Tuple}) = length(T.parameters) == 0 ? 1 : *(T.parameters...)
Base.@pure tuple_prod(T::Tuple) = prod(T)
Base.@pure tuple_minimum(T::Type{<:Tuple}) = length(T.parameters) == 0 ? 0 : minimum(tuple(T.parameters...))
Base.@pure tuple_minimum(T::Tuple) = minimum(T)

"""
    size_to_tuple(::Type{S}) where S<:Tuple
Convert a size given by `Tuple{N,M,...}` into a tuple `(N, M, ...)`
"""
Base.@pure function size_to_tuple(::Type{S}) where S<:Tuple
    return tuple(S.parameters...)
end

# +
# Something doesn't match up type wise
function check_array_parameters(Size, T, N, L)
    (!isa(Size, DataType) || (Size.name !== Tuple.name)) && throw(ArgumentError("Static Array parameter Size must be a Tuple type, got $Size"))
    !isa(T, Type) && throw(ArgumentError("Static Array parameter T must be a type, got $T"))
    !isa(N.parameters[1], Int) && throw(ArgumenError("Static Array parameter N must be an integer, got $(N.parameters[1])"))
    !isa(L.parameters[1], Int) && throw(ArgumentError("Static Array parameter L must be an integer, got $(L.parameters[1])"))
    # shouldn't reach here. Anything else should have made it to the function below
    error("Internal error. Please file a bug")
end

@generated function check_array_parameters(::Type{Size}, ::Type{T}, ::Type{Val{N}}, ::Type{Val{L}}) where {Size,T,N,L}
    if !all(x->isa(x, Int), Size.parameters)
        return :(throw(ArgumentError("Static Array parameter Size must be a tuple of Ints (e.g. `SArray{Tuple{3,3}}` or `SMatrix{3,3}`).")))
    end

    if L != tuple_prod(Size) || L < 0 || tuple_minimum(Size) < 0 || tuple_length(Size) != N
        return :(throw(ArgumentError("Size mismatch in Static Array parameters. Got size $Size, dimension $N and length $L.")))
    end

    return nothing
end

# +
"""
   TrivialView
Use to drop static dimensions to override dispatch
"""
struct TrivialView{A, T, N} <: AbstractArray{T, N}
    a::A
end

size(a::TrivialView) = size(a.a)
getindex(a::TrivialView) = getindex(a.a, inds...)
setindex!(a::TrivialView) = setindex!(a.a, inds...)
Base.IndexStyle(::Type{<:TrivialView{A}}) where {A} = IndexStyle{A}

# +
@inline drop_sdims(a::StaticArrayLike) = TrivialView(a)
@inline drop_sdims(a) = a

Base.@propagate_inbounds function invperm(p::StaticVector)
    # in difference to base, this does not check if p is a permutation (every value unique)
     ip = similar(p)
     ip[p] = 1:length(p)
     similar_type(p)(ip)
end


text:
using Weave
start_dir = pwd()
cd(@__DIR__)

weave("../examples/FIR_design.jmd",
     informat="markdown", out_path = "build/examples", doctype = "pandoc")

weave("../examples/FIR_design.jmd",
      informat="markdown", out_path = "build/examples", doctype = "md2html")

cp("../examples/FIR_design.jmd",
    "build/examples/FIR_design.jmd", force = true)

cp("build/examples/FIR_design.md",
      "build/examples/FIR_design.txt", force = true)

weave("../examples/FIR_design_plots.jl", out_path = "build/examples")

cp("../examples/FIR_design_plots.jl",
    "build/examples/FIR_design_plots.jl", force = true)

if !haskey(ENV, "TRAVIS")
    weave("../examples/FIR_design.jmd",
        informat="markdown", out_path = "build/examples", doctype = "md2pdf")
    weave("../examples/FIR_design_plots.jl", doctype = "md2pdf", out_path = "build/examples")
end

cd(start_dir)


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/11/2014

module CASCoordination

export
    AbstractCASCoord,

    initialize,
    update,

    setRecord,
    getRecord,
    getAll,

    GenericCoord

using AbstractCASCoordImpl

using GenericCoordImpl

end


text:
#= matrixlib.jl
=#

function matrixlib(::Type{T}, name::Symbol, args...) where T
  if     name == :cauchy   return matrixlib_cauchy(T, args...)
  elseif name == :fourier  return matrixlib_fourier(T, args...)
  elseif name == :hilb     return matrixlib_hilb(T, args...)
  else                     throw(ArgumentError("name"))
  end
end
matrixlib(name::Symbol, args...) = matrixlib(Float64, name, args...)

function matrixlib_cauchy(::Type{T}, x::AbstractVector, y::AbstractVector) where T
  m = length(x)
  n = length(y)
  A = Array{T}(undef, m, n)
  @inbounds for j = 1:n
    @simd for i = 1:m
      A[i,j] = 1/(x[i] - y[j])
    end
  end
  A
end

function matrixlib_fourier(::Type{T}, x::AbstractVector, y::AbstractVector) where T
  S = eltype(complex(zero(T)))
  m = length(x)
  n = length(y)
  A = Array{S}(undef, m, n)
  @inbounds for j = 1:n, i = 1:m
    A[i,j] = exp(-2im*pi*x[i]*y[j])
  end
  A
end
matrixlib_fourier(::Type{T}, m::Integer, n::Integer) where {T} =
  matrixlib_fourier(T, 0:m-1, (0:n-1)/n)
matrixlib_fourier(::Type{T}, n::Integer) where {T} = matrixlib_fourier(T, n, n)

function matrixlib_hilb(::Type{T}, m::Integer, n::Integer) where T
  m >= 0 || throw(ArgumentError("m"))
  n >= 0 || throw(ArgumentError("n"))
  A = Array{T}(undef, m, n)
  @inbounds for j = 1:n, i = 1:m
    A[i,j] = 1/(i + j - 1)
  end
  A
end
matrixlib_hilb(::Type{T}, n::Integer) where {T} = matrixlib_hilb(T, n, n)


text:
"""
    str_allele_length("myvcf.gz", chrpos)

Build a simple lookup table mapping the allele
to the length of that allele for a given STR or
indel variant in a vcf file.  Vcf must be compressed
and indexed; location should be given as 1:123456.

Dependencies: bcftools
"""
function str_allele_length(vcf, chrpos)

    lup = Dict{Int64, Int64}()

    if Sys.which("bcftools") == nothing
        error("bcftools not found. Please add bcftools to your path.")
    end

    run(`bcftools index -f  $vcf`)

    (chr, pos) = split(chrpos, ':')

    s = read(`bcftools view -H -r $chrpos $vcf`, String)

    if ( length(s) < 1 ) && ( occursin(r"_fam_", vcf) )
        error("the STR specified at $chrpos was not found!")
    end

    t = split(s, '\t')

    lup[0] = length(t[4]) #ref allele
    lup[-1] = 0           #missing data (./. recoded to -1/-1)

    x = split(t[5], ',')

    for (i, j) in enumerate(x)  #all alt alleles
        k = length(j)
        lup[i] = k
    end

    return lup

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
# Lotka-Volterra, np = 1
id = :Sahlodin2011b

pL = [2.95]
pU = [3.05]
Sahlodin2011b_x0(p::Vector{T}) where T = T[1.2*one(T); 1.1*one(T)]
function Sahlodin2011b_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    du[1] = p[1]*u[1]*(1.0 - u[2])
    du[2] = p[1]*u[2]*(u[1] - 1.0)
    return
end
tspan = (0.0, 2.0)

prob = DBB.ODERelaxProb(Sahlodin2011b_f!, tspan, Sahlodin2011b_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0168927411000316"
source = "Sahlodin2011"
desc = "Lotka-Volterra, np = 1"


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
# Writen by Alinson S. Xavier <axavier@anl.gov>

using JuMP, MathOptInterface, DataStructures
import JuMP: value, fix, set_name

# Extend some JuMP functions so that decision variables can be safely replaced by
# (constant) floating point numbers.
function value(x::Float64)
    x
end

function fix(x::Float64, v::Float64; force)
    abs(x - v) < 1e-6 || error("Value mismatch: $x != $v")
end

function set_name(x::Float64, n::String)
    # nop
end

"""
Create a JuMP model using the variables and constraints defined by
the collection of `UCComponent`s in `formulation`.

Parameters
===
  * `isf`: injection shift factors
  * `lodf`: line outage distribution factors
"""
function build_model(;
                     filename::Union{String, Nothing}=nothing,
                     instance::Union{UnitCommitmentInstance, Nothing}=nothing,
                     isf::Union{Array{Float64,2}, Nothing}=nothing,
                     lodf::Union{Array{Float64,2}, Nothing}=nothing,
                     isf_cutoff::Float64=0.005,
                     lodf_cutoff::Float64=0.001,
                     optimizer=nothing,
                     model=nothing,
                     variable_names::Bool=false,
                     formulation::Vector{UCComponent} = UnitCommitment.DefaultFormulation,
                    ) :: UnitCommitmentModel2

    if (filename == nothing) && (instance == nothing)
        error("Either filename or instance must be specified")
    end

    if filename != nothing
        @info "Reading: $(filename)"
        time_read = @elapsed begin
            instance = UnitCommitment.read(filename)
        end
        @info @sprintf("Read problem in %.2f seconds", time_read)
    end

    if length(instance.buses) == 1
        isf = zeros(0, 0)
        lodf = zeros(0, 0)
    else
        if isf == nothing
            @info "Computing injection shift factors..."
            time_isf = @elapsed begin
                isf = UnitCommitment.injection_shift_factors(lines=instance.lines,
                                                             buses=instance.buses)
            end
            @info @sprintf("Computed ISF in %.2f seconds", time_isf)

            @info "Computing line outage factors..."
            time_lodf = @elapsed begin
                lodf = UnitCommitment.line_outage_factors(lines=instance.lines,
                                                          buses=instance.buses,
                                                          isf=isf)
            end
            @info @sprintf("Computed LODF in %.2f seconds", time_lodf)

            @info @sprintf("Applying PTDF and LODF cutoffs (%.5f, %.5f)", isf_cutoff, lodf_cutoff)
            isf[abs.(isf) .< isf_cutoff] .= 0
            lodf[abs.(lodf) .< lodf_cutoff] .= 0
        end
    end

    @info "Building model..."
    time_model = @elapsed begin
        if model == nothing
            if optimizer == nothing
                mip = Model()
            else
                mip = Model(optimizer)
            end
        else
            mip = model
        end
    @info "About to build model"
        model = UnitCommitmentModel2(mip,         # JuMP.Model
                                     DotDict(),   # vars
                                     DotDict(),   # eqs
                                     DotDict(),   # exprs
                                     instance,    # UnitCommitmentInstance
                                     isf,         # injection shift factors
                                     lodf,        # line outage distribution factors
                                     AffExpr(),   # obj
                                     formulation, # formulation
                                    )

        # Prepare variables
        for var in get_required_variables(formulation)
          add_variable(mip, model, instance, UnitCommitment.var_list[var])
        end # prepare variables

        # Prepare constraints
        for constr in get_required_constraints(formulation)
          add_constraint(mip, model, instance, constr)
        end # prepare constraints

        # Prepare expressions (in this case, affine expressions that are later used as part of constraints or objective)
        #   * :startup_cost => contribution to objective of startup costs
        for field in [:startup_cost] #[:net_injection]
            setproperty!(model.exprs, field, OrderedDict())
        end

        # Add components to mip
        for c in formulation
          c.add_component(c, mip, model)
        end

        # Add objective function
        build_obj_function!(model)
    end # end timing of building model
    @info @sprintf("Built model in %.2f seconds", time_model)

    if variable_names
        set_variable_names!(model)
    end

    return model
end # build_model

"""
Add a particular variable to `model.vars`.
"""
function add_variable(mip::JuMP.Model,
                      model::UnitCommitmentModel2,
                      instance::UnitCommitmentInstance,
                      var::UCVariable)
  setproperty!(model.vars, var.name, OrderedDict())
  x = getproperty(model.vars, var.name)
  if !isnothing(var.add_variable)
    var.add_variable(var, x, mip, instance)
    return
  end

  # The following is a bit complex-looking, but the idea is ultimately straightforward
  # We want to loop over the possible index values for var,
  # for every dimension of var (e.g., looping over units and time)
  # The OrderedDict `ind_to_field` maps a UCElement to the corresponding field name within a UnitCommitmentInstance
  # NB: this can be an array of field names, such as [:x, :y], which means we want to access instance.x.y
  # Furthermore, `var` has an array `indices` of UCElement values, describing which index loops over
  # So all we want is to extract the _length_ of the corresponding field of `instance`
  # We create a Tuple so we can feed it to CartesianIndices
  fields = UnitCommitment.ind_to_field(var.indices)
  num_indices = UnitCommitment.num_indices(fields)

  # There is some really complicated logic below that one day needs to be improved
  # (we need to handle nested indices, and this is one way that hopefully works, but it is definitely not intuitive)
  loop_primitive = UnitCommitment.loop_over_indices(UnitCommitment.get_indices_tuple(instance, fields))
  indices = UnitCommitment.get_indices(loop_primitive) # returns an array of tuples? or a unit range maybe.

  for ind in indices
    # For each of the indices, check if the field corresponding to that index has a name
    # Then we will index the variable by that name instead of the integer
    curr_tuple = Tuple(ind)
    new_tuple = ()
    for i in 1:num_indices
      curr_field = UnitCommitment.get_nested_field(instance, fields, i, curr_tuple)
      if :name in propertynames(curr_field)
        new_tuple = (new_tuple..., curr_field.name)
      else
        new_tuple = (new_tuple..., curr_tuple[i])
      end
    end
    name = string(var.name, "[")
    for (i,val) in enumerate(new_tuple)
      name = string(name, val, i < num_indices ? "," : "")
    end
    name = string(name, "]")
    if num_indices == 1
      new_tuple = new_tuple[1]
    end
    x[new_tuple] = @variable(mip,
                             lower_bound=var.lb,
                             upper_bound=var.ub,
                             integer=var.integer,
                             base_name=name)
  end
  ### DEBUG
  #if var.name == :reserve_shortfall
  #  @show var.name, num_indices, loop_primitive, indices, x
  #  #@show JuMP.all_variables(mip)
  #end
  ### DEBUG
end # add_variable

"""
Add constraint to `model.eqs` (set of affine expressions represent left-hand side of constraints).
"""
function add_constraint(mip::JuMP.Model,
                        model::UnitCommitmentModel2,
                        instance::UnitCommitmentInstance,
                        constr::Symbol)
  setproperty!(model.eqs, constr, OrderedDict())
end # add_constraint

"""
Components of the objective include, summed over time:
  * production cost above minimum
  * minimum production cost if generator is on
  * startup cost
  * shutdown cost
  * cost of not meeting shortfall
  * penalty for not meeting or exceeding load (using curtai variable)
  * shutdown cost
"""
function build_obj_function!(model::UnitCommitmentModel2)
    @objective(model.mip, Min, model.obj)
end # build_obj_function

function enforce_transmission(;
                              model::UnitCommitmentModel2,
                              violation::Violation,
                              isf::Array{Float64,2},
                              lodf::Array{Float64,2})::Nothing

    instance, mip, vars = model.instance, model.mip, model.vars
    limit::Float64 = 0.0

    if violation.outage_line == nothing
        limit = violation.monitored_line.normal_flow_limit[violation.time]
        @info @sprintf("    %8.3f MW overflow in %-5s time %3d (pre-contingency)",
                       violation.amount,
                       violation.monitored_line.name,
                       violation.time)
    else
        limit = violation.monitored_line.emergency_flow_limit[violation.time]
        @info @sprintf("    %8.3f MW overflow in %-5s time %3d (outage: line %s)",
                       violation.amount,
                       violation.monitored_line.name,
                       violation.time,
                       violation.outage_line.name)
    end

    fm = violation.monitored_line.name
    t = violation.time
    flow = @variable(mip, base_name="flow[$fm,$t]")

    # |flow| <= limit + overflow
    overflow = vars.overflow[violation.monitored_line.name, violation.time]
    @constraint(mip,  flow <= limit + overflow)
    @constraint(mip, -flow <= limit + overflow)

    if violation.outage_line == nothing
        @constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] *
                                         isf[violation.monitored_line.offset, b.offset]
                                     for b in instance.buses
                                     if b.offset > 0))
    else
        @constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] * (
                                         isf[violation.monitored_line.offset, b.offset] + (
                                             lodf[violation.monitored_line.offset, violation.outage_line.offset] *
                                             isf[violation.outage_line.offset, b.offset]
                                         )
                                     )
                                     for b in instance.buses
                                     if b.offset > 0))
    end
    nothing
end # enforce_transmission

function set_variable_names!(model::UnitCommitmentModel2)
    @info "Setting variable and constraint names..."
    time_varnames = @elapsed begin
        #set_jump_names!(model.vars) # amk: already set
        set_jump_names!(model.eqs)
    end
    @info @sprintf("Set names in %.2f seconds", time_varnames)
end # set_variable_names

function set_jump_names!(dict)
    for name in keys(dict)
        for idx in keys(dict[name])
            idx_str = isa(idx, Tuple) ? join(map(string, idx), ",") : idx
            set_name(dict[name][idx], "$name[$idx_str]")
        end
    end
end # set_jump_names

function get_solution(model::UnitCommitmentModel2)
    instance, T = model.instance, model.instance.time
    function timeseries(vars, collection)
        return OrderedDict(b.name => [round(value(vars[b.name, t]), digits=5) for t in 1:T]
                           for b in collection)
    end
    function production_cost(g)
        return [value(model.vars.is_on[g.name, t]) * g.min_power_cost[t] +
                sum(Float64[value(model.vars.segprod[g.name, k, t]) * g.cost_segments[k].cost[t]
                            for k in 1:length(g.cost_segments)])
                for t in 1:T]
    end
    function production(g)
        return [value(model.vars.is_on[g.name, t]) * g.min_power[t] +
                sum(Float64[value(model.vars.segprod[g.name, k, t])
                            for k in 1:length(g.cost_segments)])
                for t in 1:T]
    end
    function startup_cost(g)
        #S = length(g.startup_categories)
        #return [sum(g.startup_categories[s].cost * value(model.vars.startup[g.name, s, t])
        #            for s in 1:S)
        #        for t in 1:T]
        return [ value.(model.exprs.startup_cost[g.name, t]) for t in 1:T ]
    end
    sol = OrderedDict()
    sol["Production (MW)"] = OrderedDict(g.name => production(g) for g in instance.units)
    sol["Production cost (\$)"] = OrderedDict(g.name => production_cost(g) for g in instance.units)
    sol["Startup cost (\$)"] = OrderedDict(g.name => startup_cost(g) for g in instance.units)
    sol["Is on"] = timeseries(model.vars.is_on, instance.units)
    sol["Switch on"] = timeseries(model.vars.switch_on, instance.units)
    sol["Switch off"] = timeseries(model.vars.switch_off, instance.units)
    sol["Reserve (MW)"] = timeseries(model.vars.reserve, instance.units)
    sol["Net injection (MW)"] = timeseries(model.vars.net_injection, instance.buses)
    sol["Load curtail (MW)"] = timeseries(model.vars.curtail, instance.buses)
    if !isempty(instance.lines)
        sol["Line overflow (MW)"] = timeseries(model.vars.overflow, instance.lines)
    end
    if !isempty(instance.price_sensitive_loads)
        sol["Price-sensitive loads (MW)"] = timeseries(model.vars.loads, instance.price_sensitive_loads)
    end
    return sol
end # get_solution

function fix!(model::UnitCommitmentModel2, solution)::Nothing
    vars, instance, T = model.vars, model.instance, model.instance.time
    for g in instance.units
        for t in 1:T
            is_on = round(solution["Is on"][g.name][t])
            production = round(solution["Production (MW)"][g.name][t], digits=5)
            reserve = round(solution["Reserve (MW)"][g.name][t], digits=5)
            JuMP.fix(vars.is_on[g.name, t], is_on, force=true)
            JuMP.fix(vars.prod_above[g.name, t], production - is_on * g.min_power[t], force=true)
            JuMP.fix(vars.reserve[g.name, t], reserve, force=true)
        end
    end
end # fix!

function set_warm_start!(model::UnitCommitmentModel2, solution)::Nothing
    vars, instance, T = model.vars, model.instance, model.instance.time
    for g in instance.units
        for t in 1:T
            JuMP.set_start_value(vars.is_on[g.name, t], solution["Is on"][g.name][t])
            JuMP.set_start_value(vars.switch_on[g.name, t], solution["Switch on"][g.name][t])
            JuMP.set_start_value(vars.switch_off[g.name, t], solution["Switch off"][g.name][t])
        end
    end
end # set_warm_start

function optimize!(model::UnitCommitmentModel2;
                   time_limit=3600,
                   gap_limit=1e-4,
                   two_phase_gap=true,
                   )::Nothing

    function set_gap(gap)
        try
            JuMP.set_optimizer_attribute(model.mip, "MIPGap", gap)
            @info @sprintf("MIP gap tolerance set to %f", gap)
        catch
            @warn "Could not change MIP gap tolerance"
        end
    end

    instance = model.instance
    initial_time = time()

    large_gap = false
    has_transmission = (length(model.isf) > 0)

    if has_transmission && two_phase_gap
        set_gap(1e-2)
        large_gap = true
    else
        set_gap(gap_limit)
    end

    while true
        time_elapsed = time() - initial_time
        time_remaining = time_limit - time_elapsed
        if time_remaining < 0
            @info "Time limit exceeded"
            break
        end

        @info @sprintf("Setting MILP time limit to %.2f seconds", time_remaining)
        JuMP.set_time_limit_sec(model.mip, time_remaining)

        @info "Solving MILP..."
        JuMP.optimize!(model.mip)

        has_transmission || break

        violations = find_violations(model)
        if isempty(violations)
            @info "No violations found"
            if large_gap
                large_gap = false
                set_gap(gap_limit)
            else
                break
            end
        else
            enforce_transmission(model, violations)
        end
    end

    nothing
end # optimize!

"""
Identify which transmission lines are violated.
See find_violations description from screening.jl.
"""
function find_violations(model::UnitCommitmentModel2)
    instance, vars = model.instance, model.vars
    length(instance.buses) > 1 || return []
    violations = []
    @info "Verifying transmission limits..."
    time_screening = @elapsed begin
        non_slack_buses = [b for b in instance.buses if b.offset > 0]
        net_injections = [value(vars.net_injection[b.name, t])
                          for b in non_slack_buses, t in 1:instance.time]
        overflow = [value(vars.overflow[lm.name, t])
                    for lm in instance.lines, t in 1:instance.time]
        violations = UnitCommitment.find_violations(instance=instance,
                                                    net_injections=net_injections,
                                                    overflow=overflow,
                                                    isf=model.isf,
                                                    lodf=model.lodf)
    end
    @info @sprintf("Verified transmission limits in %.2f seconds", time_screening)
    return violations
end # find_violations

function enforce_transmission(model::UnitCommitmentModel2, violations::Array{Violation, 1})
    for v in violations
        enforce_transmission(model=model,
                             violation=v,
                             isf=model.isf,
                             lodf=model.lodf)
    end
end # enforce_transmission

export UnitCommitmentModel2, build_model, get_solution, optimize!


text:
function augment_states(m::SmetsWouters{T}, TTT::Matrix{T}, RRR::Matrix{T},
                        CCC::Vector{T}) where {T<:AbstractFloat}
    endo = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    exo = m.exogenous_shocks

    n_endo = n_states(m)
    n_exo = n_shocks_exogenous(m)
    @assert (n_endo, n_endo) == size(TTT)
    @assert (n_endo, n_exo) == size(RRR)
    @assert n_endo == length(CCC)

    # Initialize augmented matrices
    n_addl_states = length(endo_addl)
    n_addl_eqs = n_addl_states
    TTT_aug = zeros(n_endo + n_addl_eqs, n_endo + n_addl_states)
    TTT_aug[1:n_endo, 1:n_endo] = TTT
    RRR_aug = [RRR; zeros(n_addl_eqs, n_exo)]
    CCC_aug = [CCC; zeros(n_addl_eqs)]

    ### TTT modifications

    # Track Lags
    TTT_aug[endo_addl[:y_t1], endo[:y_t]] = 1.0
    TTT_aug[endo_addl[:c_t1], endo[:c_t]] = 1.0
    TTT_aug[endo_addl[:i_t1], endo[:i_t]] = 1.0
    TTT_aug[endo_addl[:w_t1], endo[:w_t]] = 1.0
    TTT_aug[endo_addl[:π_t1], endo[:π_t]] = 1.0
    TTT_aug[endo_addl[:L_t1], endo[:L_t]]  = 1.0

    ## We construct state for expected inflation using the fact
    ##
    ##   E_t[s_{t+1}] = CCC + TTT*s_{t} = CCC + TTT * (CCC + TTT*s_{t-1} + RRR)
    ##                = (CCC+TTT*CCC) + (TTT^2)s_{t-1} + TTT*RRR
    ##
    ## So to construct the state for E_t[p_{t+1}], we need to make all these
    ## objects, and just index out the row relevant to pi_t

    T2 = TTT^2
    TR = TTT*RRR
    CTC = CCC+TTT*CCC

    TTT_aug[endo_addl[:Et_π_t],:] = [T2[endo[:π_t],:]; zeros(n_addl_states)]

    RRR_aug[endo_addl[:Et_π_t],:] = TR[endo[:π_t],:]
    CCC_aug[endo_addl[:Et_π_t],:] = CTC[endo[:π_t],:]

    return TTT_aug, RRR_aug, CCC_aug
end


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries
using Test

@testset "Mutating functions" begin
    t1 = Taylor1(6)

    # Dictionaries with calls
    @test length(TaylorSeries._dict_binary_ops) == 5
    # @test length(TaylorSeries._dict_unary_ops) == 22 # why are these tested?

    @test all([haskey(TaylorSeries._dict_binary_ops, op)
        for op in [:+, :-, :*, :/, :^]])
    @test all([haskey(TaylorSeries._dict_binary_calls, op)
        for op in [:+, :-, :*, :/, :^]])

    for kk in keys(TaylorSeries._dict_binary_ops)
        res = TaylorSeries._internalmutfunc_call(
            TaylorSeries._InternalMutFuncs(TaylorSeries._dict_binary_ops[kk]) )
        @test TaylorSeries._dict_binary_calls[kk] == res
    end
    for kk in keys(TaylorSeries._dict_unary_ops)
        res = TaylorSeries._internalmutfunc_call(
            TaylorSeries._InternalMutFuncs(TaylorSeries._dict_unary_ops[kk]) )
        @test TaylorSeries._dict_unary_calls[kk] == res
    end

    # Some examples
    t1 = Taylor1(5)
    t2 = zero(t1)
    TaylorSeries.pow!(t2, t1, 2, 2)
    @test t2[2] == 1.0
    #
    res = zero(t1)
    TaylorSeries.add!(res, t1, t2, 3)
    @test res[3] == 0.0
    TaylorSeries.add!(res, 1, t2, 3)
    @test res[3] == 0.0
    TaylorSeries.add!(res, t2, 3, 0)
    @test res[0] == 3.0
    TaylorSeries.subst!(res, t1, t2, 2)
    @test res[2] == -1.0
    TaylorSeries.subst!(res, t1, 1, 0)
    @test res[0] == -1.0
    @test res[2] == -1.0
    TaylorSeries.subst!(res, 1, t2, 2)
    @test res[2] == -1.0

    res[3] = rand()
    TaylorSeries.mul!(res, t1, t2, 3)
    @test res[3] == 1.0
    TaylorSeries.mul!(res, res, 2, 3)
    @test res[3] == 2.0
    TaylorSeries.mul!(res, 0.5, res, 3)
    @test res[3] == 1.0

    res[0] = rand()
    TaylorSeries.div!(res, t2-1, 1+t1, 0)
    res[1] = rand()
    TaylorSeries.div!(res, t2-1, 1+t1, 1)
    @test res[0] == (t1-1)[0]
    @test res[1] == (t1-1)[1]
    TaylorSeries.div!(res, res, 2, 0)
    @test res[0] == -0.5

    res = zero(t1)
    TaylorSeries.identity!(res, t1, 0)
    @test res[0] == t1[0]
    TaylorSeries.zero!(res, t1, 0)
    TaylorSeries.zero!(res, t1, 1)
    @test res[0] == zero(t1[0])
    @test res[1] == zero(t1[1])
    TaylorSeries.one!(res, t1, 0)
    TaylorSeries.one!(res, t1, 0)
    @test res[0] == one(t1[0])
    @test res[1] == zero(t1[1])

    res = zero(t1)
    TaylorSeries.abs!(res, -1-t2, 2)
    @test res[2] == 1.0
    @test_throws DomainError TaylorSeries.abs!(res, t2, 2)

    res = zero(t1)
    TaylorSeries.abs2!(res, 1-t1, 1)
    @test res[1] == -2.0
    TaylorSeries.abs2!(res, t1, 2)
    @test res[2] == 1.0

end


text:
"""
```
filter(m, data, system, s_0 = [], P_0 = []; cond_type = :none,
    include_presample = true, in_sample = true,
    outputs = [:loglh, :pred, :filt])
```

Computes and returns the filtered values of states for the state-space
system corresponding to the current parameter values of model `m`.

### Inputs

- `m::AbstractModel`: model object
- `data`: `DataFrame` or `nobs` x `hist_periods` `Matrix{S}` of data for
  observables. This should include the conditional period if `cond_type in
  [:semi, :full]`
- `system::System`: `System` object specifying state-space system matrices for
  the model
- `s_0::Vector{S}`: optional `Nz` x 1 initial state vector
- `P_0::Matrix{S}`: optional `Nz` x `Nz` initial state covariance matrix

where `S<:AbstractFloat`.

### Keyword Arguments

- `cond_type::Symbol`: conditional case. See `forecast_all` for documentation of
  all `cond_type` options
- `include_presample::Bool`: indicates whether to include presample periods in
  the returned vector of `Kalman` objects
- `in_sample::Bool`: indicates whether or not to discard out of sample rows in `df_to_matrix` call
- `outputs::Vector{Symbol}`: which Kalman filter outputs to compute and return.
  See `?kalman_filter`

### Outputs

- `kal::Kalman`: see `?Kalman`
"""
function filter{S<:AbstractFloat}(m::AbstractModel, df::DataFrame, system::System{S},
    s_0::Vector{S} = Vector{S}(), P_0::Matrix{S} = Matrix{S}(0, 0);
    cond_type::Symbol = :none, include_presample::Bool = true, in_sample::Bool = true,
    outputs::Vector{Symbol} = [:loglh, :pred, :filt])

    data = df_to_matrix(m, df; cond_type = cond_type, in_sample = in_sample)
    start_date = max(date_presample_start(m), df[1, :date])
    filter(m, data, system, s_0, P_0; start_date = start_date,
           include_presample = include_presample, outputs = outputs)
end

function filter{S<:AbstractFloat}(m::AbstractModel, data::Matrix{S}, system::System,
    s_0::Vector{S} = Vector{S}(0), P_0::Matrix{S} = Matrix{S}(0, 0);
    start_date::Date = date_presample_start(m), include_presample::Bool = true,
    outputs::Vector{Symbol} = [:loglh, :pred, :filt])

    # Partition sample into pre- and post-ZLB regimes
    # Note that the post-ZLB regime may be empty if we do not impose the ZLB
    regime_inds = zlb_regime_indices(m, data, start_date)

    # Get system matrices for each regime
    TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs = zlb_regime_matrices(m, system, start_date)

    # If s_0 and P_0 provided, check that rows and columns corresponding to
    # anticipated shocks are zero in P_0
    if !isempty(s_0) && !isempty(P_0)
        ant_state_inds = setdiff(1:n_states_augmented(m), inds_states_no_ant(m))
        @assert all(x -> x == 0, P_0[:, ant_state_inds])
        @assert all(x -> x == 0, P_0[ant_state_inds, :])
    end

    # Specify number of presample periods if we don't want to include them in
    # the final results
    Nt0 = include_presample ? 0 : n_presample_periods(m)

    # Run Kalman filter, construct Kalman object, and return
    out = kalman_filter(regime_inds, data, TTTs, RRRs, CCCs, QQs,
                        ZZs, DDs, EEs, s_0, P_0; outputs = outputs, Nt0 = Nt0)

    return Kalman(out...)
end


text:
module ImageBase

export
    # two-fold downsampling
    # originally from ImageTransformations.jl
    restrict,

    # finite difference on one-dimension
    # originally from Images.jl
    fdiff,
    fdiff!,

    # basic image statistics, from Images.jl
    minfinite,
    maxfinite,
    maxabsfinite,
    meanfinite

using Reexport

using Base.Cartesian: @nloops
@reexport using ImageCore
using ImageCore.OffsetArrays

include("diff.jl")
include("restrict.jl")
include("statistics.jl")
include("compat.jl")
include("deprecated.jl")

if VERSION >= v"1.4.2" # work around https://github.com/JuliaLang/julia/issues/34121
    include("precompile.jl")
    _precompile_()
end

end


text:
# test constructors
l = 7
vars = (x = 1:3, v = 4:6, z = 1:6, w = 7)
params = (a = 1, b = 1.0)

ps = MVector{l}(zeros(l))
p1 = Particle(ps; variables = vars, parameters = params)
p2 = Particle(eltype(ps), length(ps); variables = vars, parameters = params)

@test p1.state == p2.state == ps
@test p1.views == p2.views
@test p1.params == p2.params == params
@test p1.params === p2.params === params

# test assertions
@test_throws AssertionError Particle(ps; variables = (z = 1:6, w = 7), parameters = (a = 1, z = 2))
@test_throws AssertionError Particle(ps; variables = (z = 1:6, state = 7), parameters = (a = 1, b = 2))
@test_throws AssertionError Particle(ps; variables = (z = 1:6, w = 7), parameters = (a = 1, state = 2))

# test getproperty
s = rand(7)
p = Particle(s; variables = (x = 1:3, v = 4:6, z = 1:6, w = 7), parameters = (a = 1, b = 1.0))

# test state vector
@test p.state === s

# test views
@test p.x == p.state[1:3]
@test p.v == p.state[4:6]
@test p.z == p.state[1:6]
@test p.w == fill(p.state[7])

# test parameters
@test p.a == 1
@test p.b == 1.0


text:
function model(du, u, p::Tuple, t)
    k, x, t0, σ = p
    t0 *= 1000 # input has to be smaller for optimization purposes
    pump_strength = gauss(t, t0, σ)

    # number of reservoirs
    N = length(du)

    # reset deltas to 0
    du .= 0.0

    # main loop
    for i in 1:N, j in 1:N
        i == j && continue
        du[i] -= k[i,j] * u[i]
        du[i] += k[j,i] * u[j]
        x[i,j] ≠ 0 && (du[i] -= u[i] * x[i,j] * pump_strength) # pumping
        x[j,i] ≠ 0 && (du[i] += u[j] * x[j,i] * pump_strength) # pumping don't flip j and i
    end

    nothing
end

function parameters_to_rate_matrix!(M, A, p)
    for i in eachindex(A)
        A[i] ≠ 0 && (M[i] = p[A[i]])
    end
    nothing
end

function parameters_to_rate_matrix!(M, a::Dict{Tuple{T,T}, N}, p) where {T, N}
    for key in keys(a)
        i = a[key] # parameter index
        M[key...] = p[i] # set the value
    end
    nothing
end

time_zero_from_parameters(p) = @view p[end]

num_states(a::Dict) = maximum(maximum.(keys(a)))
num_states(a::Array{T,2}) where T = size(a,1)

function gauss(t, t0, σ)
    exp(-(t - t0)^2 / σ^2)
end

struct State
    target::Int
    source::Array{Int,1}
end
State() = State(0, Int[])

function population_to_attenuation!(a, u, states)
    for i in eachindex(states)
        if states[i].target == 0
            a[i] = 1.0
            continue
        end
        a[i] = 0.0
        for source_state in states[i].source
            a[i] += u[source_state]
        end
        u_target = u[states[i].target]
        a[i] -= u_target
    end
end

function attenuated_sfspectrum!(y, x, a, A, ω, Γ)
    a .*= A
    sfspectrum!(y, x, a, ω, Γ)
end

function solution_to_attenuated_spectra!(Y, a, wn, dl, spec_params, sol, states)
    Threads.@threads for i in eachindex(dl)
        _a = @view a[:,i]
        _u = @view sol[:,i]
        _y = @view Y[:,i]
        population_to_attenuation!(_a, _u, states)
        attenuated_sfspectrum!(_y, wn, _a, spec_params[1], spec_params[2], spec_params[3])
    end
end

function fit(data, spec_params, u0, p_initial, assignments, state_array)
    function loss(p)
            SFGReservoirs.parameters_to_rate_matrix!(k, ak, p)
            SFGReservoirs.parameters_to_rate_matrix!(x, ax, p)

            # these should be at fixed position
            t0 = p[end-1]
            σ  = p[end]

            tspan = extrema(delaytimes)
            params = (k,x,t0,σ)
            prob = ODEProblem(model, u0, tspan, params)
            sol = solve(prob; saveat=delaytimes)

            SFGReservoirs.solution_to_attenuated_spectra!(
                Y, a, wavenumbers, delaytimes, spec_params, sol, state_array
            )

            @. R = Y - Y′
            sum(abs2.(R))
    end

    # unpack data
    delaytimes, wavenumbers, Y′ = data
    ak, ax = assignments

    n_states = length(u0)
    # preallocations
    k  = @MMatrix zeros(n_states, n_states)
    x  = @MMatrix zeros(n_states, n_states)
    # caclulated values
    Y = zeros(ComplexF64, (length(wavenumbers), length(delaytimes)))
    # residuals
    R = zeros(ComplexF64, (length(wavenumbers), length(delaytimes)))
    # attenuation factors
    a = zeros(length(state_array), length(delaytimes))

    optimize(loss, p_initial, NelderMead(),
            Optim.Options(
                time_limit=120,
            );
    )
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test

# Tests for domsort

## Test that domsort doesn't mangle single-argument phis (#29262)
let m = Meta.@lower 1 + 1
    @assert Meta.isexpr(m, :thunk)
    src = m.args[1]::Core.CodeInfo
    src.code = Any[
        # block 1
        Expr(:call, :opaque),
        Expr(:gotoifnot, Core.SSAValue(1), 10),
        # block 2
        Core.PhiNode(Any[8], Any[Core.SSAValue(7)]), # <- This phi must not get replaced by %7
        Core.PhiNode(Any[2, 8], Any[true, false]),
        Expr(:gotoifnot, Core.SSAValue(1), 7),
        # block 3
        Expr(:call, :+, Core.SSAValue(3), 1),
        # block 4
        Core.PhiNode(Any[5, 6], Any[0, Core.SSAValue(6)]),
        Expr(:call, >, Core.SSAValue(7), 10),
        Expr(:gotoifnot, Core.SSAValue(8), 3),
        # block 5
        Core.PhiNode(Any[2, 8], Any[0, Core.SSAValue(7)]),
        Expr(:return, Core.SSAValue(10)),
    ]
    nstmts = length(src.code)
    src.ssavaluetypes = nstmts
    src.codelocs = fill(Int32(1), nstmts)
    src.ssaflags = fill(Int32(0), nstmts)
    ir = Core.Compiler.inflate_ir(src)
    Core.Compiler.verify_ir(ir)
    domtree = Core.Compiler.construct_domtree(ir.cfg)
    ir = Core.Compiler.domsort_ssa!(ir, domtree)
    Core.Compiler.verify_ir(ir)
    @test isa(ir.stmts[3], Core.PhiNode) && length(ir.stmts[3].edges) == 1
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using VDPTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using PmapProgressMeter
using MCTS
using Plots

@show N = 1000

@everywhere begin
    using POMDPs
    using POMDPToolbox
    using MCTS
    using VDPTag

    mdp = VDPTagMDP()
    dmdp = DiscreteVDPTagMDP()

    n = 10000
    seed = 42

    planners = Dict{String, Policy}(
        "continuous_dpw" => begin
            solver = DPWSolver(depth=20,
                               exploration_constant=40.0,
                               rng=MersenneTwister(seed),
                               estimate_value=RolloutEstimator(ToNextML(mdp)),
                               n_iterations=n,
                               k_action=8.0,
                               alpha_action=1/20,
                               k_state=4.0,
                               alpha_state=1/20)
            solve(solver, mdp)
        end,

        "discrete_dpw" => begin
            ro = translate_policy(ToNextML(mdp), mdp, dmdp, dmdp)
            solver = DPWSolver(depth=20,
                               exploration_constant=40.0,
                               rng=MersenneTwister(seed),
                               estimate_value=RolloutEstimator(ro),
                               n_iterations=n,
                               k_action=8.0,
                               alpha_action=1/20,
                               k_state=4.0,
                               alpha_state=1/20)
            translate_policy(solve(solver, dmdp), dmdp, mdp, dmdp)
        end,

        "discrete_mcts" => begin
            ro = translate_policy(ToNextML(mdp), mdp, dmdp, dmdp)
            solver = MCTSSolver(depth=20,
                                exploration_constant=40.0,
                                rng=MersenneTwister(seed),
                                estimate_value=RolloutEstimator(ro),
                                n_iterations=n)
            translate_policy(solve(solver, dmdp), dmdp, mdp, dmdp)
        end,

        "discrete_random" => translate_policy(RandomPolicy(dmdp, rng=MersenneTwister(seed)), dmdp, mdp, dmdp),

        "discrete_heur" => begin
            dheur = translate_policy(ToNextML(mdp), mdp, dmdp, dmdp)
            translate_policy(dheur, dmdp, mdp, dmdp)
        end,

        "continuous_heur" => ToNextML(mdp)
    )
end

for (k,p) in planners
    prog = Progress(N, desc="Simulating...")
    rewards = pmap(prog, 1:N) do i
        hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
        hist = simulate(hr, mdp, p)
        discounted_reward(hist)
    end
    @show k
    @show mean(rewards)
end


text:
using Test, WrappedR

X = randn(100,5)
y = 7.0 .+ X[:,1] .+ 2.0 .* X[:,2] .* X[:,3] + 3*randn(100)

@test size(X) == (100, 5)
@test size(y) == (100,)

m  = fitsprintr(X, y)
pr = predictsprintr(m, X)

@test m[:n] == 100
@test m[:p] == 5
#@test m[:a0] #changes. setseed
#@test m[:compact] #changes. setseed
#@test m[:fit]


text:
"""
A Julia package for process neuroimaing data.
"""
module Neuroimaging

using Logging,  # For user feedback
    Unitful,
    DataFrames,
    Distances,
    ProgressMeter,
    BDF,
    CSV,
    DelimitedFiles,
    DSP,
    Distributions,
    Plots,
    Images,
    BDF,
    MAT,
    Printf,
    Statistics,
    FFTW,
    DataDeps

export new_processing_key,
    find_keys_containing,
    fileparts,
    add_dataframe_static_rows,
    _find_closest_number_idx,
    # File type reading and writing
    import_biosemi,
    channelNames_biosemi_1020,
    create_channel,
    create_events,
    read_avr,
    read_bsa,
    read_dat,
    read_sfp,
    read_elp,
    write_dat,
    prepare_dat,
    write_avr,
    read_evt,
    read_rba_mat,
    # Pre-processing
    epoch_rejection,
    channel_rejection,
    highpass_filter,
    lowpass_filter,
    bandpass_filter,
    compensate_for_filter,
    remove_template,
    rereference,
    clean_triggers,
    validate_triggers,
    extra_triggers,
    downsample,
    # Reshaping of data
    extract_epochs,
    average_epochs,
    create_sweeps,
    # Statistics
    ftest,
    gfp,
    # Type - Volume Image
    VolumeImage,
    read_VolumeImage,
    plot,
    normalise,
    find_dipoles,
    std,
    mean,
    isequal,
    ==,
    # Type - SSR
    SSR,
    samplingrate,
    modulationrate,
    channelnames,
    read_SSR,
    trim_channel,
    add_triggers,
    remove_channel!,
    keep_channel!,
    add_channel,
    assr_frequency,
    save_results,
    trigger_channel,
    write_SSR,
    merge_channels,
    hcat,
    # Source analysis
    Sensor,
    Electrode,
    label,
    labels,
    x,
    y,
    z,
    show,
    match_sensors,
    EEG_64_10_20,
    EEG_Vanvooren_2014,
    EEG_Vanvooren_2014_Left,
    EEG_Vanvooren_2014_Right,
    Coordinate,
    SPM,
    BrainVision,
    Talairach,
    UnknownCoordinate,
    convert,
    conv_bv2tal,
    conv_spm_mni2tal,
    Dipole,
    best_dipole,
    Leadfield,
    match_leadfield,
    find_location,
    project,
    euclidean,
    # Plotting
    plot,
    oplot,
    plot_dat,
    plot_spectrum,
    plot_timeseries,
    plot_single_channel_timeseries,
    plot_multi_channel_timeseries,
    oplot_dipoles,
    SSR_spectrogram,
    plot_filter_response,
    plot_ftest

# Helper functions
include("miscellaneous/helper.jl")
include("datasets/datasets.jl")

# Pre-processing
include("preprocessing/data_rejection.jl")
include("preprocessing/filtering.jl")
include("preprocessing/reference.jl")
include("preprocessing/triggers.jl")

# Reshaping of data
include("reshaping/epochs.jl")
include("reshaping/sweeps.jl")

# Statistics
include("statistics/ftest.jl")
include("statistics/gfp.jl")

# Type - Coordinates
include("types/Coordinates/Coordinates.jl")

# Type - Sensors
include("types/Sensors/Sensors.jl")

# Type - SSR
include("types/SSR/SSR.jl")
include("types/SSR/Preprocessing.jl")
include("types/SSR/ReadWrite.jl")
include("types/SSR/Reshaping.jl")
include("types/SSR/Statistics.jl")

# Type - Dipole
include("types/Dipole/Dipole.jl")
include("types/Dipole/Operations.jl")

include("types/Sensors/Operations.jl")
include("types/Sensors/Sets.jl")
include("types/Coordinates/Operations.jl")

# Source analysis
include("source_analysis/dipoles.jl")

# Type - Leadfield
include("types/Leadfield/Leadfield.jl")
include("types/Leadfield/Operations.jl")

# Type - Volume Image
include("types/VolumeImage/VolumeImage.jl")
include("types/VolumeImage/ReadWrite.jl")
include("types/VolumeImage/Plotting.jl")
include("types/VolumeImage/Dipoles.jl")
include("types/VolumeImage/Operations.jl")

# Plotting functions
include("types/SSR/plotting.jl")
include("plotting/plots.jl")
include("types/Dipole/Plotting.jl")

# File type reading and writing
include("read_write/avr.jl")
include("read_write/bdf.jl")
include("read_write/bsa.jl")
include("read_write/dat.jl")
include("read_write/elp.jl")
include("read_write/evt.jl")
include("read_write/rba.jl")
include("read_write/sfp.jl")

end # module


text:
@testset begin
    config = _default_config()

    project_contents = """
        [deps]
        Example = "7876af07-990d-54b4-ab0e-23690620f79a"
    """

    project_filename = _contents_to_filename(; contents = project_contents)

    projs_and_mnfsts = [
        Project(FileContents(project_contents)),
        Project(FileName(project_filename)),
        Project(FileName(project_filename); script = PkgServerCacheWarmer._DEFAULT_SCRIPT),
    ]

    @test warmup_cache(config, projs_and_mnfsts) isa Nothing
end


text:
@testset "VarInt" begin

    io = IOBuffer()
    write(io, UInt8[0x05])
    seekstart(io)
    @test BTCParser.read_varint(io) == 0x0000_0000_0000_0005

    io = IOBuffer()
    write(io, UInt8[0xfd, 0x00, 0x05])
    seekstart(io)
    @test BTCParser.read_varint(io) == 0x0000_0000_0000_0500

    io = IOBuffer()
    write(io, UInt8[0xfe, 0x00, 0x00, 0x00, 0x05])
    seekstart(io)
    @test BTCParser.read_varint(io) == 0x0000_0000_0500_0000

    io = IOBuffer()
    write(io, UInt8[0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x05])
    seekstart(io)
    @test BTCParser.read_varint(io) == 0x0500_0000_0000_0000

end


text:
module Time

using JavaCall: JString

using ..OrekitWrapper.TimeWrapper

import Dates

export AbsoluteDate, get_tai

TimeWrapper.get_tai() = get_tai(TimeScalesFactory)

Base.show(io::IO, tai::TAIScale) = print(io, to_string(tai))

function TimeWrapper.AbsoluteDate(year::Int, month::Int, day::Int, hour::Int, minute::Int, second::Float64, scale)
    return AbsoluteDate(Int32(year), Int32(month), Int32(day),
        Int32(hour), Int32(minute), second, convert(TimeScale, scale))
end

TimeWrapper.AbsoluteDate(str::String, scale) = AbsoluteDate(JString(str), convert(TimeScale, scale))

function TimeWrapper.AbsoluteDate(dt::Dates.DateTime, scale)
    y = Dates.year(dt)
    m = Dates.month(dt)
    d = Dates.day(dt)
    h = Dates.hour(dt)
    mn = Dates.minute(dt)
    s = Dates.millisecond(dt) * 1e-3 + Dates.second(dt)
    return AbsoluteDate(y, m, d, h, mn, s, scale)
end

function Dates.DateTime(dt::AbsoluteDate, scale)
    comp = get_components(dt, convert(TimeScale, scale))
    date = Dates.Date(get_date(comp))
    time = Dates.Time(get_time(comp))
    return Dates.DateTime(date, time)
end

function Dates.Date(date::DateComponents)
    y = get_year(date)
    m = get_month(date)
    d = get_day(date)
    return Dates.Date(y, m, d)
end

function Dates.Time(time::TimeComponents)
    h = get_hour(time)
    m = get_minute(time)
    s = get_second(time)
    sec = round(Int, s)
    ms = (s - sec) * 1e3
    milli = round(Int, ms)
    μs = (ms - milli) * 1e3
    micro = round(Int, μs)
    ns = (μs - micro) * 1e3
    nano = round(Int, ns)
    return Dates.Time(h, m, sec, milli, micro, nano)
end

end


text:
@testset "TinayHanabiEnv" begin

    env = TinyHanabiEnv()

    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
function newton_LambertW(x ; rel_tol=0.0001)
    w0 = 0
    if x<10
        w1 = 1
    else
        w1 = max(1,log(x)-log(log(x))+ log(log(x))/log(x))
    end

    while abs((w0-w1)/min(abs(w0),abs(w1))) > rel_tol
        w0 = w1
        w1 = w1 - (w1 - x * exp(-w1))/(1+w1)
    end
    return w1
end

function slicing(x0_eff, σ_eff, bl_eff, n, former_λ ; n_it = 5, bounds = 5)
    a = x0_eff + n * σ_eff^2
    t_eff = bl_eff * σ_eff^2

    if (exp(a) * t_eff) < 10000
        best_λi = a - newton_LambertW(t_eff * exp(a), rel_tol=0.000001)
    else
        lt = log(t_eff) + a
        best_λi = a - lt + log(lt) - (log(lt)/lt)
    end
    max_λi = best_λi + bounds * σ_eff
    min_λi = best_λi - bounds * σ_eff
    extant_λi = max_λi - min_λi

    max_f = n * best_λi - bl_eff * exp(best_λi) - (best_λi - x0_eff)^2/(2 * σ_eff^2)
    j = 1.
    function fx(x ; bl = bl_eff, ne = n, x0 = x0_eff, s = 2*σ_eff^2)
         exp(x * ne - bl * exp(x) - (x-x0)^2/s - max_f)
     end

    λ = best_λi
    fλ = fx(λ)

    if -1e10 < max_f < 1e10
        for j in 1:n_it
            min_eff = min_λi
            max_eff = max_λi

            u = fλ * rand()
            reject = true
            if u < 1
                n_iter = 0
                while reject
                    n_iter += 1
                    λ = rand() * (max_eff - min_eff) + min_eff
                    fλ = fx(λ)
                    reject = (fλ < u)
                    if reject
                        if λ < best_λi
                            min_eff = λ
                        else
                            max_eff = λ
                        end
                    end

                end
            else
                λ = best_λi
            end
        end
    else
        λ = log(former_λ)
        #new_λ = former_λ
    end

    new_λ = exp(λ)
    if isnan(new_λ) || new_λ <= 0.
        new_λ = former_λ#rates[i + 1]
    end

    return new_λ
end

function slicing0(S, n_events, former_λ0 ; n_it = 5, bounds = 5)

    best_λ0 = log(n_events/S)
    max_λ0 = best_λ0+20
    min_λ0 = best_λ0-20
    extant_λ0 = max_λ0 - min_λ0

    max_f = n_events * (best_λ0 - 1)

    function fx(x ; ne = n_events, s = S)
         exp(x * ne - s * exp(x) - max_f)
     end

    λ = best_λ0
    fλ = fx(λ)

    if -1e10 < max_f < 1e10
        for j in 1:n_it
            min_eff = min_λ0
            max_eff = max_λ0

            u = fλ * rand()
            reject = true
            if u < Inf
                while reject
                    λ = rand() * (max_eff - min_eff) + min_eff
                    if max_eff - min_eff <= 0
                        println(" oh oh $j, $(max_eff - min_eff), $fλ, $u, $max_f, $min_eff, $max_eff, $best_λ0, $min_λ0, $max_λ0")
                        break
                    end
                    fλ = fx(λ)
                    reject = fλ < u
                    if reject
                        if λ < best_λ0
                            min_eff = λ
                        else
                            max_eff = λ
                        end
                    end

                end
            else
                λ = log(former_λ0)
            end
        end
    else
        new_λ = former_λ0
    end

    new_λ0 = exp(λ)

    if isnan(new_λ0) || new_λ0 <= 0.
        new_λ0 = former_λ0#rates[i + 1]
    end

    return new_λ0
end

function extract_S(tree::Tree ,edge_trees::Array{EdgeTreeRates2,1})

    S::Float64 = 0.
    for et in edge_trees
        S += et.stem_rate[1] * et.effective_bl
    end

    return S

end

function extract_nextinct(tree::Tree ,edge_trees::Array{EdgeTreeRates2,1})

    n::Int64 = 2 - n_tip(tree)
    i=0
    for et in edge_trees
        n += n_tip(et.tree) - et.tip_number
    end

    return n
end

function n_tip(tree::Tree ,edge_trees::Array{EdgeTreeRates2,1})
    n = n_tip(tree)
    for i in 1:length(edge_trees)
        n += n_tip(edge_trees[i].tree) - 1
    end
    return n
end

function n_tip(tree::Tree)
    return Int64((tree.n_nodes+1)/2)
end

function draw_λ0!(tree::Tree, ε::Float64, rates::Array{Float64,1},
    edge_trees::Array{EdgeTreeRates2,1})

    S = extract_S(tree, edge_trees)
    S *= 1+ε
    n = n_tip(tree, edge_trees)
    n_events = 2 * n - 2 - n_extant_tips(tree,edge_trees)
    former_λ0 = rates[1]
    S /= former_λ0
    λ0 = slicing0(S, n_events, former_λ0)

    #rates[1] = λ0
    ratio = λ0 / former_λ0
    rates .*= ratio

    update_rates!(tree, ratio)
    for i in 1:length(edge_trees)
        edge_trees[i].stem_rate[1] *= ratio
    end
end

function draw_σ(relative_rates, α ; α0 = 1., β0 = 0.01)
    n = length(relative_rates)
    m = mean(relative_rates)
    α_n = α0 + n/2
    log_α = log(α)#m #log(m)
    β_n = β0

    for rate in relative_rates
        β_n += (rate - log_α)^2 /2
    end

    if isnan(β_n)
        println("inv gam $α0 $β0; $α_n , $β_n ;    ")
        print(relative_rates)
    end
    sqrt(rand(InverseGamma(α_n,β_n)))
end

function draw_α(relative_rates, σ ; α_0 = 0., σ_0 = 1)
    # know posterior for the mean of a normal with normal prior, see for example
    # www.ams.sunysb.edu/~zhu/ams570/Bayesian_Normal.pdf

    # so here we have a normal prior on log(α) (with default parameters (α_0,σ_0)=(0,1))
    n = length(relative_rates)
    #α_n = (α_0 * σ^2 + σ_0 * sum(map(x -> x^2, relative_rates)))/(σ^2 + n * σ_0^2)
    σ_n = sqrt((σ^2 * σ_0^2)/(σ^2 + n * σ_0^2))
    α_n = (α_0 * σ^2 + σ_0^2 * sum(relative_rates))/(σ^2 + n * σ_0^2)

    a = exp(rand(Normal(α_n, σ_n)))
    if a==0
        println("$α_n , $σ_n, $α_0 ; $σ ; $relative_rates")
    end
    return a
end

function draw_ε_crown(S::Float64, n_extinct::Int64, n_cond::Int64 ; n_it = 10) where {T<:Number}

    n_extinct_eff = n_extinct + 1  # because we want a flat prior on the nonloged value
    n_cond_eff = n_cond - 1
    δ = sqrt((n_extinct_eff + n_cond_eff - S)^2 + 4*n_extinct_eff*S)
    best_ε = log(max((n_extinct_eff + n_cond_eff - S + δ)/(2*S),(n_extinct_eff + n_cond_eff - S - δ)/(2*S)))
    max_f = best_ε * n_extinct_eff +
        log(exp(best_ε)+1) * n_cond_eff -
        S * exp(best_ε)

    function fx(x ; ne = n_extinct_eff, nc = n_cond_eff, s = S)
         exp(x * ne  + log(exp(x)+1) * nc - s * exp(x) - max_f)
     end

    λ = best_ε
    fλ = fx(λ)

    for j in 1:n_it
        u = fλ * rand()
        reject = true

        min_eff = best_ε - 10
        max_eff = best_ε + 10

        while reject
            λ = rand() * (max_eff - min_eff) + min_eff
            fλ = fx(λ)

            #println("$fλ , $u , $λ, $best_ε ")

            reject = fλ < u

            if reject
                if λ < best_ε
                    min_eff = λ
                else
                    max_eff = λ
                end
            end
        end
    end

    if !isnan(λ)
        return exp(λ)
    else
        return 0.001
    end
end

function draw_ε_crown(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}, lefts::Array{Int64,1})
    S = extract_S(tree, edge_trees)
    n_extinct = extract_nextinct(tree, edge_trees)
    crown_edges = [1, 1 + lefts[1]]

    n_cond = 0
    if tree.offsprings[1].n_nodes > 1
        n_cond += 1
    end
    if tree.offsprings[2].n_nodes > 1
        n_cond += 1
    end
    bl = [tree.offsprings[1].branch_length, tree.offsprings[2].branch_length]
    for edge in 1:2
        if edge_trees[crown_edges[edge]].tree.n_nodes > 1
            n_cond -= 1
            ltt = LTT(edge_trees[crown_edges[edge]].tree)
            s = edge_trees[crown_edges[edge]].scale
            b = bl[edge]*s
            for e in 1:length(ltt[2])
                if ltt[2][e] == 1
                    n_cond += 1
                end
                if ltt[1][e] >= b
                    break
                end

            end
        end
    end

    return draw_ε_crown(S, n_extinct, n_cond)
end

function draw_ε_crown_priorUnif(S::Float64, n_extinct::Int64, n_cond::Int64 ; n_it = 10) where {T<:Number}

    n_extinct_eff = n_extinct + 1  # because we want a flat prior on the nonloged value
    n_cond_eff = n_cond - 1
    δ = sqrt((n_extinct_eff + n_cond_eff - S)^2 + 4*n_extinct_eff*S)
    best_ε = log(max((n_extinct_eff + n_cond_eff - S + δ)/(2*S),(n_extinct_eff + n_cond_eff - S - δ)/(2*S)))
    best_ε = min(best_ε,log(990))
    max_f = best_ε * n_extinct_eff +
        log(exp(best_ε)+1) * n_cond_eff -
        S * exp(best_ε)

    xmax = log(1000)

    function fx(x ; ne = n_extinct_eff, nc = n_cond_eff, s = S)
        if x > xmax
            0
        else
            exp(x * ne  + log(exp(x)+1) * nc - s * exp(x) - max_f)
        end
     end

    λ = best_ε
    fλ = fx(λ)
    for j in 1:n_it
        u = fλ * rand()
        reject = true

        min_eff = best_ε - 10
        max_eff = min(best_ε + 10, xmax)

        while reject
            λ = rand() * (max_eff - min_eff) + min_eff
            fλ = fx(λ)

            #println("$fλ , $u , $λ, $best_ε ")

            reject = fλ < u

            if reject
                if λ < best_ε
                    min_eff = λ
                else
                    max_eff = λ
                end
            end
        end
    end

    if !isnan(λ)
        return exp(λ)
    else
        return 0.001
    end
end

function draw_ε_crown_priorUnif(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}, lefts::Array{Int64,1})
    S = extract_S(tree, edge_trees)
    n_extinct = extract_nextinct(tree, edge_trees)
    crown_edges = [1, 1 + lefts[1]]

    n_cond = 0
    if tree.offsprings[1].n_nodes > 1
        n_cond += 1
    end
    if tree.offsprings[2].n_nodes > 1
        n_cond += 1
    end
    bl = [tree.offsprings[1].branch_length, tree.offsprings[2].branch_length]
    for edge in 1:2
        if edge_trees[crown_edges[edge]].tree.n_nodes > 1
            n_cond -= 1
            ltt = LTT(edge_trees[crown_edges[edge]].tree)
            s = edge_trees[crown_edges[edge]].scale
            b = bl[edge]*s
            for e in 1:length(ltt[2])
                if ltt[2][e] == 1
                    n_cond += 1
                end
                if ltt[1][e] >= b
                    break
                end

            end
        end
    end

    return draw_ε_crown_priorUnif(S, n_extinct, n_cond)
end

function draw_ε_crown_priorln(S::Float64, n_extinct::Int64, n_cond::Int64 ; n_it = 10, logε0 = 0., sd = 0.5) where {T<:Number}

    n_extinct_eff = n_extinct
    n_cond_eff = n_cond - 1
    δ = sqrt((n_extinct_eff + n_cond_eff - S)^2 + 4*n_extinct_eff*S)

    best_ε = log(max((n_extinct_eff + n_cond_eff - S + δ)/(2*S),(n_extinct_eff + n_cond_eff - S - δ)/(2*S)))
    max_f = best_ε * n_extinct_eff +
        log(exp(best_ε)+1) * n_cond_eff -
        S * exp(best_ε)

    function fx(x ; ne = n_extinct_eff, nc = n_cond_eff, s = S)
         exp(x * ne  + log(exp(x)+1) * nc - s * exp(x) - max_f - (x - logε0)^2/(2*sd^2))
     end
     function df(x;  ne = n_extinct_eff, nc = n_cond_eff, s = S)
         ne + exp(x)/(exp(x) + 1) * nc - s * exp(x) - (x - logε0)/(sd^2)
    end

    λ = best_ε
    fλ = fx(λ)

    for j in 1:n_it
        u = fλ * rand()
        reject = true

        min_eff = best_ε - 10
        max_eff = best_ε + 10

        while reject
            λ = rand() * (max_eff - min_eff) + min_eff
            fλ = fx(λ)

            #println("$fλ , $u , $λ, $best_ε ")

            reject = fλ < u

            dfλ = df(λ)
            if reject
                if dfλ > 0
                    min_eff = λ
                else
                    max_eff = λ
                end
            end
        end
    end

    if !isnan(λ)
        return exp(λ)
    else
        return 0.001
    end
end

function draw_ε_crown_priorln(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}, lefts::Array{Int64,1} ; logε0 = 1., sd = 0.5)
    S = extract_S(tree, edge_trees)
    n_extinct = extract_nextinct(tree, edge_trees)
    crown_edges = [1, 1 + lefts[1]]

    n_cond = 0
    if tree.offsprings[1].n_nodes > 1
        n_cond += 1
    end
    if tree.offsprings[2].n_nodes > 1
        n_cond += 1
    end
    bl = [tree.offsprings[1].branch_length, tree.offsprings[2].branch_length]
    for edge in 1:2
        if edge_trees[crown_edges[edge]].tree.n_nodes > 1
            n_cond -= 1
            ltt = LTT(edge_trees[crown_edges[edge]].tree)
            s = edge_trees[crown_edges[edge]].scale
            b = bl[edge]*s
            for e in 1:length(ltt[2])
                if ltt[2][e] == 1
                    n_cond += 1
                end
                if ltt[1][e] >= b
                    break
                end

            end
        end
    end

    return draw_ε_crown_priorln(S, n_extinct, n_cond,  logε0 = logε0, sd = sd)
end

function draw_λi_quad!(rates::Array{Float64,1}, edge_trees::Array{EdgeTreeRates2,1},
    σ::Float64, α::Float64, ε::Float64, tree::Tree, lefts::Array{Int64,1}; n_it = 5, bounds = 10)

    log_α = log(α)
    function aux(subtree, id, parent_rate, sampled_rates)

        if subtree.n_nodes < 2
            former_λ = subtree.attributes[1]
            x0_eff = parent_rate + log_α

            n = 0
            bl_eff = 0.
            if id>0
                bl_eff = edge_trees[id].effective_bl * (1 + ε)
                #s = sum(extract_branch_lengths(edge_trees[id].tree) .*extract_rates(edge_trees[id].tree) / edge_trees[id].scale)
                n = edge_trees[id].tree.n_nodes - edge_trees[id].tip_number
            end

            λ = slicing(x0_eff, σ, bl_eff, n, former_λ)
            pushfirst!(sampled_rates, λ)

            return  n, bl_eff*λ

        else

            tip_λ = subtree.attributes[1]
            former_λ = subtree.attributes[1]
            if id>0
                tip_λ *= edge_trees[id].tip_rate
            end
            tip_λ = log(tip_λ)

            n_right, bl_right = aux(subtree.offsprings[2], id + 1 + lefts[id+1], tip_λ, sampled_rates)
            n_left, bl_left = aux(subtree.offsprings[1], id + 1, tip_λ, sampled_rates)

            x0_eff = parent_rate + log_α
            n = n_left + n_right
            bl_eff = (bl_right + bl_left) / former_λ

            λ = 0.

            if id>0
                bl_eff += edge_trees[id].effective_bl * (1 + ε)
                n += edge_trees[id].tree.n_nodes - edge_trees[id].tip_number

                λ = slicing(x0_eff, σ, bl_eff, n, former_λ)
            else
                λ = slicing0(bl_eff, n, former_λ)
            end

            ratio = λ/former_λ

            for k in 1:(subtree.n_nodes - 1)
                sampled_rates[k] *= ratio
            end

            pushfirst!(sampled_rates, λ)
            return  n, bl_eff*λ
        end
    end

    r = Array{Float64,1}(undef,0)
    aux(tree, 0, 0., r);
    rates[1] = r[1]

    for i in 1:length(edge_trees)
        rates[i+1] = r[i+1]
        edge_trees[i].stem_rate[1] = rates[i+1]
    end

    update_rates!(tree,rates)
end

function draw_λi_quad!(rates::Array{Float64,1}, edge_trees::Array{EdgeTreeRates2,1}, σ::Float64, α::Float64, ε::Float64, tree::Tree; n_it = 5, bounds = 10)
    lefts = n_left(tree)
    draw_λi_quad!(rates, edge_trees, σ, α, ε, tree, lefts, n_it = n_it, bounds = bounds)
end

function draw_λ0_slicing!(rates::Array{Float64,1}, edge_trees::Array{EdgeTreeRates2,1}, σ::Float64, α::Float64, ε::Float64, lefts; n_it = 5, bounds = 20)

    parent_λ = NaN

    daughter_λ = map(log,[rates[2],rates[2+lefts[1]]])
    x0_eff = (-2*log(α)+daughter_λ[1]+daughter_λ[2])/2
    σ_eff = σ / sqrt(2)

    law = Normal(x0_eff, σ_eff)
    new_λ = exp(rand(law,1)[1])
    rates[1] = new_λ

    return new_λ
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# test abstract trees
let
    ID_TYPE = VtxID
    EL_TYPE = Int
    a = GraphUtils.TreeNode{EL_TYPE,ID_TYPE}(1)
    b = GraphUtils.TreeNode{EL_TYPE,ID_TYPE}(2)
    c = GraphUtils.TreeNode{EL_TYPE,ID_TYPE}(3)
    GraphUtils._id_type(a)
    set_parent!(b,a)
    @test get_parent(a) === a
    @test get_parent(b) === a
    @test haskey(get_children(a),node_id(b))
    @test validate_tree(a)
    @test validate_tree(b)
    # Check cached behavior
    update_element!(a,2)
    @test GraphUtils.cached_node_up_to_date(a)
    @test !GraphUtils.cached_node_up_to_date(b)
    update_element!(b,2)
    @test GraphUtils.cached_node_up_to_date(b)
    @test GraphUtils.get_cached_value!(b) == 2
    # make tree cyclic and check for error in validation
    set_parent!(c,b)
    set_parent!(a,c)
    @test_throws ErrorException validate_tree(b)
    # check that parent-child edge is removed when a new parent is set.
    d = GraphUtils.TreeNode{EL_TYPE,ID_TYPE}(4)
    set_parent!(d,a)
    @test haskey(get_children(a),node_id(d))
    set_parent!(d,b)
    @test haskey(get_children(b),node_id(d))
    @test !haskey(get_children(a),node_id(d))
end
# test cached trees
let

end


text:
using GLVisualize, GLAbstraction, FileIO, GeometryTypes, Reactive, GLWindow, Colors

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end

description = """
This example shows how to set up a custom camera,
animate it and display it in a second screen. The first screen
allows you to look at the camera in the scene.
"""

"""
functions to halve some rectangle
"""
xhalf(r)  = SimpleRectangle(r.x, r.y, r.w÷2, r.h)
xhalf2(r) = SimpleRectangle(r.w÷2, r.y, r.w÷2, r.h)

"""
Makes a spiral
"""
function spiral(i, start_radius, offset)
	Point3f0(sin(i), cos(i), i/10f0) * (start_radius + ((i/2pi)*offset))
end
# 2D particles
curve_data(i, N) = Point3f0[spiral(i+x/20f0, 1, (i/20)+1) for x=1:N]
# create a spiraling camera path
const camera_path  = curve_data(1f0, 360)

area_a, area_b = x_partition(window.area, 50)
# create first screen for the camera
camera_screen = Screen(
	window, name=:camera_screen,
	area = area_a
)
# create second screen to _view the scene
scene_screen = Screen(
	window, name=:scene_screen,
	area = area_b
)

# create an camera eyeposition signal, which follows the path
eyeposition = map(timesignal) do t
    len = length(camera_path)
    index = round(Int, (t*(len-1))+1) # mod1, just to be on the save side
    Vec3f0(camera_path[index])
end

# create the camera lookat and up vector
lookatposition = Signal(Vec3f0(0))
upvector = Signal(Vec3f0(0,0,1))

# create a camera from these
cam = PerspectiveCamera(camera_screen.area, eyeposition, lookatposition, upvector)

"""
Simple visualization of a camera
"""
function showcam(cam::PerspectiveCamera)
    lookvec, posvec, upvec = map(f-> getfield(cam, f), (:lookat, :eyeposition, :up))
    positions = map((a,b) -> Point3f0[a,b], lookvec, posvec)
    lines = map(lookvec, posvec, upvec) do l,p,u
        dir = p-l
        right = normalize(cross(dir,u))
        Point3f0[
            l,p,
            p, p+u,
            p, p+right
        ]
    end
    colors = RGBA{Float32}[
        RGBA{Float32}(1,0,0,1),
        RGBA{Float32}(1,0,0,1),

        RGBA{Float32}(0,1,0,1),
        RGBA{Float32}(0,1,0,1),

        RGBA{Float32}(0,0,1,1),
        RGBA{Float32}(0,0,1,1),
    ]
    poses = visualize((Sphere(Point3f0(0), 0.05f0), positions))
    lines = visualize(lines, :linesegment, color=colors)
    Context(poses, lines)
end

# add the camera to the camera screen as the perspective camera
camera_screen.cameras[:perspective] = cam

# something to look at
cat = visualize(GLNormalMesh(loadasset("cat.obj")))

# visualize the camera path
camera_points = visualize(
    (Circle(Point2f0(0), 0.03f0), camera_path),
    color=RGBA{Float32}((Vec3f0(0,206,209)/256)..., 1f0), billboard=true
)

camera_path_line = visualize(camera_path, :lines)

# _view everything on the appropriate screen.
# we need to copy the cat, because _view inserts the camera into the
# actual render object. this is sub optimal and will get changed!
# Note, that this is a shallow copy, so the actual data won't be copied,
# just the data structure that holds the camera
_view(copy(cat), camera_screen, camera=:perspective)
_view(copy(cat), scene_screen, camera=:perspective)
_view(showcam(cam), scene_screen, camera=:perspective)
_view(camera_points, scene_screen, camera=:perspective)
_view(camera_path_line, scene_screen, camera=:perspective)

if !isdefined(:runtests)
	renderloop(window)
end


text:
function krest!(r::Robot)
    for side in (Nord,West,Sud,Ost)
    move_to_GorMarker!(r,side)
    move_to_Gor!(r,inverse(side))

    end
end

function move_to_GorMarker!(r::Robot,side::HorizonSide)
    while !isborder(r,side)
        move!(r,side)
        putmarker!(r)
    end
end

function move_to_Gor!(r::Robot,side::HorizonSide)
    while ismarker(r) == true
        move!(r,side)
    end
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))


text:
import Infernal_jll

"""
    path(name)

Path to Infernal binary `name`.
"""
path(name::String) = joinpath(Infernal_jll.artifact_dir, "bin", name)

"""
    cmd(name, [args...])

Returns a command to call an Infernal binary `name`.
"""
cmd(name::String, args...) = `$(path(name)) $args`

"""
    exe(name, [args...])

Runs the Infernal binary `name` with arguments `args`.
"""
exe(name::String, args...) = run(cmd(name, args...))