text:
# Space filling designs as in
# Surrogate Modelling of Computer Experiments with Sequential Experimental Design.", Crombecq, 2011

@doc """
    MonteCarloThDesign(; dims :: Int ) <: PointIterator

Return an iterator that provides samples within a hyper-
rectangle in ``ℝ^n`` where ``n`` is the dimension provided
to the iterator as `dims`.

Optional keyword arguments are:

* `n_points=100*dims`
    The maximum number of points to generate.
* `lb=zeros(dims)`
    Lower bounds of hyper-rectangle. Must all be finite.
* `ub=ones(dims)`
    Upper bounds of hyper-rectangle. Must all be finite.
* `seeds = []`.
    An initial set of points to add samples to. Is copied and internally modified.
* `clean_seeds=true`
    Throw away seeds that violate the box constraints.
* `spawn_factor = 100`
    In each iteration `dims*spawn_factor` random points are generated.
* `max_rand_points`
    Upper bound on the number of random points.
* `p_dist_th_factor = 0.5`
    The samples must have a projected distance exceeding `2*p_dist_th_factor/N`.
"""
@with_kw struct MonteCarloThDesign{
    VT <: AbstractVector{<:AbstractFloat},
    ST <: AbstractVector{<:VT}
    } <: PointIterator

    dims :: Int # only mandatory argument

    n_points :: Union{Nothing,Int} = 100 * dims

    lb :: VT = zeros(dims)
    ub :: VT = ones(dims)

    seeds :: ST = Vector{ eltype(lb) }[]
    clean_seeds :: Bool = true

    spawn_factor :: Int = 100
    max_rand_points :: Int = typemax(Int)

    p_dist_th_factor :: Float64 = .5

    needs_scaling :: Bool = !(iszero(lb) && all(ub .== 1))

    function MonteCarloThDesign{VT,ST}( dims :: Int, n_points :: Int, lb :: VT, ub :: VT, seeds :: ST, clean_seeds :: Bool,
        spawn_factor :: Int, max_rand_points :: Int, p_dist_th_factor, needs_scaling :: Bool ) where{ VT, ST }

        @assert dims > 0 "`dims` must be positive."
        @assert isnothing(n_points) || n_points >= 0 "`n_points` must be non-negative."
        @assert length(lb) == dims "`dims` and length of variable bounds must be the same."
        @assert length(lb) == length(ub) "Lower bounds `lb` and upper bounds `up` must have same length."
        @assert all(isfinite.(lb)) "All lower bounds `lb` must be finite."
        @assert all(isfinite.(ub)) "All upper bounds `ub` must be finite."

        copied_seeds = deepcopy(seeds)
        return new{VT,ST}(dims,n_points,lb,ub,copied_seeds,clean_seeds,spawn_factor,max_rand_points,p_dist_th_factor,needs_scaling)
    end

    # constructor for when vectors with non-float entries are provided
    function MonteCarloThDesign(dims, n_points, lb :: RVec, ub :: RVec,
        seeds :: RVecArr, clean_seeds, spawn_factor, max_rand_points, p_dist_th_factor, needs_scaling )
        T = Base.promote_eltype( lb, ub, eltype(seeds), Float64 )

        # we simply go with the array type of `lb` for `new_lb,new_ub,new_seeds`
        # to avoid abstract common types, e.g. promote_type( Vector{Float64}, SVector{3,Float64}) = AbstractVector{Float64}
        new_lb = similar(lb, T)
        new_lb[:] = [convert(T, x) for x in lb]

        new_ub = similar(lb, T)
        new_ub[:] = [convert(T, x) for x in ub]

        new_seeds = similar(seeds, typeof(new_lb))
        for (i,s) in enumerate(new_seeds)
            s[:] = [ convert(T, x) for x in seeds[i] ]
        end

        return MonteCarloThDesign(dims, n_points, new_lb, new_ub, new_seeds, clean_seeds, spawn_factor, max_rand_points, p_dist_th_factor, needs_scaling )
    end
end

function clean_seeds!( des :: MonteCarloThDesign )
    if des.clean_seeds
        cleaned_seeds = copy(des.seeds[good_indices( des.seeds, des.lb, des.ub)])
        empty!(des.seeds)
        append!(des.seeds, cleaned_seeds)
    end
    return nothing
end

function scale_seeds!( des :: MonteCarloThDesign ) :: Nothing
    if des.needs_scaling
        for seed in des.seeds
            scale_to_unit_square!( seed, des.lb, des.ub )
        end
    end
    return nothing
end

function unscale_seeds!( des :: MonteCarloThDesign ) :: Nothing
    if des.needs_scaling
        for seed in des.seeds
            unscale_from_unit_square!( seed, des.lb, des.ub )
        end
    end
    return nothing
end

function unscale_result( x, des :: MonteCarloThDesign )
    if des.needs_scaling
        return unscale_from_unit_square( x, des.lb, des.ub )
    else
        return x
    end
end

# NOTE / TODO
#=
Stateful iterators were introduced after I finished this.
The external state `point_array` could easily be stored in `des`.
=#

"Return the first design point (within original variable boundaries)
and the state -- a list of points in unit hypersquare generated so far."
function Base.iterate( des :: MonteCarloThDesign{VT,ST} ) where {VT, ST}
    if des.n_points > 0
        # throw out seeds that are not in variable bounds
        clean_seeds!(des)
        # scale seeds to unit hypersquare [0,1]^dims
        scale_seeds!(des)

        # design has no point associated, initalize with zero vector or first seed
        next_point = isempty( des.seeds ) ? VT(zeros(des.dims)) : des.seeds[1]
        point_array = [ next_point, ]
        return ( unscale_result(next_point, des), point_array )
    else
        return nothing
    end
end

function _rand_candidates( VT :: Type, d, N) :: Vector{VT}
    T = eltype(VT)
    return [ VT(rand(T,d)) for _ = 1 : N  ]
end

function Base.iterate( des :: MonteCarloThDesign{VT,ST}, point_array ) where {VT, ST}

    n_points_so_far = length(point_array)
    if n_points_so_far < length(des)
        if n_points_so_far < length( des.seeds )
            # as long as we have fewer points than seeds:
            # add seed to points
            next_point = des.seeds[ n_points_so_far + 1 ]
        else
            d = des.dims
            N = n_points_so_far
            spawn_factor = des.spawn_factor

            # spawn random candidate points
            num_candidates = min(N * spawn_factor, des.max_rand_points)
            candidates = _rand_candidates(VT, d, num_candidates)

            # compute candidate scores
            intersite_factor = (( N + 1 )^( 1/d ) - 1)/2
            pdist_factor = (N+1)/2
            th = 2 * des.p_dist_th_factor / N

            function (score_fn(p :: AbstractVector)::Float64)
                return intersite_factor * minimum( distance( p, point_array ) ) +
                    pdist_factor * minimum( projected_distance_thresholded(p, point_array, th) )
            end
            scores = score_fn.(candidates) :: Vector{Float64}

            # find point that maximizes score (i.e., is good for intersite distance and
            # projected distance )
            best_index = argmax( scores ) :: Int
            next_point = candidates[ best_index ] :: VT
        end
        push!( point_array, next_point )

        return unscale_result(next_point, des), point_array
    end
    return nothing
end

function Base.IteratorSize( des :: MonteCarloThDesign )
    if isnothing(des.n_points)
        Base.IsInfinite()
    else
        Base.HasLength()
    end
end
Base.length( des:: MonteCarloThDesign ) :: Union{Float64, Int} = isnothing(des.n_points) ? Inf : des.n_points

# so that collect returns the right kind of array
Base.IteratorEltype( :: MonteCarloThDesign ) = Base.HasEltype()
Base.eltype( :: MonteCarloThDesign{VT,ST} ) where {VT,ST} = VT

## Some legacy functions
@doc """
    monte_carlo_th( n_points = 10, n_dims = 2; seeds = [], spawn_factor = 50, pdist_threshold_tolerance = 0.5 )

Return an array of length `n_points` containing real vectors
representing points in space with `n_dims` dimensions.
The points are iteratively chosen from random point
sets to maximize a space-filling criterion as described in

"Surrogate Modelling of Computer Experiments with Sequential Experimental Design.", Crombecq, 2011

The returned point set is constructed starting with the points in `seeds`.
If `seeds` is empty (default), then the singleton set containing the zero vector is used.
"""
function monte_carlo_th(
    n_points :: Integer, n_dims :: Integer;
    seeds :: RVecArr = Vector{Float64}[], spawn_factor :: Integer = 100,
    pdist_threshold_tolerance :: Float64 = 0.5
    )

    iterator = MonteCarloThDesign(;
        dims = n_dims,
        n_points = n_points,
        seeds = seeds,
        spawn_factor = spawn_factor,
        p_dist_th_factor = pdist_threshold_tolerance
    )

    return collect( iterator )
end

@doc "Scale the design returned by the unconstrained version of this
function to the box defined by `lb` and `ub`."
function monte_carlo_th( n_points :: Integer, lb = RVec, ub = RVec;
    seeds :: RVecArr = Vector{Float64}[], spawn_factor :: Integer = 100,
    pdist_threshold_tolerance :: Float64 = 0.5, clean_seeds ::Bool = true
    )

    iterator = MonteCarloThDesign(;
        dims = length(lb),
        n_points = n_points,
        lb = lb,
        ub = ub,
        seeds = seeds,
        clean_seeds = clean_seeds,
        spawn_factor = spawn_factor,
        p_dist_th_factor = pdist_threshold_tolerance
    )

    return collect( iterator )
end


text:
module Link

using Gtk.ShortNames

export
    # Types
    AbstractState,
    State,
    # Functions
    link,
    disconnect,
    get,
    set!,
    set_quietly!,
    widget,
    id,
    lLabel,
    lEntry,
    lScale

abstract type AbstractState{T} end
abstract type AbstractLinkedWidget{T,W<:Gtk.GtkWidget} end

### AbstractState objects hold values and sync to UI elements
# They have at least three fields: `value`, `widgets`, and `canvases`

Base.eltype(::Type{AbstractState{T}}) where T = T
Base.eltype(::Type{S}) where {S<:AbstractState} = eltype(super(S))

Base.get(state::AbstractState) = state.value

add_widget!(state::AbstractState, w) = push!(state.widgets, w)
add_canvas!(state::AbstractState, w) = push!(state.canvases, w)

function Base.show(io::IO, state::AbstractState)
    print(io, typeof(state).name.name, "(", state.value)
    for w in state.widgets
        n = get_gtk_property(w.widget, :name, String)
        if !isempty(n)
            print(io, ",\"", n, "\"")
        end
    end
    print(io, ")")
end

### State

mutable struct State{T} <: AbstractState{T}
    value::T
    widgets::Vector
    canvases::Vector{Canvas}
end

function State(value::T;
         widgets::Vector = Array{AbstractLinkedWidget{T}}(undef,0),
         canvases::Vector{Canvas} = Array{Canvas}(undef,0)) where T
    State{T}(value, widgets, canvases)
end

function Base.similar(s::State)
    v = get(s)
    if !isimmutable(v)
        v = copy(v)
    end
    State(v, copy(s.widgets), copy(s.canvases))
end

function set!(state::State{T}, value) where {T}
    state.value = value
    for w in state.widgets
        set_quietly!(w, state.value)
    end
    for c in state.canvases
        draw(c)
    end
    state
end

# Set the value, skipping the input widget
function set!(state::State{T}, widget::AbstractLinkedWidget) where {T}
    value = get(widget)
    state.value = value
    for w in state.widgets
        w == widget && continue
        set_quietly!(w, state.value)
    end
    for c in state.canvases
        draw(c)
    end
    state
end

### Widgets linked to AbstractState objects

mutable struct LinkedWidget{T,W,S<:AbstractState} <: AbstractLinkedWidget{T,W}
    widget::W
    id::UInt
    state::S
end

widget(w::AbstractLinkedWidget) = w.widget
id(w::AbstractLinkedWidget) = w.id

Gtk.get_gtk_property(w::AbstractLinkedWidget, s, ::Type{T}) where {T} = get_gtk_property(widget(w), s, T)
Gtk.set_gtk_property!(w::AbstractLinkedWidget, s, val) = (set_gtk_property!(widget(w), s, val); w)
Base.push!(c::Gtk.GtkRadioButtonGroup, w::LinkedWidget) = push!(c, widget(w))  # ambiguity
Base.push!(c::Gtk.GtkContainer, w::LinkedWidget) = push!(c, widget(w))

"""
`set(w, val)` sets the value of the linked widget `w` and fires
the callback, thereby updating all other linked widgets.
"""
function set!(w::AbstractLinkedWidget, val)
    _set!(w, val)  # this might fire the callback, depending on the widget
    emit(w)        # provide this method for those that need explicit firing
    w
end

set!(w::LinkedWidget, val) = set!(w.state, val)

emit(w::AbstractLinkedWidget) = nothing   # fallback method

"""
`set_quietly!(w, val)` sets the value of the linked widget `w` without
firing the callback.
"""
function set_quietly!(w::AbstractLinkedWidget, val)
    ID = id(w)
    ID != 0 && signal_handler_block(w.widget, ID)
    _set!(w, val)
    ID != 0 && signal_handler_unblock(w.widget, ID)
    w
end

