text:
# The solver type
"""
    ValueIterationSolver <: Solver

The solver type. Contains the following parameters that can be passed as keyword arguments to the constructor

    - max_iterations::Int64, the maximum number of iterations value iteration runs for (default 100)
    - belres::Float64, the Bellman residual (default 1e-3)
    - verbose::Bool, if set to true, the bellman residual and the time per iteration will be printed to STDOUT (default false)
    - include_Q::Bool, if set to true, the solver outputs the Q values in addition to the utility and the policy (default true)
    - init_util::Vector{Float64}, provides a custom initialization of the utility vector. (initializes utility to 0 by default)
"""
mutable struct ValueIterationSolver <: Solver
    max_iterations::Int64 # max number of iterations
    belres::Float64 # the Bellman Residual
    verbose::Bool
    include_Q::Bool
    init_util::Vector{Float64}
end
# Default constructor
function ValueIterationSolver(;max_iterations::Int64 = 100,
                               belres::Float64 = 1e-3,
                               verbose::Bool = false,
                               include_Q::Bool = true,
                               init_util::Vector{Float64}=Vector{Float64}(undef, 0))
    return ValueIterationSolver(max_iterations, belres, verbose, include_Q, init_util)
end

@POMDP_require solve(solver::ValueIterationSolver, mdp::Union{MDP,POMDP}) begin
    P = typeof(mdp)
    S = statetype(P)
    A = actiontype(P)
    @req discount(::P)
    @subreq ordered_states(mdp)
    @subreq ordered_actions(mdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A,::S)
    @req stateindex(::P,::S)
    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    as = actions(mdp)
    ss = states(mdp)
    @req length(::typeof(ss))
    @req length(::typeof(as))
    a = first(as)
    s = first(ss)
    dist = transition(mdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)
end

#####################################################################
# Solve runs the value iteration algorithm.
# The policy input argument is either provided by the user or
# initialized during the function call.
# Verbose is a flag that triggers text output to the command line
# Example code for running the function:
# mdp = GridWorld(10, 10) # initialize a 10x10 grid world MDP (user written code)
# solver = ValueIterationSolver(max_iterations=40, belres=1e-3)
# policy = ValueIterationPolicy(mdp)
# solve(solver, mdp, policy, verbose=true)
#####################################################################
function solve(solver::ValueIterationSolver, mdp::MDP; kwargs...)

    # deprecation warning - can be removed when Julia 1.0 is adopted
    if !isempty(kwargs)
        @warn("Keyword args for solve(::ValueIterationSolver, ::MDP) are no longer supported. For verbose output, use the verbose option in the ValueIterationSolver")
    end

    @warn_requirements solve(solver, mdp)

    # solver parameters
    max_iterations = solver.max_iterations
    belres = solver.belres
    discount_factor = discount(mdp)
    ns = length(states(mdp))
    na = length(actions(mdp))

    # intialize the utility and Q-matrix
    if !isempty(solver.init_util)
        @assert length(solver.init_util) == ns "Input utility dimension mismatch"
        util = solver.init_util
    else
        util = zeros(ns)
    end
    include_Q = solver.include_Q
    if include_Q
        qmat = zeros(ns, na)
    end
    pol = zeros(Int64, ns)

    total_time = 0.0
    iter_time = 0.0

    # create an ordered list of states for fast iteration
    state_space = ordered_states(mdp)

    # main loop
    for i = 1:max_iterations
        residual = 0.0
        iter_time = @elapsed begin
        # state loop
        for (istate,s) in enumerate(state_space)
            sub_aspace = actions(mdp, s)
            if isterminal(mdp, s)
                util[istate] = 0.0
                pol[istate] = 1
            else
                old_util = util[istate] # for residual
                max_util = -Inf
                # action loop
                # util(s) = max_a( R(s,a) + discount_factor * sum(T(s'|s,a)util(s') )
                for a in sub_aspace
                    iaction = actionindex(mdp, a)
                    dist = transition(mdp, s, a) # creates distribution over neighbors
                    u = 0.0
                    for (sp, p) in weighted_iterator(dist)
                        p == 0.0 ? continue : nothing # skip if zero prob
                        r = reward(mdp, s, a, sp)
                        isp = stateindex(mdp, sp)
                        u += p * (r + discount_factor * util[isp])
                    end
                    new_util = u
                    if new_util > max_util
                        max_util = new_util
                        pol[istate] = iaction
                    end
                    include_Q ? (qmat[istate, iaction] = new_util) : nothing
                end # action
                # update the value array
                util[istate] = max_util
                diff = abs(max_util - old_util)
                diff > residual ? (residual = diff) : nothing
            end
        end # state
        end # time
        total_time += iter_time
        solver.verbose ? @printf("[Iteration %-4d] residual: %10.3G | iteration runtime: %10.3f ms, (%10.3G s total)\n", i, residual, iter_time*1000.0, total_time) : nothing
        residual < belres ? break : nothing
    end # main
    if include_Q
        return ValueIterationPolicy(mdp, qmat, util, pol)
    else
        return ValueIterationPolicy(mdp, utility=util, policy=pol, include_Q=false)
    end
end

function solve(::ValueIterationSolver, ::POMDP)
    throw("""
           ValueIterationError: `solve(::ValueIterationSolver, ::POMDP)` is not supported,
          `ValueIterationSolver` supports MDP models only, look at QMDP.jl for a POMDP solver that assumes full observability.
           If you still wish to use the transition and reward from your POMDP model you can use the `UnderlyingMDP` wrapper from POMDPModelTools.jl as follows:
           ```
           solver = ValueIterationSolver()
           mdp = UnderlyingMDP(pomdp)
           solve(solver, mdp)
           ```
           """)
end


text:
@inline function LocalRegistry.find_package_path(package::ClonedPackage)
    return package.path
end

@inline function LocalRegistry.find_registry_path(registry::ClonedRegistry)
    return registry.path
end

ext: jl


