text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

"""
```
become!(lk::Link, bhv::Function, args1...; kwargs...)
become!(name::Symbol, ....)
```
Cause an actor to change behavior.

# Arguments
- actor `lk::Link` (or `name::Symbol` if registered),
- `bhv`: function implementing the new behavior,
- `args1...`: first arguments to `bhv` (without a possible `msg` argument),
- `kwargs...`: keyword arguments to `bhv`.
"""
become!(lk::Link, bhv::F, args::Vararg{Any, N}; kwargs...) where {F<:Function,N} =
    send!(lk, Become(Func(bhv, args...; kwargs...)))
become!(name::Symbol, args...; kwargs...) = become!(whereis(name), args...; kwargs...)

"""
```
call!(lk::Link, from::Link, args2...)
call!(lk::Link, args2...; timeout::Real=5.0)
call!(name::Symbol, ....)
```
Call an actor to execute its behavior and to send a
[`Response`](@ref) with the result.

# Arguments
- actor `lk::Link` (or `name::Symbol` if registered),
- `from::Link`: sender link; If `from` is omitted, `call!`
    blocks and returns the result.
- `args2...`: second arguments to the actor.
- `timeout::Real=5.0`: timeout in seconds.
"""
call!(lk::Link, from::Link, args...) = send!(lk, Call(args, from))
call!(lk::Link, args...; timeout::Real=5.0) = request!(lk, Call, args...; timeout=timeout)
call!(name::Symbol, args...; kwargs...) = call!(whereis(name), args...; kwargs...)

"""
```
cast!(lk::Link, args2...)
cast!(name::Symbol, args2...)
```
Cast a message to the actor `lk` (or `name` if registered) to
execute its behavior with `args2...` without sending a
response.

*Note:* you can prompt the returned value with [`query!`](@ref).
"""
cast!(lk::Link, args...) = send!(lk, Cast(args))
cast!(name::Symbol, args...) = cast!(whereis(name), args...)

"""
```
exec!(lk::Link, from::Link, f::Function, args...; kwargs...)
exec!(lk::Link, from::Link, fu::Func)
exec!(lk::Link, fu::Func; timeout::Real=5.0)
exec!(name::Symbol, ....)
```

Ask an actor `lk` (or `name` if registered) to execute an
arbitrary function and to send the returned value as
[`Response`](@ref).

# Arguments
- `lk::Link` or `name::Symbol` of the actor,
- `from::Link`: the link a `Response` should be sent to.
    If `from` is ommitted, `exec!` blocks, waits and returns
    the result. In that case there is a `timeout`.
- `f::Function, args...; kwargs...` or
- `fu::Func`: function arguments,
- `timeout::Real=5.0`: timeout in seconds. Set `timeout=Inf`
    if you don't want a timeout.
"""
exec!(lk::Link, from::Link, f::F, args...; kwargs...) where F<:Function =
    send!(lk, Exec(Func(f, args...; kwargs...), from))
exec!(lk::Link, from::Link, fu::Func) = send!(lk, Exec(fu, from))
exec!(lk::Link, fu::Func; timeout::Real=5.0) =
    request!(lk, Exec, fu; timeout=timeout)
exec!(name::Symbol, args...; kwargs...) = exec!(whereis(name), args...; kwargs...)

"""
```
exit!(lk::Link, code=0)
exit!(name::Symbol, code=0)
```
Tell an actor `lk` (or `name` if registered) to exit. If it has a [`term`](@ref _ACT)
function, it calls it with `code` as last argument.

!!! note "This behavior is not yet fully implemented!"

    It is needed for supervision.

"""
exit!(lk::Link, code=0) = send!(lk, Stop(code))
exit!(name::Symbol, code=0) = exit!(whereis(name), code)

"""
```
init!(lk::Link, f::Function, args...; kwargs...)
init!(name::Symbol, ....)
```
Tell an actor `lk` to save the function `f` with the given
arguments as an [`init`](@ref _ACT) function, to execute it
and to save the returned value as state [`sta`](@ref _ACT)
variable.

The `init` function will be called at actor restart.

!!! note "This behavior is not yet implemented!"

    It is needed for supervision.
"""
init!(lk::Link, f::F, args...; kwargs...) where F<:Function =
    send!(lk, Init(Func(f, args...; kwargs...)))
init!(name::Symbol, args...; kwargs...) = init!(whereis(name), args...; kwargs...)

"""
```
query!(lk::Link, from::Link, s::Symbol)
query!(lk::Link, s::Symbol; timeout::Real=5.0)
query!(name::Symbol, ....)
```

Ask the `lk` actor to send a [`Response`](@ref) message to
`from` with an internal state variable `s`.

If `from` is omitted, `query!` blocks and returns the response.
In that case there is a `timeout`.

- `s::Symbol` can be one of `:sta`, `:res`, `:bhv`, `:dsp`.

# Examples

```julia
julia> f(x, y; u=0, v=0) = x+y+u+v  # implement a behavior
f (generic function with 1 method)

julia> fact = Actor(f, 1)     # start an actor with it
Channel{Message}(sz_max:32,sz_curr:0)

julia> cast!(fact, 1)         # cast a second parameter to it
YAActL.Cast{Tuple{Int64}}((1,))

julia> query!(fact, :res)     # query the result
2

julia> query!(fact, :bhv)     # query the behavior
f (generic function with 1 method)

julia> set!(fact, state)      # set dispatch mode
YAActL.Update{Dispatch}(:dsp, state)

julia> query!(fact, :dsp)     # query the dispatch mode
state::Dispatch = 1

julia> update!(fact, 10)      # update the state
YAActL.Update{Int64}(:sta, 10)

julia> query!(fact)           # query the state variable
10

julia> call!(fact, 1)
11
```
"""
query!(lk::Link, from::Link, s::Symbol=:sta) = send!(lk, Query(s, from))
query!(lk::Link, s::Symbol=:sta; timeout::Real=5.0) = request!(lk, Query, s, timeout=timeout)
query!(name::Symbol, args...; kwargs...) = query!(whereis(name), args...; kwargs...)

"""
    self()

Get the [`Link`](@ref) of your actor.
"""
self() = task_local_storage("_ACT").self

"""
```
set!(lk::Link, dsp::Dispatch)
set!(name::Symbol, dsp::Dispatch)
```
Set the `lk` actor's [`Dispatch`](@ref) to `dsp`.
"""
set!(lk::Link, dsp::Dispatch) = update!(lk, dsp, s=:dsp)
set!(name::Symbol, dsp::Dispatch) = set!(whereis(name), dsp)

"""
```
term!(lk::Link, f::Function, args1...; kwargs...)
term!(name::Symbol, ....)
```
Tell an actor `lk` to execute a function `f` with the given
arguments when it terminates. `f` must accept a `code=0`
as last argument. This is added by the actor to `args1...`
when it [`exit!`](@ref)s.

!!! note "This behavior is not yet implemented!"

    It is needed for supervision.
"""
term!(lk::Link, f::F, args...; kwargs...) where F<:Function =
    send!(lk, Term(Func(f, args...; kwargs...)))
term!(name::Symbol, args...; kwargs...) = term!(whereis(name), args...; kwargs...)

"""
```
update!(lk::Link, x; s::Symbol=:sta)
update!(lk::Link, arg::Args)
update!(name::Symbol, ....)
```
Update an actor's internal state `s` with `args...`.

# Arguments
- an actor `lk::Link` or `name::Symbol` (if registered),
- `x`: value/variable to update the choosen state with,
- `arg::Args`: arguments to update,
- `s::Symbol`: can be one of `:sta`, `:dsp`, `:arg`, `:self`, `:name`.

*Note:* If you want to update the stored arguments to the
behavior function with `s=:arg`, you must pass an [`Args`](@ref)
to `arg`. If `Args` has keyword arguments, they are merged
with existing keyword arguments to the behavior function.

# Example
```julia
julia> update!(fact, 5)       # note that fact is in state dispatch
YAActL.Update{Int64}(:sta, 5)

julia> call!(fact, 5)         # call it with 5
10

julia> update!(fact, Args(0, u=5));  # update arguments

julia> call!(fact, 5)         # add the last result, 5 and u=5
20
```
"""
update!(lk::Link, x; s::Symbol=:sta) = send!(lk, Update(s, x))
update!(lk::Link, arg::Args) = send!(lk, Update(:arg, arg))
update!(name::Symbol, args...; kwargs...) = update!(whereis(name), args...; kwargs...)


text:
@recipe function f(fCCF::CCF)

    # plot configuration
    legend:= false

    titlefont:= font(8, "Courier")
    xtickfont:= font(3, "Courier")
    ytickfont:= font(3, "Courier")
    yguidefont:= font(5, "Courier")
    xguidefont:= font(5, "Courier")

    N = fCCF.N
    type = (fCCF.type == "cor") ? "Correlation" : "Covariance"
    auto_prefix = fCCF.auto ? "Auto" : "Cross"
    yguide := auto_prefix*'-'*type
    xguide := "Lag"

    # Center ticks
    lr = length(fCCF.ccf)
    gap = div(lr,11)
    mp = div(lr+1,2)
    rs = mp:gap:lr
    ls = mp-gap:-gap:1

    if fCCF.auto
        xticks := 1:gap:lr
    else
        xt = vcat(reverse(collect(ls)),collect(rs))
        xticks := (xt, xt .- xt[div(end+1,2)])
    end

    @series begin
        seriestype := :sticks
        ymirror := false
        guide_position := :left
        linewidth := 100/(fCCF.lag+1)
        fCCF.ccf
    end

    if fCCF.type == "cor"
        lg = fCCF.lag
        a1 = fCCF.alpha[1]; c1 = fCCF.ci[1]
        a2 = fCCF.alpha[2]; c2 = fCCF.ci[2]
        dc = c2-c1

        @series begin
            seriestype := :hline
            seriescolor := :black
            linealpha := 0.5
            linestyle := :dash
            [-c1,c1]
        end

        @series begin
            seriestype := :hline
            seriescolor := :black
            linealpha := 0.5
            linestyle := :dot
            ax = fCCF.auto ? lg : 2*lg + 1
            annotations :=
                [(ax,c1+dc/8,string("CI %",a1*100),font(3, "Courier")),
                 (ax,c2+dc/8,string("CI %",a2*100),font(3, "Courier"))]
            [-c2,c2]
        end
    end
end