"""
`w_linked = link(state, widget)` links the value of the user-interface
element `widget` to the value of the `AbstractState` `state`. The two
will henceforth be synchronized: calling `get(state)` or
`get(w_linked)` returns the current value, and `set!(state, val)` or
`set!(w_linked, val)` will change the value for all mutually-linked
objects.

`link(state, c)`, where `c` is a `Canvas`, makes `c` a
listener for `state`. There is no return value.
"""
function link(val::AbstractState{T}, widget::Gtk.GtkWidget) where {T}
    w = LinkedWidget{T,typeof(widget),typeof(val)}(widget, 0, val)
    _link(val, w)
end

function link(val::AbstractState{T}, c::Canvas) where {T}
    add_canvas!(val, c)
    nothing
end

function _link(val::AbstractState{T}, w::AbstractLinkedWidget) where {T}
    w.id = create_callback(val, w)
    _set!(w, get(val))
    add_widget!(val, w)
    w
end

function disconnect(val::AbstractState)
    for w in val.widgets
        signal_handler_disconnect(w.widget, w.id)
    end
end

function disconnect(val::AbstractState, c::Canvas)
    index = findfirst(val.canvases, c)
    if index != 0
        deleteat!(val.canvases, index)
    end
    index
end

function Base.show(io::IO, w::LinkedWidget)
    print(io, "Linked ", typeof(w.widget).name.name, "(", get(w), ")")
    n = get_gtk_property(w.widget, :name, String)
    if !isempty(n)
        print(io, ",\"", n, "\"")
    end
end

## Label

function create_callback(val::AbstractState{T}, w::LinkedWidget{T,W}) where {T,W<:Label}
    0
end

function Base.get(w::LinkedWidget{T,W}) where {T<:AbstractString,W<:Label}
    val = get_gtk_property(w.widget, :label, String)
    convert(T, val)::T
end

function Base.get(w::LinkedWidget{T,W}) where {T<:Number,W<:Label}
    val = get_gtk_property(w.widget, :label, String)
    parse(T, val)::T
end

function _set!(w::LinkedWidget{T,W}, value) where {T,W<:Label}
    set_gtk_property!(w.widget, :label, string(value))
end

function lLabel(state::AbstractState)
    lbl = Label(string(get(state)))
    link(state, lbl)
end

## Entry

function create_callback(val::AbstractState{T}, w::LinkedWidget{T,W}) where {T,W<:Entry}
    signal_connect(w.widget, :activate) do widget
        set!(val, get(w))
    end
end

function Base.get(w::LinkedWidget{T,W}) where {T<:AbstractString,W<:Entry}
    val = get_gtk_property(w.widget, :text, String)
    convert(T, val)::T
end

function Base.get(w::LinkedWidget{T,W}) where {T<:Number,W<:Entry}
    val = get_gtk_property(w.widget, :text, String)
    parse(T, val)::T
end

function _set!(w::LinkedWidget{T,W}, value) where {T,W<:Entry}
    set_gtk_property!(w.widget, :text, string(value))
end

emit(w::LinkedWidget{T,W}) where {T,W<:Entry} = signal_emit(w.widget, :activate, Nothing)

function lEntry(state::AbstractState; kwargs...)
    e = Entry(;kwargs...)
    link(state, e)
end

## Scale

# To avoid a segfault, we have to use a low-level approach. See Gtk issue #161
function scale_cb(scaleptr::Ptr, state)
    widget = convert(Scale, scaleptr)
    v = G_.value(widget)
    set!(state, v)
    nothing
end

function create_callback(val::AbstractState{T}, w::LinkedWidget{T,W}) where {T,W<:Scale}
    signal_connect(scale_cb, w.widget, "value-changed", Nothing, (), false, val)
end

function Base.get(w::LinkedWidget{T,W}) where {T,W<:Scale}
    adj = Adjustment(w.widget)
    val = get_gtk_property(adj, :value, Int)
    convert(T, val)::T
end

function _set!(w::LinkedWidget{T,W}, value) where {T,W<:Scale}
    adj = Adjustment(w.widget)
    set_gtk_property!(adj, :value, value)
end

function lScale(state::AbstractState, args...)
    sc = Scale(args...)
    link(state, sc)
end

end

ext: jl


text:
using MPI
MPI.Init()

import mpiFFT4jl

function foo(comm::MPI.Comm)
    const N = [32, 64, 64]    # Global shape of mesh
    const L = [2pi, 2pi, 2pi] # Real size of mesh

    FFT = mpiFFT4jl.slab.r2c(N, L, comm)

    # Create a random array and its transform
    U = rand(Float64, mpiFFT4jl.slab.real_shape(FFT))
    U_hat  = Array{Complex{Float64}}(mpiFFT4jl.slab.complex_shape(FFT))
    U_hat2 = similar(U_hat)
    U2 = similar(U)
    U3 = similar(U)

    mpiFFT4jl.slab.rfft3(FFT, U_hat, U)

    for i in 1:10
        mpiFFT4jl.slab.irfft3(FFT, U2, U_hat)
#         mpiFFT4jl.slab.irfft3(FFT, U3, U_hat, 3)
    end

    (MPI.Reduce(sumabs2(U2-U), MPI.SUM, 0, comm), FFT)
end

comm = MPI.COMM_WORLD
@time (k, FFT) = foo(comm)
if MPI.Comm_rank(comm) == 0 println("Error  = $(k)") end
MPI.Finalize()


text:
import DielectricUpdateTechique as DUTCH
using Test

nCells = [5,12,10];
nUpdates = [4,4,4]
jv,eI,a,A,efft,pfft,pifft = DUTCH._allocate_memory_for_volume_integral_equation(nCells,nUpdates)

tl = (nCells[1]+2)*(nCells[2]+1)*(nCells[3]+1) + (nCells[1]+1)*(nCells[2]+2)*(nCells[3]+1) + (nCells[1]+1)*(nCells[2]+1)*(nCells[3]+2)

@test length(jv) == sum(nUpdates)
@test abs(sum(jv)) == 0
@test length(eI) == tl
@test abs(sum(eI)) == 0
@test length(a) == tl
@test abs(sum(a)) == 0
@test length(A) == 16*32*32*3
@test abs(sum(A)) == 0
@test length(efft) == 16*32*32*3
@test abs(sum(efft)) == 0

x,p,r,rt,u,v,q,uq = DUTCH._allocate_memory_cgs(nCells)

@test length(x) == tl
@test abs(sum(x)) == 0
@test length(p) == tl
@test abs(sum(p)) == 0
@test length(r) == tl
@test abs(sum(r)) == 0
@test length(rt) == tl
@test abs(sum(rt)) == 0
@test length(u) == tl
@test abs(sum(u)) == 0
@test length(v) == tl
@test abs(sum(v)) == 0
@test length(q) == tl
@test abs(sum(q)) == 0
@test length(uq) == tl
@test abs(sum(uq)) == 0


text:
using Test
using SDPNAL

# Example page 4 of SeDuMi_Guide_105R5.pdf
@testset "Linear Programming example" begin
    blk = ["l"  4]
    c = [ 1.0
         -1.0
          0.0
          0.0]
    A = [10.0 -7.0 -1.0 0.0
          1.0  0.5  0.0 1.0]
    b = [5.0
         3.0]
    obj, X, s, y, Z1, Z2, y2, v, info, runhist = sdpnalplus(blk, [Matrix(A')], [c], b, printlevel = 0)
    tol = 1e-5
    @test obj[1] ≈ -1/8 atol=tol rtol=tol
    @test obj[2] ≈ -1/8 atol=tol rtol=tol
    @test X[1] ≈ [47/24, 25//12, 0, 0] atol=tol rtol=tol
    @test isempty(s)
    @test y ≈ [1/8, -1/4] atol=tol rtol=tol
    @test Z1[1] ≈ [0, 0, 1/8, 1/4] atol=tol rtol=tol
    @test length(Z2) == 1
    @test size(Z2[1]) == (4, 1)
    @test all(iszero, Z2[1])
    @test isempty(y2)
    @test isempty(v)
    @test info["termcode"] == 0.0
end


text:
#
# Demo for unsteady flows
#
using Revise
using Compat.@info
using JLD
using FileIO
using HDF5

include("timestep_base.jl")

size_with_colorbar = [4,3]
size_no_colorbar = [3,3]

function hook_plot(solver, sol, ti, i)
    @time begin
        figure(2)
        clf()
        hook_vorticity(solver, sol, ti, i)
        hook_quiver(solver, sol, ti, i)
        println("* Plotting done")
    end
    show()
end

function hook_quiver(solver, sol, ti, i)
    # Unpack
    U1 = sol.U1
    U2 = sol.U2
    X = solver.X
    Y = solver.Y
    dcurve = solver.dcurve
    interior = solver.interior

    @. U1[interior==false] = NaN
    @. U2[interior==false] = NaN

    k = 20*2 # Point skip in quiver plot
    quiver(X[1:k:end, 1:k:end], Y[1:k:end, 1:k:end], U1[1:k:end, 1:k:end], U2[1:k:end, 1:k:end])
    plot_boundary(solver.dcurve)
    axis("image")
    axis("off")
end

function hook_vorticity(solver, sol, ti, i)
    X = solver.X
    Y = solver.Y
    vorticity = sol.dU2dx - sol.dU1dy
    interior = solver.interior
    @. vorticity[interior==false] = NaN
    s = 1
    pcolormesh(X[1:s:end,1:s:end], Y[1:s:end,1:s:end], vorticity[1:s:end,1:s:end],
               cmap=PyPlot.cm[:coolwarm])
    clim(-3,3)
    if figure(2).o[:get_size_inches]() == size_with_colorbar
        colorbar(shrink=0.8, pad=0, label="Vorticity")
    end
    plot_boundary(solver.dcurve)
    axis("image")
    axis("off")
end

function precompute(Re = 20, dt=0.1; R=0.5, Ngrid=150, pux_R_over_R=1/4, numpanels=100, sdc_order=1)
    dt_substep = sdc_substep(dt, sdc_order)
    alpha = sqrt(Re/dt_substep)
    @info("alpha = $alpha")
    figure(1)
    circle = AnalyticDomains.starfish(n_arms=0, amplitude=0.0, radius=R)
    solver = precompute_solver(alpha=alpha,
                               pux_R=R*pux_R_over_R,
                               Ngrid=Ngrid,
                               curve = circle,
                               numpanels = numpanels
                               )
    return solver, dt, sdc_order
end

function timestep(pre, Nstep=10; velocity=1)
    solver, dt, sdc_order = pre
    # Set BCs
    dcurve = solver.dcurve
    xbdry = dcurve.points[1,:]
    ybdry = dcurve.points[2,:]
    n1 = dcurve.normals[1,:]
    n2 = dcurve.normals[2,:]
    phi = (x, y) -> atan2(y, x)
    k = 1
    ubc1 = @. -n1*sin(k*phi(xbdry, ybdry))*velocity
    ubc2 = @. -n2*sin(k*phi(xbdry, ybdry))*velocity

    compat = @. (n1*ubc1 + n2*ubc2)*dcurve.dS
    info("Compability condition=", sum(compat))
    assert(sum(compat) < 1e-10)

    # Estimate BC resolution
    panelorder = dcurve.panelorder
    L = mss.legendre_matrix(panelorder)
    bc_resolution = 0.0
    for i=1:dcurve.numpanels
        xcoeff = L*ubc1[(1:panelorder) + panelorder*(i-1)]
        ycoeff = L*ubc2[(1:panelorder) + panelorder*(i-1)]
        xres = max(abs(xcoeff[end]), abs(xcoeff[end-1]))
        yres = max(abs(ycoeff[end]), abs(ycoeff[end-1]))
        res = max(xres, yres)
        bc_resolution = max(bc_resolution, res)
    end
    info("BC resolution: ", bc_resolution)
    # Start from rest
    X = solver.X
    Y = solver.Y
    U1, U2 = zeros(X), zeros(X)
    RHS1, RHS2 = zeros(X), zeros(X)
    dU1dx, dU1dy = zeros(X), zeros(X)
    dU2dx, dU2dy = zeros(X), zeros(X)

    sol0 = Solution(U1, U2, dU1dx, dU1dy, dU2dx, dU2dy, RHS1, RHS2, [], [])
    # Run
    t0 = 0.0
    istart = 0
    hook = hook_plot
    t, sol = timestep(Nstep, dt, sdc_order, hook, solver, t0, sol0, ubc1, ubc2; istart=istart)
end

#  pre = precompute()
# timestep(pre)

function instab1()
    close("all")
    figure(2, figsize=size_no_colorbar)
    pre = precompute(50, 0.04; R=0.5, Ngrid=150*2, pux_R_over_R=1/4, numpanels=100*2, sdc_order=1)
    timestep(pre, 45, velocity=1)
end

function instab2()
    close("all")
    figure(2, figsize=size_no_colorbar)
    pre = precompute(50, 0.02; R=0.5, Ngrid=150*2, pux_R_over_R=1/4, numpanels=100*2, sdc_order=1)
    timestep(pre, 200, velocity=1)
end

function instab3()
    close("all")
    figure(2, figsize=size_no_colorbar)
    pre = precompute(100, 0.02; R=0.5, Ngrid=150*2, pux_R_over_R=1/4, numpanels=100*2, sdc_order=1)
    timestep(pre, 45, velocity=1)
end

function instab4()
    close("all")
    figure(2, figsize=size_with_colorbar)
    pre = precompute(100, 0.01; R=0.5, Ngrid=150*2, pux_R_over_R=1/4, numpanels=100*2, sdc_order=1)
    timestep(pre, 200, velocity=1)
end

function write_figure(num)
    pngname = "../docs/paper/figs/time_stability_$(num).png"
    savefig(pngname, dpi=200)
    Base.run(`mogrify -trim $pngname`)
end

function paper_scripts()
    instab1()
    write_figure(1)
    instab2()
    write_figure(2)
    instab3()
    write_figure(3)
    instab4()
    write_figure(4)
end


text:
@testset "soss compat" begin
    k = SqExponentialKernel()
    y = randn(3)
    X = randn(3, 1)
    x = [rand(1) for _ in 1:3]
    GPRegression = @model X begin
        # Priors.
        α ~ LogNormal(0.0, 0.1)
        ρ ~ LogNormal(0.0, 1.0)
        σ² ~ LogNormal(0.0, 1.0)

        # Realized covariance function
        kernel = α * transform(SqExponentialKernel(), 1 / ρ)
        f = GP(kernel)

        # Sampling Distribution.
        y ~ f(X, σ²)
    end
    # Test for matrices
    m = GPRegression(X=RowVecs(X))
    @test length(dynamicHMC(m, (y=y,), 5)) == 5
    # Test for vectors of vector
    m = GPRegression(X=x)
    @test length(dynamicHMC(m, (y=y,), 5)) == 5
end


text:
import Base: copy, copyto!, ==, llvmcall
import Random: rand, seed!
import RandomNumbers: gen_seed, union_uint, seed_type, unsafe_copyto!, unsafe_compare

"The key for AESNI."
mutable struct AESNIKey
    key1::__m128i
    key2::__m128i
    key3::__m128i
    key4::__m128i
    key5::__m128i
    key6::__m128i
    key7::__m128i
    key8::__m128i
    key9::__m128i
    key10::__m128i
    key11::__m128i
    AESNIKey() = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
end

copyto!(dest::AESNIKey, src::AESNIKey) = unsafe_copyto!(dest, src, UInt128, 11)

copy(src::AESNIKey) = copyto!(AESNIKey(), src)

==(key1::AESNIKey, key2::AESNIKey) = unsafe_compare(key1, key2, UInt128, 11)

"""
Assistant function for AES128. Compiled from the C++ source code:
```cpp
R123_STATIC_INLINE __m128i AES_128_ASSIST (__m128i temp1, __m128i temp2) {
    __m128i temp3;
    temp2 = _mm_shuffle_epi32 (temp2 ,0xff);
    temp3 = _mm_slli_si128 (temp1, 0x4);
    temp1 = _mm_xor_si128 (temp1, temp3);
    temp3 = _mm_slli_si128 (temp3, 0x4);
    temp1 = _mm_xor_si128 (temp1, temp3);
    temp3 = _mm_slli_si128 (temp3, 0x4);
    temp1 = _mm_xor_si128 (temp1, temp3);
    temp1 = _mm_xor_si128 (temp1, temp2);
    return temp1;
}
```
"""
_aes_128_assist(a::__m128i, b::__m128i) = llvmcall(
    """%3 = bitcast <2 x i64> %1 to <4 x i32>
    %4 = shufflevector <4 x i32> %3, <4 x i32> undef, <4 x i32> <i32 3, i32 3, i32 3, i32 3>
    %5 = bitcast <4 x i32> %4 to <2 x i64>
    %6 = bitcast <2 x i64> %0 to <16 x i8>
    %7 = shufflevector <16 x i8> <i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 0, i8 0, i8 0, i8 0>, <16 x i8> %6, <16 x i32> <i32 12, i32 13, i32 14, i32 15, i32 16, i32 17, i32 18, i32 19, i32 20, i32 21, i32 22, i32 23, i32 24, i32 25, i32 26, i32 27>
    %8 = bitcast <16 x i8> %7 to <2 x i64>
    %9 = xor <2 x i64> %8, %0
    %10 = shufflevector <16 x i8> <i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 0, i8 0, i8 0, i8 0>, <16 x i8> %7, <16 x i32> <i32 12, i32 13, i32 14, i32 15, i32 16, i32 17, i32 18, i32 19, i32 20, i32 21, i32 22, i32 23, i32 24, i32 25, i32 26, i32 27>
    %11 = bitcast <16 x i8> %10 to <2 x i64>
    %12 = xor <2 x i64> %9, %11
    %13 = shufflevector <16 x i8> <i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 undef, i8 0, i8 0, i8 0, i8 0>, <16 x i8> %10, <16 x i32> <i32 12, i32 13, i32 14, i32 15, i32 16, i32 17, i32 18, i32 19, i32 20, i32 21, i32 22, i32 23, i32 24, i32 25, i32 26, i32 27>
    %14 = bitcast <16 x i8> %13 to <2 x i64>
    %15 = xor <2 x i64> %12, %5
    %16 = xor <2 x i64> %15, %14
    ret <2 x i64> %16""",
    __m128i, Tuple{__m128i, __m128i},
    a, b
)

function _aesni_expand!(k::AESNIKey, rkey::__m128i)
    k.key1 = rkey
    tmp = _aes_key_gen_assist(rkey, Val(0x1))
    rkey = _aes_128_assist(rkey, tmp)
    k.key2 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x2))
    rkey = _aes_128_assist(rkey, tmp)
    k.key3 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x4))
    rkey = _aes_128_assist(rkey, tmp)
    k.key4 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x8))
    rkey = _aes_128_assist(rkey, tmp)
    k.key5 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x10))
    rkey = _aes_128_assist(rkey, tmp)
    k.key6 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x20))
    rkey = _aes_128_assist(rkey, tmp)
    k.key7 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x40))
    rkey = _aes_128_assist(rkey, tmp)
    k.key8 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x80))
    rkey = _aes_128_assist(rkey, tmp)
    k.key9 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x1b))
    rkey = _aes_128_assist(rkey, tmp)
    k.key10 = rkey

    tmp = _aes_key_gen_assist(rkey, Val(0x36))
    rkey = _aes_128_assist(rkey, tmp)
    k.key11 = rkey

    k