text:
@doc raw"""
    K, γ, info = glover_mcfarlane(G::AbstractStateSpace, γ = 1.1; W1=1, W2=1)

Design a controller for `G` that maximizes the stability margin ϵ = 1/γ with normalized coprime factor uncertainty using the method of Glover and McFarlane
```
γ = 1/ϵ = ||[K;I] inv(I-G*K)*inv(M)||∞
G = inv(M + ΔM)*(N + ΔN)
```
γ is given as a relative factor above γmin and must be greater than 1, i.e., if γ = 1.1, the controller will be designed for γ = 1.1*γmin.

We want γmin (which is always ≥ 1) as small as possible, and we usually require that γmin is less than 4, corresponding to 25% allowed coprime uncertainty.

Performance modeling is incorporated in the design by calling `glover_mcfarlane` on the shaped system `Gs = W2*G*W1` and then forming the controller as `K = W1*Ks*W2`. Using this formulation, traditional loop shaping can be done on `Gs = W2*G*W1`. The plant shaping is handled internally if keyword arguments `W1, W2` are used and the returned controller is already scaled. In this case, `Gs` and `Ks` are included in the `info` named tuple for inspection.

See also [`glover_mcfarlane_2dof`](@ref) to design a feedforward filter as well.

# Example:
Example 9.3 from the reference below.
```julia
using RobustAndOptimalControl, ControlSystems, Plots, Test
G = tf(200, [10, 1])*tf(1, [0.05, 1])^2     |> ss
Gd = tf(100, [10, 1])                       |> ss
W1 = tf([1, 2], [1, 1e-6])                  |> ss
K, γ, info = glover_mcfarlane(G, 1.1; W1)
@test info.γmin ≈ 2.34 atol=0.005
Gcl = extended_gangoffour(G, K) # Form closed-loop system

bodeplot([G, info.Gs, G*K], lab=["G" "" "G scaled" "" "Loop transfer"]) |> display
bodeplot(Gcl, lab=["S" "KS" "PS" "T"], plotphase=false) |> display # Plot gang of four

plot( step(Gd*feedback(1, info.Gs), 3), lab="Initial controller")
plot!(step(Gd*feedback(1, G*K), 3), lab="Robustified") |> display

nyquistplot([info.Gs, G*K], ylims=(-2,1), xlims=(-2, 1),
    Ms_circles = 1.5,
    lab = ["Initial controller" "Robustified"],
    title = "Loop transfers with and without robustified controller"
    ) |> display
```

Ref: Sec 9.4.1 of Skogestad, "Multivariable Feedback Control: Analysis and Design"

# Extended help

Skogestad gives the following general advice:
1. Scale the plant outputs and inputs. This is very important for most design
    procedures and is sometimes forgotten. In general, scaling improves the
    conditioning of the design problem, it enables meaningful analysis to be made
    of the robustness properties of the feedback system in the frequency domain,
    and for loop-shaping it can simplify the selection of weights. There are a variety
    of methods available including normalization with respect to the magnitude of
    the maximum or average value of the signal in question. If one is to go straight to a design the following variation has
    proved useful in practice:
    (a) The outputs are scaled such that equal magnitudes of cross-coupling into each
        of the outputs is equally undesirable.
    (b) Each input is scaled by a given percentage (say 10%) of its expected range
        of operation. That is, the inputs are scaled to reflect the relative actuator
        capabilities. An example of this type of scaling is given in the aero-engine
        case study of Chapter 12.
2. Order the inputs and outputs so that the plant is as diagonal as possible. The
    relative gain array [`rga`](@ref) can be useful here. The purpose of this pseudo-diagonalization
    is to ease the design of the pre- and post-compensators which, for simplicity, will
    be chosen to be diagonal.

    Next, we discuss the selection of weights to obtain the shaped plant $G_s = W_2 G W_1$
    where $W_1 = W_p W_a W_g$
3. Select the elements of diagonal pre- and post-compensators $W_p$ and $W_2$ so that
    the singular values of $W_2 G W_p$ are desirable. This would normally mean high
    gain at low frequencies, roll-off rates of approximately 20 dB/decade (a slope of
    about 1) at the desired bandwidth(s), with higher rates at high frequencies. Some
    trial and error is involved here. $W_2$ is usually chosen as a constant, reflecting the
    relative importance of the outputs to be controlled and the other measurements
    being fed back to the controller. For example, if there are feedback measurements
    of two outputs to be controlled and a velocity signal, then $W_2$ might be chosen
    to be `diag([1, 1, 0.1])`, where 0.1 is in the velocity signal channel. $W_p$ contains the
    dynamic shaping. Integral action, for low frequency performance; phase-advance
    for reducing the roll-off rates at crossover, and phase-lag to increase the roll-off
    rates at high frequencies should all be placed in $W_p$ if desired. The weights should
    be chosen so that no unstable hidden modes are created in $G_s$.
5. Optional: Introduce an additional gain matrix $W_g$ cascaded with $W_a$ to provide
    control over actuator usage. $W_g$ is diagonal and is adjusted so that actuator rate
    limits are not exceeded for reference demands and typical disturbances on the
    scaled plant outputs. This requires some trial and error.

6. Robustly stabilize the shaped plant $G_s = W_2 G W_1$ , where $W_1 = W_p W_a W_g$,
    using `glover_mcfarlane`. First, the maximum stability
    margin $ϵ_{max} = 1/γ_{min}$ is calculated. If the margin is too small, $ϵmax < 0.25$, then go back and modify the weights. Otherwise, a γ-suboptimal controller is synthesized. There is usually no advantage to be gained by using the optimal controller. When $ϵ_{max}$ > 0.25
    (respectively $γ_{min}$ < 4) the design is usually successful. In this case, at least
    25% coprime factor uncertainty is allowed, and we also find that the shape of the
    open-loop singular values will not have changed much after robust stabilization.
    A small value of ϵmax indicates that the chosen singular value loop-shapes are
    incompatible with robust stability requirements. That the loop-shapes do not
    change much following robust stabilization if γ is small (ϵ large), is justified
    theoretically in McFarlane and Glover (1990).

7. Analyze the design and if all the specifications are not met make further
    modifications to the weights.
8. Implement the controller. The configuration shown in below has been found
    useful when compared with the conventional set up. This is because
    the references do not directly excite the dynamics of $K$, which can result in large amounts of overshoot (classical derivative kick). The constant prefilter ensures a steady-state gain of 1 between r and y, assuming integral action in $W_1$ or $G$ (note, the K returned by this function has opposite sign compared to that of Skogestad, so we use negative feedback here).

Anti-windup can be added to $W_1$ but putting $W_1$ on Hanus form after the synthesis, see [`hanus`](@ref).

```
       ┌─────────┐      ┌────────┐      ┌────────┐
    r  │         │    us│        │  u   │        │  y
   ───►│(K*W2)(0)├──+──►│   W1   ├─────►│   G    ├────┬──►
       │         │  │-  │        │      │        │    │
       └─────────┘  │   └────────┘      └────────┘    │
                    │                                 │
                    │                                 │
                    │   ┌────────┐      ┌────────┐    │
                    │   │        │  ys  │        │    │
                    └───┤   K    │◄─────┤   W2   │◄───┘
                        │        │      │        │
                        └────────┘      └────────┘
```
Keywords: nfcsyn, coprimeunc
"""
function glover_mcfarlane(G::AbstractStateSpace{Continuous}, γ = 1.1; W1=1, W2=1)
    γ > 1 || throw(ArgumentError("γ must be greater than 1"))
    Gs = W2*G*W1
    A,B,C,D = ssdata(Gs)

    R = I + D*D'
    S = I + D'D
    Sl = lu(S)
    # arec(A, B, R, Q, S) solves A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0

    Ā = A - B*(S\D'C)
    Z,_ = MatrixEquations.arec(Ā', C', R, B*(Sl\B'))
    X,_ = MatrixEquations.arec(Ā, B, S, C'*(R\C))

    γmin = sqrt(1 + ρ(X*Z))

    γ *= γmin

    L = (1-γ^2)*I + X*Z
    F = -(Sl\(D'C + B'X))
    BK = γ^2*(L'\Z)*C'
    AK = A + B*F + BK*(C + D*F)
    CK = B'X
    DK = -D'
    Ks = -ss(AK, BK, CK, DK)
    Gcl = extended_gangoffour(Gs, Ks)
    imargin, ω = hinfnorm2(Gcl)
    K = W1*Ks*W2
    Gcl = extended_gangoffour(G, K)
    K, γ, (; Gcl, margin = inv(imargin), ω, γmin, Ks, Gs, Z, X)
end

"""
    K, γ, info = glover_mcfarlane(G::AbstractStateSpace{<:Discrete}, γ = 1.1; W1=1, W2=1, strictly_proper=false)

For discrete systems, the `info` tuple contains also feedback gains `F, L` and observer gain `Hkf` such that the controller on observer form is given by
```math
x′ = Ax + Bu + H_{kf}*(Cx - y)\\\\
u = Fx + L*(Cx - y)
```
Note, this controller is *not* strictly proper, i.e., it has a non-zero D matrix.
The controller can be transformed to observer form for the scaled plant (`info.Gs`)
by `Ko = observer_controller(info)`, in which case the following holds `G*K == info.Gs*Ko` (realizations are different).

If `strictly_proper = true`, the returned controller `K` will have `D == 0`.
This can be advantageous in implementations where computational delays are present. In this case, `info.L == 0` as well.

Ref discrete version: Iglesias, "The Strictly Proper Discrete-Time Controller for the Normalized Left-Coprime Factorization Robust Stabilization Problem"
"""
function glover_mcfarlane(G::AbstractStateSpace{<:Discrete}, γ = 1.1; W1=1, W2=1, strictly_proper=false)
    γ > 1 || throw(ArgumentError("γ must be greater than 1"))
    Gs = W2*G*W1
    A,B,C,D = ssdata(Gs)
    iszero(D) || throw(ArgumentError("System must be strictly proper (D must be 0)"))

    X,_,Flq = MatrixEquations.ared(A, B, I, C'C)
    Z,_ = MatrixEquations.ared(A', C', I, B*B')
    Q1 = I + C*Z*C'

    if strictly_proper
        T = Q1 # They appear to be the same
        T12 = sqrt(T)
        X12 = sqrt(X)
        M1 = T12\C*Z*A'X12
        Q2 = I + X12*Z*X12

        QM = [
            sqrt(Q1 + 1/4 * M1*M1') -1/2*M1
            -1/2*M1'       sqrt(Q2 + 1/4 * M1'*M1)
        ]
        γmin = ρ(QM)
        γ *= γmin

        W = (γ^2 - 1)*I - Z*X
        XW = X/W
        Hkf = -A*Z*C'/Q1

        # Flq is the feedback gain from the ARE above, ref Rowe and Maciejowski, "Tuning MPC using H∞ Loop Shaping" between eq 36-37.
        # Flq = -B'X*((I + B*B'X)\A) # NOTE: Flq might be required for inverse optimal control
        F = -γ^2*Flq/W # Original paper
        # F = -γ^2*B'XW*((I + γ^2*B*B'XW)\A) # Merl paper
        # @show [F; Flq; F0]
        Ak = A + B*F + Hkf*C
        # NOTE: It's unclear what the best realization is. Both papers talk about a "dual" realization, but they appear to disagree on the formulation, and none of the expressions for the Q-matrices, in the Merl paper or the Maciejowski paper, do not lead to the same caluclated feedback gain. I think they might have left out a similarity transform on some matrices. The equations used here are thus not any of the dual forms, even if those appear to be recommended.

        L = 0
        Ck = F
        Bk = -Hkf
        Dk = 0
    else

        γmin = sqrt(1 + ρ(X*Z))
        γ *= γmin

        W = (γ^2 - 1)*I - Z*X
        Hkf = -A*Z*C'/Q1
        Akf = A+Hkf*C

        XW = X/W
        hest = lu(I + γ^2*B*B'XW)
        γ2B = γ^2*B'
        arne = γ2B*XW
        Ak = hest\Akf
        Ck = arne*Ak
        Bk = hest\Hkf
        Dk = arne*Bk

        F0 = -γ2B*XW/hest
        F = F0*A
        L = F0*Hkf
    end

    Ks = -ss(Ak, Bk, Ck, Dk, G.timeevol)
    Gcl = extended_gangoffour(Gs, Ks)
    imargin, ω = hinfnorm2(Gcl)
    K = W1*Ks*W2
    Gcl = extended_gangoffour(G, K)
    K, γ, (; Gcl, margin = inv(imargin), ω, γmin, Ks, Gs, Z, X, F, L, Hkf, W)
end

"""
    observer_controller(glover_mcfarlane_info::NamedTuple)

Return a controller on observer form (observer of the scaled plant `info.Gs`).
The observer controller has input vector `y`.

Ref: eq (2.5)-(2.6) of Iglesias, "The Strictly Proper Discrete-Time Controller for the Normalized Left-Coprime Factorization Robust Stabilization Problem"
"""
function ControlSystems.observer_controller(info::NamedTuple)
    isdiscrete(info.Gs) || throw(ArgumentError("Observer controller can only be generated for Glover McFarlane designs on discrete systems."))
    A,B,C,D = ssdata(info.Gs)
    iszero(D) || throw(ArgumentError("observer_controller does not support non-zero D matrix")) # not sure if this is a strict limitation or the paper author simplified.
    H, L, F = info.Hkf, info.L, info.F
    Ao = A + B*(F + L*C) + H*C
    Bo = (B*L+H)
    Co = F + L*C
    Do = L
    ss(Ao, Bo, Co, Do, info.Gs.timeevol)
end

"""
    observer_predictor(glover_mcfarlane_info::NamedTuple)

Return a predictor for the scaled plant `info.Gs`.
The observer predictor has input vector `[u; y]`.

Ref: eq (2.5)-(2.6) of Iglesias, "The Strictly Proper Discrete-Time Controller for the Normalized Left-Coprime Factorization Robust Stabilization Problem"
"""
function ControlSystems.observer_predictor(info::NamedTuple)
    isdiscrete(info.Gs) || throw(ArgumentError("Observer predictor can only be generated for Glover McFarlane designs on discrete systems."))
    A,B,C,D = ssdata(info.Gs)
    H, L, F = info.Hkf, info.L, info.F
    Ao = A + H*C
    Bo = [B-H*D -H]
    Co = C
    Do = [D zeros(size(D,1), size(H, 2))]
    ss(Ao, Bo, Co, Do, info.Gs.timeevol)
end

"""
    K, γ, info = glover_mcfarlane_2dof(G::AbstractStateSpace{Continuous}, Tref::AbstractStateSpace{Continuous}, γ = 1.1, ρ = 1.1;
    W1 = 1, Wo = I, match_dc = true, kwargs...)

Joint design of feedback and feedforward compensators
```math
K = \\left\\[K_1 & K_2\\right\\]
```
```
   ┌──────┐   ┌──────┐        ┌──────┐    ┌─────┐
r  │      │   │      │        │      │    │     │
──►│  Wi  ├──►│  K1  ├───+───►│  W1  ├───►│  G  ├─┐y
   │      │   │      │   │    │      │    │     │ │
   └──────┘   └──────┘   │    └──────┘    └─────┘ │
                         │                        │
                         │    ┌──────┐            │
                         │    │      │            │
                         └────┤  K2  ◄────────────┘
                              │      │
                              └──────┘
```
Where the returned controller `K` takes the measurement vector `[r; y]` (positive feedback),
i.e., it includes all blocks `Wi, K1, K2, W1`.
If `match_dc = true`, `Wi` is automatically computed to make sure the static gain matches `Tref` exactly, otherwise `Wi` is set to `I`.
The `info` named tuple contains the feedforward filter for inspection (`info.K1 = K1*Wi`).

# Arguments:
- `G`: Plant model
- `Tref`: Reference model
- `γ`: Relative γ
- `ρ`: Design parameter, typically 1 < ρ < 3. Increase to emphasize model matching at the expense of robustness.
- `W1`: Pre-compensator for loop shaping.
- `Wo`: Output selction matrix. If there are more measurements than controlled variables, this matrix let's you select which measurements are to be controlled.
- `kwargs`: Are sent to [`hinfsynthesize`](@ref).

Ref: Sec. 9.4.3 of Skogestad, "Multivariable Feedback Control: Analysis and Design".
The reference contains valuable pointers regarding gain-scheduling implementation of the designed controller as an observer with feedback from estimated states.
In order to get anti-windup protection when `W1` contains an integrator,
transform `W1` to self-conditioned Hanus form (using [`hanus`](@ref)) and implement the controller like this
```julia
W1h = hanus(W1)             # Perform outside loop

# Each iteration
us = filter(Ks, [r; y])     # filter inputs through info.Ks (filter is a fictive function that applies the transfer function)
u  = filter(W1h, [us; ua])  # filter us and u-actual (after input saturation) through W1h
ua = clamp(u, lower, upper) # Calculate ua for next iteration as the saturated value of u
```

# Example:
```julia
P = tf([1, 5], [1, 2, 10]) # Plant
W1 = tf(1,[1, 0]) |> ss    # Loop shaping controller

Tref = tf(1, [1, 1]) |> ss # Reference model

K1dof, γ1, info1 = glover_mcfarlane(ss(P), 1.1; W1)
K2dof, γ2, info2 = glover_mcfarlane_2dof(ss(P), Tref, 1.1, 1.1; W1)

G1 = feedback(P*K1dof)
G2 = info2.Gcl

bodeplot(info2.K1, w, lab="Feedforward filter")
plot([step(G1, 15), step(G2, 15), step(Tref, 15)], lab=["1-DOF" "2-DOF" "Tref"])
```
"""
function glover_mcfarlane_2dof(G::AbstractStateSpace{Continuous}, Tref::AbstractStateSpace{Continuous}, γ = 1.1, ρ = 1.1; W1=1, Wo = I, match_dc = true, kwargs...)
    γ > 1 || throw(ArgumentError("γ must be greater than 1"))
    ρ > 1 || throw(ArgumentError("ρ must be greater than 1"))
    Gs = G*W1
    As,Bs,Cs,Ds = ssdata(Gs)
    Ar,Br,Cr,Dr = ssdata(Tref)
    nr,nr = size(Ar)
    lr,mr = size(Dr)
    ns,ns = size(As)
    ls,ms = size(Ds)
    Rs = I + Ds*Ds'
    sRs = sqrt(Rs) # Example in Skogestad uses matlab sqrt which is elementwise
    Ss = lu!(I + Ds'Ds)
    A1 = (As - Bs*(Ss\Ds'Cs))
    R1 = Cs'*(Rs\Cs)
    Q1 = Bs*(Ss\Bs')
    Zs, _ = MatrixEquations.arec(A1', R1, Q1)

    A = cat(As, Ar, dims=(1,2))
    B1 = [
        zeros(ns,mr) ((Bs*Ds')+(Zs*Cs'))/sRs
        Br zeros(nr,ls)
    ]
    B2 = [Bs; zeros(nr,ms)]
    C1 = [zeros(ms,ns+nr); Cs zeros(ls,nr); ρ*Wo*Cs -ρ^2*Cr]
    C2 = [zeros(mr,ns+nr); Cs zeros(ls,nr)]
    D11 = [zeros(ms,mr+ls); zeros(ls,mr) sRs;-ρ^2*Dr ρ*Wo*sRs]
    D12 = [I(ms);Ds;ρ*Ds]
    D21 = [ρ*I(mr) zeros(mr,ls);zeros(ls,mr) sRs]
    D22 = [zeros(mr,ms); Ds]
    P = ss(A, B1, B2, C1, C2, D11, D12, D21, D22)
    _, Ks, γopt = hinfsynthesize(P, γrel = γ; kwargs...)

    u1 = 1:ms
    K1 = Ks[:, u1]
    if match_dc
        K2 = Ks[:, ms+1:end]
        sens = output_sensitivity(Gs, -K2)*Gs*K1 # eq. 9.89
        Wi = (Wo*dcgain(sens, 1e-6))\dcgain(Tref)
        K1 = K1*Wi
        Ks.B[:,u1] .= K1.B
        Ks.D[:,u1] .= K2.D
    else
        Wi = I
    end
    K = W1*Ks
    Gcl = feedback(G*K, ss(1), W1=1:ms, U1=ms+1:2ms, pos_feedback=true)
    info = (; Gcl, Ks, Gs, Wi, K1)
    K, γopt, info
end

"Spectral radius"
function ρ(X)
    e = eigvals(X)
    abs(e[end])
end

"""
    Wh = hanus(W)

Return `Wh` on Hanus form. `Wh` has twice the number of inputs, where the second half of the inputs are "actual inputs", e.g., potentially saturated. This is used to endow `W` with anti-windup protection.
`W` must have an invertable `D` matrix and be minimum phase.

Ref: Sec 9.4.5 of Skogestad, "Multivariable Feedback Control: Analysis and Design"
"""
function hanus(W)
    A,B,C,D = ssdata(W)
    nu = W.nu
    BD = B/D
    A2 = A - BD*C
    B2 = [0*I(nu) BD]
    D2 = [D 0*I(nu)]
    ss(A2, B2, C, D2)
end

"""
    extended_gangoffour(P, C)

Returns a single statespace system that maps
- `w1` reference or measurement noise
- `w2` load disturbance
to
- `z1` control error
- `z2` control input
```
      z1          z2
      ▲  ┌─────┐  ▲      ┌─────┐
      │  │     │  │      │     │
w1──+─┴─►│  C  ├──┴───+─►│  P  ├─┐
    │    │     │      │  │     │ │
    │    └─────┘      │  └─────┘ │
    │                 w2         │
    └────────────────────────────┘
```

The returned system has the transfer-function matrix
```math
\\begin{bmatrix}
I \\\\ C
\\end{bmatrix} (I + PC)^{-1} \\begin{bmatrix}
I & P
\\end{bmatrix}
```

The gang of four can be plotted like so
```julia
Gcl = extended_gangoffour(G, C) # Form closed-loop system
bodeplot(Gcl, lab=["S" "CS" "PS" "T"], plotphase=false) |> display # Plot gang of four
```
Note, the last output of Gcl is the negative of the `CS` and `PS` transfer functions from `gangoffour2`.
See [`glover_mcfarlane`](@ref) for an extended example. See also [`ncfmargin`](@ref).
"""
function extended_gangoffour(P, C)
    ny,nu = size(P)
    S = feedback(ss(I(ny+nu), P.timeevol), [0*I(ny) P; -C 0*I(nu)], pos_feedback=true)
    Gcl = S + cat(0*I(ny), -I(nu), dims=(1,2))
    Gcl
end

"""
    m, ω = ncfmargin(P, K)

Normalized coprime factor margin, defined has the inverse H∞ norm of
```math
\\begin{bmatrix}
I \\\\ K
\\end{bmatrix} (I + PK)^{-1} \\begin{bmatrix}
I & P
\\end{bmatrix}
```
A margin ≥ 0.25-0.3 is a reasonable for robustness.

If controller `K` stabilizes `P` with margin `m`, then `K` will also stabilize `P̃` if `nugap(P, P̃) < m`.

See also [`extended_gangoffour`](@ref), [`diskmargin`](@ref).
"""
function ncfmargin(P, K)
    Gcl = extended_gangoffour(P, K)
    im, w = hinfnorm2(Gcl)
    inv(im), w
end


text:
using Test
using LinearAlgebra
using SparseArrays
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function allequal(u,v)
    return all(u .== v)
end

function allapprox(u,v)
    return all(u .≈ v)
end

function allapprox(v1, v2, tol)
    flag = length(v1) == length(v2)
    return flag && all([isapprox(v1[i],v2[i],atol=tol) for i = 1:length(v1)])
end

function boundary_classifier(x)
    if x[2] ≈ 0.0 || x[1] ≈ 0.0
        return :mixed
    else
        return :traction
    end
end

function mixed_bc_components(x)
    if x[2] ≈ 0.0
        return [0.0,1.0],[1.,0.]
    elseif x[1] ≈ 0.0
        return [1.,0.],[0.,1.]
    end
end

function mixed_bc_values(x)
    return 0.0,0.0
end

function traction_bc_data(x,xR)
    if x[1] ≈ xR
        return [1.0,0.0]
    else
        return [0.0,0.0]
    end
end

function distance_function(coords,xc)
    return coords[1,:] .- xc
end

polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
mesh = UniformMesh([0.,0.],[1.,1.],[1,1])
interface_location = 0.5
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->distance_function(x,interface_location),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
lambda1,mu1 = 1.,2.
lambda2,mu2 = 1.5,2.5
stabilization = 0.1
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda1,mu1))
D2 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda2,mu2))

cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D2,stabilization)
hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
facetohelid = hybrid_element_numbering.facetohelid
interfacehelid = hybrid_element_numbering.interfacehelid
number_of_hybrid_elements = hybrid_element_numbering.number_of_hybrid_elements
dofsperelement = ufs.dofsperelement

bc_data = Dict(:mixed_components => mixed_bc_components,
               :mixed_values => mixed_bc_values,
               :traction => x->traction_bc_data(x,1.0))

system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

# HDGElasticity.assemble_uniform_boundary_faces!(system_matrix,dgmesh,ufs,
#     cellsolvers,facetohelid,boundary_classifier,bc_data)
HDGElasticity.assemble_cut_boundary_faces!(system_matrix,dgmesh,ufs,
    cellsolvers,facetohelid,interfacehelid,boundary_classifier,bc_data)
HDGElasticity.assemble_boundary_faces!(system_rhs,dgmesh,ufs,
    facetohelid,boundary_classifier,bc_data)
HDGElasticity.assemble_incoherent_interface!(system_matrix,dgmesh,ufs,
    cellsolvers,1,4,1:4,8,5:8,dofsperelement)

ndofs = dofsperelement*number_of_hybrid_elements

K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)
sol = K\R

H1 = sol[1:16]
L1 = cellsolvers[1,1].iLLxLH*H1
S1 = -D1*reshape(L1[1:12],3,:)
U1 = reshape(L1[13:20],2,:)

H2 = sol[17:end]
L2 = cellsolvers[2,1].iLLxLH*H2
S2 = -D2*reshape(L2[1:12],3,:)
U2 = reshape(L2[13:20],2,:)

strain_xx(t,l,m) = t/((l+2m) - l^2/(l+2m))
strain_yy(exx,l,m) = -l/(l+2m)*exx

ex1 = strain_xx(1.0,lambda1,mu1)
ey1 = strain_yy(ex1,lambda1,mu1)
ex2 = strain_xx(1.0,lambda2,mu2)
ey2 = strain_yy(ex2,lambda2,mu2)

ux2 = ex2*0.5
uy2 = ey2*1.0

ux1 = ux2+ex1*0.5
uy1 = ey1*1.0

@test allapprox(0.5(U1[:,1]+U1[:,3]),[ux2,0.0],1e-10)
@test allapprox(0.5(U1[:,2]+U1[:,4]),[ux2,uy1],1e-10)
@test allapprox(0.5(U2[:,1]+U2[:,3]),[ux2,0.0],1e-10)
@test allapprox(0.5(U2[:,2]+U2[:,4]),[ux2,uy2],1e-10)

testS = zeros(3,4)
testS[1,:] .= 1.0
@test allapprox(testS,S1,1e-10)
@test allapprox(testS,S2,1e-10)


text:
using Gherkin
using Base.Test

function _test_tokens(f)
    result = tokenize(f)
    expected = split(open(readall, f * ".tokens"), '\n')[1:end-1]
    @test length(result) == length(expected)
    for line in 1:length(expected)
        if expected[line] !== "(41:1)Other:/\"\"\"/"
            # TODO escaped "
            continue
        elseif !ismatch(r"\(\d+:1\)Empty://", result[line])
            @test result[line] == expected[line]
        else
            @test ismatch(r"\(\d+:1\)(Other|Empty)://", expected[line])
        end
    end
end

const GOOD = joinpath(dirname(@__FILE__), "good")
for f in filter(x -> endswith(x, ".feature"), readdir(GOOD))
    _test_tokens(joinpath(GOOD, f))
end


text:
module TestGetIndex

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running getindex tests:")

## test for changes in DataFrame behaviour
df = convert(DataFrame, rand(8, 4))
df2 = DataFrame(a = rand(8), b = rand(8))
@test isa(df[1:2], DataFrame)
@test isa(df[1:2, 3:4], DataFrame)
@test isa(df[1, 3:4], DataFrame)
@test isa(df[1], Array{Float64,1})
@test isa(df[2:4, 1], Array{Float64,1})
@test isa(df2[1], DataArray{Float64,1})
@test isa(df2[2:4, 1], DataArray{Float64,1})
@test isa(df[2, 2], Float64)

## init test values
vals = rand(30, 4)
dats = Date[Date(2013, 7, ii) for ii=1:30]
nams = [:A, :B, :C, :D]

for t = (:Timedata, :Timenum, :Timematr)
    eval(quote

        println("   Running getindex for ", $t)

        tn = $(t)(vals, nams, dats)

        ## multiple columns
        tmp = tn[2:4]
        @test isa(tmp, $(t))

        ## multiple rows, multiple columns
        tmp = tn[3:5, 1:2]
        @test isa(tmp, $(t))

        ## single row, multiple columns
        tmp = tn[5, :]
        @test isa(tmp, $(t))

        ## single column
        tmp = tn[2]
        @test isa(tmp, $(t))

        ## multiple rows, single column
        tmp = tn[5:8, 2]
        @test isa(tmp, $(t))

        ## single row, single column
        tmp = tn[5, 3]
        @test isa(tmp, $(t))

        ## multiple rows, single column
        tmp = tn[4:10, :A]
        @test isa(tmp, $(t))

        ## logical indexing
        logicCol = [true, false, true, false]
        tmp = tn[logicCol]
        @test isa(tmp, $(t))

        ## logical indexing
        logicRow = repmat([true, false, true], 10, 1)[:]
        tmp = tn[logicRow, logicCol]
        @test isa(tmp, $(t))

        ## logical indexing rows
        tmp = tn[logicRow, :]
        @test isa(tmp, $(t))

        ## logical indexing from expression
        ## ex = :(A .> 0.5)
        ## tmp = tn[ex, :]
        ## @test isa(tmp, $(t))

        ## indexing single row by Date
        tmp = tn[Date(2013, 07, 04)]
        @test isa(tmp, $(t))

        ## indexing rows by idx
        idxToFind = Date[Date(2013, 07, ii) for ii=12:18]
        tmp = tn[idxToFind]
        @test isa(tmp, $(t))

        ## indexing rows by idx, columns
        @test isa(tn[idxToFind, 2:3], $(t))
        @test isa(tn[idxToFind, :A], $(t))
        @test isa(tn[Date(2013,01,03):Date(2013,07,12)], $(t))
        @test isa(tn[Date(2013,01,03):Date(2013,07,12), :D], $(t))
        @test isa(tn[Date(2013,01,03):Date(2013,07,12), 2:3], $(t))
        @test isa(tn[Date(2013,01,03):Date(2013,07,12),
                     [true, false, false, true]], $(t))

        ## test empty instance
        @test isa(tn[Date(2013,01,03):Date(2013,07,12),
                     [false, false, false, false]], $(t))
    end)
end

end


text:
export ScatteringTransform1d

abstract type ScatteringTransform end

"""
    ScatteringTransform1d

A structure for the wavelet scattering transform for 1d signals.

## Fields: ScatteringTransform1d
 | **Field** | **Description** |
 |:----------|:----------------|
 | :D        | number of layers |
 | :Q        | Number of wavelets per octave |
 | :J        | log2(max filter scale) |
"""
struct ScatteringTransform1d <: ScatteringTransform
    D::Int # depth/number of layers
    Q::Vector{Int} # number of filters per octave for each layer
    J::Vector{Int} # log2 maximum scale of filters for each layer
end
depth(S::ScatteringTransform1d) = S.D


text:
mutable struct Dirac1D <:  TSSM.TimeSplittingSpectralMethodComplex1D{Float64}
    m::Schroedinger1D
    epsilon::Float64
    delta::Float64
    nu::Float64
    d::Array{Float64, 1}
    Qd::Array{Float64, 1}
    Qo::Array{Float64, 1}
    function Dirac1D(nx::Integer, xmin::Real, xmax::Real; epsilon::Real=1, delta::Real=1, nu::Real=1,
                     potential::Function=TSSM.none_1D) #, potential_t::Function=TSSM.none_2D)
        m = Schroedinger1D(nx, xmin, xmax, potential=potential)
        #CHECK: maybe the order of elements in mu is different!!!
        mu = [2*pi*l/(xmax-xmin) for l =vcat(0:(div(nx,2)-1), (-div(nx,2):-1))] # nx should be even
        #mu = [2*pi*l/(xmax-xmin) for l =-div(nx,2):div(nx,2)-1] # nx should be even
        eta = [sqrt(nu^2+delta^2*epsilon^2*x^2) for x=mu]
        d = [x/(delta*epsilon^2) for x=eta]
        f = [1/sqrt(2*x*(x+nu)) for x=eta]
        Qd = f.*[x+nu for x=eta]
        Qo = f.*[delta*epsilon*x for x=mu]
        new(m, epsilon, delta, nu, d, Qd, Qo)
    end
end

mutable struct WfDirac1D <: TSSM.WaveFunctionComplex1D{Float64}
    m::Dirac1D
    psi1::WfSchroedinger1D
    psi2::WfSchroedinger1D
    function WfDirac1D(m::Dirac1D)
        psi1 = WfSchroedinger1D(m.m)
        psi2 = WfSchroedinger1D(m.m)
        new(m, psi1, psi2)
    end
end

function TSSM.wave_function(m::Dirac1D )
    WfDirac1D(m)
end

set_potential!(m::Dirac1D, V::Function) = TSSM.set_potential!(m.m, V)
set_potential_t!(m::Dirac1D, V::Function) = TSSM.set_potential_t!(m.m, V)

function TSSM.to_real_space!(psi::WfDirac1D)
    to_real_space!(psi.psi1)
    to_real_space!(psi.psi2)
end

function TSSM.to_frequency_space!(psi::WfDirac1D)
    to_frequency_space!(psi.psi1)
    to_frequency_space!(psi.psi2)
end

function set!(psi::WfDirac1D, f1::Function, f2::Function)
    TSSM.set!(psi.psi1, f1)
    TSSM.set!(psi.psi2, f2)
end

function distance(psi1::WfDirac1D, psi2::WfDirac1D)
    m = psi1.m
    if m ≠ psi2.m
        error("psi1 and psi2 must belong to the same method")
    end
    sqrt(TSSM.distance(psi1.psi1, psi2.psi1)^2+TSSM.distance(psi1.psi1, psi2.psi1)^2)
end

function TSSM.set_time!(psi::WfDirac1D, t::Number)
    set_time!(psi.psi1, t)
    set_time!(psi.psi2, t)
end

TSSM.get_time(psi::WfDirac1D) = get_time(psi.psi1)
TSSM.set_propagate_time_together_with_A!(m::Dirac1D, flag::Bool) = set_propagate_time_together_with_A!(m.m, flag)
TSSM.get_propagate_time_together_with_A(m::Dirac1D) = get_propagate_time_together_with_A(m.m)

function TSSM.copy!(psi1::WfDirac1D, psi2::WfDirac1D)
    TSSM.copy!(psi1.psi1, psi2.psi1)
    TSSM.copy!(psi1.psi2, psi2.psi2)
end

function TSSM.propagate_B!(psi::WfDirac1D, dt::Number)
    propagate_B!(psi.psi1, dt)
    propagate_B!(psi.psi2, dt)
end

function TSSM.propagate_A!(psi::WfDirac1D, dt::Number)
    to_frequency_space!(psi)
    u1 = get_data(psi.psi1, true)
    u2 = get_data(psi.psi2, true)
    m = psi.m
    for l=1:length(u1)
        ee = exp(-1im*dt*m.d[l])
        h1 = ee*(m.Qd[l]*u1[l]+m.Qo[l]*u2[l])
        h2 =   (-m.Qo[l]*u1[l]+m.Qd[l]*u2[l])/ee
        u1[l] = (m.Qd[l]*h1-m.Qo[l]*h2)
        u2[l] = (m.Qo[l]*h1+m.Qd[l]*h2)
    end
end


text:
"""
    Bounds.Ellipsoid([T=Float64], N)
    Bounds.Ellipsoid(center::AbstractVector, A::AbstractMatrix)

An `N`-dimensional ellipsoid defined by

```math
(x - center)^T A (x - center) = 1
```

where `size(center) == (N,)` and `size(A) == (N,N)`.

This implementation follows the algorithm presented in Mukherjee et al. (2006).[^2]

[^2]: Pia Mukherjee, et al., 2006, ApJ 638 L51 ["A Nested Sampling Algorithm for Cosmological Model Selection"](https://iopscience.iop.org/article/10.1086/501068)
"""
mutable struct Ellipsoid{T} <: AbstractBoundingSpace{T}
    center::Vector{T}
    A::Matrix{T}
    axes::Matrix{T}
    axlens::Vector{T}
    volume::T
end

function Ellipsoid(center::AbstractVector, A::AbstractMatrix)
    axes, axlens = decompose(A)
    Ellipsoid(center, A, axes, axlens, _volume(A))
end
Ellipsoid(ndim::Integer) = Ellipsoid(Float64, ndim)
Ellipsoid(T::Type, ndim::Integer) = Ellipsoid(zeros(T, ndim), diagm(0 => ones(T, ndim)))
Ellipsoid{T}(center::AbstractVector, A::AbstractMatrix) where {T} = Ellipsoid(T.(center), T.(A))

Base.broadcastable(e::Ellipsoid) = (e,)

Base.ndims(ell::Ellipsoid) = length(ell.center)

# Returns the volume of an ellipsoid given its axes matrix
_volume(A::AbstractMatrix{T}) where {T} = T(volume_prefactor(size(A, 1))) / sqrt(det(A))
volume(ell::Ellipsoid) = ell.volume

# Returns the principal axes
axes(ell::Ellipsoid) = ell.axes

decompose(A::AbstractMatrix) = decompose(Symmetric(A))  # ensure that eigen() always returns real values

function decompose(A::Symmetric)
    E = eigen(A)
    axlens = @. 1 / sqrt(E.values)
    axes = E.vectors * Diagonal(axlens)
    return axes, axlens
end

# axes and axlens
decompose(ell::Ellipsoid) = ell.axes, ell.axlens

# Scale to new volume
function scale!(ell::Ellipsoid, factor)
    # linear factor
    f = factor^(1 / ndims(ell))
    ell.A ./= f^2
    ell.axes .*= f
    ell.axlens .*= f
    ell.volume *= factor
    return ell
end

function endpoints(ell::Ellipsoid)
    axes, axlens = decompose(ell)
        # find major axis
    major_axis = axes[:, argmax(axlens)]
    return ell.center .- major_axis, ell.center .+ major_axis
end

function Base.in(x::AbstractVector, ell::Ellipsoid)
    d = x .- ell.center
    return dot(d, ell.A * d) ≤ 1.0
end

randoffset(rng::AbstractRNG, ell::Ellipsoid{T}) where {T} = axes(ell) * randball(rng, T, ndims(ell))
Base.rand(rng::AbstractRNG, ell::Ellipsoid) = ell.center .+ randoffset(rng, ell)

fit(E::Type{<:Ellipsoid}, x::AbstractMatrix{S}; pointvol = 0) where {S} = fit(E{float(S)}, x; pointvol = pointvol)

function fit(E::Type{<:Ellipsoid{R}}, x::AbstractMatrix{S}; pointvol = 0) where {R,S}
    T = float(promote_type(R, S))
    x = T.(x)
    ndim, npoints = size(x)

    # single element is an n-sphere with pointvol volume
    if npoints == 1
        pointvol > 0 || error("Cannot compute bounding ellipsoid with one point without a valid pointvol (got $pointvol)")
        d = log(pointvol) - log(volume_prefactor(ndim))
        r = exp(d / ndim)
        A = diagm(0 => fill(1 / r^2, ndim))
        return Ellipsoid(vec(x), A)
    end
    # get estimators
    center, cov = mean_and_cov(x, 2)
    delta = x .- center
    # Covariance is smaller than r^2 by a factor of 1/(n+2)
    cov .*= ndim + 2
    # Ensure cov is nonsingular
    targetprod = (npoints * pointvol / volume_prefactor(ndim))^2
    make_eigvals_positive!(cov, targetprod)

    # get transformation matrix. Note: use pinv to avoid error when cov is all zeros
    A = pinv(cov)

    # calculate expansion factor necessary to bound each points
    f = diag(delta' * (A * delta))
    fmax = maximum(f)

    # try to avoid round-off errors s.t. furthest point obeys
    # x^T A x < 1 - √eps
    flex = 1 - sqrt(eps(T))
    if fmax > flex
        A .*= flex / fmax
    end

    ell = E(vec(center), A)

    if pointvol > 0
        minvol = npoints * pointvol
        vol = volume(ell)
        vol < minvol && scale!(ell, minvol / vol)
    end

    return ell
end

# ---------------------------------------------
# Helper functions

"""
    volume_prefactor(::Integer)

Volume constant for an n-dimensional sphere:

for n even:      (2pi)^(n    /2) / (2 * 4 * ... * n)
for n odd :  2 * (2pi)^((n-1)/2) / (1 * 3 * ... * n)
"""
function volume_prefactor(n::Integer)
    f, range = iseven(n) ? (1.0, 2:2:n) : (2.0, 3:2:n)
    for i in range
        f *= 2π / i
    end
    return f
end

# sample N samples from unit D-dimensional ball
randball(T::Type, D::Integer, N::Integer) = randball(GLOBAL_RNG, T, D, N)
function randball(rng::AbstractRNG, T::Type, D::Integer, N::Integer)
    z = randn(rng, T, D, N)
    z .*= rand(rng, T, 1, N).^(1 ./ D) ./ sqrt.(sum(p->p^2, z, dims = 1))
    return z
end

# sample from unit D-dimensional ball
randball(T::Type, D::Integer) = randball(GLOBAL_RNG, T, D)
function randball(rng::AbstractRNG, T::Type, D::Integer)
    z = randn(rng, T, D)
    z .*= rand(rng)^(1 / D) / sqrt(sum(p->p^2, z))
    return z
end

function make_eigvals_positive!(cov::AbstractMatrix, targetprod)
    E = eigen(cov)
    mask = E.values .< 1e-10
    if any(mask)
        nzprod = prod(E.values[.!mask])
        nzeros = count(mask)
        E.values[mask] .= (targetprod / nzprod)^(1 / nzeros)
        cov .= E.vectors * Diagonal(E.values) / E.vectors
    end
    return cov
end

make_eigvals_positive(cov::AbstractMatrix, targetprod) = make_eigvals_positive!(copy(cov), targetprod)

max_stars_repo_name: mileslucas/NestedSamplers

max_issues_repo_name: mileslucas/NestedSamplers

max_forks_repo_name: mileslucas/NestedSamplers


text:
edgeTable = [
0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0]

triTable = [
[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
[3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
[3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
[3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
[9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
[9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
[2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
[8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
[9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
[4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
[3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
[1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
[4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
[4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
[9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
[5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1],
[2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1],
[9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
[0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
[2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1],
[10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1],
[4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1],
[5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1],
[5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1],
[9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1],
[0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1],
[1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1],
[10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1],
[8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1],
[2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1],
[7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1],
[9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1],
[2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1],
[11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1],
[9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1],
[5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1],
[11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1],
[11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
[1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1],
[9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1],
[5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1],
[2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
[0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
[5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1],
[6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1],
[3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1],
[6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1],
[5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1],
[1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
[10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1],
[6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1],
[8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1],
[7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1],
[3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
[5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1],
[0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1],
[9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1],
[8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1],
[5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1],
[0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1],
[6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1],
[10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1],
[10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1],
[8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1],
[1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1],
[3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1],
[0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1],
[10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1],
[3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1],
[6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1],
[9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1],
[8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1],
[3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1],
[6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1],
[0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1],
[10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1],
[10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1],
[2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1],
[7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1],
[7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1],
[2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1],
[1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1],
[11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1],
[8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1],
[0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1],
[7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
[10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
[2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
[6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1],
[7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1],
[2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1],
[1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1],
[10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1],
[10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1],
[0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1],
[7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1],
[6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1],
[8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1],
[9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1],
[6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1],
[4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1],
[10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1],
[8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1],
[0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1],
[1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1],
[8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1],
[10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1],
[4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1],
[10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
[5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
[11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1],
[9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
[6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1],
[7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1],
[3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1],
[7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1],
[9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1],
[3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1],
[6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1],
[9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1],
[1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1],
[4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1],
[7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1],
[6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1],
[3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1],
[0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1],
[6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1],
[0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1],
[11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1],
[6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1],
[5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1],
[9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1],
[1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1],
[1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1],
[10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1],
[0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1],
[5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1],
[10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1],
[11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1],
[9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1],
[7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1],
[2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1],
[8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1],
[9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1],
[9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1],
[1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1],
[9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1],
[9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1],
[5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1],
[0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1],
[10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1],
[2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1],
[0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1],
[0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1],
[9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1],
[5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1],
[3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1],
[5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1],
[8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1],
[0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1],
[9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1],
[0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1],
[1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1],
[3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1],
[4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1],
[9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1],
[11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1],
[11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1],
[2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1],
[9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1],
[3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1],
[1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1],
[4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1],
[4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1],
[0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1],
[3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1],
[3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1],
[0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1],
[9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1],
[1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]
]


text:
global μ₀=4pi*1e-7

mutable struct UsovParticle
    n::Int
    νs::Array{Float64,2}
    rs2ds::Array{Float64,2}
    ns::Array{Float64,3}
    ψs::Array{Float64,1}
    ϕs::Array{Float64,1}
    λ::Float64
    hprev::Float64
    hsat::Float64

    function UsovParticle(r::Float64,R::Float64,n::Int,η::Float64,K_M2::Float64)
        ps = random_sphere_pack(r, R, n)
        νs = get_points_cartesian_random(n)
        ϕs = acos.(νs[:,3])
        ψs = atan2.(νs[:,2], νs[:,1])
        rs2ds = zeros(n, n)
        ns = zeros(3, n, n)
        for i = 1:n
            for j = 1:n
                ds = ps[j,:] -  ps[i,:]
                rs2ds[i,j] = (r/norm(ds))^3
                ns[:,i,j] = ds/norm(ds)
            end
        end
        λ = μ₀/(6*K_M2)
        return new(n, νs, rs2ds, ns, ψs, ϕs, λ, 0.0, 4.0)
    end
end

function total_energy(n, νs, ns, λ, ϕs, ψs, rs2ds, h)
    es = [ sin.(ϕs).*cos.(ψs) sin.(ϕs).*sin.(ψs)  cos.(ϕs)]

    #sum(es.*nus,2) - scallar multiplication of all particle vectors es and nus
    energy_anis_and_zeeman = sum( 0.5(1-sum(es.*νs, 2).^2) - h*es[:,3]) #external sum - by particles
    energy_dip_dip = 0
    for i = 1:n
        # external sum by j-th neighbouring particles
        #v1 = λ*rs2ds[i+1:n,i]
        #v21 = (es[i,:]'*ns[:,i+1:n,i])
        #v22 = sum( es[i+1:n,:].*ns[:,i+1:n,i]',2 )
        #v2 = v21'.*v22
        #v3 = es[i+1:n,:]*es[i,:]
        # energy_dip_dip -= sum(v1.*(3v2.-v3))
        energy_dip_dip -= sum(λ*rs2ds[i+1:n,i].* (3*(es[i,:]'*ns[:,i+1:n,i])'.*sum( es[i+1:n,:].*ns[:,i+1:n,i]',2 ) - es[i+1:n,:]*es[i,:]) )
    end

    return energy_anis_and_zeeman + energy_dip_dip
end

function total_energy_for_loop(n, νs, ns, λ, ϕs, ψs, rs2ds, h)
    es = @. [ sin(ϕs)*cos(ψs) sin(ϕs)*sin(ψs)  cos(ϕs)]

    anisotropy = 0.0
    zeeman = 0.0
    dip_dip = 0.0
    for i = 1:n
        anisotropy += 0.5(1-sum(es[i,:].*νs[i,:]).^2)
        zeeman += h*es[i,3]
        for j = i+1:n
            dip_dip += λ*rs2ds[j,i]* (3*sum(es[i,:].*ns[:,j,i])*sum(es[j,:].*ns[:,j,i]) - sum(es[j,:].*es[i,:]))
        end
    end
    return anisotropy - zeeman - dip_dip
end

function apply_field_by_step!(p::UsovParticle, h::Float64)
    energy(x) = total_energy(p.n, p.νs, p.ns, p.λ, x[1:p.n], x[p.n+1:2p.n], p.rs2ds, h);
    optimization_result = optimize(energy, [p.ϕs; p.ψs], BFGS());
    angles = Optim.minimizer(optimization_result);
    p.ϕs[1:end] = angles[1:p.n]
    p.ψs[1:end] = angles[p.n+1:2p.n]
    p.hprev = h
end

function apply_field!(p::UsovParticle, applied_field::Float64)
    hstep = 0.01
    if applied_field>p.hprev
        field_values = p.hprev:hstep:applied_field
    else
        field_values = p.hprev:-hstep:applied_field
    end

    if field_values[end] != applied_field
        push!(field_values,applied_field)
    end

    for h in field_values
        apply_field_by_step!(p, h)
    end

end

function create_for_constant_concentration(r::Float64,n::Int,η::Float64,K::Float64,M::Float64)
    R = r*(n/η)^(1/3)
    K_M2 = K/(M*M)
    p = UsovParticle(r, R, n, η, K_M2)
end

function draw_representation(p::UsovParticle)
    hstep = 0.01
    hmax = p.hsat
    h = [collect(hmax:-hstep:-hmax); collect(-hmax:hstep:hmax)]
    len = length(h)
    println(len)
    m = zeros(len)

    for i=1:len
        apply_field!(p, h[i])
        m[i] = sum(cos.(p.ϕs))/p.n
    end
    plot(h, m, xlim = (-p.hsat, p.hsat), ylim = (-2,2));
end

function get_magnetization(p::UsovParticle)
    return sum(cos.(p.ϕs))/p.n
end

function magnetize_particle(p::UsovParticle, h::Array{Float64, 1})
    len = length(h)

    m = zeros(len)
    for i = 1:len
        apply_field!(p, h[i])
        m[i] = get_magnetization(p)
    end

    return m
end


text:
function sumnondiagonal(x::Matrix{T}) where {T}
  @assert size(x,1) == size(x,2)
  s = zero(T)
  @inbounds for j in 1:size(x,2)
    for i in 1:size(x,1)
      s += (i == j) ? zero(T) : x[i,j]
    end
  end
  s
end

function backsumnondiagonal(x::AbstractMatrix,Δ::T) where {T}
  @assert size(x,1) == size(x,2)
  s = zero(x)
  @inbounds for j in 1:size(x,2)
    for i in 1:size(x,1)
      s[i,j] = (i == j) ? zero(T) : Δ
    end
  end
  s
end

sumnondiagonal(x::Flux.Tracker.TrackedMatrix) = Flux.Tracker.track(sumnondiagonal,x)
Flux.Tracker.@grad function sumnondiagonal(x)
  return(sumnondiagonal(Flux.data(x)), Δ -> (backsumnondiagonal(Flux.data(x) , Δ),))
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
module AbstractType

"""
    AbstractAtoms

All types defined for different representations of atomic systems
(different [Atoms](@ref atoms) types) must be a subtype of AbstractAtoms.
"""
abstract type AbstractAtoms end

"""
    AbstractFeatureDescriptor

All [FeatureDescriptors](@ref fd) defined for different types of features must be a
subtype of AbstractFeatureDescriptor.
"""
abstract type AbstractFeatureDescriptor end

"""
    AbstractFeaturization

All types defined for different [Featurizations](@ref fzn) must be a subtype
of AbstractFeaturization.
"""
abstract type AbstractFeaturization end

"""
    AbstractCodec

All [Codecs](@ref codecs) defined for different encoding-decoding schemes
must be a subtype of AbstractCodec.
"""
abstract type AbstractCodec end

end


text:
"""
    ∇²(i, j, k, grid, c)

Calculates the Laplacian of c via

    1/V * [δxᶜᵃᵃ(Ax * δxᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * δyᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * δzᵃᵃᶠ(c))]

which will end up at the location `ccc`.
"""
@inline function ∇²(i, j, k, grid, c)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, δᴶxᶠᵃᵃ, c) +
                                    δyᵃᶜᵃ(i, j, k, grid, δᴶyᵃᶠᵃ, c) +
                                    δzᵃᵃᶜ(i, j, k, grid, δᴶzᵃᵃᶠ, c))
end


text:
#Contains all the evaluation stuff
using JLD2
using PyPlot
using FileIO
using Statistics

a= load("Results/10_26_evaluate_pend/1_all.jld2")["RESULT"]
b= load("Results/10_26_evaluate_pend/2_all.jld2")["RESULT"]
c= load("Results/10_26_evaluate_pend/3-4_all.jld2")["RESULT"]

a_all = 4
iters = 500

RE  = zeros(a_all,iters)
A2C = zeros(a_all,iters)
TNPG= zeros(a_all,iters)
TRPO= zeros(a_all,iters)
PPO = zeros(a_all,iters)
TD3 = zeros(a_all,iters)

RE[1,:]   = a[1:6:end]
A2C[1,:]  = a[2:6:end]
TNPG[1,:] = a[3:6:end]
TRPO[1,:] = a[4:6:end]
PPO[1,:]  = a[5:6:end]
TD3[1,:]  = a[6:6:end]

RE[2,:]   = b[1:6:end]
A2C[2,:]  = b[2:6:end]
TNPG[2,:] = b[3:6:end]
TRPO[2,:] = b[4:6:end]
PPO[2,:]  = b[5:6:end]
TD3[2,:]  = b[6:6:end]

RE[3,:]   = c[1:6:3000]
A2C[3,:]  = c[2:6:3000]
TNPG[3,:] = c[3:6:3000]
TRPO[3,:] = c[4:6:3000]
PPO[3,:]  = c[5:6:3000]
TD3[3,:]  = c[6:6:3000]

RE[4,:]   = c[3001:6:end]
A2C[4,:]  = c[3002:6:end]
TNPG[4,:] = c[3003:6:end]
TRPO[4,:] = c[3004:6:end]
PPO[4,:]  = c[3005:6:end]
TD3[4,:]  = c[3006:6:end]

scale = 1
vis   = 0.5
offset = 0

function f_filt(x)
    l = length(x)
    tau = 0.1
    y = zeros(l)
    y[1] = x[1]

    for i = 2:l
        y[i] = (1-tau)*y[i-1] + tau *x[i]
    end

    return y
end

re_mean = f_filt(vec(mean(RE,dims = 1)) .+ offset)
re_std  = f_filt(scale*vec(std(RE,dims=1)))

a2c_mean = f_filt(vec(mean(A2C,dims = 1)) .+ offset)
a2c_std  = f_filt(scale*vec(std(A2C,dims=1)))

tnpg_mean = f_filt(vec(mean(TNPG,dims = 1)) .+ offset)
tnpg_std  = f_filt(scale*vec(std(TNPG,dims=1)))

trpo_mean = f_filt(vec(mean(TRPO,dims = 1)) .+ offset)
trpo_std  = f_filt(scale*vec(std(TRPO,dims=1)))

ppo_mean = f_filt(vec(mean(PPO,dims = 1)) .+ offset)
ppo_std  = f_filt(scale*vec(std(PPO,dims=1)))

td3_mean = f_filt(vec(mean(TD3,dims = 1)) .+ offset)
td3_std  = f_filt(scale*vec(std(TD3,dims=1)))

figure("pend_v0")
clf()
grid(alpha=0.25)
fill_between(1:iters,re_mean .+ re_std, re_mean .- re_std, alpha=vis)
plot(re_mean, label="REINFORCE")
fill_between(1:iters,a2c_mean .+ a2c_std, a2c_mean .- a2c_std, alpha=vis)
plot(a2c_mean, label="A2C")
fill_between(1:iters,tnpg_mean .+ tnpg_std, tnpg_mean .- tnpg_std, alpha=vis)
plot(tnpg_mean, label="TNPG")
fill_between(1:iters,trpo_mean .+ trpo_std, trpo_mean .- trpo_std, alpha=vis)
plot(trpo_mean, label="TRPO")
fill_between(1:iters,ppo_mean .+ ppo_std, ppo_mean .- ppo_std, alpha=vis)
plot(ppo_mean, label="PPO")
fill_between(1:iters,td3_mean .+ td3_std, td3_mean .- td3_std, alpha=vis)
plot(td3_mean, label="TD3")
xlabel("Iterations", fontsize = 16)
ylabel("Cost", fontsize = 16)
#legend()

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
"""
    generate_conditioned_dist(covar::CovarianceAtDate, conditioning_draws::Array{Symbol,T}) where T<:Real
Given some subset of known stochastic integral values this generates a conditional multivariate normal distribution.
"""
function generate_conditioned_distribution(covar::CovarianceAtDate, conditioning_draws::Dict{Symbol,T}) where T<:Real
    # https://stats.stackexchange.com/questions/30588/deriving-the-conditional-distributions-of-a-multivariate-normal-distribution
    what_conditioned_on = collect(keys(conditioning_draws))
    len = length(what_conditioned_on)
    if len < 1 error("Nothing that was input was in the covariance matrix. So nothing to condition on.") end
    conditioning_indices = Array{typeof(Integer(1)),1}()
    for i in 1:len
        append!(conditioning_indices, findall(what_conditioned_on[i] .==  covar.covariance_labels_))
    end
    other_indices = setdiff(1:length(covar.covariance_labels_), conditioning_indices)
    labels = covar.covariance_labels_[other_indices]
    # Segmenting the covariance matrix.
    sigma11 = covar.covariance_[other_indices,other_indices]
    sigma12 = covar.covariance_[other_indices,conditioning_indices]
    sigma21 = covar.covariance_[conditioning_indices,other_indices]
    sigma22 = covar.covariance_[conditioning_indices,conditioning_indices]
    mu1     = zeros(length(other_indices))
    mu2     = zeros(length(conditioning_indices))
    conditioned_values = map( x -> conditioning_draws[x], what_conditioned_on)
    sigma12_invsigma22 = sigma12 * inv(sigma22)
    conditional_mu = mu1 + sigma12_invsigma22 * (conditioned_values - mu2)
    conditional_sigma = sigma11 -  sigma12_invsigma22 * sigma21
    return conditional_mu, conditional_sigma, labels
end


text:
"""
    set_spins(N, index::Int, σ::SparseMatrixCSC)

Return a list of Pauli matrices.

# Examples
```julia-repl
julia> set_spins(3, 2, σᶻ) == [σ⁰, σᶻ, σ⁰]
true
```
"""
function set_spins(N, index::Int, σ::SparseMatrixCSC)
    index = pbc(index, N)
    list_mats = fill(σ⁰, N)
    list_mats[index] = σ
    return list_mats
end

"""
    set_spins(N, siteindeces::Vector{Int}, σs)

Return a list of Pauli matrices.

# Examples
```julia-repl
julia> set_spins(3, [2], [σᶻ]) == [σ⁰, σᶻ, σ⁰]
true

julia> set_spins(4, [2], [σᶻ]) == [σ⁰, σᶻ, σ⁰, σ⁰]
true

julia> set_spins(4, [2, 3], [σᶻ, σˣ]) == [σ⁰, σᶻ, σˣ, σ⁰]
true
```
"""
function set_spins(N, siteindeces::Vector{Int}, σs)
    siteindeces = pbc(siteindeces, N)
    list_mats = fill(σ⁰, N)
    for (i, σ) in zip(siteindeces, σs)
            list_mats[i] = σ
    end
    return list_mats
end

"""
    set_spins(N, siteindeces::CartesianIndex, σs)

Return a list of Pauli matrices, where siteindeces is the row and column number
of an adjacency matrix.

# Examples
```julia-repl
julia> set_spins(4, CartesianIndex(2, 3), [σᶻ, σˣ]) == [σ⁰, σᶻ, σˣ, σ⁰]
true
```
"""
function set_spins(N, siteindeces::CartesianIndex, σs)
    @assert length(siteindeces) == 2 "len(siteindeces) must be 2."
    @assert length(σs) == 2 "len(σs) must be 2."
    list_mats = fill(σ⁰, N)
    list_mats[siteindeces[1]] = σs[1]
    list_mats[siteindeces[2]] = σs[2]
    return list_mats
end


text:
# William Kahan's method, accuracy is within 2 ulp

function csd_ad_minus_bc{T<:AbstractFloat}(a::T, b::T, c::T, d::T)
    bcHi  = b*c
    bcLo  = fma(-b, c,  bcHi)  # exact
    hi    = fma( a, d, -bcHi)  # close
    hi - bcLo                  # compensated
end

csd_ab_minus_cd{T<:AbstractFloat}(a::T, b::T, c::T, d::T) = ad_minus_bc(a,c,d,b)

function csdDet2x2{T<:AbstractFloat}(m2x2::Matrix{T})
   if size(k) != (2,2)
       throw(DomainError())
   end
   csd_ad_minus_bc(reshape(m2x2,4,1)...)
end

function csdCross3D{T<:AbstractFloat}(a::Vector{T}, b::Vector{T})
    a1,a2,a3 = a[1],a[2],a[3]
    b1,b2,b3 = b[1],b[2],b[3]

    x = ad_minus_bc(a2, a3, b2, b3)
    y = ad_minus_bc(a3, a1, b3, b1)
    z = ad_minus_bc(a1, a2, b1, b2)

    [x,y,z]
end


text:
function test_subsets(members, f, itr, order=:descending)
    alreadySeen = Dict()
    lastValue = order == :descending ? Inf : -Inf
    totalCount = 0
    for s in itr
        if !isempty(s)
            if order == :descending
                @test f(s) <= lastValue
            else
                @test f(s) >= lastValue
            end

            lastValue = f(s)
        end

        alreadySeen[s] = get(alreadySeen, s, 0) + 1
        @test alreadySeen[s] == 1
        totalCount += 1
    end
    @test totalCount == 2^length(members)
end

# test orderedsubsets using a minimum function
members = [4.19897,0.0835427,0.0452934,1.77052,0.0945746,4.88804,0.135961]
test_subsets(members, minimum, ESValues.orderedsubsets(members, minimum))

# test orderedsubsets using a value that rejects subsets larger than a given value
P = 7
function subsetValue(x)

    s = length(x)
    if s == 0
        return 1e12
    elseif s > P/2
        return -1
    end

    w = (P-1)/(binomial(P,s)*s*(P-s))
    minimum(x)*w
end
test_subsets(members, subsetValue, ESValues.orderedsubsets(members, subsetValue))

# test orderedsubsets ascending using a maximum function
members = [4.19897,0.0835427,0.0452934,1.77052,0.0945746,4.88804,0.135961]
test_subsets(members, maximum, ESValues.orderedsubsets(members, maximum, :ascending), :ascending)

# test orderedsubsets ascending using a sum function
members = [4.19897,0.0835427,0.0452934,1.77052,0.0945746,4.88804,0.135961]
test_subsets(members, sum, ESValues.orderedsubsets(members, sum, :ascending), :ascending)

# test eskernelsubsets using a value that rejects subsets larger than a given value
members = collect(1:P)
variances = [4.19897,0.0835427,0.0452934,1.77052,0.0945746,4.88804,0.135961]
function subsetValue2(x)
    x = convert(Vector{Int64}, x[1])
    s = length(x)
    if s == P || s == 0
        return 1e12
    elseif s > P/2
        return -1
    elseif s == 1
        return 1e10 + minimum(variances[x])
    end

    w = (P-1)/(binomial(P,s)*s*(P-s))
    minimum(variances[x])*w
end
test_subsets(members, subsetValue2, ESValues.eskernelsubsets(members, variances))


text:
#
# finds local minimum between st1 and st3 along contour
# works by taking minimum over uniform sampling of points and then iterating once.
# returns t value of minimum along contour: tmin
# also returns integrand value at minimum: fmin
#
function locate_minimum(st1::Complex{Float64}, st3::Complex{Float64}, a::Float64, b::Float64, x::Float64, y::Float64, derivs::Bool)
    rst1,rst3 = real(st1),real(st3)
    t1,t2 = rst1,rst3
    tmin,fmin = 0.0,0.0
    # number of loops
    for k=1:n_its
        linspace!(lm_ts,t1,t2,n_t)
        # Precalculations
        for i=1:n_t
            temp1c=2.0(lm_ts[i] - rst3) + W
            temp2c=-4.0(lm_ts[i] - rst1) + V
            f = M_F_1 * atan(temp1c) - M_F_2
            g = M_G_1 * atan(temp2c) - M_G_2
            lm_gam[i] = complex(lm_ts[i],f*g)
        end
        integrand!(lm_integ,lm_integx,lm_integy,lm_gam,a,b,x,y,derivs,n_t)
        # function to be minimized: square of two norm of u, ux, uy
        for i=1:n_t
            test_min[i] = abs2(lm_integ[i])+abs2(lm_integx[i])+abs2(lm_integy[i])
        end

        fmin,idx = findmin(test_min)
        if idx == 1 || idx == n_t || k == n_its
            tmin = lm_ts[idx]
        else
            t1,t2 = lm_ts[idx-1],lm_ts[idx+1]
        end
    end
    return tmin,fmin
end


text:
import TrajectoryOptimization.Dynamics: forces, moments, inertia, inertia_inv, mass_matrix
import TrajectoryOptimization: Rotation, control_dim

export
    orientation,
    linear_velocity,
    angular_velocity

@inline Base.size(model::RigidBody{<:UnitQuaternion}) = 13, control_dim(model)
@inline Base.size(model::RigidBody) = 12, control_dim(model)

function Base.rand(model::RigidBody{D}) where {D}
    n,m = size(model)
    r = @SVector rand(3)
    q = rand(D)
    v = @SVector rand(3)
    ω = @SVector rand(3)
    x = build_state(model, r, q, v, ω)
    u = @SVector rand(m)  # NOTE: this is type unstable
    return x,u
end

function Base.zeros(model::RigidBody{D}) where D
    n,m = size(model)
    r = @SVector zeros(3)
    q = zero(D)
    v = @SVector zeros(3)
    ω = @SVector zeros(3)
    x = build_state(model, r, q, v, ω)
    u = @SVector rand(m)  # NOTE: this is type unstable
    return x,u
end

@inline rotation_type(::RigidBody{D}) where D = D
@inline rotation_type(model::InfeasibleModel) where D = rotation_type(model.model)

@inline Base.position(model::RigidBody, x) = SVector{3}(x[1],x[2],x[3])
orientation(model::RigidBody{R}, x::AbstractVector{T}, renorm=false) where {R,T} =
    R(T,x[4],x[5],x[6])
@inline linear_velocity(model::RigidBody, x) = SVector{3}(x[7],x[8],x[9])
@inline angular_velocity(model::RigidBody, x) = SVector{3}(x[10],x[11],x[12])

function orientation(model::RigidBody{UnitQuaternion{T,D}}, x::AbstractVector{T2},
        renorm=false) where {T,D,T2}
    q = UnitQuaternion{T2,D}(x[4],x[5],x[6],x[7])
    if renorm
        q = normalize(q)
    end
    return q
end
@inline linear_velocity(model::RigidBody{<:UnitQuaternion}, x) = SVector{3}(x[8],x[9],x[10])
@inline angular_velocity(model::RigidBody{<:UnitQuaternion}, x) = SVector{3}(x[11],x[12],x[13])

function flipquat(model::RigidBody{<:UnitQuaternion}, x)
    return @SVector [x[1], x[2], x[3], -x[4], -x[5], -x[6], -x[7],
        x[8], x[9], x[10], x[11], x[12], x[13]]
end

function parse_state(model::RigidBody, x, renorm=false)
    r = position(model, x)
    p = orientation(model, x, renorm)
    v = linear_velocity(model, x)
    ω = angular_velocity(model, x)
    return r, p, v, ω
end

function build_state(model::RigidBody{R}, x, q::Rotation, v, ω) where R <: Rotation
    q = SVector(R(q))
    build_state(model, x, q, v, ω)
end

function build_state(model::RigidBody{R}, x, q::SVector{4}, v, ω) where R <: Rotation
    @SVector [x[1], x[2], x[3],
              q[1], q[2], q[3], q[4],
              v[1], v[2], v[3],
              ω[1], ω[2], ω[3]]
end

function build_state(model::RigidBody{R}, x, q::SVector{3}, v, ω) where R <: Rotation
    @SVector [x[1], x[2], x[3],
              q[1], q[2], q[3],
              v[1], v[2], v[3],
              ω[1], ω[2], ω[3]]
end

function build_state(model::RigidBody{R}, x, q::Vector, v, ω) where R <: Rotation
    @SVector [x[1], x[2], x[3],
              q[1], q[2], q[3],
              v[1], v[2], v[3],
              ω[1], ω[2], ω[3]]
end

function build_state(model::RigidBody{<:UnitQuaternion}, x, q::Vector, v, ω) where R <: Rotation
    if length(q) == 3
        push!(q,q[1])
    end
    @SVector [x[1], x[2], x[3],
              q[1], q[2], q[3], q[4],
              v[1], v[2], v[3],
              ω[1], ω[2], ω[3]]
end

function fill_state(model::RigidBody{<:UnitQuaternion}, x::Real, q::Real, v::Real, ω::Real)
    @SVector [x,x,x, q,q,q,q, v,v,v, ω,ω,ω]
end

function fill_state(model::RigidBody, x::Real, q::Real, v::Real, ω::Real)
    @SVector [x,x,x, q,q,q, v,v,v, ω,ω,ω]
end

function fill_error_state(model::RigidBody, x::Real, q::Real, v::Real, ω::Real)
    @SVector [x,x,x, q,q,q, v,v,v, ω,ω,ω]
end

function fill_error_state(model::RigidBody{UnitQuaternion{T,IdentityMap}},
        x::Real, q::Real, v::Real, ω::Real) where T
    @SVector [x,x,x, q,q,q,q, v,v,v, ω,ω,ω]
end

############################################################################################
#                                DYNAMICS
############################################################################################
function dynamics(model::RigidBody{D}, x, u) where D

    r,q,v,ω = parse_state(model, x)

    F,τ = wrenches(model, x, u)
    M = mass_matrix(model, x, u)
    J = inertia(model, x, u)
    Jinv = inertia_inv(model, x, u)

    xdot = v
    qdot = kinematics(q,ω)
    vdot = M\F
    ωdot = Jinv*(τ - ω × (J*ω))

    build_state(model, xdot, qdot, vdot, ωdot)
end

function wrenches(model::RigidBody, x::SVector, u::SVector)
    F = forces(model, x, u)
    M = moments(model, x, u)
    return F,M
end

@inline mass_matrix(::RigidBody, x, u) = throw(ErrorException("Not Implemented"))
@inline forces(::RigidBody, x, u)::SVector{3} = throw(ErrorException("Not implemented"))
@inline moments(::RigidBody, x, u)::SVector{3} = throw(ErrorException("Not implemented"))
@inline inertia(::RigidBody, x, u)::SMatrix{3,3} = throw(ErrorException("Not implemented"))
@inline inertia_inv(::RigidBody, x, u)::SMatrix{3,3} = throw(ErrorException("Not implemented"))

############################################################################################
#                          STATE DIFFERENTIAL METHODS
############################################################################################
function state_diff(model::RigidBody, x::SVector{N}, x0::SVector{N}) where {N}
    r,q,v,ω = parse_state(model, x)
    r0,q0,v0,ω0 = parse_state(model, x0)
    δr = r - r0
    δq = q ⊖ q0
    δv = v - v0
    δω = ω - ω0
    build_state(model, δr, δq, δv, δω)
end

@generated function state_diff_jacobian(model::RigidBody{UnitQuaternion{T,D}},
        x0::SVector{N,T}) where {N,T,D}
    if D == IdentityMap
        :(I)
    else
        quote
            q0 = orientation(model, x0)
            G = TrajectoryOptimization.∇differential(q0)
            I1 = @SMatrix [1 0 0 0 0 0 0 0 0 0 0 0;
                           0 1 0 0 0 0 0 0 0 0 0 0;
                           0 0 1 0 0 0 0 0 0 0 0 0;
                           0 0 0 G[1] G[5] G[ 9] 0 0 0 0 0 0;
                           0 0 0 G[2] G[6] G[10] 0 0 0 0 0 0;
                           0 0 0 G[3] G[7] G[11] 0 0 0 0 0 0;
                           0 0 0 G[4] G[8] G[12] 0 0 0 0 0 0;
                           0 0 0 0 0 0 1 0 0 0 0 0;
                           0 0 0 0 0 0 0 1 0 0 0 0;
                           0 0 0 0 0 0 0 0 1 0 0 0;
                           0 0 0 0 0 0 0 0 0 1 0 0;
                           0 0 0 0 0 0 0 0 0 0 1 0;
                           0 0 0 0 0 0 0 0 0 0 0 1.]
        end
    end
end

function state_diff_jacobian(model::RigidBody{<:Rotation}, x0::SVector)
    q0 = orientation(model, x0)
    G = TrajectoryOptimization.∇differential(q0)
    return @SMatrix [
        1 0 0 0 0 0 0 0 0 0 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0 0;
        0 0 0 G[1] G[4] G[7] 0 0 0 0 0 0;
        0 0 0 G[2] G[5] G[8] 0 0 0 0 0 0;
        0 0 0 G[3] G[6] G[9] 0 0 0 0 0 0;
        0 0 0 0 0 0 1 0 0 0 0 0;
        0 0 0 0 0 0 0 1 0 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 0 0 0 0 0 0 0 0 0 1;
    ]
    # return I # I1 = Diagonal(@SVector ones(N))
end

state_diff_size(::RigidBody) = 12
state_diff_size(::RigidBody{UnitQuaternion{T,Union{IdentityMap}}}) where T = 13
state_diff_size(::RigidBody{UnitQuaternion{T,Union{ReNorm}}}) where T = 13

function TrajectoryOptimization.∇²differential(model::RigidBody,
        x::SVector, dx::AbstractVector)
      q = orientation(model, x)
      dq = SVector(orientation(model, dx, false))
      G2 = TrajectoryOptimization.∇²differential(q, dq)
      return @SMatrix [
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 G2[1] G2[4] G2[7] 0 0 0 0 0 0;
            0 0 0 G2[2] G2[5] G2[8] 0 0 0 0 0 0;
            0 0 0 G2[3] G2[6] G2[9] 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
      ]
end

function TrajectoryOptimization.inverse_map_jacobian(model::RigidBody{<:UnitQuaternion},
        x::SVector)
    q = orientation(model, x)
    G = TrajectoryOptimization.inverse_map_jacobian(q)
    return @SMatrix [
            1 0 0 0 0 0 0 0 0 0 0 0 0;
            0 1 0 0 0 0 0 0 0 0 0 0 0;
            0 0 1 0 0 0 0 0 0 0 0 0 0;
            0 0 0 G[1] G[4] G[7] G[10] 0 0 0 0 0 0;
            0 0 0 G[2] G[5] G[8] G[11] 0 0 0 0 0 0;
            0 0 0 G[3] G[6] G[9] G[12] 0 0 0 0 0 0;
            0 0 0 0 0 0 0 1 0 0 0 0 0;
            0 0 0 0 0 0 0 0 1 0 0 0 0;
            0 0 0 0 0 0 0 0 0 1 0 0 0;
            0 0 0 0 0 0 0 0 0 0 1 0 0;
            0 0 0 0 0 0 0 0 0 0 0 1 0;
            0 0 0 0 0 0 0 0 0 0 0 0 1;
    ]
end

function TrajectoryOptimization.∇²differential(model::RigidBody{UnitQuaternion{T,IdentityMap}},
        x::SVector, dx::SVector) where T
    return I*0
end

function TrajectoryOptimization.inverse_map_jacobian(model::RigidBody, x::SVector)
    return I
end

function TrajectoryOptimization.inverse_map_∇jacobian(model::RigidBody{<:UnitQuaternion},
        x::SVector, b::SVector)
    q = orientation(model, x)
    bq = @SVector [b[4], b[5], b[6]]
    ∇G = TrajectoryOptimization.inverse_map_∇jacobian(q, bq)
    return @SMatrix [
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 ∇G[1] ∇G[5] ∇G[ 9] ∇G[13] 0 0 0 0 0 0;
        0 0 0 ∇G[2] ∇G[6] ∇G[10] ∇G[14] 0 0 0 0 0 0;
        0 0 0 ∇G[3] ∇G[7] ∇G[11] ∇G[15] 0 0 0 0 0 0;
        0 0 0 ∇G[4] ∇G[8] ∇G[12] ∇G[16] 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
    ]

end

function TrajectoryOptimization.inverse_map_∇jacobian(model::RigidBody,
        x::SVector, b::SVector)
    return I*0
end

function TrajectoryOptimization.discrete_dynamics(::Type{RK3},
        model::RigidBody{UnitQuaternion{T,ReNorm}},
        z::KnotPoint) where T
    x = discrete_dynamics(RK3, model, state(z), control(z), z.t, z.dt)
    r,q,v,ω = parse_state(model, x)
    q = normalize(q)
    build_state(model, r,q,v,ω)
end


text:
fname = :TRKdense
c = Combi(hessian_dense,PDataK,solve_modelKTR,preprocessKTR,decreaseKTR,TparamsKTR())
include("Template.jl")


text:
function perimeter_2!(r)

    dir = []
    set(Nord, West, dir)

    for i in [Sud, Ost, Nord, West]
        go_forward(i)
    end

    while length(dir) > 0
        move!(r, inverse_side(pop!(dir)))
    end

end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
using Test
using CryptoGroups: Hash, PRG, RO

h = Hash("sha256")

s = hex2bytes("000102030405060708090a0b0c0d0e0f101112131415161718191a1b1c1d1e1f")
prg = PRG(h, s)

@test bytes2hex(prg[1:128]) == "70f4003d52b6eb03da852e93256b5986b5d4883098bb7973bc5318cc66637a8404a6950a06d3e3308ad7d3606ef810eb124e3943404ca746a12c51c7bf7768390f8d842ac9cb62349779a7537a78327d545aaeb33b2d42c7d1dc3680a4b23628627e9db8ad47bfe76dbe653d03d2c0a35999ed28a5023924150d72508668d244"

# As can be found on page 36 in:
# Wikstrom, “How To Implement A Stand-Alone Veriﬁer for the Veriﬁcatum Mix-Net.”
ro = RO(h, 65)
d = hex2bytes("000102030405060708090a0b0c0d0e0f101112131415161718191a1b1c1d1e1f")
@test bytes2hex(ro(d)) == "001a8d6b6f65899ba5"

ro = RO(h, 261)
d = hex2bytes("000102030405060708090a0b0c0d0e0f101112131415161718191a1b1c1d1e1f")
@test bytes2hex(ro(d)) == "1c04f57d5f5856824bca3af0ca466e283593bfc556ae2e9f4829c7ba8eb76db878"


text:
#=

This script extracts information about firm sizes in Eindhoven from the Excel files
provided by UCSC, scales them to a population of 10_000 and writes the results
back to a CSV file as a list of individual employer sizes (one employer per row).

It was tricky to find a way to do a reasonable scaling while preserving the jobs
ratio for the general population. The logic can be inferred from the code, but
here is a quick rundown (I'm happy to add more explanations on request):

- We read the number of employees by company size (but not the number of firms)
  for both public and private sector.
- We aggregate those two datasets together, using the more precise bin sizes
  from the private sector data. The public sector jobs are distributed in those
  "sub-bins" at the same pro-rata that private sector jobs are. This relies on
  the assumption that no private sector bin overlaps with two public sector
  bins, but that is the case in our data.
- For each bin, find the "optimal" employer size, which is more or less the size
  that minimizes the number of orphan jobs. See the `find_optimal_employer_size`
  for details.
- Generate individual company sizes from numbers computed at the previous step.
- Create one more company size to capture the orphan jobs.

The result is a list of individual company sizes written to a csv file that we
can directly use to initialize the model.

The script was written for Julia 0.6.4.

-- Nicolas 2018-08-31
-- edited Simone 2019-09-20

=#
using DataFrames, ExcelReaders

function make_df(xl_file, sheet, headers, jobs)
  function ranges(hs)
    parse_range(a) = @. parse(Int, a) + [0, 0]
    rs = @. parse_range(matchall(r"\d+", string(hs)))
    rs[end][2] = typemax(Int)
    [r[1]:r[2] for r ∈ rs]
  end
  DataFrame(
    ranges = vec(ranges(readxl(xl_file, "$(sheet)!$(headers)"))),
    jobs = vec(readxl(xl_file, "$(sheet)!$(jobs)"))
  )
end

function distribute_jobs(df_private, df_public)
  falls_within(inner, outer) = inner.start >= outer.start && inner.stop <= outer.stop
  vcat(map(eachrow(df_public)) do row
    df = df_private[falls_within.(df_private[:range], [row[:range]]), :]
    df[:jobs] = df[:jobs] + df[:jobs] / sum(df[:jobs]) * row[:jobs]
    df
  end...)
end

function aggregate_big_employers(df_total)
  # identify the ranges where the number of jobs is
  # not enough to meet the lower bound of the range
  df_total[:enough] = df_total[:jobs] .>= getfield.(df_total[:ranges], :start)
  (by(df_total, :enough) do df
    if df[:enough][1] df[[:ranges, :jobs]] else
      # and aggregate all of them in a single range
      # with the lowest lower bound and the max upper bound
      DataFrame(
        ranges = [df[:ranges][1].start:df[:ranges][end].stop],
        jobs = sum(df[:jobs])
      )
	end
	end)[[:ranges, :jobs]]
end

function find_optimal_employer_size(df_total)
  by(df_total, names(df_total)) do df
	if floor(Int, df[:jobs][1]) == 38
		ranges = df[:ranges][1].start:min(df[:ranges][1].stop, 801)
	else
		ranges = df[:ranges][1].start:min(df[:ranges][1].stop, floor(Int, df[:jobs][1]))
	end
	println(ranges)
	# debug
	#println([n for n ∈ ranges])
    # for each range, find the employer size that minimizes the number of "orphaned" jobs
    # and the reciprocal of the number of companies (doubly weighted). The latter is
    # to prevent the last category to contain one big employer, which would result in
    # way too many links in the complete job network
    scores = [rem(df[:jobs][1], n) + 2 / div(df[:jobs][1], n) for n ∈ ranges]
	# debug
	#println(scores)
	# debug
	#println(ranges[indmin(scores)])
    DataFrame(employer_size = ranges[indmin(scores)])
  end
end

xl_file    = openxl("/Users/digitaldust/Documents/clients/LABSS-ISTC-CNR/PROTON-OC/inputs/eindhoven/raw/data Eindhoven Proton.xlsx")
#df_private = make_df(xl_file, "Private sector", "C3:P3", "C7:P7")
#df_public  = make_df(xl_file, "Public sector", "C19:H19", "C22:H22")
df_total   = make_df(xl_file, "firms by size", "C1:M1", "C3:M3")

eindhoven_pop = readxl(xl_file, "1b!F24")
sim_pop     = 10_000

println(df_total)

df_total[:jobs] = df_total[:jobs] * (sim_pop / eindhoven_pop)

println(df_total)

df_total = aggregate_big_employers(df_total)

println(df_total)

df_total = find_optimal_employer_size(df_total)

println(df_total)

make_sizes(row) = fill(row[:employer_size], Int(row[:jobs] ÷ row[:employer_size]))
sizes = map(make_sizes, eachrow(df_total)) |> Iterators.flatten |> collect
orphans = round(Int, sum(rem.(df_total[:jobs], df_total[:employer_size])))
push!(sizes, orphans)

open(io -> writedlm(io, sizes), "/Users/digitaldust/Documents/clients/LABSS-ISTC-CNR/PROTON-OC/inputs/eindhoven/data/employer_sizes.csv", "w")

max_stars_repo_name: CherryKitty/PROTON-OC


text:
N = 50000
x = randn(N)
y = randn(N)
z = randn(N)
t1 = rand.(Gamma.(10,1*exp.(x-0.3y)))
t2 = rand(Gamma(15,1),N)

W = [(t2[i]>t1[i]) ? Event(t1[i], false) : Event(t2[i], true) for i = 1:N]

df = DataFrame(x = x, y = y, z = z, a = W)

res = aft(@formula(a ~ 1 + x +y + x*y+ z), df, PGamma(); tol = 1e-3)

true_res = [log(10), log(10), 1., -0.3, 0., 0.]
@test res.coefmat.cols[1] ≈ true_res atol = 1e-1
println("Aft works:")
println("$(res.coefmat.cols[1])~$(true_res)")
t = @time res = aft(@formula(a ~ 1 + x +y + x*y+ z), df, PGamma(); tol = 1e-3)
println(t)


text:
__precompile__()

module ImageMethodReverb

using Random, LinearAlgebra

include("rim_solver.jl")
include("utils.jl")

end

max_stars_repo_name: nantonel/JuRIM

max_issues_repo_name: nantonel/JuRIM

max_forks_repo_name: nantonel/JuRIM


text:
using JLD, DataFrames

include(joinpath(dirname(@__FILE__),"../src/dataframe.jl"))

r  = load(joinpath(dirname(@__FILE__),"resp.jld"), "resp")
df = dataframe(r)

facts("dataframe works on Request object") do

  context("there are three rows") do
      @fact  size(df,1) --> 30
      @fact  size(df,1) --> 30
  end

  context("there are six columns") do
      @fact  size(df,2) --> 13
      @fact  size(df,2) --> 13
  end

  context("correct values") do
      @fact  df[:Close][1]  --> 115.01
      @fact  df[:Close][30] --> 109.95
  end

  context("oldest date last") do
      @fact df[1,1] < df[30,1] --> true
  end
end


text:
include("kaiming_init.jl")
using Flux
using Flux:@treelike
using NNlib: @fix, σ_stable

################################### NAC ########################################

struct NAC{F,S,T}
  Ŵ::S
  M̂::F
  b::T
end

NAC(in::Integer, out::Integer;
               initŴ = kaiming_uniform, initM̂ = kaiming_uniform, initb = zeros) =
    NAC(param(initŴ(out, in; a=√5.0)), param(initM̂(out, in; a=√5.0)), param(initb(out)))

@treelike NAC

function (nac::NAC)(x)
  Ŵ, M̂, b = nac.Ŵ, nac.M̂, nac.b
  W = tanh.(Ŵ) .* σ_stable.(M̂)
  @fix W*x .+ b
end

Base.show(io::IO, l::NAC) = print(io, "NAC(", size(l.Ŵ, 2), ", ", size(l.Ŵ, 1), ")")

################################## NALU ########################################

struct NALU{N,S,T,E}
  nac::N
  G::S
  b::T
  ϵ::E
end

NALU(in::Integer, out::Integer;
               initĜ = kaiming_uniform, initb = zeros) =
  NALU(NAC(in, out), param(initĜ(out, in; a=√5.0)), param(initb(out)), 1e-10)

@treelike NALU

function (nalu::NALU)(x)
  nac, G, b, ϵ = nalu.nac, nalu.G, nalu.b, nalu.ϵ
  a = nac(x)
  log_inp = log.(abs.(x) .+ ϵ)
  m = exp.(nac(log_inp))
  g = σ_stable.(G*x .+ b)
  @fix g .* a .+ (1.0 .- g) .* m
end

Base.show(io::IO, l::NALU) = print(io, "NALU(", size(l.G, 2), ", ", size(l.G, 1), ")")


text:
LogisticClassifier = @load LogisticClassifier pkg="ScikitLearn"
RandomForestClassifier = @load RandomForestClassifier pkg="ScikitLearn"

"""
    linear_classifier(X, Y)

    Returns a pipeline for the linear classifier used in our analyses
"""
function linear_classifier(X, Y)
    pipe = @pipeline FeatureSelector FillImputer Standardizer OneHotEncoder LogisticClassifier(max_iter=1000)
    return machine(pipe, remove_lowvar(X), Y)
end

"""
    rf_classifier(X, Y)

    Returns a pipeline for the random forest classifier used in our analyses
"""
function rf_classifier(X, Y)
    pipe = @pipeline FeatureSelector FillImputer Standardizer OneHotEncoder RandomForestClassifier(n_estimators=100)
    return machine(pipe, remove_lowvar(X), Y)
end