@recipe function f(fSTL::STL)

    # subplots configuration
    legend:= false
    grid:= false
    layout := @layout [Data
                       Seasonal
                       Trend
                       Remainder]

    xtickfont:= font(3, "Courier")
    ytickfont:= font(3, "Courier")
    yguidefont:= font(5, "Courier")

    # reference bar
    a,b = extrema(skipmissing(values(fSTL.ta[:Seasonal].+
                                     fSTL.ta[:Trend].+
                                     fSTL.ta[:Remainder]))); hd = b-a
    a,b = extrema(skipmissing(values(fSTL.ta[:Seasonal]))); hs = b-a
    a,b = extrema(skipmissing(values(fSTL.ta[:Trend]))); ht = b-a
    a,b = extrema(skipmissing(values(fSTL.ta[:Remainder]))); hr = b-a
    mh = min(hd,hs,ht,hr)

    inset_subplots := [(1, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (2, bbox(-0.012, 0, 0.01, 1.0, :right)),
                       (3, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (4, bbox(-0.012, 0, 0.01, 1.0, :right))]

    @series begin
        subplot := 1
        yguide := "Data"
        xaxis := nothing
        bottom_margin := -7Plots.mm
        fSTL.ta[:Seasonal] .+ fSTL.ta[:Trend] .+ fSTL.ta[:Remainder]
    end

    @series begin
        subplot := 2
        yguide := "Trend"
        xaxis := nothing
        ymirror:=true
        guide_position:=:left
        bottom_margin := -7Plots.mm
        fSTL.ta[:Trend]
    end

    @series begin
        subplot := 3
        yguide := "Seasonal"
        xaxis := nothing
        seriestype := :sticks
        bottom_margin := -7Plots.mm
        fSTL.ta[:Seasonal]
    end

    @series begin
        subplot := 4
        yguide := "Remainder"
        seriestype := :sticks
        #bottom_margin:=0Plots.mm
        ymirror:=true
        guide_position:=:left
        fSTL.ta[:Remainder]
    end

    @series begin
        subplot := 5
        background_color_inside := nothing
        framestyle := :none
        seriestype := :bar
        ylims := (0,hd)
        mh:mh
    end

    @series begin
        subplot := 6
        background_color_inside := nothing
        framestyle := :none
        seriestype := :bar
        ylims := (0,ht)
        mh:mh
    end

    @series begin
        subplot := 7
        framestyle := :none
        seriestype := :bar
        ylims := (0,hs)
        mh:mh
    end

    @series begin
        subplot := 8
        framestyle := :none
        seriestype := :bar
        ylims := (0,hr)
        mh:mh
    end

end

ext: jl


text:
mutable struct SelfLoopCumulativeResNet <: Block
    self::Union{Block,Nothing}
    degree::Int
    function SelfLoopCumulativeResNet(self::Block, n::Int)
        new(self, n)
    end
    function SelfLoopCumulativeResNet(n::Int)
        new(nothing, n)
    end
end

function forward(m::SelfLoopCumulativeResNet, x::Variable)
    for i = 1:m.degree
        x = forward(m.self, x) + x
    end
    return x
end

function predict(m::SelfLoopCumulativeResNet,  x::AbstractArray)
    for i = 1:m.degree
        x = predict(m.self, x) + x
    end
    return x
end

function clone(this::SelfLoopCumulativeResNet; type::Type=Array{Float32})
    cloned = SelfLoopCumulativeResNet(this.degree)
    cloned.self = clone(this.self, type=type)
    return cloned
end

function unbiasedof(m::SelfLoopCumulativeResNet)
    return unbiasedof(m.self)
end

function weightsof(m::SelfLoopCumulativeResNet)
    return weightsof(m.self)
end

function gradsof(m::SelfLoopCumulativeResNet)
    return gradsof(m.self)
end

function zerograds!(m::SelfLoopCumulativeResNet)
    return zerograds!(m.self)
end

function paramsof(m::SelfLoopCumulativeResNet)
    return paramsof(m.self)
end

function xparamsof(m::SelfLoopCumulativeResNet)
    return xparamsof(m.self)
end

function nparamsof(m::SelfLoopCumulativeResNet)
    return nparamsof(m.self)
end

function bytesof(m::SelfLoopCumulativeResNet)
    return bytesof(m.self)
end

function to(type::Type, m::SelfLoopCumulativeResNet)
    m.self = to(type, m.self)
    return m
end

function to!(type::Type, m::SelfLoopCumulativeResNet)
    m.self = to(type, m.self)
    return nothing
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
maxintfloat(::Type{ArbFloat{P}}) where {P} = ArbFloat{P}(BigInt(2)^(P))
maxintfloat(::Type{ArbReal{P}}) where {P} = ArbReal{P}(BigInt(2)^(P))


text:
import Base: show

export xml

####### Show s-expression ###########
function show(io::IO, doc::Document)
    print(io, "(Document")
    if hasroot(doc)
        show(io, root(doc), 1)
    else
        print(io, ")")
    end
end

function show(io::IO, p::Node, depth::Integer = 0)
    println(io)
    print(io, "  "^depth)
    print(io, "(", typeof(p), " \"")
    print(io, istext(p) ? nodecontent(p) : nodename(p))
    print(io, "\"")
    children = nodes(p)
    if isempty(children)
        print(io, ")")
    else
        # println(io)
        for n in children
            show(io, n, depth + 1)
        end
        # print(io, "  "^depth)
        print(io, ")")
    end
end

####### XML ###########
xml(node) = xml(stdout, node)

function xml(io::IO, doc::Document)
    if hasroot(doc)
        xml(io, root(doc), 1)
    else
        print(io, "Document()")
    end
end

function xml(io::IO, n::TextNode)
    print(io, n.content)
end

function xml(io::IO, n::Node, depth::Integer = 0)
    print(io, "Unknown node type")
end

function xml(io::IO, n::AttributeNode)
   print(io, n.name, "=\"", n.value,"\"")
end

function xml(io::IO, n::TextNode, depth::Integer)
    print(io, "  "^depth)
    println(io, n.content)
end

function xml(io::IO, parent::ElementNode, depth::Integer = 0)
    print(io, "  "^depth)

    tag = nodename(parent)
    print(io, "<$tag")
    attrs = map(x -> x.name * "=\"$(x.value)\"", attributes(parent))
    attr_str = join(attrs, " ")

    if !isempty(attr_str)
        print(io, " ", attr_str)
    end

    children = nodes(parent)
    len = length(children)

    if len == 0
        println(io, "/>")
    elseif len == 1 && istext(first(children))
        print(io, ">")
        for n in children xml(io, n) end
    else
        println(io, ">")
        for n in children
            xml(io, n, depth + 1)
        end
        print(io, "  "^depth)
    end

    if len != 0
        println(io, "</$tag>")
    end
end


text:
using ITensors
using Printf
#using Logging
##Set global logging level to Debug to show timers
#global_logger(SimpleLogger(stdout, Logging.Debug))

function main()
  N = 100
  sites = spinOneSites(N)

  ampo = AutoMPO(sites)
  for j=1:N-1
    add!(ampo,"Sz",j,"Sz",j+1)
    add!(ampo,0.5,"S+",j,"S-",j+1)
    add!(ampo,0.5,"S-",j,"S+",j+1)
  end
  H = toMPO(ampo)

  psi = randomMPS(sites)

  sweeps = Sweeps(5)
  maxdim!(sweeps,10,20,100,100,200)
  cutoff!(sweeps,1E-11)
  @show sweeps

  energy,psi = @time dmrg(H,psi,sweeps,maxiter=2)
  @printf "Final energy = %.12f\n" energy
end
main()


text:
#=
ellipse.jl
=#

#using ImagePhantoms #: Object, Object2d

export Ellipse
export Circle
export phantom, radon, spectrum

"""
    Ellipse <: AbstractShape2
"""
struct Ellipse <: AbstractShape2 end

# constructors

"""
    Ellipse(cx, cy, rx=1, ry=rx, ϕ=0, value::Number=1)
    Ellipse(center::NTuple{2,RealU}, radii::NTuple{2,RealU}=(1,1), ϕ::RealU=0, v=1)
    Ellipse([6-vector])
Construct `Ellipse` object from parameters.
"""
function Ellipse(
    cx::RealU,
    cy::RealU,
    rx::RealU = oneunit(cx),
    ry::RealU = rx,
    ϕ::RealU = 0,
    value::Number = 1,
)
    (cx, cy, rx, ry) = promote(cx, cy, rx, ry)
    Object(Ellipse(), (cx,cy), (rx,ry), ϕ, value)
end

function Ellipse(
    center::NTuple{2,RealU},
    radii::NTuple{2,RealU} = (1,1) .* oneunit(center[1]),
    ϕ::RealU = 0,
    value::Number = 1,
)
    Ellipse(center..., radii..., ϕ, value)
end

function Ellipse(v::AbstractVector{<:Number})
    length(v) == 6 || throw(ArgumentError("$v wrong length"))
    Ellipse(v...)
end

# circles as a special case

"""
    Circle(x,y,r,v=1) (circle of radius `r` centered at `(x,y)`)
    Circle((x,y), r=1, v=1) ditto
    Circle([4-vector]) ditto
    Circle(r, v=1) centered at origin
Construct circle objects as special cases of `Ellipse` objects.
"""
Circle(cx::RealU, cy::RealU, r::RealU, v::Number = 1) =
    Ellipse(cx, cy,  r, r, 0, v)
Circle(center::NTuple{2,RealU}, r::RealU = oneunit(center[1]), v::Number = 1) =
    Circle(center..., r, v)
Circle(r::RealU, v::Number = 1) = Circle((zero(r), zero(r)), r, v)

function Circle(v::AbstractVector{<:Number})
    length(v) == 4 || throw(ArgumentError("$v wrong length"))
    Circle(v...)
end

# helper

# methods

"""
    phantom(ob::Object2d{Ellipse})
Returns function of `(x,y)` for making image.
"""
phantom(ob::Object2d{Ellipse}) =
    (x,y) -> (sum(abs2, coords(ob, x, y)) ≤ 1) * ob.value

"""
    radon_ellipse(r, ϕ, cx, cy, rx, ry, θ)
Radon transform of ellipse at point `(r,ϕ)`.
"""
function radon_ellipse(r, ϕ, cx, cy, rx, ry, θ)
    (sinϕ, cosϕ) = sincos(ϕ)
    r -= cx * cosϕ + cy * sinϕ # Radon translation property
    (sinϕ, cosϕ) = sincos(ϕ - θ) # Radon rotation property
    rp2 = abs2(rx * cosϕ) + abs2(ry * sinϕ) # square of projected radius
    return 2rx*ry / rp2 * sqrt(max(rp2 - abs2(r), 0*oneunit(r)^2))
end

"""
    radon(ob::Object2d{Ellipse})
Returns function of `(r,ϕ)` for making a sinogram.
"""
radon(ob::Object2d{Ellipse}) = (r,ϕ) -> ob.value *
    radon_ellipse(r, ϕ, ob.center..., ob.width..., ob.angle[1])

function spectrum_ellipse(fx, fy, cx, cy, rx, ry, θ)
    (kx, ky) = rotate2d(fx, fy, θ) # rect is rotated first, then translated
    return 2rx * exp(-2im*π*fx*cx) * # width=diameter=2*radius
           2ry * exp(-2im*π*fy*cy) * jinc(2sqrt(abs2(kx * rx) + abs2(ky * ry)))
end

"""
    spectrum(ob::Object2d{Ellipse})
Returns function of ``(f_x,f_y)`` for the spectrum (2D Fourier transform).
"""
spectrum(ob::Object2d{Ellipse}) = (fx,fy) -> ob.value *
    spectrum_ellipse(fx, fy, ob.center..., ob.width..., ob.angle[1])


text:
# Wrapping effect example from Rump paper (2015)

using AffineArithmetic
using ValidatedNumerics
using LinearAlgebra

A = 0.5 * [1 2; -1 1]
eigvals(A)
maximum(abs.(eigvals(A)))  # spectral radius

reset_affine_index()
X = Affine(-1..1)
Y = Affine(-1..1)

range(X)

XX = [X, Y]  # vector

for i in 1:10
    XX = A * XX
    println(i)
    println(XX)
    println(range.(XX))
    println()
end

vertices = []
using Iterators
for (ε1, ε2, ε3, ε4) in product( ([-1, 1] for i in 1:4)... )
    @show (ε1, ε2, ε3, ε4)
    x = 2 + ε1 - 2ε2 + 3ε3 - ε4
    y = 1 + 3ε1 - ε3 + 2ε4

    @show x, y

    push!(vertices, [x, y])
end

vertices
xs = [p[1] for p in vertices]
ys = [p[2] for p in vertices]

using Plots
gr()

scatter(xs, ys, aspect_ratio=:equal)
# Need the convex hull of these points for the zonotope


text:
using MATLAB

x = range(-10.0, stop=10.0, length=500)
mat"plot($x, sin($x))"  # evaluate a MATLAB function

y = range(2.0, stop=3.0, length=500)
mat"""
    $u = $x + $y
	$v = $x - $y
"""
@show u v

max_stars_repo_name: niladridas/UKF

max_issues_repo_name: niladridas/UKF

max_forks_repo_name: niladridas/UKF


text:
using PyPlot

xkcd()
fig = figure(figsize=(10,5))
ax = axes()

x = linspace(1,10,101)
y = sin(3x + cos(5x))
p = plot(x,y)

ax[:set_xlim]([0.0,6])
annotate("A little\nsketchy",xy=[0.98,.001],
arrowprops=["arrowstyle"=>"->"],xytext=[1.3,-0.3])

xticks([])
yticks([])
xlabel("TIME")
ylabel("PRODUCTIVITY")
title("An xkcd-style plot in Julia")

ax[:spines]["top"][:set_color]("none")
ax[:spines]["right"][:set_color]("none")


text:
"""Implementation of the Chebychev propagator."""

using SpecialFunctions
using LinearAlgebra

"""Calculate Chebychev coefficients.

Return an array of coefficiencts larger than `limit`.

# Arguments

* `Δ`: the spectral radius of the underlying operator
* `dt`: the time step
"""
function cheby_coeffs(Δ, dt; limit=1e-12)
    α = abs(0.5 * Δ * dt)
    coeffs = Float64[]
    a = besselj(0, α)
    append!(coeffs, a)
    ϵ = abs(a)
    i = 1
    while ϵ > limit
        a = 2 * besselj(i,α)
        append!(coeffs, a)
        ϵ = abs(a)
        i += 1
    end
    return coeffs
end

"""Calculate Chebychev coefficients in-place.

```julia
n = cheby_coeffs!(coeffs, Δ, dt, limit=1e-12)
```

overwrites the first `n` values in `coeffs` with new coefficients larger than
`limit` for the given new spectral radius `Δ` and time step `dt`. The `coeffs`
array will be resized if necessary, and may length > `n` on exit.
"""
function cheby_coeffs!(coeffs, Δ, dt, limit=1e-12)
    α = abs(0.5 * Δ * dt)
    a = besselj(0, α)
    coeffs[1] = a
    N = length(coeffs)
    ϵ = abs(a)
    n = 1
    while ϵ > limit
        n += 1
        a = 2 * besselj(n-1,α)
        coeffs[n] = a
        ϵ = abs(a)
        if (n >= N)
            N *= 2
            resize!(coeffs, N)
        end
    end
    return n
end

"""
Workspace for the Chebychev propagation routine.

```julia
    ChebyWrk(Ψ, Δ, E_min, dt; limit=1e-12)
```

initializes the workspace for the propagation of a state similar to Ψ under a
Hamiltonian with eigenvalues between `E_min` and `E_min + Δ`, and a time step
dt. Chebychev coefficients smaller than the given `limit` are discarded.
"""
mutable struct ChebyWrk{T, CFS, FT<:AbstractFloat}
    v0 :: T
    v1 :: T
    v2 :: T
    coeffs :: CFS
    n_coeffs :: Int64
    Δ :: FT
    E_min :: FT
    dt :: FT
    limit :: FT
    function ChebyWrk(Ψ::T, Δ::FT, E_min::FT, dt::FT;
                      limit::FT=1e-12) where {T, FT}
        v0::T = similar(Ψ)
        v1::T = similar(Ψ)
        v2::T = similar(Ψ)
        coeffs = cheby_coeffs(Δ, dt; limit=limit)
        n_coeffs = length(coeffs)
        new{T, typeof(coeffs), FT}(
            v0, v1, v2, coeffs, n_coeffs, Δ, E_min, dt, limit
        )
    end
end

"""Evaluate `Ψ = exp(-i H dt) Ψ` in-place.

```julia
cheby!(Ψ, H, dt, wrk; E_min=nothing, check_normalization=false)
```

Args:

* `Ψ`: on input, initial vector. Will be overwritten with result.
* `H`: Hermitian operator
* `dt`: time step
* `wrk`: internal workspace
* `E_min`: minimum eigenvalue of H, to be used instead of the `E_min` from the
   initialization of `wrk`. The same `wrk` may be used for different values
   `E_min`, as long as the spectra radius `Δ` and the time step `dt` are the
   same as those used for the initialization of `wrk`.
* `check_normalizataion`: perform checks that the H does not exceed the
  spectral radius for which the the workspace was initialized.

The routine will not allocate any internal storage. This implementation
requires `copyto!` `lmul!`, and `axpy!` to be implemented for `Ψ`, and the
three-argument `mul!` for `Ψ` and `H`.
"""
function cheby!(Ψ, H, dt, wrk; kwargs...)

    E_min = get(kwargs, :E_min, nothing)
    check_normalization = get(kwargs, :check_normalization, false)

    Δ = wrk.Δ
    β::typeof(wrk.E_min) = (Δ / 2) + wrk.E_min  # "normfactor"
    if E_min ≠ nothing
        β = (Δ / 2) + E_min
    end
    @assert dt ≈ wrk.dt "wrk was initialized for dt=$(wrk.dt), not dt=$dt"
    if dt > 0
        c = -2im / Δ
    else
        c = 2im / Δ
    end
    a = wrk.coeffs
    ϵ = wrk.limit
    @assert length(a) > 1 "Need at least 2 Chebychev coefficients"
    v0 = wrk.v0
    v1 = wrk.v1
    v2 = wrk.v2

    # v0 ↔ Ψ; Ψ = a[1] * v0
    copyto!(v0, Ψ)
    lmul!(a[1], Ψ)

    # v1 = -i * H_norm * v0 = c * (H * v0 - β * v0)
    mul!(v1, H, v0)
    axpy!(-β, v0, v1)
    lmul!(c, v1)

    # Ψ += a[2] * v1
    axpy!(a[2], v1, Ψ)

    c *= 2

    for i = 3 : wrk.n_coeffs

        # v2 = -2i * H_norm * v1 + v0 = c * (H * v1 - β * v1) + v0
        mul!(v2, H, v1)
        axpy!(-β, v1, v2)
        lmul!(c, v2)
        if check_normalization
            map_norm = abs(dot(v1, v2)) / (2 * norm(v1)^2)
            @assert(
                map_norm <= (1.0 + ϵ), "Incorrect normalization (E_min, wrk.Δ)"
            )
        end
        v2 .+= v0

        # Ψ += a[i] * v2
        axpy!(a[i], v2, Ψ)

        v0, v1, v2 = v1, v2, v0  # switch w/o copying

    end

    lmul!(exp(-im * β * dt), Ψ)

end

"""Evaluate `Ψ = exp(i- H dt) Ψ.

```julia
Ψ_out = cheby(Ψ, H, dt, wrk; E_min=nothing, check_normalization=false)
```

acts like [`cheby!`](@ref) but does not modify `Ψ` in-place.
"""
function cheby(Ψ, H, dt, wrk; kwargs...)

    E_min = get(kwargs, :E_min, nothing)
    check_normalization = get(kwargs, :check_normalization, false)

    Δ = wrk.Δ
    β::typeof(wrk.E_min) = (Δ / 2) + wrk.E_min  # "normfactor"
    if E_min ≠ nothing
        β = (Δ / 2) + E_min
    end
    @assert dt ≈ wrk.dt "wrk was initialized for dt=$(wrk.dt), not dt=$dt"
    if dt > 0
        c = -2im / Δ
    else
        c = 2im / Δ
    end
    a = wrk.coeffs
    ϵ = wrk.limit
    @assert length(a) > 1 "Need at least 2 Chebychev coefficients"
    v0 = wrk.v0
    v1 = wrk.v1
    v2 = wrk.v2

    v0 = Ψ
    Ψ = a[1] * v0

    v1 = c * (H * v0 - β * v0)
    Ψ += a[2] * v1

    c *= 2

    for i = 3 : wrk.n_coeffs

        v2 = H * v1
        v2 += - v1 * β
        v2 = c * v2
        if check_normalization
            map_norm = abs(dot(v1, v2)) / (2 * norm(v1)^2)
            @assert(
                map_norm <= (1.0 + ϵ), "Incorrect normalization (E_min, wrk.Δ)"
            )
        end
        v2 += v0

        Ψ += a[i] * v2

        v0, v1, v2 = v1, v2, v0  # switch w/o copying

    end

    return exp(-im * β * dt) * Ψ

end


text:
using SQLite
using Negotiations

if (!("test.sqlite" in readdir()))

    db = initialize_db("test")

    SQLite.execute(
        db,
        """
        INSERT INTO party
        VALUES
            (1, \"P1\", \"PARTY_1\"),
            (2, \"P2\", \"PARTY_2\"),
            (3, \"P3\", \"PARTY_3\");
        """
    )

    SQLite.execute(
        db,
        """
        INSERT INTO statement
        VALUES
            (1, \"Test statement 1\", \"This is the first test statement.\"),
            (2, \"Test statement 2\", \"This is the second test statement.\"),
            (3, \"Test statement 3\", \"This is the third test statement.\");
        """
    )

    SQLite.execute(
        db,
        """
        INSERT INTO opinion
        VALUES
            (1, 1, 1, \"We hold this position because of A.\"),
            (1, 2, -1, \"We hold this position because of B.\"),
            (1, 3, 0, \"We hold this position because of C.\"),
            (2, 1, -1, \"We hold this position because of D.\"),
            (2, 2, 0, \"We hold this position because of E.\"),
            (2, 3, -1, \"We hold this position because of F.\"),
            (3, 1, 1, \"We hold this position because of G.\"),
            (3, 2, 1, \"We hold this position because of H.\"),
            (3, 3, 0, \"We hold this position because of I.\");
        """
    )

end


text:
# test for tree2Matrix
# Claudia October 2015

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

T = readTopology("(A,(B,(C,D)));");
S = ["A","B","C","D"]
M = tree2Matrix(T,S)


text:
struct NoMethodError <: Exception
    func
    types::Tuple
end
NoMethodError(f, types::Type) = NoMethodError(f, (types.parameters...))
function Base.showerror(io::IO, e::NoMethodError)
    args = join(map(t->"::$t", e.types), ", ")
    print(io, "$(e.func)($args) couldn't be found")
end

# Give string based stages a type
struct CodeLLVM
    source::String
end
function CodeLLVM(f, types)
    src = sprint() do io
        code_llvm(io, f, types)
    end
    CodeLLVM(src)
end
struct CodeNative
    source::String
end
function CodeNative(f, types)
    src = sprint() do io
        code_native(io, f, types)
    end
    CodeNative(src)
end

# deal with all variances in base that should really be tuples but are something else
to_tuple(x) = (x,)
to_tuple(x::Core.SimpleVector) = tuple(x...)
to_tuple(x::Tuple) = x
to_tuple(x::AbstractVector) = tuple(x...)
to_tuple(x::Type{Tuple}) = tuple(x.parameters...)

# typeof working with concrete and types at the same time
function _typeof(x::Type{T}) where T
	Type{T}
end
function _typeof(x::T) where T
	T
end

function juliabasepath(file)
    srcdir = joinpath(Sys.BINDIR,"..","..","base")
    releasedir = joinpath(Sys.BINDIR,"..","share","julia","base")
    normpath(joinpath(isdir(srcdir) ? srcdir : releasedir, file))
end

function get_source_file(path::AbstractString, ln::Int)
    isfile(path) && return path
    # if not a file, it might be in julia base
    file = juliabasepath(path)
    if !isfile(file)
        throw(LoadError(path, ln, ErrorException("file $path not found")))
    end
    file
end

function get_lambda(pass, ftype::DataType, types)
    world = typemax(UInt)
    if !isclosure(ftype)
        ftype = Type{ftype}
    end
    tt = Tuple{ftype, to_tuple(types)...}
    (ti, env, meth) = Base._methods_by_ftype(tt, 1, world)[1]
    meth = Base.func_for_method_checked(meth, tt)
    params = Core.Inference.InferenceParams(world)
    (_, code, ty) = Core.Inference.typeinf_code(meth, ti, env, false, false, params)
    code
end

function get_ast(pass, f, types)
    lambda = get_lambda(pass, f, types)
    get_ast(lambda)
end
function get_ast(li::Core.CodeInfo)
    ast = li.code
    if isa(ast, Vector{UInt8})
        return Base.uncompressed_ast(li)
    end
    ast
end

function get_static_parameters(f, types)
    m = get_method(f, types)
    mi = m.specializations.func
    spnames = map(x->x.name, to_tuple(m.tvars))
    sptypes = to_tuple(mi.sparam_vals)
    spnames, sptypes
end

function get_lambda(pass, f, types, optimize = false)
    if isintrinsic(f)
        error("$f is an intrinsic function")
    end
    lambda = try
        if pass == code_typed
            pass(f, types, optimize = optimize)
        else
            pass(f, types)
        end
    catch e
        error("Couldn't get lambda for $f $types:\n$e")
    end
    if isa(lambda, Vector)
        if isempty(lambda)
            throw(NoMethodError(f, types))
        elseif length(lambda) == 1
            lam = first(lambda)
            if isa(lam, Pair)
                return first(lam)
            else
                return lam
            end
        else
            error("
                More than one method found for signature $f $types.
                Please use more specific types!
            ")
        end
    else
        isa(lambda, LambdaInfo) && return lambda
        error("Not sure what's up with returntype of $pass. Returned: $lambda")
    end
end
function get_method(f, types::Type)
    get_method(f, (types.parameters...))
end
function get_method(ftype::DataType, types::Tuple)
    world = typemax(UInt)
    if !isclosure(ftype)
        ftype = Type{ftype}
    end
    tt = Tuple{ftype, to_tuple(types)...}
    (ti, env, meth) = Base._methods_by_ftype(tt, 1, world)[1]
    Base.func_for_method_checked(meth, tt)
end
function get_method(f, types::Tuple)
    if !all(isleaftype, types)
        error("Not all types are concrete: $types")
    end
    # make sure there is a specialization with precompile
    # TODO, figure out a better way, since this method is not very reliable.
    # (I think, e.g. anonymous functions don't work)
    precompile(f, (types...))
    x = methods(f, types)
    if isempty(x)
        throw(NoMethodError(f, types))
    elseif length(x) != 1
        error("
            More than one method found for signature $f $types.
            Please use more specific types!
        ")
    end
    first(x)
end

"""
Looks up the source of `method` in the file path found in `method`.
Returns the AST and source string, might throw an LoadError if file not found.
"""
function get_source_at(file, linestart::Int)
    file = get_source_file(file, linestart)
    code, str = open(file) do io
        line = ""
        for i=1:linestart-1
            line = readline(io)
        end
        try # lines can be one off, which will result in a parse error
            parse(line)
        catch e
            line = readline(io)
        end
        while !eof(io)
            line = line * "\n" * readline(io)
            e = Base.parse_input_line(line; filename=file)
            if !(isa(e,Expr) && e.head === :incomplete)
                return e, line
            end
        end
    end
    code, str
end

function get_source(method::Method)
    file = string(method.file)
    linestart::Int = method.line
    code, str = get_source_at(file, linestart)
    # for consistency, we always return the `function f(args...) end` form
    long = MacroTools.longdef(code)
    # and return only the body
    long.args[2], str
end

function slot_vector(lam_typed)
    if isa(lam_typed.slotnames, Void) || isa(lam_typed.slottypes, Void)
        return [(SlotNumber(-1), ("", Void))]
    end
    slotnames = copy(lam_typed.slotnames)
    slottypes = copy(lam_typed.slottypes)

    slots = [(SlotNumber(i), (string(name), slottypes[i])) for (i, name) in enumerate(slotnames)]
    ssaslot = [(SSAValue(i-1), ("ssa_$(i-1)", t)) for (i,t) in enumerate(lam_typed.ssavaluetypes)]
    vcat(slots, ssaslot)
end

function slot_dictionary(lam_typed)
    slots = slot_vector(lam_typed)
    s_dict = Dict()
    for (k, (name, T)) in slots
        if !isa(k, SSAValue)
            s_dict[k] = T
        end
        s_dict[Symbol(name)] = T
    end
    s_dict
end

"""
Returns an AST most similar to what you would get from a macro
"""
function macro_form(f, types)
    method = get_method(f, types)
    local code::Expr; local str::String
    try
        code, str = get_source(method)
    catch e
        if isa(e, LoadError)
            # file for method not found, we need to get the source the low level way
            # (this is the case when e.g. code comes from sysimg without source install or REPL/IJulia)
            # allthough for the REPL/IJUlia we might be able to figure something out, since it's theoretically saved somewhere
            # .. and what about fused functions generated for broadcast?
            code = code_lowered_clean(f, types)
            str  = sprint() do io
                show(io, code)
            end
        else
            rethrow(e)
        end
    end
end

"""
Returns only one return type. If multiple are applicable, returns a Union.
"""
function return_type(f, types)
    x = Base.return_types(f, types)
    if isempty(x)
        throw(NoMethodError(f, types))
    elseif length(x) == 1
        x[1]
    else
        Union{x...}
    end
end

function clean_form(f, types)
    ast, str = macro_form(f, types)
    static_params = get_static_parameters(f, types)
    ast = filter_expr(x-> !is_linenumber(x), ast)
    remove_static_params(ast, static_params)
end

function clean_typed(f, types)
    ast = clean_form(f, types)
    sdict = slot_dictionary(get_lambda(code_typed, f, types))
    insert_types(ast, sdict)
end

const _source_map = Dict{Function, Expr}()

macro preserve_source(expr)
    # TODO this doesn't seem to work for curly
    result = @match expr begin
        (
            f_(args__) = body_ |
            function f_(args__) body_ end |
            f_{sargs__}(args__) = body_ |
            function f_{sargs__}(args__) body_ end
        ) => (f, args, body)
    end
    if result == nothing
        error("Expr $expr doesn't declare a function")
    end
    func_expr = MacroTools.longdef(expr) # always have the function be in long form (function x(args...) end)
    quote
        # evaluate function
        $expr
        # insert into source map
        _source_map[$f] = $(Expr(:quote, func_expr))
    end
end


text:
#
# This example shows how to write a basic calculator with variables.
#
# Comparison with Lark
#

using Lerche

calc_grammar = """
    ?start: sum
          | NAME "=" sum    -> assign_var

    ?sum: product
        | sum "+" product   -> add
        | sum "-" product   -> sub

    ?product: atom
        | product "*" atom  -> mul
        | product "/" atom  -> div

    ?atom: NUMBER           -> number
         | "-" atom         -> neg
         | NAME             -> var
         | "(" sum ")"

    %import common.CNAME -> NAME
    %import common.NUMBER
    %import common.WS_INLINE

    %ignore WS_INLINE
"""

struct CalculateTree <: Transformer
    vars::Dict
end

CalculateTree() = CalculateTree(Dict())

@inline_rule assign_var(s::CalculateTree, name, value) = begin
    s.vars[name] = value
    return value
end

@inline_rule var(s::CalculateTree, name) =  s.vars[name]
@rule add(s::CalculateTree, nums) = nums[1]+nums[2]
@rule sub(s::CalculateTree, nums) = nums[1]-nums[2]
@rule mul(s::CalculateTree, nums) = nums[1]*nums[2]
@rule div(s::CalculateTree, nums) = nums[1]/nums[2]
@rule number(s::CalculateTree, nums) = Base.parse(Float64,nums[1])
@rule neg(s::CalculateTree, nums) = -1*nums[1]

calc_parser = Lark(calc_grammar, parser="lalr", transformer=CalculateTree())
calc(a::String) = Lerche.parse(calc_parser,a)

test() = begin
    println(calc("a = 1+2"))
    println(calc("1+a*-3"))
end

main() = begin
    while true
        print("> ")
        s = ""
        s = readline(stdin)
        if s == "" break end
        println(calc(s))
    end
end

main()
#test()


text:
using Documenter, FortranFiles

makedocs(
   modules  = [FortranFiles],
   sitename = "FortranFiles.jl",
   format = Documenter.HTML(prettyurls = false),
   pages = [
      "Home" => "index.md",
      "files.md",
      "types.md",
      "read.md",
      "write.md",
      "exceptions.md",
      "tests.md",
      "Index" => "theindex.md"
   ]
)

deploydocs(
   repo   = "github.com/traktofon/FortranFiles.jl.git",
   branch = "gh-pages"
)


text:
module ImageBinarization

using ImageContrastAdjustment
using ColorTypes
using ColorVectorSpace
using LinearAlgebra
using HistogramThresholding
using Polynomials

abstract type BinarizationAlgorithm end
struct Otsu <: BinarizationAlgorithm end
struct Polysegment <: BinarizationAlgorithm end
struct UnimodalRosin <: BinarizationAlgorithm end
struct MinimumIntermodes <: BinarizationAlgorithm end
struct Intermodes <: BinarizationAlgorithm end
struct MinimumError <: BinarizationAlgorithm end
struct Balanced <: BinarizationAlgorithm end
struct Yen <: BinarizationAlgorithm end

include("balanced.jl")
include("otsu.jl")
include("polysegment.jl")
include("unimodal.jl")
include("minimum.jl")
include("intermodes.jl")
include("minimum_error.jl")
include("yen.jl")

export
	# main functions
    binarize,
    Otsu,
    Balanced,
    Yen,
    Polysegment,
    MinimumIntermodes,
    Intermodes,
    MinimumError,
    UnimodalRosin
end # module


text:
"""
    Map <: AbstractDict{MOI.VariableIndex, AbstractBridge}

Mapping between bridged constraints and the bridge that bridged the constraint.
"""
struct Map <: AbstractDict{MOI.ConstraintIndex, AbstractBridge}
    # Constraint Index of bridged constraint -> Bridge.
    # It is set to `nothing` when the constraint is deleted.
    bridges::Vector{Union{Nothing, AbstractBridge}}
    # Constraint Index of bridged constraint -> Constraint type.
    constraint_types::Vector{Tuple{DataType, DataType}}
    # The order of the keys is used in `keys_of_type` which is used by
    # `ListOfConstraintIndices`. Therefore they need to be in the order
    # of creation so we need `OrderedDict` and not `Dict`.
    # For `SingleVariable` constraints: (variable, set type) -> bridge
    single_variable_constraints::OrderedDict{Tuple{Int64, DataType}, AbstractBridge}
    # For `VectorVariable` constraints: (variable, set type) -> bridge
    vector_of_variables_constraints::OrderedDict{Tuple{Int64, DataType}, AbstractBridge}
end
function Map()
    return Map(Union{Nothing, AbstractBridge}[],
               Tuple{DataType, DataType}[],
               OrderedDict{Tuple{Int64, DataType}, AbstractBridge}(),
               OrderedDict{Tuple{Int64, DataType}, AbstractBridge}())
end

# Implementation of `AbstractDict` interface.

function Base.isempty(map::Map)
    return all(bridge -> bridge === nothing, map.bridges) &&
        isempty(map.single_variable_constraints) &&
        isempty(map.vector_of_variables_constraints)
end
function Base.empty!(map::Map)
    empty!(map.bridges)
    empty!(map.constraint_types)
    empty!(map.single_variable_constraints)
    empty!(map.vector_of_variables_constraints)
    return map
end
function Base.length(map::Map)
    return count(bridge -> bridge !== nothing, map.bridges) +
        length(map.single_variable_constraints) +
        length(map.vector_of_variables_constraints)
end
function Base.haskey(map::Map, ci::MOI.ConstraintIndex{F, S}) where {F, S}
    return 1 ≤ ci.value ≤ length(map.bridges) &&
        map.bridges[ci.value] !== nothing &&
        (F, S) == map.constraint_types[ci.value]
end
function Base.getindex(map::Map,
                       ci::MOI.ConstraintIndex)
    return map.bridges[ci.value]
end
function Base.haskey(map::Map, ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S
    return haskey(map.single_variable_constraints, (ci.value, S))
end
function Base.getindex(map::Map,
                       ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S
    return map.single_variable_constraints[(ci.value, S)]
end
function Base.haskey(map::Map, ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S
    return haskey(map.vector_of_variables_constraints, (ci.value, S))
end
function Base.getindex(map::Map,
                       ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S
    return map.vector_of_variables_constraints[(ci.value, S)]
end
function Base.delete!(map::Map, ci::MOI.ConstraintIndex)
    map.bridges[ci.value] = nothing
    return map
end
function Base.delete!(map::Map, ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S
    delete!(map.single_variable_constraints, (ci.value, S))
    return map
end
function Base.delete!(map::Map, ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S
    delete!(map.vector_of_variables_constraints, (ci.value, S))
    return map
end
function Base.values(map::Map)
    return Base.Iterators.flatten((
        # See comment in `values(::Variable.Map)`.
        Base.Iterators.Filter(bridge -> bridge !== nothing, map.bridges),
        values(map.single_variable_constraints),
        values(map.vector_of_variables_constraints)
    ))
end

# Implementation of iterate: it should combine non-variablewise constraints,
# `SingleVariable` constraints and `VectorOfVariables` constraints.
function _iterate_vov(map::Map, elem_state=iterate(map.vector_of_variables_constraints))
    if elem_state === nothing
        return nothing
    else
        i, S = elem_state[1].first
        bridge = elem_state[1].second
        ci = MOI.ConstraintIndex{MOI.VectorOfVariables, S}(i)
        return ci => bridge, (3, elem_state[2])
    end
end
function _iterate_sv(map::Map, elem_state=iterate(map.single_variable_constraints))
    if elem_state === nothing
        return _iterate_vov(map)
    else
        i, S = elem_state[1].first
        bridge = elem_state[1].second
        ci = MOI.ConstraintIndex{MOI.SingleVariable, S}(i)
        return ci => bridge, (2, elem_state[2])
    end
end
function _iterate(map::Map, state=1)
    while state ≤ length(map.bridges) && map.bridges[state] === nothing
        state += 1
    end
    if state > length(map.bridges)
        return _iterate_sv(map)
    else
        F, S = map.constraint_types[state]
        return MOI.ConstraintIndex{F, S}(state) => map.bridges[state], (1, state + 1)
    end
end
Base.iterate(map::Map) = _iterate(map)
function Base.iterate(map::Map, state)
    if state[1] == 1
        return _iterate(map, state[2])
    elseif state[1] == 2
        return _iterate_sv(map, iterate(map.single_variable_constraints, state[2]))
    else
        return _iterate_vov(map, iterate(map.vector_of_variables_constraints, state[2]))
    end
end

# Custom interface for information needed by `AbstractBridgeOptimizer`s that is
# not part of the `AbstractDict` interface.

"""
    number_of_type(map::Map, C::Type{<:MOI.ConstraintIndex})

Return the number of keys of type `C` in `map`.
"""
function number_of_type end

"""
    keys_of_type(map::Map, C::Type{<:MOI.ConstraintIndex})

Return a list of all the keys of type `C` in `map` in order order in which they
were created with `add_key_for_bridge`.
"""
function keys_of_type end

function number_of_type(map::Map, C::Type{MOI.ConstraintIndex{F, S}}) where {F, S}
    return count(i -> haskey(map, C(i)), eachindex(map.bridges))
end
function keys_of_type(map::Map, C::Type{MOI.ConstraintIndex{F, S}}) where {F, S}
    return Base.Iterators.Filter(
        ci -> haskey(map, ci),
        MOIU.LazyMap{C}(
            i -> C(i), eachindex(map.bridges)))
end
function number_of_type(map::Map, C::Type{MOI.ConstraintIndex{MOI.SingleVariable, S}}) where S
    return count(key -> key[2] == S, keys(map.single_variable_constraints))
end
function keys_of_type(map::Map, C::Type{MOI.ConstraintIndex{MOI.SingleVariable, S}}) where S
    return MOIU.LazyMap{C}(
        key -> C(key[1]),
        Base.Iterators.Filter(key -> key[2] == S, keys(map.single_variable_constraints))
    )
end
function number_of_type(map::Map, C::Type{MOI.ConstraintIndex{MOI.VectorOfVariables, S}}) where S
    return count(key -> key[2] == S, keys(map.vector_of_variables_constraints))
end
function keys_of_type(map::Map, C::Type{MOI.ConstraintIndex{MOI.VectorOfVariables, S}}) where S
    return MOIU.LazyMap{C}(
        key -> C(key[1]),
        Base.Iterators.Filter(key -> key[2] == S, keys(map.vector_of_variables_constraints))
    )
end

"""
    list_of_key_types(map::Map)

Return a list of all the different concrete type of keys in `map`.
"""
function list_of_key_types(map::Map)
    list = Set{Tuple{DataType, DataType}}()
    for i in eachindex(map.bridges)
        if map.bridges[i] !== nothing
            push!(list, map.constraint_types[i])
        end
    end
    for key in keys(map.single_variable_constraints)
        push!(list, (MOI.SingleVariable, key[2]))
    end
    for key in keys(map.vector_of_variables_constraints)
        push!(list, (MOI.VectorOfVariables, key[2]))
    end
    return list
end

"""
    variable_constraints(map::Map, vi::MOI.VariableIndex)

Return the list of all keys corresponding to [`MathOptInterface.SingleVariable`](@ref)
constraints on the variable `vi`.
"""
function variable_constraints(map::Map, vi::MOI.VariableIndex)
    cis = MOI.ConstraintIndex{MOI.SingleVariable}[]
    for key in keys(map.single_variable_constraints)
        if key[1] == vi.value
            push!(cis, MOI.ConstraintIndex{MOI.SingleVariable, key[2]}(vi.value))
        end
    end
    return cis
end

"""
    vector_of_variables_constraints(map::Map)

Return the list of all keys that correspond to
[`MathOptInterface.VectorOfVariables`](@ref) constraints.
"""
function vector_of_variables_constraints(map::Map)
    return MOIU.LazyMap{MOI.ConstraintIndex{MOI.VectorOfVariables}}(
        key -> MOI.ConstraintIndex{MOI.VectorOfVariables, key[2]}(key[1]),
        keys(map.vector_of_variables_constraints)
    )
end

"""
    has_bridges(map::Map)::Bool

Return a `Bool` indicating whether any bridge was added yet. Note that it
returns `false` even if all bridges were deleted while `isempty` would return
`true`. It is computed in `O(1)` while `isempty` needs `O(n)` hence it is used
by [`MathOptInterface.Bridges.AbstractBridgeOptimizer`](@ref) to shortcut
operations in case variable bridges are not used.
"""
function has_bridges(map::Map)
    return !isempty(map.bridges) ||
        !isempty(map.single_variable_constraints) ||
        !isempty(map.vector_of_variables_constraints)
end

"""
    add_key_for_bridge(map::Map, bridge::AbstractBridge,
                       func::MOI.AbstractFunction, set::MOI.AbstractSet)

Return a new constraint index `ci` and stores the mapping `ci => bridge`.
"""
function add_key_for_bridge end

function add_key_for_bridge(map::Map, bridge::AbstractBridge,
                            func::MOI.AbstractFunction, set::MOI.AbstractSet)
    push!(map.bridges, bridge)
    push!(map.constraint_types, (typeof(func), typeof(set)))
    return MOI.ConstraintIndex{typeof(func), typeof(set)}(length(map.bridges))
end
function add_key_for_bridge(map::Map, bridge::AbstractBridge,
                            func::MOI.SingleVariable, set::MOI.AbstractScalarSet)
    map.single_variable_constraints[(func.variable.value, typeof(set))] = bridge
    return MOI.ConstraintIndex{MOI.SingleVariable, typeof(set)}(func.variable.value)
end
function add_key_for_bridge(map::Map, bridge::AbstractBridge,
                            func::MOI.VectorOfVariables, set::MOI.AbstractVectorSet)
    index = first(func.variables).value
    map.vector_of_variables_constraints[(index, typeof(set))] = bridge
    return MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(index)
end

"""
    EmptyMap <: AbstractDict{MOI.ConstraintIndex, AbstractBridge}

Empty version of [`Map`](@ref). It is used by
[`MathOptInterface.Bridges.Variable.SingleBridgeOptimizer`](@ref) as it does
not bridge any constraint.
"""
struct EmptyMap <: AbstractDict{MOI.VariableIndex, AbstractBridge} end
Base.isempty(::EmptyMap) = true
function Base.empty!(::EmptyMap) end
Base.keys(::EmptyMap) = MOIU.EmptyVector{MOI.VariableIndex}()
Base.values(::EmptyMap) = MOIU.EmptyVector{AbstractBridge}()
has_bridges(::EmptyMap) = false
number_of_type(::EmptyMap, ::Type{<:MOI.ConstraintIndex}) = 0


text:
using Test

@testset "Decapodes" begin
  include("Decapodes.jl")
end

@testset "PetriNets" begin
  include("PetriNets.jl")
end

@testset "Examples" begin
  include("Examples.jl")
end


text:
stack = Int[]           # []
@show push!(stack, 1)   # [1]
@show push!(stack, 2)   # [1, 2]
@show push!(stack, 3)   # [1, 2, 3]
@show pop!(stack)       # 3
@show length(stack)     # 2
@show empty!(stack)     # []
@show isempty(stack)    # true

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo2Real(a))
end

@test begin
    a = 1
    !isreal(Exo2Real(a, a, a, a))
end

@test begin
    z = random(Exo2Real{Int})
    z == +(z)
end

@test begin
    a = 1
    real(Exo2Real(a)) == a
end

@test begin
    a = 1
    b = 2
    c = 3
    l = unreal(Exo2Real(0, a, b, c))
    r = [a, b, c]
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2Real(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, a: 2, b: 3, ab: 4]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = Exo2Real(Exo1Real(a, b))
    r = Exo2Real(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = Exo2Real(a, b, c)
    r = Exo2Real(a, b, c, 0)
    l == r
end

@test begin
    x = random(Exo2Real{BigFloat})
    y = random(Exo2Real{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo2Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Exo1Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo2Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Exo1Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = associator(x, y, z)
    r = zero(Exo2Real{BigInt})
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo2Real{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo2Real{Rational{BigInt}})
    y = random(Exo2Real{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo2Real{Rational{BigInt}})
    y = random(Exo2Real{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo2Real{Rational{BigInt}})
    y = random(Exo2Real{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo2Real{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Exo1Real{BigInt})
    x = random(Exo2Real{BigInt})
    l = a * x
    r = x * a
    l != r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    z = random(Exo2Real{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo2Real{BigInt})
    quadrance(x) > big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(Exo2Real{BigInt})
    y = random(Exo2Real{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo2Real(0, 0.0))
end

@test_throws ErrorException begin
    inv(Exo2Real(0, 1, 1, 1))
end

@test_throws ErrorException begin
    inv(Exo2Real(0, 0, 2, 1))
end

@test_throws ErrorException begin
    random(Exo2Real{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo2Real{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo2Real{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo2Real{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo2Real{BigInt})
    z == z + zero(Exo2Real{BigInt})
end

@test begin
    z = random(Exo2Real{BigInt})
    z == z * one(Exo2Real{BigInt})
end

@test begin
    z = random(Exo2Real{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo2Real{BigInt})
    z == z * one(z)
end

@test begin
    x = random(Exo2Real{BigInt})
    cloak(cloak(x)) == x
end

@test begin
    x = random(Exo2Real{BigInt})
    dagger(dagger(x)) == x
end

@test begin
    x = random(Exo2Real{BigInt})
    star(star(x)) == dagger(x)
end

@test begin
    x = random(Exo2Real{Rational{BigInt}})
    selfstar(selfstar(x)) == selfstar(x)
end

@test begin
    x = random(Exo2Real{Rational{BigInt}})
    anti_selfstar(anti_selfstar(x)) == anti_selfstar(x)
end


text:
# low-level wrappers of the CUFFT library

cufftGetVersion() = ccall((:cufftGetVersion,libcufft), Cint, ())

function cufftGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  @check ccall((:cufftGetProperty, libcufft), cufftStatus_t,
               (Cint, Ptr{Cint}),
               property, value_ref)
  value_ref[]
end

cufftDestroy(plan) = ccall((:cufftDestroy,libcufft), Nothing, (cufftHandle_t,), plan)

function cufftPlan1d(plan, nx, type, batch)
    @check ccall((:cufftPlan1d,libcufft),cufftStatus_t,
                 (Ptr{cufftHandle_t}, Cint, cufftType, Cint),
                 plan, nx, type, batch)
end

function cufftPlan2d(plan, nx, ny, type)
    @check ccall((:cufftPlan2d,libcufft),cufftStatus_t,
                 (Ptr{cufftHandle_t}, Cint, Cint, cufftType),
                 plan, nx, ny, type)
end

function cufftPlan3d(plan, nx, ny, nz, type)
    @check ccall((:cufftPlan3d,libcufft),cufftStatus_t,
                 (Ptr{cufftHandle_t}, Cint, Cint, Cint, cufftType),
                 plan, nx, ny, nz, type)
end

function cufftPlanMany(plan, rank, n, inembed, istride, idist, onembed, ostride, odist, type, batch)
    @check ccall((:cufftPlanMany,libcufft),cufftStatus_t,
                 (Ptr{cufftHandle_t}, Cint, Ptr{Cint},
                  Ptr{Cint}, Cint, Cint,
                  Ptr{Cint}, Cint, Cint,
                  cufftType, Cint),
                 plan, rank, n, inembed, istride, idist, onembed, ostride, odist, type, batch)
end

function cufftExecC2C(plan, idata, odata, direction)
    @check ccall((:cufftExecC2C,libcufft), cufftStatus_t,
                 (cufftHandle_t, CuPtr{cufftComplex}, CuPtr{cufftComplex}, Cint),
                 plan, idata, odata, direction)
end

function cufftExecC2R(plan, idata, odata)
    @check ccall((:cufftExecC2R,libcufft), cufftStatus_t,
                 (cufftHandle_t, CuPtr{cufftComplex}, CuPtr{cufftComplex}),
                 plan, idata, odata)
end

function cufftExecR2C(plan, idata, odata)
    @check ccall((:cufftExecR2C,libcufft), cufftStatus_t,
                 (cufftHandle_t, CuPtr{cufftReal}, CuPtr{cufftComplex}),
                 plan, idata, odata)
end

function cufftExecZ2Z(plan, idata, odata, direction)
    @check ccall((:cufftExecZ2Z,libcufft), cufftStatus_t,
                 (cufftHandle_t, CuPtr{cufftDoubleComplex}, CuPtr{cufftDoubleComplex},
                  Cint),
                 plan, idata, odata, direction)
end

function cufftExecZ2D(plan, idata, odata)
    @check ccall((:cufftExecZ2D,libcufft), cufftStatus_t,
                 (cufftHandle_t, CuPtr{cufftDoubleComplex}, CuPtr{cufftDoubleComplex}),
                 plan, idata, odata)
end

function cufftExecD2Z(plan, idata, odata)
    @check ccall((:cufftExecD2Z,libcufft), cufftStatus_t,
                 (cufftHandle_t, CuPtr{cufftDoubleReal}, CuPtr{cufftDoubleComplex}),
                 plan, idata, odata)
end


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-
cd(dirname(@__FILE__))
include("../HANTS.jl")
using .HANTS, Plots

y = [
    0.055952, 0.037081, 0.062657, 0.120110, 0.178219, 0.244443, 0.213190, 0.494648,
    0.328767, 0.561776, 0.380300, 0.641233, 0.532486, 0.827757, 0.654052, 0.816695,
    0.611424, 0.661557, 0.398067, 0.339881, 0.113957, 0.086790, 0.058600
]

nbase = 365
nfreq = 3
tseries = 1:16:365
validrange = (-1, 1)
fet = 0.02
dod = 1
δ = 0.1

A, φ, yrec = hants(
    y, fet, dod, δ; nbase=nbase, nfreq=nfreq, validrange=validrange, tseries=tseries, outlier=nothing
)
y_reconstruct = reconstruct(A, φ, nbase)

A_Hi, φ_Hi, yrec_Hi = hants(
    y, fet, dod, δ; nbase=nbase, nfreq=nfreq, validrange=validrange, tseries=tseries, outlier=:Hi
)
y_reconstruct_Hi = reconstruct(A_Hi, φ_Hi, nbase)

A_Lo, φ_Lo, yrec_Lo = hants(
    y, fet, dod, δ; nbase=nbase, nfreq=nfreq, validrange=validrange, tseries=tseries, outlier=:Lo
)
y_reconstruct_Lo = reconstruct(A_Lo, φ_Lo, nbase)

fig1 = plot(tseries, y, lc=:black, label="Original Data", shape=:circle, mc=:black)
plot!(tseries, yrec, lc=:blue, label="HANTS - None", shape=:circle, mc=:blue)
plot!(tseries, yrec_Hi, lc=:red, label="HANTS - Hi", shape=:circle, mc=:red)
plot!(tseries, yrec_Lo, lc=:green, label="HANTS - Lo", shape=:circle, mc=:green)
title!("Testing HANTS Algorithm")

fig2 = plot(tseries, y, lc=:black, label="Original Data", shape=:circle, mc=:black)
plot!(1:365, y_reconstruct, lc=:blue, label="HANTS - None")
plot!(1:365, y_reconstruct_Hi, lc=:red, label="HANTS - Hi")
plot!(1:365, y_reconstruct_Lo, lc=:green, label="HANTS - Lo")
title!("Testing HANTS Reconstruct")

plot(fig1, fig2; layout=(1, 2), size=(800,300))

max_stars_repo_name: shenrq/JuliaHANTS


text:
#####################################################################################
# Regions required to build the Hamiltonians in hamiltonian.jl
#####################################################################################

"New method to Quantica.combine that accepts Regions"
Quantica.combine(regtup::Tuple{Quantica.RegionPresets.Region,Quantica.RegionPresets.Region}) =
    Quantica.combine(regtup[1], regtup[2])
Quantica.combine(reg1::Quantica.RegionPresets.Region,
    reg2::Quantica.RegionPresets.Region) =
 Quantica.RegionPresets.Region{2}(r-> reg1.f(r) && reg2.f(r))

"Bravais vectors for the periodic hamiltonian"
function nanoribbonbravais(p)
    @unpack edges = p
    if edges == :armchair_zigzag
        br1 = (1, 1)
        br2 = (1, -1)
        br3 = (1, 0)
    elseif edges == :zigzag_armchair
        br1 = (1, 1)
        br2 = (1, -1)
        br3 = (0, 1)
    else
        nothing
    end
    return br1, br2, br3
end

"x and y boundaries for all types of edges :armchair_zigzag, :zigzag_armchair used in hamiltonians.jl"
function xyregs(p)
    @unpack a0, Ln, Ls, W, scale, edges = p
    ϵ = a0/1000
    a0 = a0 * scale
    Ln, Ls, W = edgecondition(Ln, Ls, W, a0, edges)
    if edges == :armchair_zigzag
        xreg = Quantica.RegionPresets.Region{2}(r -> -Ln/2 -Ls < r[1] < Ln/2 + Ls)
        yreg =  Quantica.RegionPresets.Region{2}(r -> ifelse(sign(r[3]) < 0,
            - ϵ <= r[2] <= W + a0/(2*√3)  + ϵ,
            -a0/(2*√3) - ϵ <= r[2] <= W + a0/(2*√3)  + ϵ))
    elseif edges == :zigzag_armchair
        xreg = Quantica.RegionPresets.Region{2}(r -> ifelse(sign(r[3]) < 0,
            -Ln/2 -Ls -a0/(2*√3) - ϵ <= r[1] <= Ln/2 + a0/(2*√3) + Ls + ϵ,
            -Ln/2 - Ls -a0/(2*√3) - ϵ <= r[1] <= Ln/2 + a0/(2*√3) + Ls + ϵ))
        yreg = Quantica.RegionPresets.Region{2}(r ->  0 < r[2] < W -a0/2 + ϵ)
    end
    return xreg, yreg
end

"regions used in model.jl"
function modelregs(p)
    @unpack a0, dinter, dintra, Ls, Ln, W, scale, edges = p

    a0 = a0 * scale
    Ln, Ls, W = edgecondition(Ln, Ls, W, a0, edges)
    W = ifelse(edges == :armchair_zigzag, W + a0/sqrt(3), W-a0/2)
    ϵ = a0*1e-3
    edgereg(r) = ifelse(edges == :zigzag_armchair,
    -Ln/2 -a0/(2*√3) - ϵ <= r[1] <=  -Ln/2 - ϵ  || Ln/2 <= r[1] <=  Ln/2 + a0/(2*√3) + ϵ,
    -Ln/2 < r[1] < -Ln/2 + 2a0/(2*√3) || Ln/2 - 2a0/(2*√3) < r[1] < Ln/2)

    Redgereg(r) = ifelse(edges == :zigzag_armchair,
    Ln/2 <= r[1] <=  Ln/2 + a0/(2*√3) + ϵ,
    Ln/2 - 2a0/(2*√3) < r[1] < Ln/2)

    normalreg(r) = ifelse(edges == :zigzag_armchair,
    -Ln/2 -a0/(2*√3) - ϵ <= r[1] <= Ln/2 + a0/(2*√3) + ϵ,
    -Ln/2 < r[1] < Ln/2)
    #- Ln/2 -a0/(2*sqrt(3)) <= r[1] <= Ln/2 +a0/(2*sqrt(3))
    scRreg(r) = r[1] > ifelse(edges == :zigzag_armchair, Ln/2 + a0/(2*√3) + ϵ, Ln/2)
    scLreg(r) = r[1] < ifelse(edges == :zigzag_armchair,-Ln/2 -a0/(2*√3) - ϵ, -Ln/2)
    scRreg2(r) = r[1] > ifelse(edges == :zigzag_armchair, Ln/2 + a0/(2*√3) + ϵ, Ln/2 - a0/(2*√3))
    scLreg2(r) = r[1] < ifelse(edges == :zigzag_armchair,-Ln/2 -a0/(2*√3) - ϵ, -Ln/2 + a0/(2*√3))
    screg(r) = scRreg(r) || scLreg(r)
    screg2(r) = scRreg2(r) || scLreg2(r)
    return a0, Ln, Ls, W, (edgereg, Redgereg, normalreg, scRreg, scLreg, screg, screg2)
end

"condition to prevent dangling bonds in both the zigzag and armchair direction"
edgecondition(Ln, Ls, W, a0, edges) = ifelse(edges == :zigzag_armchair,
    (zigzagcond(Ln, a0), zigzagcond(Ls, a0), armchaircond(W, a0)),
    (armchaircond(Ln, a0), armchaircond(Ls, a0), zigzagcond(W, a0)))
"Condition to prevent dangling bonds in the zigzag direction"
zigzagcond(length, a0) = 6a0/sqrt(3) * round((length)/(6a0/sqrt(3)))
"Condition to prevent dangling bonds in the armchair direction"
armchaircond(length, a0) = 4a0 * round(length/(4a0))


text:
using Sundials

## f routine. Compute function f(t,y).

function f(t, y, ydot)
    ydot[1] = -0.04 * y[1] + 1.0e4 * y[2] * y[3]
    ydot[3] = 3.0e7 * y[2] * y[2]
    ydot[2] = -ydot[1] - ydot[3]
    return Sundials.CV_SUCCESS
end

t = [0.0; 4 * exp10.(range(-1.0, stop = 7.0, length = 9))]
y0 = [1.0, 0.0, 0.0]
res = Sundials.cvode(f, y0, t)


text:
struct ZeroPotential <: Potential end

UU(phi,  ::ZeroPotential) = 0
UUp(phi, ::ZeroPotential) = 0

parameters(potential::ZeroPotential) = ()

"""
* `oophiM2` : 1/ϕM^2
* `oophiQ`  : 1/ϕQ

Note that, in the notation used in the paper XXXX.XXXX,
```math
λ₄ = 1/(4 ϕM^2), λ₆ = 1/ϕQ
```
"""
Base.@kwdef struct Phi8Potential{T} <: Potential
    oophiM2 :: T = 0.0
    oophiQ  :: T = 0.0
end

parameters(potential::Phi8Potential) = (oophiM2=potential.oophiM2, oophiQ=potential.oophiQ)

@doc raw"""
```math
U = -\frac{1}{3} + \frac{ϕ^2}{3 ϕM^2} + \frac{ϕ^2}{2 ϕM^4} - \frac{2 ϕ^2}{ϕQ}
 + \frac{4 ϕ^4}{3 ϕQ} + \frac{6 ϕ^4}{ϕM^2 ϕQ} - \frac{ϕ^4}{12 ϕM^4}
 + \frac{18 ϕ^6}{ϕQ^2} - \frac{2 ϕ^6}{3 ϕM^2 ϕQ}
 - \frac{4ϕ^8}{3 ϕQ^2}
```
"""
function UU(phi, potential::Phi8Potential)
    phi2 = phi  * phi
    phi4 = phi2 * phi2
    phi6 = phi2 * phi4
    phi8 = phi4 * phi4

    oophiM2 = potential.oophiM2
    oophiM4 = oophiM2 * oophiM2
    oophiQ  = potential.oophiQ
    oophiQ2 = oophiQ * oophiQ

    -1/3 + phi2 * (oophiM4 / 2 + oophiM2 / 3 - 2 * oophiQ) +
        phi4 * ( -oophiM4 / 12 + 4 * oophiQ / 3 + 6 * oophiM2 * oophiQ) +
        phi6 * (18 * oophiQ2 - 2 * oophiM2 * oophiQ / 3) -
        phi8 * 4 * oophiQ2 / 3
end

@doc raw"""
```math
Up = \frac{dU}{dϕ}
```
"""
function UUp(phi, potential::Phi8Potential)
    phi2 = phi  * phi
    phi3 = phi  * phi2
    phi5 = phi3 * phi2
    phi7 = phi5 * phi2

    oophiM2 = potential.oophiM2
    oophiM4 = oophiM2 * oophiM2
    oophiQ  = potential.oophiQ
    oophiQ2 = oophiQ * oophiQ

    phi * (oophiM4 + 2 * oophiM2 / 3 - 4 * oophiQ) +
        phi3 * ( -oophiM4 / 3 + 16 * oophiQ / 3 + 24 * oophiM2 * oophiQ) +
        phi5 * (108 * oophiQ2 - 4 * oophiM2 * oophiQ) -
        phi7 * 32 * oophiQ2 / 3
end

#New potential coming from a superpotential with extra term alpha*phi^8

Base.@kwdef struct PhiAlphaBetaPotential{T} <: Potential
    oophiM2 :: T   = 0.0
    oophiQ  :: T   = 0.0
    alpha   :: T   = 0.0
    beta    :: Int = 1
end

parameters(potential::PhiAlphaBetaPotential) = (oophiM2=potential.oophiM2, oophiQ=potential.oophiQ, alpha=potential.alpha)

function UU(phi, potential::PhiAlphaBetaPotential)
    phi2  = phi  * phi
    phi4  = phi2 * phi2
    phi6  = phi2 * phi4
    phi8  = phi4 * phi4

    oophiM2 = potential.oophiM2
    oophiM4 = oophiM2 * oophiM2
    oophiQ  = potential.oophiQ
    oophiQ2 = oophiQ * oophiQ
    alpha   = potential.alpha
    alpha2  = alpha * alpha
    beta    = potential.beta
    beta2   = beta * beta

    return -1/3+phi2*(1/2*oophiM4+1/3*oophiM2-2*oophiQ)+
            phi4*(4/3*oophiQ+6*oophiM2*oophiQ-1/12*oophiM4)+
            phi6*(18*oophiQ2-2/3*oophiM2*oophiQ)-4/3*phi8*oophiQ2-
            phi^beta*(2/3*alpha*oophiM2+6*alpha*beta*oophiQ)+
            +1/2*alpha2*beta2*phi^(2*beta-6)+
            phi^(beta-4)*(4*alpha-alpha*beta)+
            4/3*phi^(2*beta-4)*alpha2+
            phi^(beta-2)*(4/3*alpha+alpha*beta*oophiM2)-
            8/3*alpha*oophiQ*phi^(2+beta)

end

@doc raw"""
```math
Up = \frac{dU}{dϕ}
```
"""
function UUp(phi, potential::PhiAlphaBetaPotential)
    phi2  = phi  * phi
    phi3  = phi  * phi2
    phi5  = phi3 * phi2
    phi7  = phi5 * phi2

    oophiM2 = potential.oophiM2
    oophiM4 = oophiM2 * oophiM2
    oophiQ  = potential.oophiQ
    oophiQ2 = oophiQ * oophiQ
    alpha   = potential.alpha
    alpha2  = alpha * alpha
    beta    = potential.beta
    beta2   = beta * beta

    return 2*phi*(1/2*oophiM4+1/3*oophiM2-2*oophiQ)+
            4*phi3*(4/3*oophiQ+6*oophiM2*oophiQ-1/12*oophiM4)+
            6*phi5*(18*oophiQ2-2/3*oophiM2*oophiQ)-32/3*phi7*oophiQ2-
            beta*phi^(beta-1)*(2/3*alpha*oophiM2+6*alpha*beta*oophiQ)+
            +1/2*alpha2*beta2*(2*beta-6)*phi^(2*beta-7)+
            (beta-4)*phi^(beta-5)*(4*alpha-alpha*beta)+
            4/3*(2*beta-4)*phi^(2*beta-5)*alpha2+
            (beta-2)*phi^(beta-3)*(4/3*alpha+alpha*beta*oophiM2)-
            8/3*alpha*oophiQ*(2+beta)*phi^(1+beta)

end


text:
import Requests: get, json

base_url = "https://pypi.python.org/pypi/"
package_list_file = "packages.txt"

make_pypi_url = p::String -> string(base_url, p, "/json")

function parse_json(j::Dict)::Tuple{Int, Int, Int}
    dls = j["info"]["downloads"]
    dls["last_day"], dls["last_week"], dls["last_month"]
end

function get_stats()
    packages = strip.(read_package_list())
    @printf("%-15s: %7s %7s %7s\n", "package", "day", "week", "month")

    @sync begin         # will wait for all enclosed async tasks to complete
        for package in packages
            @async begin
                j = json(get(make_pypi_url(package)))
                d, w, m = parse_json(j)
                o = @sprintf("%-15s: %7d %7d %7d", package, d, w, m)
                println(o)          # have to sprintf to a string and then println to avoid intermingled output
            end
        end
    end
end

function read_package_list()::Array{String}
    try
        readlines(open(package_list_file))
    catch e
        println("error opening packages.txt")
        quit()
    end
end

function main()
    get_stats()
end

main()


text:
module Resolve

using DocStringExtensions

import ..ApFilterTrials
using ..Configuration
import ..Delayed
import ..GMM
import ..Graph
import ..IndexObjects
import ..LogGMM
import ..Resolve
import ..Samplers
import ..Trials

include("posterior.jl")
include("nodes.jl")

end # module Resolve


text:
baremodule StopWatches

export @stopwatch

using StopWatchesBase: @stopwatch, Time, append_times!, clear

function summary end
function time_dataframe end
function time_table end

module Implementations

using DataFrames: DataFrame
using StopWatchesBase: Time, append_times!, clear
using StructArrays: StructVector
using Tables: Tables
using Unitful: @u_str, uconvert

using ..StopWatches: StopWatches

include("utils.jl")
include("conversions.jl")

end  # module Implementations

Implementations.define_docstrings()

end  # baremodule StopWatches


text:
function Variance()
    return Variance(())
end

function Variance(arg0::SecondMoment)
    return Variance((SecondMoment,), arg0)
end

function Variance(arg0::Variance)
    return Variance((Variance,), arg0)
end

function Variance(arg0::jboolean)
    return Variance((jboolean,), arg0)
end

function Variance(arg0::jboolean, arg1::SecondMoment)
    return Variance((jboolean, SecondMoment), arg0, arg1)
end

function aggregate(obj::Variance, arg0::Variance)
    return jcall(obj, "aggregate", void, (Variance,), arg0)
end

function clear(obj::Variance)
    return jcall(obj, "clear", void, ())
end

function copy(obj::Variance)
    return jcall(obj, "copy", Variance, ())
end

function evaluate(obj::Variance, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function evaluate(obj::Variance, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble, arg3::jint, arg4::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jdouble, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function evaluate(obj::Variance, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jint, jint), arg0, arg1, arg2, arg3)
end

function evaluate(obj::Variance, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble), arg0, arg1)
end

function evaluate(obj::Variance, arg0::Vector{jdouble}, arg1::jdouble, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function evaluate(obj::Variance, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::Variance)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::Variance)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::Variance, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end

function is_bias_corrected(obj::Variance)
    return jcall(obj, "isBiasCorrected", jboolean, ())
end

function with_bias_correction(obj::Variance, arg0::jboolean)
    return jcall(obj, "withBiasCorrection", Variance, (jboolean,), arg0)
end


text:
export potential,CorrugatedMorse,calc_VG

using SpecialFunctions

struct CorrugatedMorse{T <: Abstract2DLattice}
    D
    κ
    h
    lattice::T
end

function get_w(z,ch::ScatteringChannels,ki,conf::CorrugatedMorse)
    get_V(z,ch,conf) + Diagonal(ch.kz²)
end

function get_V(z,ch::ScatteringChannels,conf::CorrugatedMorse)
    @unpack gi = ch
    l = length(gi)
    V = Array{Float64,2}(undef,l,l)
    for i in 1:l, j in 1:l
        V[i,j] = get_VG(z,gi[i]-gi[j],conf)
    end
    V
end

function get_VG(z,g,conf::CorrugatedMorse{HexGrid})
    @unpack κ,h,D = conf
    if g[1]==0 && g[2]==0
        return D*(exp(-2*κ*z)+2*exp(-κ*z))
    end

    V0 = calc_VG_factor([0,0],conf,1)
    VG = calc_VG_factor(g,conf,V0)

    VG *= exp(-2*κ*z)*D

    return VG
end

function calc_VG_factor(g,conf::CorrugatedMorse{HexGrid},div)
    @unpack κ,h,D = conf
    @unpack a = conf.lattice

    kmax = 10
    α = 2*κ*h*a/3
    VG = 0.0
    for k in -kmax:kmax
        VG += besseli(k,α)*(besseli(k+g[2],α)*besseli(k+g[1],α) + besseli(k-g[2],α)*besseli(k-g[1],α))
    end
    VG / div
end

function calc_VG(g,l::Generic2DLattice,conf::CorrugatedMorse,z,steps=100)
    VG = 0+0im
    rng = range(0.0,length=steps,1)
    for i in rng , j in rng
        r1 = l.a1.*i
        r2 = l.a2.*j
        r = r1+r2
        VG += potential(conf,r[1],r[2],z)*exp(-1im*(g'*r))
    end
    VG *= 2π/(steps^2)
end

function potential(conf::CorrugatedMorse,x,y,z)
    @unpack D,κ,h,lattice = conf
    @unpack a1,a2 = lattice
    R = [x,y]

    # TODO: only for 60 deg(a1,a2)
    ξ = h*(cos(2π*(a1'*R)/norm(a1)^2)+
           cos(2π*(a2'*R)/norm(a2)^2)+
           cos(2π*((a1-a2)'*R)/norm(a1-a2)^2))
    @. D*(exp(-2*(z-ξ)/κ)-2*exp(-z/κ))

end


text:
# A hypothetical env where consecutive days off are DEFINITELY possible
test_consecutive_env = {
    "coverage"  => Array[
        [2*ones(Int,8)],
        [2*ones(Int,8)],
        [2*ones(Int,8)],
        [2*ones(Int,8)],
    ],
}

test_consecutive_employees = {
    "a" => {
        "hours_min" => 15,
        "hours_max" => 17,
        "shift_count_min" => 2,
        "shift_count_max" => 2,
        "availability"  => Array[
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
        ],
        "longest_availability" => 8*ones(Int, 8),
    },
    "b" => {
        "hours_min" => 15,
        "hours_max" => 17,
        "shift_count_min" => 2,
        "shift_count_max" => 2,
        "availability"  => Array[
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
        ],
        "longest_availability" => 8*ones(Int, 8),
    },
    "c" => {
        "hours_min" => 15,
        "hours_max" => 17,
        "shift_count_min" => 2,
        "shift_count_max" => 2,
        "availability"  => Array[
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
        ],
        "longest_availability" => 8*ones(Int, 8),
    },
    "d" => {
        "hours_min" => 15,
        "hours_max" => 17,
        "shift_count_min" => 2,
        "shift_count_max" => 2,
        "availability"  => Array[
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
            [ones(Int,8)],
        ],
        "longest_availability" => 8*ones(Int, 8),
    },
}

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
import Geometry

"""
    struct JacobianMatrix{T} <: AbstractMatrix{T}
Jacobian matrix. Block-diagonal local matrix represents intraelement operations,
off-diagonal represents connectivity between elements.

`diag[:,iS1,:,iS2,iK]` is a square matrix d(f(u[:,iS1,iK))/d(u[:,iS2,iK])
for states iS1, iS2, element iK. Designed so that `diag[:,:,:,:,iK]` can be
in-place reshaped for matrix multiplication.

`offdiag[iFN2,iS2,iN1,iS1,iF,iK1]` is a rectangular matrix
d(f(u[iN1,iS1,iK1])/d(u[iFN2,iS2,mesh.e2e[iF,iK1]]), or effect of neighbor iF on
element iK1. Again, can in-place reshape for matrix multiplication. Note that
for upper triangular portion of global matrix (`mesh.e2e[iF,iK1]<iK1`),
the matrix `offdiag[:,:,:,:,iF,iK1]` actually yields the transpose of
operator needed. That is, the lower triangular part represents the effect of
neighbor on face nodes of iK1, and the upper triangular part represents the
transpose of the effect of neighbor's face nodes on iK1. This format therefore
only needs one array to store all off-diagonals.
"""
struct JacobianMatrix{T} <: AbstractMatrix{T}
    diag::Array{T,5}
    offdiag::Array{T,6}
    mesh::Geometry.Mesh

    function JacobianMatrix(diag::Array{T,5}, offdiag::Array{T,6}, mesh::Geometry.Mesh) where {T}
        @assert size(diag,1) == size(diag,3) "Diagonal of Jacobian must be square!"
        @assert size(diag,2) == size(diag,4) "Number of Jacobian states must be consistent!"
        @assert size(offdiag,5) == Geometry.N_FACES "Number of Jacobian neighbors does not match geometry!"
        new{T}(diag,offdiag,mesh)
    end
end

"""
    size(A::JacobianMatrix)
Returns a tuple containing the dimensions of actual operator A
"""
function Base.size(A::JacobianMatrix)
    (m,ns,n,ns2,ne) = size(A.diag)
    return (m*ns*ne, n*ns2*ne)
end
"""
    getindex(A::JacobianMatrix, i::Integer)
Linear scalar indexing into Jacobian matrix
"""
@inline function Base.getindex(A::JacobianMatrix, i::Integer)
    (sz1,sz2) = size(A)
    j = rem1(i, sz1); i = div(i, sz2, RoundUp)
    return getindex(A,i,j)
end
"""
    getindex(A::JacobianMatrix{T}, i::Integer, j::Integer)
2-dimensional scalar indexing into JacobianMatrix. Not efficient.
"""
@inline function Base.getindex(A::JacobianMatrix{T}, i::Integer, j::Integer) where {T}
    @boundscheck checkbounds(Bool, A,i,j)
    (m,ns,n,ns2,ne) = size(A.diag)
    im = rem1(i ,m ); i2 = div(i ,m , RoundUp)
    is = rem1(i2,ns); i3 = div(i2,ns, RoundUp)
    ie = rem1(i3,ne)
    jn = rem1(j ,n ); j2 = div(j ,n , RoundUp)
    js = rem1(j2,ns); j3 = div(j2,ns, RoundUp)
    je = rem1(j3,ne)
    if ie==je
        # On diagonal
        if checkbounds(Bool,A.diag, im,is,jn,js,ie)
            return @inbounds A.diag[im,is,jn,js,ie]
        end
    elseif je ∈ A.mesh.e2e[:,ie]
        iF = findfirst(x->x==je,A.mesh.e2e[:,ie])
        if je < ie
            # Lower triangular part stored directly
            iFN2 = findfirst(x->x==jn,A.mesh.ef2n[:,iF])
            iFN = 0
            if !isnothing(iFN2)
                # We need the original face node
                iFN = mesh.nfn2fn[iFN2]
            else
                # Node is not on the face => no Jacobian dependence
                return zero(T)
            end
            return A.offdiag[iFN,is,jn,js,iF,ie]
        else
            # Upper triangular part stored transposed
            iFN = findfirst(x->x==im,A.mesh.ef2n[:,iF])
            iFN2 = 0
            if !isnothing(iFN)
                # We need the neighbor's face node
                iFN2 = mesh.nfn2fn[iFN]
            else
                # Node is not on the face => no Jacobian dependence
                return zero(T)
            end
            return A.offdiag[iFN2,js,im,is,iF,ie]
        end
    else
        # Non-neighboring elements => no Jacobian dependence
        return zero(T)
    end
end
"""
    setindex!(A::JacobianMatrix, v, i::Integer)
Linear scalar indexed assignment of `JacobianMatrix`
"""
@inline function Base.setindex!(A::JacobianMatrix, v, i::Integer)
    (sz1,sz2) = size(A)
    j = rem1(i, sz1); i = div(i, sz2, RoundUp)
    setindex!(A, v, i, j)
end
"""
    setindex!(A::JacobianMatrix{T}, v, i::Integer, j::Integer)
2-dimensional scalar assignment of `JacobianMatrix`. Not efficient.
"""
@inline function Base.setindex!(A::JacobianMatrix{T}, v, i::Integer, j::Integer) where {T}
    @boundscheck checkbounds(Bool, A,i,j)
    (m,ns,n,ns2,ne) = size(A.diag)
    im = rem1(i ,m ); i2 = div(i ,m , RoundUp)
    is = rem1(i2,ns); i3 = div(i2,ns, RoundUp)
    ie = rem1(i3,ne)
    jn = rem1(j ,n ); j2 = div(j ,n , RoundUp)
    js = rem1(j2,ns); j3 = div(j2,ns, RoundUp)
    je = rem1(j3,ne)
    if ie==je
        # On diagonal
        if checkbounds(Bool,A.diag, im,is,jn,js,ie)
            @inbounds setindex!(A.diag, v, im,is,jn,js,ie)
        end
    elseif je ∈ A.mesh.e2e[:,ie]
        iF = findfirst(x->x==je,A.mesh.e2e[:,ie])
        if je < ie
            # Lower triangular part stored directly
            iFN2 = findfirst(x->x==jn,A.mesh.ef2n[:,iF])
            iFN = 0
            if !isnothing(iFN2)
                # We need the original face node
                iFN = mesh.nfn2fn[iFN2]
            else
                if !iszero(v)
                    # Node is not on the face => no Jacobian dependence
                    throw(ArgumentError("Cannot set off-diagonal entry ($i,$j) to non-zero value $v"))
                else
                    return v
                end
            end
            return setindex!(A.offdiag, v, iFN,is,jn,js,iF,ie)
        else
            # Upper triangular part stored transposed
            iFN = findfirst(x->x==im,A.mesh.ef2n[:,iF])
            iFN2 = 0
            if !isnothing(iFN)
                # We need the neighbor's face node
                iFN2 = mesh.nfn2fn[iFN]
            else
                if !iszero(v)
                    # Node is not on the face => no Jacobian dependence
                    throw(ArgumentError("Cannot set off-diagonal entry ($i,$j) to non-zero value $v"))
                else
                    return v
                end
            end
            return setindex!(A.offdiag, v, iFN2,js,im,is,iF,ie)
        end
    else
        # Non-neighboring elements => no Jacobian dependence
        if !iszero(v)
            throw(ArgumentError("Cannot set off-diagonal entry ($i,$j) to non-zero value $v"))
        else
            return v
        end
    end
end

(==)(A::JacobianMatrix, B::JacobianMatrix) = (A.diag == B.diag && A.offdiag == B.offdiag)
(-)(A::JacobianMatrix) = JacobianMatrix(-A.diag, -A.offdiag, A.mesh)
function (+)(A::JacobianMatrix, B::JacobianMatrix)
    @assert A.mesh === B.mesh
    JacobianMatrix(A.diag+B.diag, A.offdiag+B.offdiag, A.mesh)
end
function (-)(A::JacobianMatrix, B::JacobianMatrix)
    @assert A.mesh === B.mesh
    JacobianMatrix(A.diag-B.diag,A.offdiag-B.offdiag,mesh)
end
# TODO: Matrix addition with LocalMatrix as well?
(*)(c::Number, A::JacobianMatrix) = JacobianMatrix(c*A.diag, c*A.offdiag, A.mesh)
(*)(A::JacobianMatrix, c::Number) = JacobianMatrix(c*A.diag, c*A.offdiag, A.mesh)
(/)(A::JacobianMatrix, c::Number) = JacobianMatrix(A.diag / c, A.offdiag / c, A.mesh)
(*)(A::JacobianMatrix, B::JacobianMatrix) = error("Too expensive an operation!")

" Store `faces` on `volume` array according to given mask: vol[mask] += face "
function distribute_face_to_nodes!(volume, faces, mask)
    volume_shaped = reshape(volume, size(volume,1), prod(size(volume)[2:end])))
    faces_shaped = reshape(faces, size(faces,1), prod(size(faces)[2:end]))
    for j = 1:size(volume_shaped,2)
        for i = 1:size(mask,1)
            volume_shaped[mask[i],j] += faces_shaped[i,j]
        end
    end
end
" Collect `volume` into `faces` array according to given mask: face = vol[mask] "
function collect_face_from_nodes!(faces, volume, mask)
    volume_shaped = reshape(volume, size(volume,1), prod(size(volume)[2:end])))
    faces_shaped = reshape(faces, size(faces,1), prod(size(faces)[2:end]))
    for j = 1:size(volume_shaped,2)
        for i = 1:size(mask,1)
            faces_shaped[i,j] = volume_shaped[mask[i],j]
        end
    end
end

function (*)(A::JacobianMatrix{T1}, x::AbstractVector{T2}) where {T1, T2}
    (m,ns,n,ns2,ne) = size(A.diag)
    b = zeros(typeof(zero(T1)*zero(T2)), m*ns*ne)
    b_shaped = reshape(b, m,ns,ne)
    b_mult   = reshape(b, m*ns,ne)
    x_shaped = reshape(x, n,ns2,ne)
    x_mult   = reshape(x, n*ns2,ne)
    faces_mult = zeros(typeof(zero(T1)*zero(T2)), A.mesh.n_face_nodes*ns)
    faces_shaped = reshape(faces, A.mesh.n_face_nodes, ns)
    for iK = 1:ne
        # Diagonal of A
        @views Adiag = reshape(A.diag[:,:,:,:,iK], m*ns, n*ns2)
        @views b_mult[:,iK] .+= Adiag*x_mult[:,iK]
        # Off-diagonal of A
        for iF = 1:Geometry.N_FACES
            nK = A.mesh.e2e[iF,iK]
            if nK > 0
                Aoff = reshape(@view A.offdiag[:,:,:,:,iF,iK], A.mesh.n_faces_nodes*ns,m*ns)
                nF = A.mesh.e2f[iF,iK]
                if nK < iK
                    # Lower triangular part stored directly
                    @views faces_mult .= Aoff*x_mult[:,iK]
                    @views distribute_face_to_nodes!(b_shaped[:,:,iK], faces_shaped, A.mesh.ef2n[:,iF])
                else
                    # Upper triangular part stored transposed
                    @views collect_face_from_nodes!(faces_shaped, x_shaped[:,:,nK], A.mesh.ef2n[:,nF])
                    @views b_mult[:,iK] .+= Aoff'*faces_mult
                end
            end
        end
    end
    return b
end
function (*)(A::JacobianMatrix{T1}, x::SolutionVector{T2}) where {T1, T2}
    (m,ns,n,ns2,ne) = size(A.diag)
    b_shaped = zeros(typeof(zero(T1)*zero(T2)), m,ns,ne)
    b_mult = reshape(b, m*ns,ne)
    x_mult = reshape(x.data, n*ns2,ne)
    faces_mult = zeros(typeof(zero(T1)*zero(T2)), A.mesh.n_face_nodes*ns)
    faces_shaped = reshape(faces, A.mesh.n_face_nodes, ns)
    @views for iK = 1:ne
        # Diagonal of A
        Adiag = reshape(A.diag[:,:,:,:,iK], m*ns, n*ns2)
        b_mult[:,iK] .+= Adiag*x_mult[:,iK]
        # Off-diagonal of A
        for iF = 1:Geometry.N_FACES
            nK = A.mesh.e2e[iF,iK]
            if nK > 0
                Aoff = reshape(A.offdiag[:,:,:,:,iF,iK], A.mesh.n_faces_nodes*ns,m*ns)
                nF = A.mesh.e2f[iF,iK]
                if nK < iK
                    # Lower triangular part stored directly
                    faces_mult .= Aoff*x_mult[:,iK]
                    distribute_face_to_nodes!(b_shaped[:,:,iK], faces_shaped, A.mesh.ef2n[:,iF])
                else
                    # Upper triangular part stored transposed
                    collect_face_from_nodes!(faces_shaped, x.data[:,:,nK], A.mesh.ef2n[:,nF])
                    b_mult[:,iK] .+= Aoff'*faces_mult
                end
            end
        end
    end
    return SolutionVector(b_shaped)
end
"""
    mul!(y::SolutionVector, A::JacobianMatrix, x::SolutionVector, α::Number, β::Number)
Return A x α + y β by overwriting y
"""
function LinearAlgebra.mul!(y::SolutionVector{T1}, A::JacobianMatrix, x::SolutionVector, α::Number, β::Number) where {T1}
    (m,ns,n,ns2,ne) = size(A.diag)
    (nx,ns2x,nex) = true_size(x)
    (my,nsy,ney) = true_size(y)
    @assert m == my && ns  == nsy  && ne == ney "JacobianMatrix $(size(A.diag)), SolutionVector $(true_size(y)) mismatched size"
    @assert n == nx && ns2 == ns2x && ne == nex "JacobianMatrix $(size(A.diag)), SolutionVector $(true_size(x)) mismatched size"
    y_shaped = reshape(y.data, m*ns*ne)
    y_mult   = reshape(y.data, m*ns,ne)
    x_shaped = reshape(x.data, n,ns2,ne)
    x_mult   = reshape(x.data, n*ns2,ne)
    faces_mult = zeros(T1, A.mesh.n_face_nodes*ns)
    faces_shaped = reshape(faces, A.mesh.n_face_nodes, ns)
    @views for iK = 1:ne
        # Diagonal of A
        Adiag = reshape(A.diag[:,:,:,:,iK], m*ns, n*ns2)
        # y[iK] = Adiag*x[iK]*α + y[iK]*β
        mul!(y_mult[:,iK], Adiag, x_mult[:,iK], α, β)
        # Off-diagonal of A
        for iF = 1:Geometry.N_FACES
            nK = A.mesh.e2e[iF,iK]
            if nK > 0
                Aoff = reshape(A.offdiag[:,:,:,:,iF,iK], A.mesh.n_faces_nodes*ns,m*ns)
                nF = A.mesh.e2f[iF,iK]
                if nK < iK
                    # Lower triangular part stored directly
                    faces_mult .= Aoff*x_mult[:,iK]
                    # y[iK] += Aoff*faces[iK]*α
                    distribute_face_to_nodes!(y_shaped[:,:,iK], faces_shaped, A.mesh.ef2n[:,iF])
                else
                    # Upper triangular part stored transposed
                    collect_face_from_nodes!(faces_shaped, x_shaped[:,:,nK], A.mesh.ef2n[:,nF])
                    # y[iK] += Aoff'*faces[nK]*α
                    mul!(y_mult[:,iK], Aoff', faces_mult, α, 1.0)
                end
            end
        end
    end
end

"""
    ldiv!(A::JacobianMatrix, x::SolutionVector)
In-place linear solve A\\x using numerical iterative solver based on
preconditioned GMRES.
TODO: Decide actual algorithm and create production-level version of code.
"""
function LinearAlgebra.ldiv!(A::JacobianMatrix, x::SolutionVector)
    println("my ldiv!")
    (m,ns,n,ns2,ne) = size(A.diag)
    (nx,ns2x,nex) = true_size(x)
    @assert nx == m == n && ns2x == ns2 == ns && ne == nex "JacobianMatrix A $(size(A.diag)) cannot \\ x $(true_size(x))"
    # TODO: program the solver here
    error("Not yet programmed!")
end
function \(A::LocalMatrix, x::SolutionVector)
    println("my \\")
    (m,ns,n,ns2,ne) = size(A.diag)
    (nx,ns2x,nex) = true_size(x)
    @assert nx == m == n && ns2x == ns2 == ns && ne == nex "JacobianMatrix A $(size(A.diag)) cannot \\ x $(true_size(x))"
    b = deepcopy(x)
    LinearAlgebra.ldiv!(A, b)
end

function L2_error!(u::SolutionVector, u_true::SolutionVector, M::JacobianMatrix)
    u_true -= u
    return L2_norm(u_true, M)
end
function L2_norm(u::SolutionVector, M::JacobianMatrix)
    Mu = M*u
    (n,ns,ne) = true_size(Mu)
    l2_norms = zeros(ns)
    for iK = 1:ne
        for iS = 1:ns
            @views l2_norms[iS] += u.data[:,iS,iK]'*Mu.data[:,iS,iK]
        end
    end
    return l2_norms
end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
using JuMP: AbstractJuMPScalar

struct ContactResult{T, M}
    β::Vector{T}
    λ::T
    c_n::T
    # Δr::FreeVector3D{SVector{3, T}}
    point::Point3D{SVector{3, M}}
    obs::Obstacle{M}
    scaling::M
end

Base.show(io::IO, cr::ContactResult) = print(io, """ContactResult(β=$(cr.β), λ=$(cr.λ), c_n=$(cr.c_n), ...)""")

struct JointLimitResult{T, M}
    λ::Vector{T}
    scaling::M
end

struct LCPUpdate{T, M, U, Tdt}
    Δt::Tdt
    state::StateRecord{T, M}
    input::Vector{U}
    contacts::Dict{RigidBody{M}, Vector{ContactResult{T, M}}}
    joint_contacts::Vector{JointLimitResult{T, M}}
end

LCPUpdate(Δt::Tdt, state::StateRecord{T}, input::AbstractVector{U}, contacts::Associative{<:RigidBody{M}}, joint_contacts::AbstractVector) where {T, M, U, Tdt} =
    LCPUpdate{T, M, U, Tdt}(Δt, state, input, contacts, joint_contacts)

_getvalue(x::Variable) = JuMP._getValue(x)
_getvalue(x::Number) = x
_getvalue(f::FreeVector3D) = FreeVector3D(f.frame, _getvalue.(f.v))
function _getvalue(a::AffExpr)
    ret = a.constant
    for it in 1:length(a.vars)
        ret += a.coeffs[it] * _getvalue(a.vars[it])
    end
    ret
end
function _getvalue(a::QuadExpr)
    ret = _getvalue(a.aff)
    for it in 1:length(a.qvars1)
        ret += a.qcoeffs[it] * _getvalue(a.qvars1[it]) * _getvalue(a.qvars2[it])
    end
    return ret
end

JuMP.getvalue(r::JointLimitResult) = JointLimitResult(_getvalue.(r.λ), r.scaling)
JuMP.getvalue(r::StateRecord{<:AbstractJuMPScalar}) = StateRecord(r.mechanism, _getvalue.(r.state))
JuMP.getvalue(up::LCPUpdate) =
    LCPUpdate(_getvalue(up.Δt),
              getvalue(up.state),
              _getvalue.(up.input),
              Dict([k => getvalue.(v) for (k, v) in up.contacts]),
              getvalue.(up.joint_contacts))
JuMP.getvalue(c::ContactResult) = ContactResult(_getvalue.(c.β),
                                                _getvalue(c.λ),
                                                _getvalue(c.c_n),
                                                # _getvalue(c.Δr),
                                                c.point,
                                                c.obs,
                                                c.scaling)

function JuMP.setvalue(r::JointLimitResult{<:AbstractJuMPScalar}, seed::JointLimitResult{<:Number})
    setvalue.(r.λ, seed.λ)
end

function JuMP.setvalue(up::LCPUpdate{<:AbstractJuMPScalar}, seed::LCPUpdate{<:Number})
    if typeof(up.Δt) <: AbstractJuMPScalar
        setvalue(up.Δt, seed.Δt)
    end
    setvalue(up.state, seed.state)
    if eltype(up.input) <: AbstractJuMPScalar
        for i in eachindex(up.input)
            var = up.input[i]
            if var.m.colCat[var.col] != :Fixed
                setvalue(var, seed.input[i])
            end
        end
    end
    setvalue.(up.contacts, seed.contacts)
    setvalue.(up.joint_contacts, seed.joint_contacts)
end

JuMP.setvalue(f::FreeVector3D{<:AbstractVector{S}}, seed::FreeVector3D{<:AbstractVector{N}}) where {S <: AbstractJuMPScalar, N <: Number} =
    setvalue.(f.v, seed.v)

function JuMP.setvalue(contact::ContactResult{<:AbstractJuMPScalar}, seed::ContactResult{<:Number})
    @assert contact.obs == seed.obs
    setvalue(contact.β, seed.β)
    setvalue(contact.λ, seed.λ)
    setvalue(contact.c_n, seed.c_n)
    # setvalue(contact.Δr, seed.Δr)
end

function JuMP.setvalue(contacts::Dict{<:RigidBody, <:Vector{<:ContactResult}}, seeds::Dict)
    for body in keys(contacts)
        setvalue.(contacts[body], seeds[body])
    end
end

function JuMP.setvalue(state::StateRecord{<:AbstractJuMPScalar}, seed::StateRecord{<:Number})
    setvalue.(state.state, seed.state)
end


text:
using VertexSafeGraphs

using Test
import Random

import LightGraphs
const LG = LightGraphs

@testset "Graph construction and basic interface" begin
    nv = 20
    g1 = VSafeGraph(nv)
    @test LG.nv(g1) == nv
    @test LG.nv(g1.g) == nv

    g2_inner = LG.CompleteGraph(nv)
    g2 = VSafeGraph(g2_inner)
    @test LG.nv(g2) == LG.nv(g2_inner)
    @test LG.ne(g2) == LG.ne(g2_inner)

    @test all(sort(collect(LG.vertices(g2))) .== sort(collect(LG.vertices(g2_inner))))

    g3 = VSafeGraph(LG.CompleteDiGraph(30))
    @test LG.is_directed(g3)
    @test !LG.is_directed(g2)
end

@testset "Vertex deletion" begin
    Random.seed!(33)
    nv = 45
    inner = LG.CompleteGraph(nv)
    g = VSafeGraph(inner)
    @test LG.ne(inner) == LG.ne(g)
    @test LG.nv(inner) == LG.nv(g)
    nrm = 0
    for _ in 1:15
        removed_ok = LG.rem_vertex!(g, rand(1:nv))
        if !removed_ok
            continue
        end
        nrm += 1
        @test LG.nv(inner) == nv
        @test LG.nv(g) == nv - nrm
        @test length(g.deleted_vertices) == nrm

        @test LG.ne(inner) == LG.ne(g)
    end
end

@testset "Generic LG algorithms work" begin
    nv = 45
    inner = LG.CompleteGraph(nv)
    g = VSafeGraph(inner)
    removed_ok = LG.rem_vertex!(g, rand(1:nv))
    @test removed_ok
    # LG broken here
    @test_throws BoundsError LG.pagerank(g)
    @test_throws BoundsError LG.kruskal_mst(g)
end


text:
module TestPhylogenetics
using Test

using Phylo
using Diversity

@testset "Phylo" begin
    species = ["Dog", "Human", "Cat"]
    abund = [0.4, 0.3, 0.3]
    nt = RootedTree(species)
    n = createnode!(nt)
    createbranch!(nt, n, species[1], 1.0)
    createbranch!(nt, n, species[2], 1.0)
    r = createnode!(nt)
    createbranch!(nt, r, n, 1.0)
    createbranch!(nt, r, species[3], 2.0)
    ph = PhyloBranches(nt)
    leafnames = gettypenames(ph, true)
    @test species == gettypenames(ph, true)
    metaphylo = Metacommunity(abund, ph)
    @test gettypenames(metaphylo, true) == species
    @test getabundance(metaphylo, true) ≈ abund
    @test getabundance(metaphylo) ≈ [0.2, 0.2, 0.15, 0.15, 0.3]
    @test getordinariness!(metaphylo) ≈ [0.4, 0.7, 0.3, 0.7, 0.3]
    @test meta_gamma(metaphylo, 0)[!,:diversity] == [2.5]
    @test sub_gamma(metaphylo, 0)[!,:diversity] == [2.5]
end

end


text:
using Stochastics
using Test

@testset "Distributions.jl" begin
    X = Stochastics.NormalDistribution(0, 1)
    @test abs(Stochastics.normal_cumulative_distribution_function(X, 1.96) - 0.975) < 1e-4
    @test abs(Stochastics.gamma_function(1) - 1) < 1e-4
    @test abs(Stochastics.gamma_function(1/2) - √(π)) < 1e-4
    @test_throws ErrorException("α must be in ℝ > 0") Stochastics.gamma_function(0)
    X_chi = Stochastics.ChiSquaredDistribution(5)
    @test abs(Stochastics.chi_squared_cumulative_distribution_function(X_chi, 11.070) - 0.95) < 1e-4
    @test abs(Stochastics.quantile_finder(X_chi, 0.9995) - 0.1581) < 1e-4
    @test abs(Stochastics.quantile_finder(X, 0.025) - 1.96) < 1e-4
    X_t = Stochastics.tDistribution(1)
    @test abs(Stochastics.quantile_finder(X_t, 0.1) - 3.0777) < 1e-4
end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
function insertat(b::Buffer, c::Char, pos)
    pr = b.text[pos[1]]
    b.text[pos[1]] = string(unirange(pr, 1:pos[2]-1),
                            c,
                            unirange(pr, pos[2]:length(pr)),
                           )
end

function deleteback(b::Buffer)
    b.cursor.pos[2] -= 1
    if b.cursor.pos[2] >= 1
        deleteat(b, pos(b), 1)
    elseif y(b) > 1
        b.cursor.pos[1] -= 1
        b.cursor.pos[2] = width(b) + 1
        join_lines(b, (y(b), y(b) + 1), "")
    end
end
function splitp(b::Buffer)
    splitlines(b, b.cursor.pos)
    b.cursor.pos[1] += 1
    b.cursor.pos[2] = 1
end
function escape_insert(b)
    b.cursor.pos[2] -= 1
    escape(b)
end
function insert_fn(c)
    function insert(b::Buffer)
        if c=='\e'
            escape_insert(b)
        elseif c=='\x7f'
            deleteback(b)
            clamp!(b, true)
        elseif c=='\r'
            splitp(b)
        else
            insertat(b, c, pos(b))
            b.cursor.pos[2] += 1
        end
    end
end

insert_mode = Mode("insert", Action(insert_fn, x->true))