end

AESNIKey(key::UInt128) = _aesni_expand!(AESNIKey(), __m128i(key))

"""
```julia
AESNI1x <: AbstractAESNI1x
AESNI1x([seed])
```

AESNI1x is one kind of AESNI Counter-Based RNGs. It generates one `UInt128` number at a time.

`seed` is an `Integer` which will be automatically converted to `UInt128`.

Only available when [`R123_USE_AESNI`](@ref).
"""
mutable struct AESNI1x <: AbstractAESNI1x
    x::__m128i
    ctr::__m128i
    key::AESNIKey
end

function AESNI1x(seed::Integer=gen_seed(UInt128))
    r = AESNI1x(0, 0, AESNIKey())
    seed!(r, seed)
    r
end

function seed!(r::AESNI1x, seed::Integer=gen_seed(UInt128))
    r.x = zero(__m128i)
    r.ctr = zero(__m128i)
    _aesni_expand!(r.key, __m128i(seed % UInt128))
    random123_r(r)
    r
end

seed_type(::Type{AESNI1x}) = UInt128

function copyto!(dest::AESNI1x, src::AESNI1x)
    dest.x = src.x
    dest.ctr = src.ctr
    copyto!(dest.key, src.key)
    dest
end

copy(src::AESNI1x) = copyto!(AESNI1x(), src)

==(r1::AESNI1x, r2::AESNI1x) = r1.x == r2.x && r1.key == r2.key && r1.ctr == r2.ctr

"""
```julia
AESNI4x <: AbstractAESNI4x
AESNI4x([seed])
```

AESNI4x is one kind of AESNI Counter-Based RNGs. It generates four `UInt32` numbers at a time.

`seed` is a `Tuple` of four `Integer`s which will all be automatically converted to `UInt32`.

Only available when [`R123_USE_AESNI`](@ref).
"""
mutable struct AESNI4x <: AbstractAESNI4x
    x::__m128i
    ctr1::__m128i
    key::AESNIKey
    p::Int
end

function AESNI4x(seed::NTuple{4, Integer}=gen_seed(UInt32, 4))
    r = AESNI4x(zero(__m128i), zero(__m128i), AESNIKey(), 0)
    seed!(r, seed)
    r
end

function seed!(r::AESNI4x, seed::NTuple{4, Integer}=gen_seed(UInt32, 4))
    key = union_uint(Tuple(x % UInt32 for x in seed))
    r.ctr1 = 0
    _aesni_expand!(r.key, __m128i(key))
    r.p = 0
    random123_r(r)
    r
end

seed_type(::Type{AESNI4x}) = NTuple{4, UInt32}

function copyto!(dest::AESNI4x, src::AESNI4x)
    unsafe_copyto!(dest, src, UInt128, 2)
    copyto!(dest.key, src.key)
    dest.p = src.p
    dest
end

copy(src::AESNI4x) = copyto!(AESNI4x(), src)
==(r1::AESNI4x, r2::AESNI4x) = unsafe_compare(r1, r2, UInt128, 2) &&
    r1.key == r2.key && r1.p == r2.p

function aesni1xm128i(input::__m128i, key::AESNIKey)
    x = key.key1 ⊻ input
    x = _aes_enc(x, key.key2)
    x = _aes_enc(x, key.key3)
    x = _aes_enc(x, key.key4)
    x = _aes_enc(x, key.key5)
    x = _aes_enc(x, key.key6)
    x = _aes_enc(x, key.key7)
    x = _aes_enc(x, key.key8)
    x = _aes_enc(x, key.key9)
    x = _aes_enc(x, key.key10)
    x = _aes_enc_last(x, key.key11)
end

@inline function random123_r(r::AESNI1x)
    r.x = aesni1xm128i(r.ctr, r.key)
    (UInt128(r.x),)
end

@inline function random123_r(r::AESNI4x)
    r.x = aesni1xm128i(r.ctr1, r.key)
    split_uint(UInt128(r.x), UInt32)
end


text:
module EntityResolution

include("entity/entity.jl")

end # module


text:
include("tensor_cache.jl")

using AbstractTrees
using QXContexts

export build_compute_graph

"""
build_compute_graph(tnc::TensorNetworkCircuit,
                   plan::Vector{<:Tuple}),
                   bond_groups::Union{Nothing, Vector{<:Vector{<:Index}}}=nothing)

Function to build a compute tree from a tensor network and plan
"""
function build_compute_graph(tnc::TensorNetworkCircuit,
                            plan::Vector{<:Tuple},
                            bond_groups::Union{Nothing, Vector{<:Vector{<:Index}}}=nothing)
    tnc = copy(tnc)
    tn = tnc.tn
    # create load and outputs nodes
    nodes = Dict{Symbol, ComputeNode}()
    tc = TensorCache()

    # add a load node for each tensor
    for t in keys(tnc)
        data_symbol = push!(tc, tensor_data(tnc, t))
        op = LoadCommand(t, data_symbol, collect(size(tensor_data(tnc, t))))
        nodes[t] = ComputeNode{LoadCommand}(op)
    end

    # overwrite tensors corresponding to output tensors with output comamnds
    for (i, o) in enumerate(output_tensors(tnc))
        op = OutputCommand(o, i, size(tnc[o])[1])
        nodes[o] = ComputeNode{OutputCommand}(op)
    end

    changed_ids = Dict{Symbol, Symbol}()
    # add view commands to slice tensors
    if bond_groups !== nothing
        for (i, bg) in enumerate(bond_groups)
            related_tensors = union([tnc[b] for b in bg]...)
            slice_sym = Symbol("v$(i)")
            slice_dim = QXTns.dim(bg[1])
            for t in related_tensors
                new_sym = Symbol("$(t)_s")
                indices = hyperindices(tnc, t, all_indices=true)
                # find which group overlaps with slice_bonds
                index_position = findfirst(x -> length(intersect(bg, x)) > 0, indices)
                op = ViewCommand(new_sym, t, slice_sym, index_position, slice_dim)
                replace_tensor_symbol!(tn, t, new_sym)
                changed_ids[t] = new_sym
                node = ComputeNode(op)
                push!(node.children, nodes[t])
                nodes[t].parent = node
                nodes[new_sym] = node
            end
        end
    end

    # now add contraction nodes following plan
    for c in plan
        A_sym, B_sym, C_sym = c
        while haskey(changed_ids, A_sym) A_sym = changed_ids[A_sym] end
        while haskey(changed_ids, B_sym) B_sym = changed_ids[B_sym] end
        r = contraction_indices(tn, A_sym, B_sym)
        op = ContractCommand(C_sym, r.c_labels, A_sym, r.a_labels, B_sym, r.b_labels)
        node = ComputeNode(op)
        for s in [A_sym, B_sym]
            nodes[s].parent = node
            push!(node.children, nodes[s])
        end
        nodes[C_sym] = node
        contract_pair!(tn, A_sym, B_sym, C_sym; mock=true)
    end

    # if there are multiple disconnected trees, contract these so there is a single tree s
    parentless = collect(keys(filter(x -> !isdefined(x[2], :parent), nodes)))
    reduce(parentless[2:end], init=parentless[1]) do x, y
        r = contraction_indices(tn, x, y)
        op = ContractCommand(:dummy, r.c_labels, x, r.a_labels, y, r.b_labels)
        sym = contract_pair!(tn, x, y, mock=true)
        op.output_name = sym
        node = ComputeNode{ContractCommand}(op)
        nodes[sym] = node
        for s in [x, y]
            push!(node.children, nodes[s])
            nodes[s].parent = node
        end
        sym # return sym to contract with next parentless node
    end
    parentless = collect(keys(filter(x -> !isdefined(x[2], :parent), nodes)))
    @assert length(parentless) == 1 "Only root node should have no parent"
    root = parentless[1]
    node = ComputeNode(SaveCommand(:output, root))
    push!(node.children, nodes[root])
    nodes[root].parent = node
    ComputeGraph(node, convert(Dict, tc))
end


text:
# # Creating supercells with pymatgen
#
# The [Pymatgen](https://pymatgen.org/) python library allows to setup
# solid-state calculations using a flexible set of classes as well as an API
# to an online data base of structures. Its `Structure` and `Lattice`
# objects are directly supported by the DFTK `load_atoms` and `load_lattice`
# functions, such that DFTK may be readily used to run calculation on systems
# defined in pymatgen. Using the `pymatgen_structure` function a conversion
# from DFTK to pymatgen structures is also possible. In the following we
# use this to create a silicon supercell and find its LDA ground state
# using direct minimisation.

# First we setup the silicon lattice in DFTK.
using DFTK

a = 10.263141334305942  # Lattice constant in Bohr
lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]];

# Next we make a `[2, 2, 2]` supercell using pymatgen
pystruct = pymatgen_structure(lattice, atoms)
pystruct.make_supercell([2, 2, 2])
lattice = load_lattice(pystruct)
atoms = [Si => [s.frac_coords for s in pystruct.sites]];

# Setup an LDA model and discretize using
# a single k-point and a small `Ecut` of 5 Hartree.
model = model_LDA(lattice, atoms)
basis = PlaneWaveBasis(model; Ecut=5, kgrid=(1, 1, 1))

# Find the ground state using direct minimisation (always using SCF is boring ...)
scfres = direct_minimization(basis, tol=1e-5);
#-
scfres.energies


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

name = "libcxxwrap_julia"
version = v"0.7.0"

const is_yggdrasil = haskey(ENV, "BUILD_BUILDNUMBER")
git_repo = is_yggdrasil ? "https://github.com/JuliaInterop/libcxxwrap-julia.git" : joinpath(ENV["HOME"], "src/julia/libcxxwrap-julia/")
unpack_target = is_yggdrasil ? "" : "libcxxwrap-julia"

# Collection of sources required to complete build
sources = [
    GitSource(git_repo, "7959f09947f95e67ee094a2ee3fb8b7f9be696f8", unpack_target=unpack_target),
    ArchiveSource("https://julialang-s3.julialang.org/bin/linux/armv7l/1.3/julia-1.3.1-linux-armv7l.tar.gz", "965c8fab2214f8ce1b3d449d088561a6de61be42543b48c3bbadaed5b02bf824"; unpack_target="julia-arm-linux-gnueabihf"),
    ArchiveSource("https://julialang-s3.julialang.org/bin/linux/x64/1.3/julia-1.3.1-linux-x86_64.tar.gz", "faa707c8343780a6fe5eaf13490355e8190acf8e2c189b9e7ecbddb0fa2643ad"; unpack_target="julia-x86_64-linux-gnu"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-apple-darwin14.tar.gz", "f2e5359f03314656c06e2a0a28a497f62e78f027dbe7f5155a5710b4914439b1"; unpack_target="julia-x86_64-apple-darwin14"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-w64-mingw32.tar.gz", "c7b2db68156150d0e882e98e39269301d7bf56660f4fc2e38ed2734a7a8d1551"; unpack_target="julia-x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""

case "$target" in
	arm-linux-gnueabihf|x86_64-linux-gnu)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/julia-1.3.1
        ;;
    x86_64-apple-darwin14|x86_64-w64-mingw32)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/juliabin
        ;;
esac

mkdir build
cd build
cmake -DJulia_PREFIX=$Julia_PREFIX -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ../libcxxwrap-julia/
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/libcxxwrap-julia*/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
    Windows(:x86_64),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libcxxwrap_julia", :libcxxwrap_julia; dlopen_flags=[:RTLD_GLOBAL]),
    LibraryProduct("libcxxwrap_julia_stl", :libcxxwrap_julia_stl; dlopen_flags=[:RTLD_GLOBAL]),
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

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

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
"""
    ispangram(input)

Return `true` if `input` contains every alphabetic character (case insensitive).

"""
ispangram(input) =  all(∈(lowercase(input)), 'a':'z')

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
module Bubu

using Genie.Renderer.Html, Genie.Renderer, Stipple, StippleUI

@reactive mutable struct Model <: ReactiveModel
    process::R{Bool} = false
    expanded::R{Vector{Bool}} = [true, false]
    output::R{String} = ""
    input::R{String} = ""
    name::R{Vector{String}} = ["fliegner", "Burmeister"]
    given::R{Vector{String}} = ["michael", "irmgard"]
    age::R{Int} = 0
    warin::R{Bool} = true
end

function handlers(model)
    on(model.process) do _
        if (model.process[])
            model.output[] = model.input[] |> reverse
            model.given[1] = model.name[2] |> reverse
            model.given[2] = model.name[1] |> reverse
            println(
                "person" *
                model.name[1] *
                " " *
                model.name[2] *
                model.given[1] *
                " " *
                model.given[2],
            )
            model.process[] = false
        end
    end
    on(model.isready) do _
        push!(model)
    end
    model
end

function ui(model)
    pinputs = ""
    for i = 1:length(model.name[])
        pinputs =
            pinputs * p(
                [
                    "DynamicName Person $(i)"
                    input(
                        "",
                        @bind(Symbol("name[$(i-1)]")),
                        @on("keyup.enter", "process = true")
                    )
                ],
            )
        pinputs =
            pinputs * p(
                [
                    "Dynamic Given Person $(i)"
                    input(
                        "",
                        @bind(Symbol("given[$(i-1)]")),
                        @on("keyup.enter", "process = true")
                    )
                ],
            )
    end
    println("pinputs")
    println(pinputs)
    page(
        model,
        class = "container",
        [
            pinputs
            ul(
                expansionitem(
                    expandseparator = true,
                    icon = "signal_wifi_off",
                    label = "Wifi settings",
                    [
                        template(
                            p(
                                [
                                    "Name {{index}} = {{name[index]}}"
                                    input(
                                        "",
                                        @bind(Symbol("name[index]")),
                                        @on("keyup.enter", "process = true")
                                    )
                                    "Given {{index}} = {{given[index]}}"
                                    input(
                                        "",
                                        @bind(Symbol("given[index]")),
                                        @on("keyup.enter", "process = true")
                                    )
                                ],
                            ),
                            @recur(:"(todo,index) in name")
                        ),
                        ul(
                            expansionitem(
                                expandseparator = true,
                                icon = "signal_wifi_on",
                                label = "Wifi settings",
                                [
                                    card([
                                        cardsection(
                                            "Lorem ipsum dolor sit amet, consectetur adipisicing elit. Quidem, eius reprehenderit eos corrupti
                                commodi magni quaerat ex numquam, dolorum officiis modi facere maiores architecto suscipit iste
                                eveniet doloribus ullam aliquid.",
                                        ),
                                    ]),
                                ],
                            ),
                        ),
                    ],
                ),
            )
            StippleUI.form(
                action = "/sub",
                method = "POST",
                [
                    textfield(
                        "What's your name *",
                        Symbol("name[0]"),
                        name = "name",
                        @iif(:warin),
                        :filled,
                        hint = "Name and surname",
                        "lazy-rules",
                        rules = "[val => val && val.length > 0 || 'Please type something']",
                    ),
                    numberfield(
                        "Your age *",
                        :age,
                        name = "age",
                        "filled",
                        :lazy__rules,
                        rules = "[val => val !== null && val !== '' || 'Please type your age',
                          val => val > 0 && val < 100 || 'Please type a real age']",
                    ),
                    btn("submit", type = "submit", color = "primary"),
                ],
            )
        ],
    )
end

model = handlers(Stipple.init(Model; transport = Genie.WebThreads))

route("/reactive") do
    html(ui(model), context = @__MODULE__)
end

route("/sub", method = POST) do
    age = params(:age, "99")
    name = params(:name, "Schnipperdey")
    println("PARAMS " * name * "  " * string(age))
    model.name = ["a", "b", "Schnipperdey"]
    redirect("/reactive")
end

isrunning(:webserver) || up()
end #module

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
# Random fields

@doc raw"""
    struct RandomField
        names::Set{String}
        factors::Dict{Vector{String},Function}
        evidence::Dict{String,Any}
    end

A representation of a random field by a collection of factors:
``\log p(x) = \sum_{f \in \mathcal F} \log \psi_f(x_f)``, where
``\mathcal F`` is the set of (log) factors and ``x_f`` is the set of variables
incident on that (log) factor. `factors` should be properly normalized *log* mass or density
functions. There is no restriction on the state space of the variables involved as long
as the factor functions can evaluate the log probability of the variables.
For example,
```
factor_ab(x) = logpdf(MvNormal(PDMat([1.0 0.5; 0.5 1.0])), x)
factor_bc(x) = logpdf(MvNormal([2.0, 2.0], PDMat([2.0 -1.0; -1.0 2.0])), x)
```
are two valid (log) factor functions -- the first corresponding to the factor ``\log \psi_{a,b}(x_a,x_b)``
and the second corresponding to the factor ``\log \psi_{b,c}(x_b, x_c)``.
Posting evidence is done using a dictionary mapping an address to value,
e.g., `evidence = Dict("b" => 3.0)`.

Calling a random field corresponds to evaluating its log probability with the passed argument,
e.g.,
```
my_rf = RandomField(...)
x = Dict("a" => 1.0, "b" => -2.1)
my_lp = my_rf(x)  # corresponds to logprob(my_rf, x)
```
"""
struct RandomField
    names::Set{String}
    factors::Dict{Vector{String},Function}
    evidence::Dict{String,Any}
end

@doc raw"""
    struct GenerativeField{T}
        field::RandomField
        proposal::T
        val::Dict
        ...
    end

A generative field is a three-tuple of a random field `field`, a proposal kernel `proposal`,
and an initialization value `val` used to start MH sampling.
"""
struct GenerativeField{T}
    field::RandomField
    proposal::T
    val::Dict
    function GenerativeField{T}(field::RandomField, proposal::T, val::Dict) where T
        for (k,v) in field.evidence
            val[k] = v
        end
        new(field, proposal, val)
    end
end
GenerativeField(field::RandomField, proposal::T, val::Dict) where T = GenerativeField{T}(field::RandomField, proposal::T, val::Dict)
export GenerativeField

@doc raw"""
    function Distributions.rand(gf::GenerativeField; num_iterations=2500)

Sample from a generative fied using `num_iterations` of MH sampling. The value
generated by MH at the end of `num_iterations` of sampling is returned.
To customize burn-in, you can just overload this function, e.g.
`rand(gf::GenerativeField) = rand(gf; num_iterations=10000)`.
"""
function Distributions.rand(gf::GenerativeField; num_iterations=2500)
    results = mh(
        gf.field, [gf.proposal], gf.val;
        burn = num_iterations, thin = 1, num_iterations = num_iterations + 1
    )
    sample = Dict(k => v[1] for (k,v) in results.values)
    for (k, v) in sample
        gf.val[k] = v
    end
    sample
end

@doc raw"""
    function logprob(rf::RandomField, x::Dict)

Evaluates the log probability of a set of values against the density described by the
random field `rf`. The values `x` should have the format `address => value`.
"""
function logprob(rf::RandomField, x::Dict)
    lp = 0.0
    evidence_vars = keys(rf.evidence)
    for (variables, factor) in rf.factors
        to_score = []
        for v in variables
            # safeguard proposing to observed components
            if v in evidence_vars
                push!(to_score, rf.evidence[v])
            else
                push!(to_score, x[v])
            end
        end
        lp += rf.factors[variables](to_score)
    end
    lp
end
(rf::RandomField)(x::Dict) = logprob(rf, x)
logprob(gf::GenerativeField, x::Dict) = logprob(gf.field, x)

@doc raw"""
    function RandomField(factors::Dict{Vector{String},Function})

Outer constuctor for `RandomField` that requires only a dict of factors.
"""
function RandomField(factors::Dict{Vector{String},Function})
    evidence = Dict{String,Any}()
    RandomField(factors, evidence)
end

@doc raw"""
    function RandomField(factors::Dict{Vector{String},Function}, evidence::Dict)

Outer constructor for `RandomField` that requires a dict of factors and allows posting
evidence when the field is created (instead of manually doing so later).
"""
function RandomField(factors::Dict{Vector{String},Function}, evidence::Dict)
    names = Set{String}()
    for k in keys(factors)
        union!(names, k)
    end
    RandomField(names, factors, evidence)
end

export RandomField, logprob

accept(x, prob_x, x_prime, prob_x_prime, log_a::Float64) = log(rand()) < log_a ? (x_prime, prob_x_prime) : (x, prob_x)
export accept


text:
using Random
using BayesianIntegral
using Distributions
using LinearAlgebra
using Statistics

twister = MersenneTwister(10)
X = rand(twister, 20,2)
y = X[:,1] .- (X[:,2] .^ 2) .+ 12
noise = 0.04

hyper = GaussianKernelHyperparameters(1.0, [0.1,0.1])
modl = KrigingModel(X, y, hyper)

# Testing it.
output = modl(X[1,:])
output.standard_error < 1e-6
output = modl(X[2,:])
output.standard_error < 1e-6
output = modl([0.5,0.5])
output.standard_error > 0.001
output = modl([0.75,0.25])
output.standard_error > 0.001

################################################################################

PointToExamine = [0.44, 0.44]
fmin = minimum(modl.f)

EI = expected_improvement(modl, fmin, PointToExamine)
EI2 = expected_improvement(modl, fmin + 0.1, PointToExamine)
EI < EI2

xx = get_next_query_point_through_expected_improvement(modl)

hyper = GaussianKernelHyperparameters(1.0, [0.1,0.1])
batch_size = Integer(floor(size(X)[1]/2))
step_multiple = 0.05
noise = 0.0001
steps = 1000
calibrated_parameters = calibrate_by_ML_with_SGD(X, y, hyper, steps, batch_size, step_multiple, noise)
modl = KrigingModel(X, y, calibrated_parameters; noise  = noise)
xx = get_next_query_point_through_expected_improvement(modl)

# Testing integral
dist = MultivariateNormal([0.5,0.5], Hermitian([1.0 0.0; 0.0 1.0]))
int = integrate(modl, [0.5,0.5], Hermitian([1.0 0.0; 0.0 1.0])).expectation
samples = 100000
numerical_int = zeros(samples)
for i in 1:samples
    x_ = rand(twister, dist, 1)
    numerical_int[i] =  modl(x_[:,1]).value
end
numerical_int = Statistics.mean(numerical_int)


text:
# see readme for module documentation
module PyAMG

export RugeStubenSolver,
       SmoothedAggregationSolver,
       solve, solve!, aspreconditioner,
       set_kwargs!

using LinearAlgebra
using PyCall
using SparseArrays

# scipy_sparse = pyimport("scipy.sparse")
# pyamg = pyimport("pyamg")
pyamg = PyNULL()
scipy_sparse = PyNULL()
function __init__()
    copy!(pyamg, pyimport_conda("pyamg", "pyamg"))
    copy!(scipy_sparse, pyimport_conda("scipy.sparse", "scipy"))
end

"""
`py_csc(A::SparseMatrixCSC) -> PyObject`

Takes a Julia CSC matrix and converts it into `PyObject`, which stores a
`scipy.sparse.csc_matrix`.
"""
function py_csc(A::SparseMatrixCSC)
   # create an empty sparse matrix in Python
   Apy = scipy_sparse[:csc_matrix](size(A))
   # write the values and indices
   Apy[:data] = copy(A.nzval)
   Apy[:indices] = A.rowval .- 1
   Apy[:indptr] = A.colptr .- 1
   return Apy
end

"""
`py_csr(A::SparseMatrixCSC) -> PyObject`

Takes a Julia CSC matrix and converts it into `PyObject`, which stored a
`scipy.sparse.csr_matrix`. (Note: it first converts it to `csc_matrix`, then
calls `tocsr()`)

TODO: this seems extremely inefficient and should at some point be fixed.
"""
py_csr(A::SparseMatrixCSC) = py_csc(A)[:tocsr]()

"""
`struct AMGSolver{T}`

Encapsulates the AMG solver types implemented in PyAMG.

Initialise using
```
using PyAMG
amg_rs = RugeStubenSolver(A)
amg_sa = SmoothedAggregationSolver(A)
```
To solve Ax = b:
```
x = solve(amg, b, tol=1e-6, accel="cg")
```
"""
mutable struct AMGSolver{T}
    po::PyObject
    id::T
    kwargs::Vector
    A::SparseMatrixCSC
end

struct RugeStuben end
struct SmoothedAggregation end
global const RugeStubenSolver = AMGSolver{RugeStuben}
global const SmoothedAggregationSolver = AMGSolver{SmoothedAggregation}

function set_kwargs!(amg::AMGSolver; kwargs...)
   amg.kwargs = collect(kwargs)
end

"""
`RugeStubenSolver(A::SparseMatrixCSC, kwargs...)`:

Create a Ruge Stuben instance of `AMGSolver`; wraps
`pyamg.ruge_stuben_solver`; see `pyamg.ruge_stuben_solver?`
for keyword arguments.  See `?AMGSolver` for usage.
"""
RugeStubenSolver(A::SparseMatrixCSC; kwargs...) =
   AMGSolver(pyamg[:ruge_stuben_solver](py_csr(A)),
             RugeStuben(), collect(Any, kwargs), A)

"""
`SmoothedAggregationSolver(A::SparseMatrixCSC, kwargs...)`

Wrapper for `pyamg.smoothed_aggregation_solver`; see `pyamg.ruge_stuben_solver?`
for keyword arguments. See `?AMGSolver` for usage.
"""
SmoothedAggregationSolver(A::SparseMatrixCSC; kwargs...) =
   AMGSolver(pyamg[:smoothed_aggregation_solver](py_csr(A)),
             SmoothedAggregation(), collect(Any, kwargs), A)

"""
`solve(A::SparseMatrixCSC, b::Vector; kwargs...)`:

PyAMG's 'blackbox' solver. See `pyamg.solve?` for `kwargs`.
"""
function solve(A::SparseMatrixCSC, b::Vector; kwargs...)
   # If kwargs contains :residuals, we need to do some conversions, since
   # Python cannot append to Julia arrays (i.e. numpy arrays).
   for (n, (key, rj)) in enumerate(kwargs)
      if key == :residuals
         rp = PyVector(Float64[])
         kwargs[n] = (:residuals, rp)
         try
            x = pyamg[:solve]( py_csr(A), b; kwargs... )
            append!(rj, collect(rp))
            return x::Vector{Float64}
         catch
            error("Something went wrong. Probably, your version of pyamg probably does not support the `residuals` keyword; please update `pyamg` (see https://github.com/pyamg/pyamg) or call `PyAMG.solve` without a `residuals` keyword.")
         end
      end
   end
   # If we are here, then we just solve and return
   return pyamg[:solve]( py_csr(A), b; kwargs... )::Vector{Float64}
end

"""
`solve(amg::AMGSolver, b, kwargs...)`

Returns a `Vector` with the result of the AMG solver. The keyword
arguments can either be passed directly, or can be stored in
`amd` via `set_kwargs!`.

### `kwargs`  (copy-pasted from Python docs)

* `x0` : Initial guess.
* `tol` : Stopping criteria: relative residual r[k]/r[0] tolerance.
* `maxiter` : Stopping criteria: maximum number of allowable iterations.
* `cycle` : {"V","W","F","AMLI"}
    Type of multigrid cycle to perform in each iteration.
* `accel` : Defines acceleration method.  Can be a string such as "cg"
    or "gmres" which is the name of an iterative solver in
    `pyamg.krylov` (preferred) or scipy.sparse.linalg.isolve.
    If accel is not a string, it will be treated like a function
    with the same interface provided by the iterative solvers in SciPy.
         (the function version is not tested in Julia!)
* `callback` : User-defined function called after each iteration.  It is
    called as callback(xk) where xk is the k-th iterate vector.
         (also not tested in Julia!)
* `residuals` : List to contain residual norms at each iteration.
"""
function solve(amg::AMGSolver, b::Vector; kwargs...)
   # If kwargs contains :residuals, we need to do some conversions, since
   # Python cannot append to Julia arrays (i.e. numpy arrays).
   for (n, (key, rj)) in enumerate(kwargs)
      if key == :residuals
         rp = PyVector(Float64[])
         kwargs[n] = (:residuals, rp)
         x = amg.po[:solve](b; amg.kwargs..., kwargs...)
         append!(rj, collect(rp))
         return x::Vector{Float64}
      end
   end
   # If we are here, then we just solve and return
   return amg.po[:solve](b; amg.kwargs..., kwargs...)::Vector{Float64}
end

# function solve(amg::AMGSolver, b; history=false, kwargs...)
#    if history
#       r = PyVector(Float64[])
#       x = amg.po[:solve](b; amg.kwargs..., kwargs..., residuals=r)
#       return x, collect(r)
#    else
#       return amg.po[:solve](b; amg.kwargs..., kwargs...)
#    end
# end

######### Capability to use PyAMG.jl as a preconditioner for
######### nonlinear optimisation, sampling, etc

# TODO: the following 4 methods still need to be tested

import Base.\, Base.*
\(amg::AMGSolver, b::Vector) = solve(amg, b; amg.kwargs...)
*(amg::AMGSolver, x::Vector) = amg.A * x

LinearAlgebra.ldiv!(x, amg::AMGSolver, b) = copyto!(x, amg \ b)
LinearAlgebra.mul!(b, amg::AMGSolver, x) = mul!(b, amg.A, x)

##############################################################################
######### Capability to use PyAMG.jl as a preconditioner for
######### iterative linear algebra

"""
`struct AMGPreconditioner`

returned by `aspreconditioner(amg)`, when `amg` is of type `AMGSolver`.
This stores `PyObject` that acts as a linear operator. This type
should be used when `PyAMG` is used as a preconditioner for
iterative linear algebra.

Overloaded methods that can be used with an `AMGPreconditioner` are
`\`, `*`
"""
struct AMGPreconditioner
  po::PyObject
  A::SparseMatrixCSC
end

"""
`aspreconditioner(amg::AMGSolver; kwargs=...)`

returns an `M::AMGPreconditioner` object that is suitable for usage
as a preconditioner for iterative linear algebra.

If `x` is a vector, then `M \\ x` denotes application of the
preconditioner (i.e. 1 MG cycle), while `M * x` denotes
multiplication with the original matrix from which `amg` was constructed.

### kwargs:
cycle : {'V','W','F','AMLI'}
    Type of multigrid cycle to perform in each iteration.
"""
aspreconditioner(amg::AMGSolver; kwargs...) =
      AMGPreconditioner(amg.po[:aspreconditioner](kwargs...), amg.A)

\(amg::AMGPreconditioner, b::Vector) = amg.po[:matvec](b)
*(amg::AMGPreconditioner, x::Vector) = amg.A * x

LinearAlgebra.ldiv!(x, amg::AMGPreconditioner, b) = copyto!(x, amg \ b)
LinearAlgebra.mul!(b, amg::AMGPreconditioner, x) = mul!(b, amg.A, x)

"""
`diagnostics(A)`:

Wrapper for `solver_diagnostics`, which is part of PyAMG-Examples. To use it
clone [PyAMG-Examples](https://github.com/pyamg/pyamg-examples), then make the
call to `diagnostics` from a directory where `solver_diagnostics.py` is located
(e.g., from `pyamg-examples/solver_diagnostics`, but the file could be copied
anywhere.

If succesful, `diagnostics(A)` will try a variety of parameter combinations for
`smoothed_aggregation_solver`, and write two files 'solver_diagnostic.txt' and
'solver_diagnostic.py' which contain information how to generate the best
solver.

This implementation is a really poor hack, and suggestions how to improve it
would be highly appreciated.
"""
function diagnostics(A::SparseMatrixCSC; kwargs...)
    # try to import solver_diagnostics
    try
        unshift!(PyVector(pyimport("sys")["path"]), "");
        solver_diagnostics = pyimport("solver_diagnostics")
    catch
        error("""
              I tried to pyimport solver_diagnostics, but it fails.
              This is probably because `solver_diagnostics.py` is not
              in the current directory. Please see
              `?PyAMG.diagonistics` on how to use this function.
              """)
    end

    # import has worked against expectations; call the diagnostics
    solver_diagnostics.solver_diagnostics(A; kwargs...)
end

end


text:
export CRRLearner

"""
    CRRLearner(;kwargs)

See paper: [Critic Regularized Regression](https://arxiv.org/abs/2006.15134).

# Keyword arguments

- `approximator`::[`ActorCritic`](@ref): used to get Q-values (Critic) and logits (Actor) of a state.
- `target_approximator`::[`ActorCritic`](@ref): similar to `approximator`, but used to estimate the target.
- `γ::Float32`, reward discount rate.
- `batch_size::Int=32`
- `policy_improvement_mode::Symbol=:exp`, type of the weight function f. Possible values: :binary/:exp.
- `ratio_upper_bound::Float32`, when `policy_improvement_mode` is ":exp", the value of the exp function is upper-bounded by this parameter.
- `beta::Float32`,  when `policy_improvement_mode` is ":exp", this is the denominator of the exp function.
- `advantage_estimator::Symbol=:mean`, type of the advantage estimate \\hat{A}. Possible values: :mean/:max.
- `m::Int=4`, when `continuous=true`, sample `m` action to estimate \\hat{A}.
- `update_freq::Int`: the frequency of updating the `approximator`.
- `update_step::Int=0`
- `target_update_freq::Int`: the frequency of syncing `target_approximator`.
- `continuous::Bool`: type of action space.
- `rng = Random.GLOBAL_RNG`
"""
mutable struct CRRLearner{
    Aq<:ActorCritic,
    At<:ActorCritic,
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Aq
    target_approximator::At
    γ::Float32
    batch_size::Int
    policy_improvement_mode::Symbol
    ratio_upper_bound::Float32
    beta::Float32
    advantage_estimator::Symbol
    m::Int
    update_freq::Int
    update_step::Int
    target_update_freq::Int
    continuous::Bool
    rng::R
    # for logging
    actor_loss::Float32
    critic_loss::Float32
end

function CRRLearner(;
    approximator::Aq,
    target_approximator::At,
    γ::Float32 = 0.99f0,
    batch_size::Int = 32,
    policy_improvement_mode::Symbol = :exp,
    ratio_upper_bound::Float32 = 20.0f0,
    beta::Float32 = 1.0f0,
    advantage_estimator::Symbol = :mean,
    m::Int = 4,
    update_freq::Int = 10,
    update_step::Int = 0,
    target_update_freq::Int = 100,
    continuous::Bool,
    rng = Random.GLOBAL_RNG,
) where {Aq<:ActorCritic, At<:ActorCritic}
    copyto!(approximator, target_approximator)
    CRRLearner(
        approximator,
        target_approximator,
        γ,
        batch_size,
        policy_improvement_mode,
        ratio_upper_bound,
        beta,
        advantage_estimator,
        m,
        update_freq,
        update_step,
        target_update_freq,
        continuous,
        rng,
        0.0f0,
        0.0f0,
    )
end

Flux.functor(x::CRRLearner) = (Q = x.approximator, Qₜ = x.target_approximator),
y -> begin
    x = @set x.approximator = y.Q
    x = @set x.target_approximator = y.Qₜ
    x
end

function (learner::CRRLearner)(env)
    s = state(env)
    s = Flux.unsqueeze(s, ndims(s) + 1)
    s = send_to_device(device(learner), s)
    if learner.continuous
        learner.approximator.actor(s; is_sampling=true) |> vec |> send_to_host
    else
        learner.approximator.actor(s) |> vec |> send_to_host
    end
end

function RLBase.update!(learner::CRRLearner, batch::NamedTuple)
    if learner.continuous
        continuous_update!(learner, batch)
    else
        discrete_update!(learner, batch)
    end
end

function continuous_update!(learner::CRRLearner, batch::NamedTuple)
    AC = learner.approximator
    target_AC = learner.target_approximator
    γ = learner.γ
    beta = learner.beta
    batch_size = learner.batch_size
    policy_improvement_mode = learner.policy_improvement_mode
    ratio_upper_bound = learner.ratio_upper_bound
    advantage_estimator = learner.advantage_estimator
    D = device(AC)

    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = reshape(a, :, batch_size)
    r = reshape(r, :, batch_size)
    t = reshape(t, :, batch_size)

    target_a_t = target_AC.actor(s′; is_sampling=true)
    target_q_input = vcat(s′, target_a_t)
    expected_target_q = target_AC.critic(target_q_input)

    target = r .+ γ .* (1 .- t) .* expected_target_q

    q_t = Matrix{Float32}(undef, learner.m, batch_size)
    for i in 1:learner.m
        a_sample = AC.actor(s; is_sampling=true)
        q_t[i, :] = AC.critic(vcat(s, a_sample))
    end

    ps = Flux.params(AC)
    gs = gradient(ps) do
        # Critic loss
        qa_t = AC.critic(vcat(s, a))
        critic_loss = Flux.Losses.mse(qa_t, target)

        # Actor loss
        log_π = AC.actor.model(s, a)

        if advantage_estimator == :max
            advantage = qa_t .- maximum(q_t, dims=1)
        elseif advantage_estimator == :mean
            advantage = qa_t .- mean(q_t, dims=1)
        else
            error("Wrong parameter.")
        end

        if policy_improvement_mode == :binary
            actor_loss_coef = (advantage .> 0.0f0)
        elseif policy_improvement_mode == :exp
            actor_loss_coef = clamp.(exp.(advantage ./ beta), 0, ratio_upper_bound)
        else
            error("Wrong parameter.")
        end

        actor_loss = mean(-log_π .* actor_loss_coef)

        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
        end

        actor_loss + critic_loss
    end

    update!(AC, gs)
end

function discrete_update!(learner::CRRLearner, batch::NamedTuple)
    AC = learner.approximator
    target_AC = learner.target_approximator
    γ = learner.γ
    beta = learner.beta
    batch_size = learner.batch_size
    policy_improvement_mode = learner.policy_improvement_mode
    ratio_upper_bound = learner.ratio_upper_bound
    advantage_estimator = learner.advantage_estimator
    D = device(AC)

    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = CartesianIndex.(a, 1:batch_size)
    r = reshape(r, :, batch_size)
    t = reshape(t, :, batch_size)

    target_a_t = softmax(target_AC.actor(s′))
    target_q_t = target_AC.critic(s′)
    expected_target_q = sum(target_a_t .* target_q_t, dims=1)

    target = r .+ γ .* (1 .- t) .* expected_target_q

    ps = Flux.params(AC)
    gs = gradient(ps) do
        # Critic loss
        q_t = AC.critic(s)
        qa_t = reshape(q_t[a], :, batch_size)
        critic_loss = Flux.Losses.mse(qa_t, target)

        # Actor loss
        a_t = softmax(AC.actor(s))

        if advantage_estimator == :max
            advantage = qa_t .- maximum(q_t, dims=1)
        elseif advantage_estimator == :mean
            advantage = qa_t .- mean(q_t, dims=1)
        else
            error("Wrong parameter.")
        end

        if policy_improvement_mode == :binary
            actor_loss_coef = (advantage .> 0.0f0)
        elseif policy_improvement_mode == :exp
            actor_loss_coef = clamp.(exp.(advantage ./ beta), 0, ratio_upper_bound)
        else
            error("Wrong parameter.")
        end

        actor_loss = mean(-log.(a_t[a]) .* actor_loss_coef)

        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
        end

        actor_loss + critic_loss
    end

    update!(AC, gs)
end


text:
module PTimersTests

include("mpiexec.jl")
run_mpi_driver(procs=4,file="driver_p_timers.jl")

end # module


text:
const Scheme = Dict{String, Vector{Tuple{Vector{String}, Float64}}}

function parse_schemes(io)
    results = Dict{String, Scheme}()

    calgo, cpair = nothing, nothing
    for line in eachline(io) @when !isempty(line)
        if (m = match(r"\*\*\*(.*)\*\*\*", line)) != nothing
            results[m[1]] = Scheme()
            calgo = results[m[1]]
        elseif (m = match(r"(.*) -> (.*) :", line)) != nothing
            if m[1] != m[2]
                calgo["$(m[1]) $(m[2])"] = []
                cpair = calgo["$(m[1]) $(m[2])"]
            else
                cpair = nothing
            end
        elseif cpair != nothing
            path, weight = split(line, '@')
            weight = parse(f64, weight)
            path = [join(edge.args, ' ') for edge in parse(path).args]
            push!(cpair, (path, weight))
        end
    end

    results
end

function renormalize(scheme)
    map(scheme) do p
        pair, paths = p
        tw = sum(cadr, paths)
        pair => [(path, weight / tw) for (path, weight) in paths]
    end
end

function renormalize!(scheme, pair)
    paths = scheme[pair]
    tw = sum(cadr, paths)
    scheme[pair] = [(path, weight / tw) for (path, weight) in paths]
end

function select_hard_nop(nodes, edges, scheme, budget)
    function meet_budget()
        d = Dict(n => 0 for n in nodes)
        for (pair, paths) in scheme, (path, weight) in paths, edge in path
            n = split(edge, ' ') |> car
            if car(n) == 's'
                d[n] += 1
            end
        end
        all(x <= budget[n] for (n, x) in d)
    end

    function cut_down(K)
        Dict(pair => length(paths) > K ? sort(paths, by=cadr, rev=true)[1:K] : paths for (pair, paths) in scheme)
    end

    K = maximum(map(length, collect(values(scheme))))

    while !meet_budget() && K > 0
        scheme = cut_down(K)
        K -= 1
    end

    renormalize(scheme)
end

function select_greedy(nodes, edges, scheme, budget)
    scheme = deepcopy(scheme)

    function dispute_one()
        d = Dict(n => [] for n in nodes)
        for (pair, paths) in scheme, (i, (path, weight)) in enumerate(paths), edge in path
            n = split(edge, ' ') |> car
            if car(n) == 's'
                push!(d[n], (pair, i, weight))
            end
        end

        overloads = Set()
        for (n, v) in d @when length(v) > budget[n]
            push!(overloads, v...)
        end

        if length(overloads) > 0
            pair, i, w = min(overloads..., by=i"3")
            deleteat!(scheme[pair], i)
            renormalize!(scheme, pair)
            true
        else
            false
        end
    end

    while dispute_one() end
    scheme
end

function select_program(nodes, edges, scheme, budget)
    # model 1: find K
    m = grb.Model("K")
    m[:setParam]("OutputFlag", false)
    K = m[:addVar](vtype=grb.GRB[:INTEGER], name="K")

    pathvar = Dict()
    for (pair, paths) in scheme
        p1, p2 = split(pair, ' ')
        pathvar[pair] = map(1:length(paths)) do i
            m[:addVar](vtype=grb.GRB[:BINARY], name="s_$(p1)_$(p2)_$(i)")
        end
    end
    m[:update]()

    # constraint 1: at least K paths
    for (pair, var) in pathvar
        p1, p2 = split(pair, ' ')
        m[:addConstr](py"sum($var) >= $K", name="d_$(p1)_$(p2)")
    end

    # constraint 2: flow table budget
    for n in nodes
        passed = []
        for (pair, paths) in scheme, (i, path) in enumerate(paths)
            for edge in car(path) @when n in split(edge, ' ')
                push!(passed, pathvar[pair][i])
                break
            end
        end
        if length(passed) > 0
            m[:addConstr](py"sum($passed) <= $(budget[n])", name="n_$n")
        end
    end

    m[:setObjective](K, sense=grb.GRB[:MAXIMIZE])
    m[:optimize]()

    K = K[:X]

    # model 2: maximize total number of paths
    m = grb.Model("totalnop")
    m[:setParam]("OutputFlag", false)

    pathvar = Dict()
    for (pair, paths) in scheme
        p1, p2 = split(pair, ' ')
        pathvar[pair] = map(1:length(paths)) do i
            m[:addVar](vtype=grb.GRB[:BINARY], name="s_$(p1)_$(p2)_$(i)")
        end
    end
    m[:update]()

    # constraint 1: at least K paths
    for (pair, var) in pathvar
        p1, p2 = split(pair, ' ')
        m[:addConstr](py"sum($var) >= $K", name="d_$(p1)_$(p2)")
    end

    # constraint 2: flow table budget
    for n in nodes
        passed = []
        for (pair, paths) in scheme, (i, path) in enumerate(paths)
            for edge in car(path) @when n in split(edge, ' ')
                push!(passed, pathvar[pair][i])
                break
            end
        end
        if length(passed) > 0
            m[:addConstr](py"sum($passed) <= $(budget[n])", name="n_$n")
        end
    end

    all_paths = [v for l in values(pathvar) for v in l]

    m[:setObjective](py"sum($all_paths)", sense=grb.GRB[:MAXIMIZE])
    m[:optimize]()

    Dict(pair => [path for (path, choice) in zip(paths, pathvar[pair]) if choice[:X] > .5]
         for (pair, paths) in scheme)
end

max_stars_repo_name: ylxdzsw/TEexp

max_issues_repo_name: ylxdzsw/TEexp

max_forks_repo_name: ylxdzsw/TEexp


text:
# Test file for 3-d Chebyshev polynomials

using ChebyshevApprox

n1 = 10
n2 = 15
n3 = 7

dom_1 = [2.0, -3.0]
dom_2 = [1.5,  0.5]
dom_3 = [0.5, -0.5]
dom = [dom_1 dom_2 dom_3]

nodesr_1  = chebyshev_nodes(n1,dom_1)
nodesr_2  = chebyshev_nodes(n2,dom_2)
nodesr_3  = chebyshev_nodes(n3,dom_3)
nodesex_1 = chebyshev_extended(n1,dom_1)
nodesex_2 = chebyshev_extended(n2,dom_2)
nodesex_3 = chebyshev_extended(n3,dom_3)

yr  = zeros(n1,n2,n3)
yex = zeros(n1,n2,n3)

for i = 1:n1
  for j = 1:n2
    for k = 1:n3

      yr[i,j,k]  = (nodesr_1[i]+4.0)^0.5+nodesr_1[i]*sqrt(nodesr_2[j])+exp(nodesr_3[k])
      yex[i,j,k] = (nodesex_1[i]+4.0)^0.5+nodesex_1[i]*sqrt(nodesex_2[j])+exp(nodesex_3[k])

    end
  end
end

order_tensor = [7, 6, 6]
order_complete = 6

wr_tensor    = chebyshev_weights(yr,[nodesr_1,nodesr_2,nodesr_3],order_tensor,dom)
wr_complete  = chebyshev_weights(yr,(nodesr_1,nodesr_2,nodesr_3),order_complete,dom)
wex_tensor   = chebyshev_weights_extended(yex,(nodesex_1,nodesex_2,nodesex_3),order_tensor,dom)
wex_complete = chebyshev_weights_extended(yex,(nodesex_1,nodesex_2,nodesex_3),order_complete,dom)

point = [1.748, 0.753, 0.119]

yr_tensor     = chebyshev_evaluate(wr_tensor,point,order_tensor,dom)
yr_complete   = chebyshev_evaluate(wr_complete,point,order_complete,dom)
yex_tensor    = chebyshev_evaluate(wex_tensor,point,order_tensor,dom)
yex_complete  = chebyshev_evaluate(wex_complete,point,order_tensor,dom)

y_actual = (point[1]+4.0)^0.5+point[1]*sqrt(point[2])+exp(point[3])

println([y_actual yr_tensor yr_complete yex_tensor yex_complete])

yr_grad_tensor       = chebyshev_gradient(wr_tensor,point,order_tensor,dom)
yr_grad_complete     = chebyshev_gradient(wr_complete,point,order_complete,dom)
yr_deriv_tensor_1    = chebyshev_derivative(wr_tensor,point,1,order_tensor,dom)
yr_deriv_complete_1  = chebyshev_derivative(wr_tensor,point,1,order_complete,dom)
yex_grad_tensor      = chebyshev_gradient(wex_tensor,point,order_tensor,dom)
yex_grad_complete    = chebyshev_gradient(wex_complete,point,order_complete,dom)
yex_deriv_tensor_1   = chebyshev_derivative(wex_tensor,point,1,order_tensor,dom)
yex_deriv_complete_1 = chebyshev_derivative(wex_complete,point,1,order_complete,dom)

println(yr_grad_tensor, yr_grad_complete, yr_deriv_tensor_1, yr_deriv_complete_1)
println(yex_grad_tensor, yex_grad_complete, yex_deriv_tensor_1, yex_deriv_complete_1)

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
# routines for creating interpolatable forms of grid data

export interpolated_history

"""
    interpolated_history(v::History,tr)

Returns a function that interpolates the history `v`, using the time specified
by time range `tr`. Uses cubic spline interpolation. If `v` is a periodic history,
then the resulting function is also periodic.
"""
function interpolated_history(v::History{T,PeriodicHistory},tr::AbstractRange) where {T}
  length(tr) == length(v)+1 || error("Incommensurate lengths of history and range")
  return CubicSplineInterpolation(tr,Base.unsafe_view(v,1:length(tr)),extrapolation_bc=Periodic())
end

function interpolated_history(v::History{T,RegularHistory},tr::AbstractRange) where {T}
  length(tr) == length(v) || error("Incommensurate lengths of history and range")
  return CubicSplineInterpolation(tr,v,extrapolation_bc=Flat())
end

"""
    interpolated_history(v::History,xr,yr,tr)

Returns a function that interpolates the history `v` both spatially and temporally,
using the x and y ranges specified by `xr` and `yr` and time specified by time range
`tr`. Uses cubic spline interpolation in all directions. If `v` is a periodic history,
then the resulting function is also periodic in time.
"""
function interpolated_history(v::History{T,H},xr::AbstractRange,yr::AbstractRange,tr::AbstractRange) where {T,H<:HistoryType}
  # v is a PeriodicHistory, so need to add first element to last here
  nx, ny = size(v[1])

  extralen = (H == PeriodicHistory) ? 1 : 0

  (length(xr) == nx && length(yr) == ny && length(tr) == length(v)+extralen) || error("Incommensurate lengths of history and range")

  vhist = zeros(nx,ny,length(v)+extralen)
  for i in 1:size(vhist,3)
    vhist[:,:,i] = v[i]
  end

  bctype = (H == PeriodicHistory) ? Periodic() : Flat()

  return CubicSplineInterpolation((xr, yr ,tr), vhist, extrapolation_bc = (Flat(),Flat(),bctype))

end


text:
using DataFrames
using LinearAlgebra
using MixedModels
using Random
using SparseArrays
using StatsModels
using Test

using MixedModels: dataset, levels, modelcols, nlevs

const LMM = LinearMixedModel

@testset "scalarReMat" begin
    ds = dataset("dyestuff")
    f1 = @formula(yield ~ 1 + (1|batch))
    y1, Xs1 = modelcols(apply_schema(f1, schema(ds), LMM), ds)
    sf = Xs1[2]
    psts = dataset("pastes")
    f2 = @formula(strength ~ 1 + (1|batch/cask))
    y2, Xs2 = modelcols(apply_schema(f2, schema(psts), LMM), psts)
    sf1 = Xs2[2]
    sf2 = Xs2[3]

    @testset "size" begin
        @test size(sf) == (30, 6)
        @test size(sf,1) == 30
        @test size(sf,2) == 6
        @test size(sf,3) == 1
        @test size(sf1) == (60, 10)
        @test size(sf2) == (60, 30)
    end

    @testset "utilities" begin
        @test levels(sf) == string.('A':'F')
        @test refpool(sf) == levels(sf)
        @test refarray(sf) == repeat(1:6, inner=5)
        @test refvalue(sf, 3) == "C"
        @test nlevs(sf) == 6
        @test eltype(sf) == Float64
        @test sparse(sf) == sparse(1:30, sf.refs, ones(30))
        fsf = Matrix(sf)
        @test size(fsf) == (30, 6)
        @test count(!iszero, fsf) == 30
        @test sort!(unique(fsf)) == [0.0, 1.0]
        @test cond(sf) == 1.0
        @test MixedModels.nθ(sf) == 1
        @test MixedModels.getθ(sf) == ones(1)
        @test MixedModels.getθ!(Vector{Float64}(undef, 1), sf) == ones(1)
        @test lowerbd(sf) == zeros(1)
        @test MixedModels.getθ(setθ!(sf, [0.5])) == [0.5]
        MixedModels.unscaledre!(Vector{Float64}(undef, 30), sf)
        @test_throws DimensionMismatch MixedModels.getθ!(Float64[], sf)
        @test_throws DimensionMismatch setθ!(sf, ones(2))
    end

    @testset "products" begin
        @test ones(30, 1)'sf == fill(5.0, (1, 6))
        @test mul!(Array{Float64}(undef, (size(sf1, 2), size(sf2, 2))), sf1', sf2) == Array(sf1'sf2)

        crp = sf'sf
        @test isa(crp, Diagonal{Float64})
        crp1 = copy(crp)
        @test crp1 == crp
        @test crp[2,6] == 0
        @test crp[6,6] == 5
        @test size(crp) == (6,6)
        @test crp.diag == fill(5.,6)
        rhs = y1'sf
        @test rhs == reshape([7525.0,7640.0,7820.0,7490.0,8000.0,7350.0], (1, 6))
        @test ldiv!(crp, copy(rhs)') == [1505.,1528.,1564.,1498.,1600.,1470.]

        @test isa(sf1'sf1, Diagonal{Float64})
        @test isa(sf2'sf2, Diagonal{Float64})
        @test isa(sf2'sf1,SparseMatrixCSC{Float64})

        @test MixedModels.lmulΛ!(sf', ones(6)) == fill(0.5, 6)
        @test MixedModels.rmulΛ!(ones(6, 6), sf) == fill(0.5, (6, 6))
    end

    @testset "reweight!" begin
        wts = rand(MersenneTwister(1234321), size(sf, 1))
        @test isapprox(vec(MixedModels.reweight!(sf, wts).wtz), wts)
    end
end

@testset "RandomEffectsTerm" begin
    slp = dataset("sleepstudy")
    contrasts =  Dict{Symbol,Any}()

    @testset "Detect same variable as blocking and experimental" begin
        f = @formula(reaction ~ 1 + (1 + subj|subj))
        @test_throws ArgumentError apply_schema(f, schema(f, slp, contrasts), LMM)
    end

    @testset "Detect both blocking and experimental variables" begin
        # note that U is not in the fixed effects because we want to make square
        # that we're detecting all the variables in the random effects
        f = @formula(reaction ~ 1 + (1 + days|subj))
        form = apply_schema(f, schema(f, slp, contrasts), LMM)
        @test StatsModels.termvars(form.rhs) == [:days, :subj]
    end

    @testset "Runtime construction of random effects terms" begin
        # operator precedence and basic terms:
        @test term(:a) | term(:b) isa RandomEffectsTerm
        @test term(1) + term(:a) | term(:b) isa RandomEffectsTerm
        @test term(1) + term(:a) + term(:a) & term(:c) | term(:b) isa RandomEffectsTerm

        # sleep study data:
        r, d, s, one = term.((:reaction, :days, :subj, 1))

        f1 = @formula(reaction ~ 1 + (1 + days | subj))
        f2 = r ~ one + (one + d | s)
        @test f2.rhs[end] isa RandomEffectsTerm
        ff1 = apply_schema(f1, schema(slp), LMM)
        ff2 = apply_schema(f2, schema(slp), LMM)
        # equality of RE terms not defined so check that they generate same modelcols
        @test modelcols(ff1.rhs[end], slp) == modelcols(ff2.rhs[end], slp)

        m1 = fit(LMM, f1, slp; progress=false)
        m2 = fit(LMM, f2, slp; progress=false)
        @test all(m1.λ .== m2.λ)

        @test StatsModels.terms(f2.rhs[end]) == [one, d, s]
        @test StatsModels.termvars(f2.rhs[end]) == [d.sym, s.sym]
    end

    @testset "Runtime construction of ZeroCorr" begin
        r, d, s, one = term.((:reaction, :days, :subj, 1))

        f1 = @formula(reaction ~ 1 + zerocorr(1 + days | subj))
        f2 = r ~ one + zerocorr(one + d | s)
        @test f2.rhs[end] isa MixedModels.ZeroCorr
        ff1 = apply_schema(f1, schema(slp), LMM)
        ff2 = apply_schema(f2, schema(slp), LMM)
        # equality of RE terms not defined so check that they generate same modelcols
        mc1 = modelcols(ff1.rhs[end], slp)
        mc2 = modelcols(ff2.rhs[end], slp)

        # test that zerocorr actually worked
        @test mc1.inds == mc2.inds == [1, 4]

        m1 = fit(LMM, f1, slp; progress=false)
        m2 = fit(LMM, f2, slp; progress=false)
        @test all(m1.λ .== m2.λ)

        @test StatsModels.terms(f2.rhs[end]) == [one, d, s]
        @test StatsModels.termvars(f2.rhs[end]) == [d.sym, s.sym]
    end

    @testset "ZeroCorr delegation" begin
        r, d, s, one = term.((:reaction, :days, :subj, 1))

        f = @formula(0 ~ 1 + days).rhs
        zc = zerocorr(one + d | s)
        @test f == zc.lhs
        @test zc.rhs.sym == :subj
    end

end

@testset "Categorical Blocking Variable" begin
    # deepcopy because we're going to modify it.  Don't need the copy if dataset returns an Arrow.Table
    #slp = deepcopy(DataFrame(dataset("sleepstudy")))
    slp = DataFrame(dataset("sleepstudy"))
    contrasts =  Dict{Symbol,Any}()
    f = @formula(reaction ~ 1 + (1|subj))

    # String blocking-variables work fine because StatsModels is smart enough to
    # treat strings as Categorical. Note however that this is a
    # far less efficient to store the original dataframe, although it doesn't
    # matter for the contrast matrix
    slp[!,:subj] = convert.(String, slp[!, :subj])
    # @test_throws ArgumentError LinearMixedModel(f, slp)
    slp.subj = parse.(Int, getindex.(slp.subj, Ref(2:4)))
    @test_throws ArgumentError LinearMixedModel(f, slp)
end

@testset "random effects term syntax" begin

    dat = (y = rand(18),
           g = string.(repeat('a':'f', inner=3)),
           f = string.(repeat('A':'C', outer=6)))

    @testset "fulldummy" begin
        @test_throws ArgumentError fulldummy(1)

        f = @formula(y ~ 1 + fulldummy(f))
        f1 = apply_schema(f, schema(dat))
        @test typeof(last(f1.rhs.terms)) <: FunctionTerm{typeof(fulldummy)}
        @test_throws ArgumentError modelcols(f1, dat)

        f2 = apply_schema(f, schema(dat), MixedModel)
        @test typeof(last(f2.rhs.terms)) <: CategoricalTerm{<:StatsModels.FullDummyCoding}
        @test modelcols(f2.rhs, dat)[1:3, :] == [1 1 0 0
                                                 1 0 1 0
                                                 1 0 0 1]

        # implict intercept
        ff = apply_schema(@formula(y ~ 1 + (f | g)), schema(dat), MixedModel)
        rem = modelcols(last(ff.rhs), dat)
        @test size(rem) == (18, 18)
        @test rem[1:3, 1:4] == [1 0 0 0
                                1 1 0 0
                                1 0 1 0]

        # explicit intercept
        ff = apply_schema(@formula(y ~ 1 + (1+f | g)), schema(dat), MixedModel)
        rem = modelcols(last(ff.rhs), dat)
        @test size(rem) == (18, 18)
        @test rem[1:3, 1:4] == [1 0 0 0
                                1 1 0 0
                                1 0 1 0]

        # explicit intercept + full dummy
        ff = apply_schema(@formula(y ~ 1 + (1+fulldummy(f) | g)), schema(dat), MixedModel)
        rem = modelcols(last(ff.rhs), dat)
        @test size(rem) == (18, 24)
        @test rem[1:3, 1:4] == [1 1 0 0
                                1 0 1 0
                                1 0 0 1]

        # explicit dropped intercept (implicit full dummy)
        ff = apply_schema(@formula(y ~ 1 + (0+f | g)), schema(dat), MixedModel)
        rem = modelcols(last(ff.rhs), dat)
        @test size(rem) == (18, 18)
        @test rem[1:3, 1:4] == [1 0 0 0
                                0 1 0 0
                                0 0 1 0]
    end

    @testset "nesting" begin
        ff = apply_schema(@formula(y ~ 1 + (1|g/f)), schema(dat), MixedModel)
        @test modelcols(last(ff.rhs), dat) == float(Matrix(I, 18, 18))

        # in fixed effects:
        d2 = (a = rand(20), b = repeat([:X, :Y], outer=10), c = repeat([:S,:T],outer=10))
        f2 = apply_schema(@formula(0 ~ 1 + b/a), schema(d2), MixedModel)
        @test modelcols(f2.rhs, d2) == [ones(20) d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
        @test coefnames(f2.rhs) == ["(Intercept)", "b: Y", "b: X & a", "b: Y & a"]

        # runtime
        fr2 = term(1) + term(:b) / term(:a)
        @test length(fr2) == 3
        @test fr2[1] isa ConstantTerm
        @test fr2[2] isa Term
        @test fr2[3] isa InteractionTerm
        frf2 = apply_schema(term(0) ~ fr2, schema(d2), MixedModel)
        @test modelcols(frf2.rhs, d2) == [ones(20) d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
        @test coefnames(frf2.rhs) == ["(Intercept)", "b: Y", "b: X & a", "b: Y & a"]

        # check promotion
        f3 = apply_schema(@formula(0 ~ 0 + b/a), schema(d2), MixedModel)
        @test modelcols(f3.rhs, d2) == [d2.b .== :X d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
        @test coefnames(f3.rhs) == ["b: X", "b: Y", "b: X & a", "b: Y & a"]

        # runtime:
        fr3 = term(0) + term(:b) / term(:a)
        @test length(fr3) == 3
        @test fr3[1] isa ConstantTerm
        @test fr3[2] isa Term
        @test fr3[3] isa InteractionTerm
        ffr3 = apply_schema(term(0) ~ fr3, schema(d2), MixedModel)
        @test modelcols(ffr3.rhs, d2) == [d2.b .== :X d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
        @test coefnames(ffr3.rhs) == ["b: X", "b: Y", "b: X & a", "b: Y & a"]

        # errors for continuous grouping
        @test_throws ArgumentError apply_schema(@formula(0 ~ 1 + a/b), schema(d2), MixedModel)

        # errors for too much nesting
        @test_throws ArgumentError apply_schema(@formula(0 ~ 1 + b/c/a), schema(d2), MixedModel)

        # fitted model to test amalgamate and fnames, and equivalence with other formulations
        psts = dataset("pastes")
        m = fit(MixedModel, @formula(strength ~ 1 + (1|batch/cask)), psts; progress=false)
        m2 = fit(MixedModel, @formula(strength ~ 1 + (1|batch) + (1|batch&cask)), psts; progress=false)
        mr = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)/term(:cask)), psts; progress=false)
        m2r = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)) + (term(1)|term(:batch)&term(:cask)), psts; progress=false)

        @test fnames(m) == fnames(m2) == fnames(mr) == fnames(m2r) == (Symbol("batch & cask"), :batch)
        @test m.λ == m2.λ == mr.λ == m2r.λ
        @test deviance(m) == deviance(m2) == deviance(mr) == deviance(m2r)
    end

    @testset "multiple terms with same grouping" begin
        dat = MixedModels.dataset(:kb07)
        sch = schema(dat)
        f1 = @formula(rt_trunc ~ 1 + (1 + prec + load | spkr))
        ff1 = apply_schema(f1, sch, MixedModel)

        retrm = last(ff1.rhs)
        @test last(retrm.lhs.terms).contrasts.contrasts isa DummyCoding

        f2 = @formula(rt_trunc ~ 1 + (1 + prec | spkr) + (0 + load | spkr))
        ff2 = apply_schema(f2, sch, MixedModel)

        retrm2 = last(ff2.rhs)
        @test last(retrm2.lhs.terms).contrasts.contrasts isa DummyCoding
    end
end


text:
"""
The basic voting function when every agent's got an equal voting power
and the candidate is being compared to the registry mean
"""
function simpleVote(benchmark, candidate, agents)
    rejectAgents = []
    proAgents = []
    quorum = length(agents)
    for agent in agents
        push!(evaluateCandidateByAgent(candidate, agent) >= benchmark ? proAgents : rejectAgents, agent)
    end
    [length(proAgents) > quorum/2, proAgents, rejectAgents]
end

"""
The voting function when every agent's got an equal voting power
but only agents who hold some tokens are eligible to participate
"""
function binaryTokenVote(benchmark, candidate, agents)
    rejectAgents = []
    proAgents = []
    quorum = 0
    for agent in agents
        if agent.balance > 0
            if (evaluateCandidateByAgent(candidate, agent) <= benchmark)
                push!(rejectAgents, agent)
            else
                push!(proAgents, agent)
            end
            quorum += 1
        end
    end
    [length(proAgents) > quorum/2, proAgents, rejectAgents]
end

"""
The voting function when every agent's got a voting power
proportional to the amount of tokens under possesion
"""
function proRataTokenVote(benchmark, candidate, agents)
    rejectAgents = []
    proAgents = []
    proTokens = 0
    quorum = 0
    for agent in agents
        if agent.balance > 0
            if (evaluateCandidateByAgent(candidate, agent) <= benchmark)
                push!(rejectAgents, agent)
            else
                proTokens += agent.balance
                push!(proAgents, agent)
            end
            quorum += agent.balance
        end
    end
    [proTokens > quorum/2, proAgents, rejectAgents]
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
### Right-hand side and solution vectors

const SaddleVector = ArrayPartition

"""
    SaddleVector(u,f)

Construct a vector of a state part `u` and constraint part `f` of a
saddle-point vector, to be associated with a [`SaddleSystem`](@ref).
"""
function SaddleVector end

"""
    solvector([;state=][,constraint=][,aux_state=])

Build a solution vector for a constrained system. This takes three optional keyword
arguments: `state`, `constraint`, and `aux_state`. If only a state is supplied,
then the constraint is set to an empty vector and the system is assumed to
correspond to an unconstrained system. (`aux_state` is ignored in this situation.)
"""
solvector(;state=nothing,constraint=nothing,aux_state=nothing) = _solvector(state,constraint,aux_state)
_solvector(::Nothing,::Nothing,::Nothing) = nothing
_solvector(s,::Nothing,::Nothing) = ArrayPartition(s,_empty(s))
_solvector(s,::Nothing,aux) = ArrayPartition(s,_empty(s))
_solvector(s,c,::Nothing) = ArrayPartition(s,c)
#_solvector(s,c,aux) = ArrayPartition(_solvector(s,c,nothing),aux)
_solvector(s,c,aux) = ArrayPartition(s,c,aux)

mainvector(u) = u
mainvector(u::ArrayPartition) = ArrayPartition(u.x[1],u.x[2])

_empty(s) = Vector{eltype(s)}()

"""
    state(x::SaddleVector)

Provide the state part of the given saddle vector `x`
"""
state(u::ArrayPartition) = mainvector(u).x[1]
state(u) = u

"""
    constraint(x::SaddleVector)

Provide the constraint part of the given saddle vector `x`
"""
constraint(u::ArrayPartition) = mainvector(u).x[2]
constraint(u) = eltype(u)[]

"""
    aux_state(x)

Provide the auxiliary state part of the given vector `x`
"""
aux_state(u) = nothing # Array{eltype(u)}(undef,0,0)
aux_state(u::ArrayPartition{T,Tuple{F1,F2,F3}}) where {T,F1,F2,F3} = u.x[3]

for f in (:state,:constraint,:aux_state)
  @eval $f(a::AbstractArray{T}) where {T<:ArrayPartition} = map($f,a)
end

#SaddleVector(u::TU,f::TF) where {TU,TF} = ArrayPartition(u,f)

"""
    r1vector([;state_r1=][,aux_r1=])

Build a vector of the `r1` functions for the state ODEs and auxiliary state ODEs.
"""
r1vector(;state_r1=nothing,aux_r1=nothing) = _r1vector(state_r1,aux_r1)
_r1vector(::Nothing,::Nothing) = nothing
_r1vector(s,::Nothing) = s
_r1vector(::Nothing,a) = nothing
_r1vector(s,a) = ArrayPartition((s,a))

hasaux(r1) = false
hasaux(r1::ArrayPartition) = true

state_r1(r1) = r1
state_r1(r1::ArrayPartition) = r1.x[1]

aux_r1(r1) = nothing
aux_r1(r1::ArrayPartition) = r1.x[2]


text:
using Test
using DirectedAcyclicGraphs

module TestDAGs

    using Test
    using ..DirectedAcyclicGraphs

    mutable struct TestINode <: DAG
        id::Int
        children::Vector{DAG}
        TestINode(i,c) = new(i,c)
    end

    mutable struct TestLNode <: DAG
        id::Int
        TestLNode(i) = new(i)
    end

    DirectedAcyclicGraphs.NodeType(::Type{<:TestINode}) = Inner()
    DirectedAcyclicGraphs.NodeType(::Type{<:TestLNode}) = Leaf()
    DirectedAcyclicGraphs.children(n::TestINode) = n.children

    @testset "DAG utils for TestNodes" begin

        l1 = TestLNode(1)
        l2 = TestLNode(2)

        @test !has_children(l1)
        @test num_children(l1) == 0
        @test isleaf(l1)
        @test !isinner(l1)

        i1 = TestINode(1,[l1])
        i2 = TestINode(2,[l2])
        i12 = TestINode(3,[l1,l2])

        @test has_children(i1)
        @test has_children(i12)

        @test num_children(i1) == 1
        @test num_children(i12) == 2

        j1 = TestINode(1,[i1,i12])
        j2 = TestINode(2,[i2])
        j12 = TestINode(3,[i1,i2])

        r = TestINode(5,[j1,j2,j12])

        @test has_children(r)
        @test num_children(r) == 3

        foreach(r) do n
            n.id += 1
        end
        @test l1.id == 2
        @test l2.id == 3
        @test i12.id == 4
        @test j2.id == 3
        @test r.id == 6

        foreach(r, l -> l.id += 1, i -> i.id -= 1)
        @test l1.id == 2+1
        @test l2.id == 3+1
        @test i12.id == 4-1
        @test j2.id == 3-1
        @test r.id == 6-1

        @test filter(n -> iseven(n.id), r) == [l2,i2,j2]

        lastvisited = nothing
        foreach(n -> lastvisited=n,r)
        @test lastvisited === r

        lastvisited = nothing
        foreach_down(n -> lastvisited=n,r)
        @test isleaf(lastvisited)

        @test num_nodes(r) == 9
        @test num_edges(r) == 12

        @test isempty(innernodes(l1))
        @test leafnodes(l1) == [l1]

        @test issetequal(innernodes(r), [i1,i2,i12,j1,j2,j12,r])
        @test issetequal(innernodes(r), [i1,i2,i12,j1,j2,j12,r])
        @test issetequal(leafnodes(r), [l1,l2])

        @test tree_num_edges(r) == 14 # unverified

        @test linearize(r)[end] == r
        @test linearize(r)[1] == l1 || linearize(r)[1] == l2
        @test linearize(l2) == [l2]
        @test length(linearize(i12)) == 3

        @test eltype(linearize(r)) == DAG
        @test eltype(linearize(l1)) == TestLNode
        @test eltype(linearize(r, Any)) == Any

        @test left_most_descendent(r) == l1
        @test right_most_descendent(r) == l2

        np = num_parents(r)
        @test np[r] == 0
        @test np[l1] == 2
        @test np[j1] == 1

        f_l(n) = 1
        f_i1(n,cs) = sum(cs) + 1
        f_i2(n,call) = mapreduce(call, +, children(n)) + 1

        @test foldup(r, f_l, f_i2, Int) == foldup_aggregate(r, f_l, f_i1, Int)

        @test tree_num_nodes(r) == 15
        @test tree_num_edges(r) == 14
        @test num_leafnodes(r) == 2
        @test num_innernodes(r) == 7

        @test j1 ∈ r
        @test TestLNode(24) ∉ r

        @test left_most_descendent(r) == l1
        @test right_most_descendent(r) == l2

        @test label_nodes(r)[r] == 9

        @test !isempty(node_stats(r))
        @test !isempty(parent_stats(r))

        layers, num_layers = @test_nowarn feedforward_layers(r)
        @test num_layers == 4
        @test layers[r] == 4
        @test layers[l2] == 1
        @test layers[i12] == 2
        @test layers[j1] == 3

        @test maximum(num_children, r) == 3
        @test minimum(num_children, r) == 0

    end

end


text:
#=
ramachandranComp:
- Julia version:
- Author: Dan
- Date: 2021-11-28
=#

push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools

pdb_filepath = "1AKE.pdb"

comp = System(pdb_filepath)

ramachandranangles(comp)

elapsed = @elapsed ramachandranangles(comp)

println(elapsed)
@show @benchmark ramachandranangles(comp)