text:
abstract type AbstractExogenous end

abstract type EmptyProcess <: AbstractExogenous end
abstract type AbstractProcess <: AbstractExogenous end
abstract type DiscreteProcess <: AbstractProcess end
abstract type ContinuousProcess <: AbstractProcess end
abstract type IIDExogenous <: ContinuousProcess end

abstract type AbstractDiscretizedProcess <: DiscreteProcess end

# this is a bit crude but not performance critical, for now
function node(::Type{Point}, dp::DiscreteProcess, i::Int)
    v = node(dp, i)
    SVector{length(v)}(v)
end
function inode(::Type{Point}, dp::DiscreteProcess, i::Int, j::Int)
    v = inode(dp, i, j)
    SVector{length(v)}(v)
end

struct ConstantProcess <: ContinuousProcess
    μ::Array{Float64,1}
end

ConstantProcess(;μ=zeros(1)) = ConstantProcess(μ)
ConstantProcess(;μ=zeros(1)) = ConstantProcess(μ)

###
### Discretized process
###

# date-t grid has a known structure
mutable struct GDP{TG} <: AbstractDiscretizedProcess
    grid::TG
    integration_nodes::Array{Matrix{Float64},1}
    integration_weights::Array{Vector{Float64},1}
end

n_nodes(dp::GDP) = n_nodes(dp.grid)
node(dp::GDP, i) = node(dp.grid, i)
n_inodes(dp::GDP, i::Int) = size(dp.integration_nodes[i], 1)
inode(dp::GDP, i::Int, j::Int) = dp.integration_nodes[i][j, :]
iweight(dp::GDP, i::Int, j::Int) = dp.integration_weights[i][j]

# function node(::Type{Point}, dp::DiscreteProcess, i::Int) where n_x
#     SVector(dp.grid.nodes[i]...)
# end
# function inode(::Type{Point}, dp::DiscreteProcess, i::Int, j::Int) where n_x
#     v = inode(dp, i, j)
#     SVector(dp.grid.integration_nodes[i][j,:]...)
# end

function node(::Type{Point{n_x}}, dp::GDP, i::Int) where n_x
    SVector{n_x}(dp.grid.nodes[i]...)
end
function inode(::Type{Point{n_x}}, dp::GDP, i::Int, j::Int) where n_x
    SVector{n_x, Float64}(dp.integration_nodes[i][j,:]...)
end

function Product(gdp1::GDP, gdp2::GDP)
  In = [ gridmake(gdp1.integration_nodes[i], gdp2.integration_nodes[j] ) for i = 1:n_nodes(gdp1)  for j = 1:n_nodes(gdp2) ]
  Iw =[kron(gdp1.integration_weights[i], gdp2.integration_weights[j] ) for i = 1:n_nodes(gdp1)  for j = 1:n_nodes(gdp2) ]

  i_grid  = ProductGrid(gdp1.grid, gdp2.grid)

  return GDP(i_grid, In, Iw)
end

# date-t grid is unstructured
mutable struct DiscreteMarkovProcess <: AbstractDiscretizedProcess
    grid::UnstructuredGrid
    transitions::Matrix{Float64}
    values::Matrix{Float64}
    i0::Int
end

ndims(dmp::DiscreteMarkovProcess) = ndims(dmp.grid)

discretize(::Type{DiscreteMarkovProcess}, mp::DiscreteMarkovProcess) = mp
discretize(mp::DiscreteMarkovProcess) = mp

DiscreteMarkovProcess(transitions::Matrix{Float64}, values::Matrix{Float64}) =
    DiscreteMarkovProcess(UnstructuredGrid{size(values,2)}(values), transitions, values)

DiscreteMarkovProcess(grid::UnstructuredGrid, transitions::Matrix{Float64}, values::Matrix{Float64}) =
    DiscreteMarkovProcess(grid, transitions, values, 1)

n_nodes(dp::DiscreteMarkovProcess) = size(dp.values, 1)
n_inodes(dp::DiscreteMarkovProcess, i::Int) = size(dp.values, 1)
iweight(dp::DiscreteMarkovProcess, i::Int, j::Int) = dp.transitions[i, j]
default_index(dp::DiscreteMarkovProcess) = dp.i0
inode(dp::DiscreteMarkovProcess, i::Int, j::Int) = dp.values[j, :]
node(dp::DiscreteMarkovProcess, i) = dp.values[i, :]
inode(::Type{Point{d}}, dp::DiscreteMarkovProcess, i::Int, j::Int) where d = SVector{d}(dp.values[j, :]...)
node(::Type{Point{d}}, dp::DiscreteMarkovProcess, i) where d = SVector{d}( dp.values[i, :] ...)

function MarkovProduct(mc1::DiscreteMarkovProcess, mc2::DiscreteMarkovProcess)
    Q = gridmake(mc1.values, mc2.values)
    P = fkron(mc1.transitions, mc2.transitions)
    return DiscreteMarkovProcess(P, Q)
end

function MarkovProduct(mcs::DiscreteMarkovProcess...)
    reduce(MarkovProduct, mcs)
end

# date-t grid is empty

mutable struct DiscretizedIIDProcess <: AbstractDiscretizedProcess
    # nodes::Matrix{Float64}
    grid::EmptyGrid
    integration_nodes::Matrix{Float64}
    integration_weights::Vector{Float64}
end

DiscretizedIIDProcess(x, w) = DiscretizedIIDProcess(EmptyGrid{size(x,2)}(), x, w)

n_nodes(dp::DiscretizedIIDProcess) = 0
n_inodes(dp::DiscretizedIIDProcess, i::Int) = size(dp.integration_nodes, 1)
inode(dp::DiscretizedIIDProcess, i::Int, j::Int) = dp.integration_nodes[j, :]
iweight(dp::DiscretizedIIDProcess, i::Int, j::Int) = dp.integration_weights[j]
node(dip::DiscretizedIIDProcess, i::Int) = fill(NaN, n_inodes(dip, 1))

node(::Type{Point{d}}, dip::DiscretizedIIDProcess, i::Int) where d = fill(0, SVector{d,Float64})
# node(::Type{Point{d}}, dip::DiscretizedIIDProcess, i::Int) where d = fill(NaN, SVector{d,Float64}) # TODO: this should be the correct version !!

inode(::Type{Point{d}}, dip::DiscretizedIIDProcess, i::Int, j::Int) where d = SVector{d}(dip.integration_nodes[j, :]...)

# Normal law

struct MvNormal <: IIDExogenous
    μ::Vector{Float64}
    Σ::Matrix{Float64}
end

MvNormal(Σ::Matrix{Float64}) = MvNormal(zeros(size(Σ, 1)), Σ)
MvNormal(σ::Float64) = MvNormal(reshape([σ^2], 1, 1))

Normal(;Σ=zeros(1,1)) = MvNormal(Σ)

UNormal(;σ=0.0) = MvNormal(reshape([σ^2], 1, 1))

function discretize(mvn::MvNormal; n=5::Union{Int, Vector{Int}})
    x, w = QE.qnwnorm(n, mvn.μ, mvn.Σ)
    DiscretizedIIDProcess(x, w)
end

function Base.rand(mvn::MvNormal, args...)
    dist = Distributions.MvNormal(mvn.μ, mvn.Σ)
    return rand(dist, args...)
end

function simulate(mvn::MvNormal, e0::Vector{Float64}; N::Integer, T::Int, stochastic::Bool=true)
    dist = QE.MVNSampler(mvn.μ, mvn.Σ)
    d = length(mvn.μ)
    out = zeros(d, N, T)
    for i in 1:N
        out[:, i, 1] = e0
    end
    if stochastic
        for t in 2:T
            out[:, :, t] = rand(dist, N)
        end
    end
    AxisArray(out, Axis{:V}(1:d), Axis{:N}(1:N), Axis{:T}(1:T))
end

function simulate(mvn::MvNormal, N::Integer, T::Int, e0::Vector{Float64}; stochastic::Bool=true)
    simulate(mvn, e0; N=N, T=T, stochastic=stochastic)
end

function simulate(mvn::MvNormal; N::Integer, T::Int, stochastic::Bool=true)
    simulate(mvn, mvn.μ; N=N, T=T, stochastic=stochastic)
end

function simulate(mvn::MvNormal, N::Integer, T::Int; stochastic::Bool=true)
    simulate(mvn, N, T, mvn.μ; stochastic=stochastic)
end

function response(mvn::MvNormal, e1::AbstractVector; T::Int=40)
    d = length(mvn.μ)
    out = zeros(d, T)
    out[:, 2] = e1
    return out
end

function simulate(process::DiscreteMarkovProcess; N::Int, T::Int, i0::Int)
    mc_qe = QE.MarkovChain(process.transitions)
    inds = Array{Int}(undef, T, N)
    QE.simulate_indices!(inds, mc_qe, init=i0)
    AxisArray(inds, Axis{:T}(1:T),  Axis{:N}(1:N))
end

function simulate(process::DiscreteMarkovProcess, N::Int, T::Int, i0::Int)
    simulate(process; N=N, T=T, i0=i0)
end

function simulate(process::DiscreteMarkovProcess, m0::AbstractVector{Float64}; N::Int, T::Int)
    # try to find index in process.values. IF we do, then simulate using
    # the corresponding row index. If we don't, then throw an error
    for i in 1:size(process.values, 1)
        if process.values[i, :] == m0
            return simulate(process; N=N, T=T, i0=i)
        end
    end
    error("Couldn't find the vector `m0` in process.values. "*
          "Try passing an integer `i0` instead")
end

function simulate(process::DiscreteMarkovProcess, N::Int, T::Int, m0::AbstractVector{Float64})
    simulate(process, m0; N=N, T=T)
end

function simulate_values(process::DiscreteMarkovProcess; N::Int, T::Int, i0::Int)
    inds = simulate(process; N=N, T=T, i0=i0)
    n_values = size(process.values, 2)
    out_values = Array{Float64}(n_values, N, T)
    for n in 1:N, t in 1:T
        out_values[:, n, t] = process.values[inds[t, n], :]
    end
    AxisArray(out_values, Axis{:n}(1:n_values), Axis{:T}(1:T), Axis{:N}(1:N))
end

function simulate_values(process::DiscreteMarkovProcess, N::Int, T::Int, i0::Int)
    simulate_values(process; N=N, T=T, i0=i0)
end

# VAR 1

mutable struct VAR1 <: ContinuousProcess
    μ::Array{Float64,1}
    R::Array{Float64,2}
    Σ::Array{Float64,2}
end

VAR1(;ρ::Float64=0.0, Σ=ones(1,1)) = VAR1(ρ, Σ)

function VAR1(R::Array{Float64,2}, Σ::Array{Float64,2})
    M = zeros(size(R, 1))
    @assert size(R)==size(Σ)
    return VAR1(M, R, Σ)
end

function VAR1(ρ::Array{Float64,1}, Σ::Array{Float64,2})
    R = diagm(ρ)
    return VAR1(R, Σ)
end

function VAR1(ρ::Float64, Σ::Array{Float64,2})
    p = size(Σ,1)
    R = Matrix(ρ*I, p, p)
    return VAR1(R, Σ)
end

function discretize(::Type{GDP}, var::VAR1; n::Union{Int, Vector{Int}}=5, n_i::Union{Int, Vector{Int}}=5,  n_std::Int=2)
    R = var.R
    M = var.μ
    Σ = var.Σ
    S = QE.solve_discrete_lyapunov(R,Σ)  # asymptotic variance
    sig = diag(S)
    min = var.μ - n_std*sqrt.(sig)
    max = var.μ + n_std*sqrt.(sig)
    if n isa Int
        N = fill(n, length(min))  ::Vector{Int}
    else
        N = n ::Vector{Int}
    end
    grid = UCGrid{length(min)}(min,max,N)
    # discretize innovations
    x,w = QE.qnwnorm(n_i, zeros(size(var.Σ,1)), var.Σ)
    integration_nodes = [ cat([(M + R*(node(grid, i)-M) + x[j,:])' for j=1:size(x,1)]...; dims=1) for i in 1:n_nodes(grid)]
    integration_weights = [w for i in 1:n_nodes(grid)]
    return GDP(grid, integration_nodes, integration_weights)
end

discretize(var::VAR1; args...) = discretize(DiscreteMarkovProcess, var; args...)

function discretize(::Type{DiscreteMarkovProcess}, var::VAR1; n::Union{Int, Vector{Int}}=3)

    # it would be good to have a special type of VAR1 process
    # which has a scalar autoregressive term
    R =  var.R
    d = size(R, 1)
    ρ = R[1, 1]
    @assert maximum(abs, R-Matrix(ρ*I,d,d))<1e-16

    σ = var.Σ

    if size(var.Σ, 1) == 1
        mc_qe = QE.rouwenhorst(n, ρ, sqrt(σ[1]))
        return DiscreteMarkovProcess(mc_qe.p, appenddim(collect(mc_qe.state_values)))
    end

    L = chol(σ) # σ = L'*L
    if n isa Int
        NN = fill(n, d) # default: same number of nodes in each dimension
    else
        NN = n
    end

    components = [QE.rouwenhorst(N, ρ, 1.0) for N in NN]
    mc_components = [DiscreteMarkovProcess(mc.p, appenddim(collect(mc.state_values))) for mc in components]
    mc_prod = MarkovProduct(mc_components...)
    mc_prod.values = mc_prod.values*L'
    return mc_prod
end

"""
```julia
simulate(var::VAR1, N::Int, T::Int, x0::Vector{Float64}
         stochastic::Bool=true, irf::Bool=false, e0::Vector{Float64}=zeros(0))
```

Construct `N` simulated paths of length `T` the vector autoregressive process
`var`, starting each simulated path from `x0`.

If `irf` is `true`, then `e0` must be a vector of the same length as `var.μ`
and `e0` will be used as the initial shocks.

If `stochastic` is `true`, construct a sequence of shocks by drawing `N` paths
of length `T` from the distribution `Normal(0, var.Σ)` (where `0` is a
vector of zeros of the appropriate length). Otherwise, if `stochastic` is
false, set all shocks to 0 (unless `irf` is true -- see above).

The output is an `AxisArray` contining variables, paths, and time on the 3 axes
"""
function simulate(var::VAR1, x0::Vector{Float64}; N::Int, T::Int,
    stochastic::Bool=true, irf::Bool=false, e0::Vector{Float64}=zeros(0))

    n = size(var.μ, 1)
    XN = zeros(n, N, T)
    E = zeros(n, N, T)

    if stochastic
        dist = QE.MVNSampler(zeros(n), var.Σ)
        for jj in 1:N
            E[:, jj, :] = rand(dist, T)
        end
    end

    if irf
        E[:, :, 1] = repeat(e0,N,1)
    end

    # Initial conditions
    for i in 1:N
        XN[:, i, 1] = x0
        for  ii in 1:T-1
            XN[:, i, ii+1] = var.μ+var.R*(XN[:, i, ii]-var.μ)+E[:, i, ii]
        end
    end
    AxisArray(XN, Axis{:V}(1:n), Axis{:N}(1:N), Axis{:T}(1:T))
end

function simulate(var::VAR1, N::Int, T::Int, x0::Vector{Float64};
    stochastic::Bool=true, irf::Bool=false, e0::Vector{Float64}=zeros(0))

    simulate(var, x0; N=N, T=T, stochastic=stochastic, irf=irf, e0=e0)
end

function simulate(var::VAR1; N::Int, T::Int, kwargs...)
    return simulate(var, var.μ; N=N, T=T, kwargs...)
end

function simulate(var::VAR1, N::Int, T::Int; kwargs...)
    return simulate(var; N=N, T=T, kwargs...)
end

function response(var::VAR1, x0::AbstractVector,
                  e1::AbstractVector; T::Int=40)
    simulate(var, x0; N=1, T=T, stochastic=false, irf=true, e0=e1)[1, :, :]
end

function response(var::VAR1, e1::AbstractVector; T::Int=40)
    simulate(var; N=1, T=T, stochastic=false, irf=true, e0=e1)[1, :, :]
end

function response(var::VAR1, x0::AbstractVector, index_s::Int; T::Int=40)
    e1 = zeros(size(var.μ, 1))
    Impulse = sqrt.(diag(var.Σ)[index_s])
    e1[index_s] = Impulse
    simulate(var, x0; N=1, T=T, stochastic=false, irf=true, e0=e1)[1, :, :]
end

function response(var::VAR1; T::Int=40)
    e1 = sqrt.(diag(var.Σ))
    simulate(var; N=1, T=T, stochastic=false, irf=true, e0=e1)[1, :, :]
end

function ErgodDist(var::VAR1, N::Int, T::Int)

    # Simulate for T period, N times
    #  #    - simulate a VAR1 process for T periods N times#
    n = size(var.μ, 1)
    # srand(123) # Setting the seed
    d = MvNormal(var.Σ)
    VAR_process = simulate(var, N, T)
    E = VAR_process[2]
    XN = VAR_process[1]
    # Computing moments
    Mean_sim = mean(squeeze(mean(VAR_process, 3), 2), 2)
    #Computing the std of simulted Processes (Covariance matrix)
    E1_d = (E[:, :, 1] - repeat(mean(E[:, :, 1], 2), inner=[1, T]))
    Σ_sim = mean(diag(cov(E1_d[:, :])))  # mean of covariances across simulations

    # Autocorrelation matrix
    X_d = zeros(N, T, 2)
    X_d0 = zeros(N, T-1, 2)
    X_d1 = zeros(N, T-1, 2)
    R_sim = zeros(2, 2)

    for ii in 1:2
        X_d[:, :, ii] = (XN[:, :, 2] - repeat(mean(XN[:, :, ii], 2), inner=[1, T]))
        X_d0[:, :, ii] = X_d[:, 1:end-1, ii]
        X_d1[:, :, ii] = (XN[:, 2:end, ii] - repeat(mean(XN[:, 2:end, ii], 2), inner=[1, T-1]))
        R_sim[ii, ii] = mean(diag( cov(X_d0[:, :, ii], X_d1[:, :, ii] )/sqrt.(var(X_d0[:, :, ii]))/sqrt.(var(X_d1[:, :, ii]))))
        if ii == 2
            R_sim[ii-1, ii] = mean(diag(cov(X_d0[:, :, ii-1], X_d1[:, :, ii])/sqrt.(var(X_d0[:, :, ii-1]))/sqrt.(var(X_d1[:, :, ii]))))
            R_sim[ii, ii-1] = mean(diag(cov(X_d0[:, :, ii], X_d1[:, :, ii-1])/sqrt.(var(X_d0[:, :, ii]))/sqrt.(var(X_d1[:, :, ii-1]))))
        end
    end

    return Mean_sim, Σ_sim, R_sim
end

#### ProductProcess

mutable struct ProductProcess{P1<:AbstractProcess,P2<:AbstractProcess} <: AbstractProcess
    process_1::P1
    process_2::P2
end

ProductProcess(p) = p

function discretize(pp::ProductProcess{ConstantProcess, <:IIDExogenous}; opt=Dict())
    diidp = discretize(pp.process_2)
    inodes = diidp.integration_nodes
    iit = hcat(
        vcat([pp.process_1.μ' for i=1:size(inodes,1)]...),
        inodes
    )
    return DiscretizedIIDProcess(diidp.grid, iit, diidp.integration_weights)
end

function discretize(::Type{DiscreteMarkovProcess}, pp::ProductProcess; opt1=Dict(), opt2=Dict())
    p1 = discretize(DiscreteMarkovProcess, pp.process_1; opt1...)
    p2 = discretize(DiscreteMarkovProcess, pp.process_2; opt2...)
    return MarkovProduct(p1,p2)
end

function discretize(::Type{GDP}, pp::ProductProcess; opt1=Dict(), opt2=Dict())
  p1 = discretize(GDP, pp.process_1; opt1...)
  p2 = discretize(GDP, pp.process_2; opt2...)
  return Product(p1,p2)
end

function discretize(pp::ProductProcess; kwargs...)
    return discretize(DiscreteMarkovProcess, pp; kwargs...)
end

## special processes (could be implemented as types later on)

function DeathProcess(μ::Float64)
    values = [0.0 1.0;]'
    transitions = [(1-μ) μ; 0 1]
    DiscreteMarkovProcess(transitions, values)
end
function PoissonProcess(μ::Float64, K::Int)
    values = (0:K)[:,:]*1.0
    transitions = zeros(K+1, K+1)
    for i=1:K
        transitions[i,i] = 1-μ
        transitions[i,i+1] = μ
    end
    transitions[K+1,K+1] = 1
    DiscreteMarkovProcess(transitions, values)
end
function AgingProcess(μ::Float64, K::Int)
    values = zeros(K+1,2)
    values[:,1] = 0:K
    values[1,2] = 1
    transitions = zeros(K+1, K+1)
    transitions[1,1] = 1
    for i=2:K
        transitions[i,i+1] = (1-μ)
        transitions[i,1] = μ
    end
    transitions[end,1] = 1
    dp = DiscreteMarkovProcess(transitions, values)
    dp.i0 = 2
    dp
end

get_integration_nodes(dprocess::Dolo.AbstractDiscretizedProcess, i::Int)= [(iweight(dprocess,i,j), inode(dprocess,i,j), j) for j in 1:n_inodes(dprocess,i) if iweight(dprocess,i,j)!=0]

get_integration_nodes(::typeof(Point), dprocess::Dolo.AbstractDiscretizedProcess, i::Int)=[(iweight(dprocess,i,j), inode(Point,dprocess,i,j), j) for j in 1:n_inodes(dprocess,i) if iweight(dprocess,i,j)!=0]

# type unstable
#
# get_integration_nodes(dprocess::Dolo.AbstractDiscretizedProcess, i::Int)=Iterators.filter( x -> (x[1]!=0), ((iweight(dprocess,i,j), inode(dprocess,i,j), j) for j in 1:n_inodes(dprocess,i)) )
#
# get_integration_nodes(::typeof(Point), dprocess::Dolo.AbstractDiscretizedProcess, i::Int)=Iterators.filter( x -> (x[1]!=0), ((iweight(dprocess,i,j), inode(Point,dprocess,i,j), j) for j in 1:n_inodes(dprocess,i)) )

# compatibility names
const AR1 = VAR1
const MarkovChain = DiscreteMarkovProcess
# const DiscretizedProcess

MarkovChain(;transitions=ones(1,1), values=[range(1,size(transitions,1))...]) = MarkovChain(transitions, values)

#### Domains

get_domain(mv::MvNormal) = EmptyDomain{length(mv.μ)}()

function get_domain(dmp::DiscreteMarkovProcess)
    points = dmp.grid.nodes
    d = ndims(dmp)
    return DiscreteDomain{d}(points)
end

function get_domain(cp::ConstantProcess)
    d = length(cp.μ)
    min  = fill(-Inf, d) ### not absolutely clear this is the right default
    max  = fill(Inf, d)
    return CartesianDomain(min, max)
end

function get_domain(iid::DiscretizedIIDProcess)
    return EmptyDomain()
end

function get_domain(var::VAR1)
    d = length(var.μ)
    return CartesianDomain(fill(-Inf, d), fill(Inf, d))
end

function get_domain(pp::ProductProcess)
    dom1 = get_domain(pp.process_1)
    dom2 = get_domain(pp.process_2)
    return ProductDomain(dom1, dom2)
end

# function discretize(pp::ProductProcess{ConstantProcess, VAR1}; opts...)
#     cp = pp.process_1
#     gdp = discretize(pp.process_2; opts...)
#     v = SVector(cp.μ...)
#     grid = [SVector(v..., e...) for e in gdp.grid.nodes]
#     integration_nodes = [cat(fill(cp.μ, size(e,1)), e; dims=2) for e in gdp.integration_nodes]
#     return DiscretizedProcess(grid, integration_nodes, gdp.integration_weights)

# end

function discretize(::Type{DiscreteMarkovProcess},pp::ProductProcess{ConstantProcess, VAR1}; opts...)
    mc = discretize(DiscreteMarkovProcess, pp.process_2; opts...)
    cp = pp.process_1
    d1 = length(cp.μ)
    d2 = length(node(mc.grid,1))
    v = SVector(cp.μ...)
    grid = [SVector(v..., e...) for e in mc.grid.nodes]
    transitions = mc.transitions
    values = copy(from_LOP(grid))
    return DiscreteMarkovProcess(transitions, values)
end

function discretize(::Type{GDP}, p::ConstantProcess)
    μ = p.μ
    grid = PointGrid(μ)
    integration_nodes = [repeat(μ',1)]
    integration_weights = [[1.0]]
    return GDP(grid, integration_nodes, integration_weights)
end


text:
module test_bukdu_plugs_csrf_protection

using Test
using Bukdu
import .Bukdu: Deps

req = Deps.Request()
conn = Conn(req)
@test applicable(plug, Plug.CSRF.Protection, conn)

end # module test_bukdu_plugs_csrf_protection

ext: jl

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
InputFolder = './Images/NucleiCytoo/';
OutputFolder = './Results/Images/NucleiCytoo/';

@iA = '*C00*.tif';

@fxg_mBlockOtsuThr [iA] > [M];
params.OpenRad = 2;
params.TopHatRad = 50;
params.BlckSize = 64;
params.BlckShft = 8;
params.MinRatio = 5;
params.Lvl = 0.5;
params.Bck = 0.1875;
params.AbsBck = 0.1;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
function sobolindices(
    models::Array{<:UQModel},
    inputs::Array{<:UQInput},
    output::Symbol,
    sim::AbstractMonteCarlo,
)
    sim_double_samples = @set sim.n = 2 * sim.n

    samples = sample(inputs, sim_double_samples)
    samples = samples[shuffle(1:size(samples, 1)), :]

    random_names = names(filter(i -> isa(i, RandomUQInput), inputs))

    evaluate!(models, samples)

    A = samples[1:(sim.n), :]
    B = samples[(sim.n + 1):end, :]

    fA = A[:, output]
    fB = B[:, output]

    fA .-= mean(fA)
    fB .-= mean(fB)

    VY = var([fA; fB])

    Si = zeros(length(random_names), 2)
    STi = zeros(length(random_names), 2)

    for (i, name) in enumerate(random_names)
        ABi = select(A, Not(name))
        ABi[:, name] = B[:, name]

        for m in models
            evaluate!(m, ABi)
        end

        ABi[:, output] .-= mean(ABi[:, output])

        first_order = x -> mean(fB .* (x .- fA)) / VY # Saltelli 2009
        total_effect = x -> (1 / (2 * sim.n)) * sum((fA .- x) .^ 2) / VY # Saltelli 2009

        # First order effects
        Si[i, 1] = first_order(ABi[:, output])
        bs = bootstrap(first_order, ABi[:, output], BasicSampling(1000))
        Si[i, 2] = stderror(bs)[1]

        # Total effects
        STi[i, 1] = total_effect(ABi[:, output])
        bs = bootstrap(total_effect, ABi[:, output], BasicSampling(1000))
        STi[i, 2] = stderror(bs)[1]
    end

    indices = DataFrame()
    indices.Variables = random_names
    indices.FirstOrder = Si[:, 1]
    indices.FirstOrderStdError = Si[:, 2]
    indices.TotalEffect = STi[:, 1]
    indices.TotalEffectStdError = STi[:, 2]

    return indices
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using Distributed
addprocs(4)

@everywhere begin
    using Pkg
    Pkg.activate(".")
end

sleep(0.1)

@everywhere begin
    using CSV
    using POMDPs
    using POMDPPolicies
    using POMDPSimulators
    using MultiAgentPOMDPs
    using FactoredValueMCTS
    using MultiAgentSysAdmin
    using MultiAgentSysAdmin: StarSysAdmin
    using Random
    using Statistics
    using Distributions
    include("pipeline.jl")
end

const MAX_STEPS = 100

function get_maxplus_solver(d, n, c, k, agent_utils, node_exp, edge_exp)
    solver = FVMCTSSolver(
        depth=d,
        n_iterations=n,
        exploration_constant=c,
        rng=MersenneTwister(8),
        coordination_strategy = MaxPlus(
            message_iters=k,
            message_norm=true,
            use_agent_utils=agent_utils,
            node_exploration=node_exp,
            edge_exploration=edge_exp,
        )
    )
    return solver
end

function get_varel_solver(d, n, c)
    solver = FVMCTSSolver(depth=d,
                             n_iterations=n,
                             exploration_constant=c,
                             rng=MersenneTwister(8),
                             coordination_strategy=VarEl())
    return solver
end

const PSET = Dict(
    "ring_global" => BiSysAdmin{true},
    "ring_local" => BiSysAdmin{false},
    "star_global" => StarSysAdmin{true},
    "star_local" => StarSysAdmin{false},
    "ringofring_global" => RingofRingSysAdmin{true},
    "ringofring_local" => RingofRingSysAdmin{false},
)

function run_experiment(prob, nagents, solver_name, nevals, d, n, c, k, logdir)
    if startswith(prob, "ringofring")
        # we build rings of rings where each ring is of size 3!
        @assert nagents % 3 == 0
        mdp = PSET[prob](nagents=nagents//3, nagents_per_ring=3)
    else
        mdp = PSET[prob](nagents=nagents)
    end

    exp_name = "SysAdmin_$(prob)_nagents=$(nagents)_$(solver_name)"
    if solver_name == :random
        solver = RandomSolver(rng=MersenneTwister(1999))
        @info exp_name
        df = online_evaluate(Dict(solver_name=>solver), mdp, nevals, MAX_STEPS)
        CSV.write(joinpath(logdir, "$(exp_name).csv"), df)
        for k in unique(df.policy)
            fdf = filter(r->r[:policy] == k, df)
            vals = fdf.discret
            uvals = fdf.undiscret
            @info "discounted" k mean=mean(vals) std=std(vals)
            @info "undiscounted" k mean=mean(uvals) std=std(uvals)
        end
        return nothing
    end

    if solver_name == :maxplus
        exp_name *= "_d=$(d),n=$(n),c=$(c),k=$(k)"
        solver = get_maxplus_solver(d, n, c, k, true, true, false)
    elseif solver_name == :varel
        exp_name *= "_d=$(d),n=$(n),c=$(c)"
        solver = get_varel_solver(d, n, c)
    end
    @info exp_name
    df = online_evaluate(Dict(solver_name=>solver), mdp, nevals, MAX_STEPS)
    CSV.write(joinpath(logdir, "$(exp_name).csv"), df)
    for k in unique(df.policy)
        fdf = filter(r->r[:policy] == k, df)
        vals = fdf.discret
        uvals = fdf.undiscret
        @info "discounted" k mean=mean(vals) std=std(vals)
        @info "undiscounted" k mean=mean(uvals) std=std(uvals)
    end
end

function main(args)
    prob = args[1]

    nagents = parse(Int, args[2])

    solver = Symbol(args[3])

    nevals = parse(Int, args[4])

    logdir = args[5]

    !isdir(logdir) && mkpath(logdir)

    if solver == :random
        d = 0
        n = 0
        c = 0
        k = 0
    else
        d = parse(Int, args[6])

        n = parse(Int, args[7])

        c = parse(Int, args[8])

        k = parse(Int, args[9])
    end
    run_experiment(prob, nagents, solver, nevals, d, n, c, k, logdir)
end

(abspath(PROGRAM_FILE) == @__FILE__) && main(ARGS)


text:
using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays

function benchmark_results()
    f = [] # Example functions from Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application
    push!(f, x->((x + sin(x)) * exp(-x^2)))
    push!(f, x->(x^4 - 10x^3 + 35x^2 - 50x + 24))
    push!(f, x->((log(x + 1.25) - 0.84x) ^ 2))
    push!(f, x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)))
    push!(f, x->(exp(x^2)))
    push!(f, x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)))
    push!(f, x->(x^6 - 15x^4 + 27x^2 + 250))
    push!(f, x->(atan(cos(tan(x)))))
    push!(f, x->(asin(cos(acos(sin(x))))))

    s = interval(0.75, 1.75)
    df = DataFrame()

    df[:Method] = ["Automatic Differentiation", "Slope Expansion"]
    for n in 1:length(f)

        t1 = ForwardDiff.derivative(f[n], s)
        t2 = slope(f[n], s, mid(s))
        df[Symbol("f" * "$n")] = [t1, t2]
    end
    a = []
    for i in 1:length(f)
        push!(a, Symbol("f" * "$i"))
    end
    df1 = stack(df, a)
    dfnew = unstack(df1, :variable, :Method, :value)
    dfnew = rename(dfnew, :variable => :Function)
    println(dfnew)
    dfnew
end

function benchmark_time()
    f = []
    push!(f, x->((x + sin(x)) * exp(-x^2)))
    push!(f, x->(x^4 - 10x^3 + 35x^2 - 50x + 24))
    push!(f, x->((log(x + 1.25) - 0.84x) ^ 2))
    push!(f, x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)))
    push!(f, x->(exp(x^2)))
    push!(f, x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)))
    push!(f, x->(x^6 - 15x^4 + 27x^2 + 250))
    push!(f, x->(atan(cos(tan(x)))))
    push!(f, x->(asin(cos(acos(sin(x))))))

    s = interval(0.75, 1.75)
    df = DataFrame()
    df[:Method] = ["Automatic Differentiation", "Slope Expansion"]
    for n in 1:length(f)

        t1 = @belapsed ForwardDiff.derivative($f[$n], $s)
        t2 = @belapsed slope($f[$n], $s, mid($s))
        df[Symbol("f" * "$n")] = [t1, t2]
    end
    a = []
    for i in 1:length(f)
        push!(a, Symbol("f" * "$i"))
    end
    df1 = stack(df, a)
    dfnew = unstack(df1, :variable, :Method, :value)
    dfnew = rename(dfnew, :variable => :Function)
    println(dfnew)
    dfnew
end


text:
@testset "MOI Tests" begin
    @testset "Supports and SolverName" begin
        optimizer = CSTestOptimizer(; branch_strategy = :ABS)
        @test MOI.get(optimizer, MOI.SolverName()) == "ConstraintSolver"
        @test MOI.supports_constraint(
            optimizer,
            MOI.VectorOfVariables,
            CS.CPE.AllDifferent,
        )
        @test MOI.supports_constraint(optimizer, MOI.VectorOfVariables, CS.CPE.AllEqual)
        @test MOI.supports_constraint(optimizer, MOI.VectorOfVariables, CS.TableSetInternal)
        @test MOI.supports_constraint(optimizer, MOI.SingleVariable, MOI.ZeroOne)
        @test MOI.supports_constraint(optimizer, MOI.SingleVariable, MOI.Integer)
        @test MOI.supports_constraint(optimizer, MOI.SingleVariable, MOI.LessThan{Float64})
        @test MOI.supports_constraint(optimizer, MOI.SingleVariable, CS.Integers)
        @test MOI.supports_constraint(
            optimizer,
            MOI.SingleVariable,
            MOI.GreaterThan{Float64},
        )
        @test MOI.supports_constraint(optimizer, MOI.SingleVariable, MOI.EqualTo{Float64})
        @test MOI.supports_constraint(optimizer, MOI.SingleVariable, MOI.Interval{Float64})

        @test MOI.supports_constraint(
            optimizer,
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        )
        @test MOI.supports_constraint(
            optimizer,
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        )
        @test MOI.supports_constraint(
            optimizer,
            MOI.ScalarAffineFunction{Float64},
            CPE.DifferentFrom{Float64},
        )

        f = MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, MOI.VariableIndex(1))),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, MOI.VariableIndex(2))),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, MOI.VariableIndex(3))),
            ],
            [0.0, 0.0],
        )
        indicator_set = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(9.0))
        @test MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))
        indicator_set = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(9.0))
        @test MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))
        indicator_set = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(9.0))
        @test MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))

        indicator_set = CS.IndicatorSet{MOI.ACTIVATE_ON_ONE, typeof(f)}(CS.CPE.AllDifferent(2))
        @test !MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))
        indicator_set = CS.IndicatorSet{MOI.ACTIVATE_ON_ZERO, typeof(f)}(CS.TableSetInternal(2, [1 2; ]))
        @test !MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))

        @test MOI.supports_constraint(
            optimizer,
            typeof(f),
            CS.ReifiedSet{MOI.ACTIVATE_ON_ONE,MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}},
        )

        @test MOI.supports_constraint(
            optimizer,
            MOI.VectorOfVariables,
            CS.ReifiedSet{MOI.ACTIVATE_ON_ZERO,MOI.VectorOfVariables,CS.CPE.AllDifferent},
        )

        # TimeLimit
        @test MOI.supports(optimizer, MOI.TimeLimitSec())

        # objective
        @test MOI.supports(optimizer, MOI.ObjectiveSense())
        @test MOI.supports(optimizer, MOI.ObjectiveFunction{MOI.SingleVariable}())
        @test MOI.supports(
            optimizer,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        )

        @test MOI.supports(optimizer, MOI.RawParameter("backtrack"))

        @test MOI.get(optimizer, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    end

    @testset "Small MOI tests" begin
        optimizer = CSTestOptimizer()
        @assert MOI.supports_constraint(
            optimizer,
            MOI.VectorOfVariables,
            CS.CPE.AllDifferent,
        )

        x1 = MOI.add_constrained_variable(optimizer, MOI.ZeroOne())
        x2 = MOI.add_constrained_variable(optimizer, MOI.Interval(1.0, 2.0))
        x3 = MOI.add_constrained_variable(optimizer, MOI.Interval(2.0, 3.0))

        # [1] to get the index
        MOI.add_constraint(optimizer, x2[1], MOI.Integer())
        MOI.add_constraint(optimizer, x3[1], MOI.Integer())

        affine_terms = Vector{MOI.ScalarAffineTerm{Float64}}()
        push!(affine_terms, MOI.ScalarAffineTerm{Float64}(1, MOI.VariableIndex(1)))
        push!(affine_terms, MOI.ScalarAffineTerm{Float64}(1, MOI.VariableIndex(2)))
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(affine_terms, 0.0),
            MOI.EqualTo(3.0),
        )
        MOI.add_constraint(
            optimizer,
            MOI.VectorOfVariables([x1[1], x2[1], x3[1]]),
            CS.CPE.AllDifferent(3),
        )

        MOI.optimize!(optimizer)
        @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(optimizer, MOI.VariablePrimal(), x1[1]) == 1
        @test MOI.get(optimizer, MOI.VariablePrimal(), x2[1]) == 2
        @test MOI.get(optimizer, MOI.VariablePrimal(), x3[1]) == 3
    end

    @testset "ErrorHandling" begin
        optimizer = CSTestOptimizer()
        @test_throws ErrorException MOI.add_constrained_variable(
            optimizer,
            MOI.Interval(NaN, 2.0),
        )
        @test_throws ErrorException MOI.add_constrained_variable(
            optimizer,
            MOI.Interval(1.0, NaN),
        )

        optimizer = CSTestOptimizer()
        x = MOI.add_constrained_variable(optimizer, MOI.Interval(1.0, 2.0))
        @test_logs (:error, r"Upper bound .* exists") MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(MOI.VariableIndex(1)),
            MOI.LessThan(2.0),
        )
        @test_logs (:error, r"Lower bound .* exists") MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(MOI.VariableIndex(1)),
            MOI.GreaterThan(1.0),
        )

        optimizer = CSTestOptimizer()
        x1 = MOI.add_constrained_variable(optimizer, MOI.Interval(1.0, 2.0))
        x2 = MOI.add_constrained_variable(optimizer, MOI.Interval(1.0, 2.0))
        # All should be integer
        @test_throws ErrorException MOI.optimize!(optimizer)

        # non existing option
        @test_logs (:error, r"option abc doesn't exist") Model(optimizer_with_attributes(
            CS.Optimizer,
            "abc" => 1,
        ))
        @test_logs (:error, r"option abc doesn't exist") model = CS.Optimizer(abc = 1)
    end

end


text:
module TestDeprecations

using IterativeSolvers
using Test

# Currently nothing here.
@test true

end # module


text:
"""
    module Loggables

Defines [`Loggables.Loggable`](#) and its subtypes.
"""
module Loggables

"""
    abstract type Loggable

Abstract type for data that [`LoggerBackend`](#)s can log.
See `subtypes(FluxTraining.Loggables.Loggable)` and [`LoggerBackend`](#)

"""
abstract type Loggable end

struct Text <: Loggable data end
struct Image <: Loggable data end
struct Audio <: Loggable data end
struct Histogram <: Loggable data end
struct Value <: Loggable data end
struct Graph <: Loggable data end
struct File <: Loggable
    file
    data
end

end  # module


text:
#### Code to create figures
#### "Slow organisms exhibit sudden population disappearances in a reddened world" by Greyson-Gaito, Gellner, & McCann.

include("packages.jl")
include("slowfast_commoncode.jl")

# Figure 2 (primer of different trajectories)
let
    sol_axial = RozMac_pert(0.01, 0.48, 1, 0.8, 9, 2500.0, 0.0:2.0:2500.0)
    sol_qc = RozMac_pert(0.01, 0.48, 1, 0.8, 6, 3600.0, 2500.0:1.0:3100.0)
    sol_equil = RozMac_pert(0.01, 0.48, 1, 0.0, 1234, 5000.0, 0.0:2.0:5000.0)
    figure3 = figure(figsize=(3,10))
    subplot(5,1,1)
    iso_plot(0.0:0.1:3.0, RozMacPar(e = 0.48))
    plot(sol_qc[1, :], sol_qc[2, :], color = "#440154FF")
    plot(sol_axial[1, :], sol_axial[2, :], color = "#73D055FF")
    xlabel("Resource")
    ylabel("Consumer")
    ylim(0,2.5)
    xlim(0,3.0)
    subplot(5,1,2)
    iso_plot(0.0:0.1:3.0, RozMacPar(e = 0.48))
    plot(sol_equil[1, :], sol_equil[2, :], color = "#FDE725FF")
    xlabel("Resource")
    ylabel("Consumer")
    ylim(0,2.5)
    xlim(0,3.0)
    subplot(5,1,3)
    plot(sol_qc.t, sol_qc.u)
    ylim(0,3.0)
    xlabel("Time")
    ylabel("Biomass")
    tick_params(axis="x", which="both", bottom=False, labelbottom=False)
    subplot(5,1,4)
    plot(sol_axial.t, sol_axial.u)
    ylim(0,3.0)
    tick_params(axis="x", which="both", bottom=False, labelbottom=False)
    xlabel("Time")
    ylabel("Biomass")
    subplot(5,1,5)
    plot(sol_equil.t, sol_equil.u)
    tick_params(axis="x", which="both", bottom=False, labelbottom=False)
    xlabel("Time")
    ylabel("Biomass")
    ylim(0,3.0)
    tight_layout()
    # return figure2
    savefig(joinpath(abpath(), "figs/phase_timeseries_examples.pdf"))
end

# Figure 3 (Proportion Real with small delta ε & ACF plots of quasi-cycles)
RCdividedata = findRCdivide_epx_data()
let
    prop_real_smalldelep = figure(figsize=(4,3))
    plot(RCdividedata[:,1], RCdividedata[:,2], color = "black")
    # fill_between(data[:,1], data[:,2], color = "blue", alpha=0.3)
    # fill_between(data[:,1], fill(1.0, length(data[:,2])), data[:,2], color = "orange", alpha=0.3)
    xlim(0,10)
    ylim(0.0,1.0)
    xticks(fontsize=12)
    yticks(fontsize=12)
    ylabel("Proportion Real", fontsize = 12)
    xlabel("1/ε", fontsize = 12)
    # return prop_real_smalldelep
    savefig(joinpath(abpath(), "figs/epsilonxaxis_propReal_smalldelep.pdf"))
end

let
    lrange = 0:1:40
    acffast = quasicycle_data(1.0, 1000)
    acfslow = quasicycle_data(0.6, 1000)
    acffigure = figure(figsize = (6,2.5))
    subplot(1,2,1)
    plot(collect(lrange), acffast, color = "black")
    hlines(0.0,0.0,40.0, linestyles=`dashed`, linewidths=0.5)
    ylim(-1,1)
    xlim(0,40)
    xlabel("Lag")
    ylabel("Average ACF")
    subplot(1,2,2)
    plot(collect(lrange), acfslow, color = "black")
    hlines(0.0,0.0,40.0, linestyles=`dashed`, linewidths=0.5)
    ylim(-1,1)
    xlim(0,40)
    xlabel("Lag")
    ylabel("Average ACF")
    tight_layout()
    # return acffigure
    savefig(joinpath(abpath(), "figs/quasicycles_ACF.pdf"))
end

# Figure 4 (Proportion real with large delta ε & quasi-canard proportions with white noise )
wn_data05_short_RozMac = CSV.read(joinpath(abpath(),"data/wn_eff05_short_RozMac.csv"), DataFrame)
wn_data06_short_RozMac = CSV.read(joinpath(abpath(),"data/wn_eff06_short_RozMac.csv"), DataFrame)
wn_data07_short_RozMac = CSV.read(joinpath(abpath(),"data/wn_eff07_short_RozMac.csv"), DataFrame)
wn_data05_long_RozMac = CSV.read(joinpath(abpath(),"data/wn_eff05_long_RozMac.csv"), DataFrame)
wn_data06_long_RozMac = CSV.read(joinpath(abpath(),"data/wn_eff06_long_RozMac.csv"), DataFrame)
wn_data07_long_RozMac = CSV.read(joinpath(abpath(),"data/wn_eff07_long_RozMac.csv"), DataFrame)

let
    prop_real_largedelep = figure(figsize=(4,3))
    plot(RCdividedata[:,1], RCdividedata[:,2], color = "black")
    hlines(converteff_prop_RCdivide(0.5), 0.0, 1000.0, linestyles=`dashed`, linewidths=0.5, color = "black")
    hlines(converteff_prop_RCdivide(0.65), 0.0, 1000.0, linestyles=`dashed`, linewidths=0.5, color = "black")
    hlines(converteff_prop_RCdivide(0.7), 0.0, 1000.0, linestyles=`dashed`, linewidths=0.5, color = "black")
    xlim(-10,1000)
    ylim(0.0,1.0)
    xticks(fontsize=12)
    yticks(fontsize=12)
    ylabel("Proportion Real", fontsize = 15)
    xlabel("1/ε", fontsize = 15)
    # return prop_real_largedelep
    savefig(joinpath(abpath(), "figs/epsilonxaxis_propReal_largedelep.pdf"))
end

let
    figure4bc = figure(figsize = (7,2.5))
    subplot(1,3,1)
    plot(1 ./ wn_data05_short_RozMac.xrange, wn_data05_short_RozMac.canard, color = "black", linestyle = "dashed")
    plot(1 ./ wn_data05_long_RozMac.xrange, wn_data05_long_RozMac.canard, color = "black", linestyle = "dotted")
    ylabel("Proportion with\nquasi-canard", fontsize = 15)
    xlabel("1/ϵ", fontsize = 15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(1,3,2)
    plot(1 ./ wn_data06_short_RozMac.xrange, wn_data06_short_RozMac.canard, color = "black", linestyle = "dashed")
    plot(1 ./ wn_data06_long_RozMac.xrange, wn_data06_long_RozMac.canard, color = "black", linestyle = "dotted")
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    xlabel("1/ϵ", fontsize = 15)
    subplot(1,3,3)
    plot(1 ./ wn_data07_short_RozMac.xrange, wn_data07_short_RozMac.canard, color = "black", linestyle = "dashed")
    plot(1 ./ wn_data07_long_RozMac.xrange, wn_data07_long_RozMac.canard, color = "black", linestyle = "dotted")
    xlabel("1/ϵ", fontsize = 15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    tight_layout()
    # return figure4bc
    savefig(joinpath(abpath(), "figs/canard_whitenoise_prop.pdf"))
end

# Figure 5 (Proportion quasicanard/axial solution with red noise)

rn_data_ep0004eff05_RozMac = CSV.read(joinpath(abpath(),"data/rn_ep0004eff05_RozMac.csv"), DataFrame)
rn_data_ep0004eff06_RozMac = CSV.read(joinpath(abpath(),"data/rn_ep0004eff06_RozMac.csv"), DataFrame)
rn_data_ep0004eff07_RozMac = CSV.read(joinpath(abpath(),"data/rn_ep0004eff07_RozMac.csv"), DataFrame)
rn_data_ep0079eff05_RozMac = CSV.read(joinpath(abpath(),"data/rn_ep0079eff05_RozMac.csv"), DataFrame)
rn_data_ep0079eff06_RozMac = CSV.read(joinpath(abpath(),"data/rn_ep0079eff06_RozMac.csv"), DataFrame)
rn_data_ep0079eff07_RozMac = CSV.read(joinpath(abpath(),"data/rn_ep0079eff07_RozMac.csv"), DataFrame)

let
    figure5 = figure(figsize = (8,4))
    subplot(2,3,1)
    fill_between(rn_data_ep0079eff05_RozMac.xrange, rn_data_ep0079eff05_RozMac.canard, color="#440154FF")
    fill_between(rn_data_ep0079eff05_RozMac.xrange, rn_data_ep0079eff05_RozMac.canard, rn_data_ep0079eff05_RozMac.canard_plus_axial, color="#73D055FF")
    fill_between(rn_data_ep0079eff05_RozMac.xrange, fill(1.0,91), rn_data_ep0079eff05_RozMac.canard_plus_axial, color="#FDE725FF")
    ylabel("Proportion", fontsize=15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,2)
    fill_between(rn_data_ep0079eff06_RozMac.xrange, rn_data_ep0079eff06_RozMac.canard, color="#440154FF")
    fill_between(rn_data_ep0079eff06_RozMac.xrange, rn_data_ep0079eff06_RozMac.canard, rn_data_ep0079eff06_RozMac.canard_plus_axial, color="#73D055FF")
    fill_between(rn_data_ep0079eff06_RozMac.xrange, fill(1.0,91), rn_data_ep0079eff06_RozMac.canard_plus_axial, color="#FDE725FF")
    xlabel("Noise correlation", fontsize=15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,3)
    fill_between(rn_data_ep0079eff07_RozMac.xrange, rn_data_ep0079eff07_RozMac.canard, color="#440154FF")
    fill_between(rn_data_ep0079eff07_RozMac.xrange, rn_data_ep0079eff07_RozMac.canard, rn_data_ep0079eff07_RozMac.canard_plus_axial, color="#73D055FF")
    fill_between(rn_data_ep0079eff07_RozMac.xrange, fill(1.0,91), rn_data_ep0079eff07_RozMac.canard_plus_axial, color="#FDE725FF")
    ylim(0,1)

    subplot(2,3,4)
    fill_between(rn_data_ep0004eff05_RozMac.xrange, rn_data_ep0004eff05_RozMac.canard, color="#440154FF")
    fill_between(rn_data_ep0004eff05_RozMac.xrange, rn_data_ep0004eff05_RozMac.canard, rn_data_ep0004eff05_RozMac.canard_plus_axial, color="#73D055FF")
    fill_between(rn_data_ep0004eff05_RozMac.xrange, fill(1.0,91), rn_data_ep0004eff05_RozMac.canard_plus_axial, color="#FDE725FF")
    ylabel("Proportion", fontsize=15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,5)
    fill_between(rn_data_ep0004eff06_RozMac.xrange, rn_data_ep0004eff06_RozMac.canard, color="#440154FF")
    fill_between(rn_data_ep0004eff06_RozMac.xrange, rn_data_ep0004eff06_RozMac.canard, rn_data_ep0004eff06_RozMac.canard_plus_axial, color="#73D055FF")
    fill_between(rn_data_ep0004eff06_RozMac.xrange, fill(1.0,91), rn_data_ep0004eff06_RozMac.canard_plus_axial, color="#FDE725FF")
    xlabel("Noise correlation", fontsize=15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,6)
    fill_between(rn_data_ep0004eff07_RozMac.xrange, rn_data_ep0004eff07_RozMac.canard, color="#440154FF")
    fill_between(rn_data_ep0004eff07_RozMac.xrange, rn_data_ep0004eff07_RozMac.canard, rn_data_ep0004eff07_RozMac.canard_plus_axial, color="#73D055FF")
    fill_between(rn_data_ep0004eff07_RozMac.xrange, fill(1.0,91), rn_data_ep0004eff07_RozMac.canard_plus_axial, color="#FDE725FF")
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    tight_layout()
    annotate("a)", (20, 280), xycoords = "figure points", fontsize = 15)
    annotate("b)", (205, 280), xycoords = "figure points", fontsize = 15)
    annotate("c)", (385, 280), xycoords = "figure points", fontsize = 15)
    annotate("d)", (20, 140), xycoords = "figure points", fontsize = 15)
    annotate("e)", (205, 140), xycoords = "figure points", fontsize = 15)
    annotate("f)", (385, 140), xycoords = "figure points", fontsize = 15)
    # return figure5
    savefig(joinpath(abpath(), "figs/canard_rednoise_prop.pdf"))
end

# Box 2 Figure 1
let
    figure6 = figure(figsize = (3,4))
    subplot(2,1,1)
    roz_mac_plot(0.01, 0.48, 2, 0.3)
    subplot(2,1,2)
    roz_mac_plot(0.01, 0.48, 2, 0.3)
    tight_layout()
    # return figure6
    savefig(joinpath(abpath(), "figs/ratchet_schematic.pdf"))
end

## Supporting Information

#SI Figure 1 (illustration of quasi-canard finder algorithm)
let
    sol_qc = RozMac_pert(0.01, 0.55, 1, 0.7, 45, 1500.0, 0.0:1.0:1500.0)
    figure2 = figure(figsize=(7,5))
    iso_plot(0.0:0.1:3.0, RozMacPar(e = 0.55))
    plot(sol_qc[1, :], sol_qc[2, :], color = "#440154FF", linewidth=5)
    xlabel("Resource", fontsize=15)
    ylabel("Consumer", fontsize=15)
    ylim(0, 2.5)
    xlim(0, 3)
    # return figure2
    savefig(joinpath(abpath(), "figs/quasicanardfinder.pdf"))
end

#SI Figure 2  (Isoclines for Rosenzweig-MacArthur consumer-resource model)
let
    isoclines = figure(figsize = (6,2.5))
    subplot(1,3,1)
    iso_plot(0.0:0.1:3.0, RozMacPar(e = 0.5))
    title("Non-excitable\n(Real λ)", fontsize=15)
    ylim(0,2.5)
    xlim(0.0,3.0)
    xticks(fontsize=12)
    yticks(fontsize=12)
    ylabel("Consumer", fontsize=15)
    xlabel("Resource", fontsize=15)
    subplot(1,3,2)
    iso_plot(0.0:0.1:3.0, RozMacPar(e = 0.65))
    title("Excitable\n(Complex λ)", fontsize=15)
    xlabel("Resource", fontsize=15)
    ylim(0,2.5)
    xlim(0.0,3.0)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(1,3,3)
    iso_plot(0.0:0.1:3.0, RozMacPar(e = 0.7))
    title("Excitable\n(Near Hopf)", fontsize=15)
    xlabel("Resource", fontsize=15)
    ylim(0,2.5)
    xlim(0.0,3.0)
    xticks(fontsize=12)
    yticks(fontsize=12)
    tight_layout()
    # return isoclines
    savefig(joinpath(abpath(), "figs/isoclinesSI.pdf"))
end

# SI Figure 3 (Yodiz & Innes model with white noise)
wn_R12_short_YodInn = CSV.read(joinpath(abpath(),"data/wn_R12_short_YodInn.csv"), DataFrame)
wn_R10_short_YodInn = CSV.read(joinpath(abpath(),"data/wn_R10_short_YodInn.csv"), DataFrame)
wn_R08_short_YodInn = CSV.read(joinpath(abpath(),"data/wn_R08_short_YodInn.csv"), DataFrame)
wn_R12_long_YodInn = CSV.read(joinpath(abpath(),"data/wn_R12_long_YodInn.csv"), DataFrame)
wn_R10_long_YodInn = CSV.read(joinpath(abpath(),"data/wn_R10_long_YodInn.csv"), DataFrame)
wn_R08_long_YodInn = CSV.read(joinpath(abpath(),"data/wn_R08_long_YodInn.csv"), DataFrame)

let
    SIfigure3 = figure(figsize = (10,6))
    subplot(2,3,1)
    iso_plot_YodInn(0.0:0.1:3.0, YodInnScalePar(R₀ = 1.2))
    ylim(0,0.5)
    xlim(0.0,3.0)
    ylabel("Consumer", fontsize=15)
    xlabel("Resource", fontsize=15)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,2)
    iso_plot_YodInn(0.0:0.1:3.0, YodInnScalePar(R₀ = 1.0))
    ylim(0,0.5)
    xlim(0.0,3.0)
    ylabel("Consumer", fontsize=15)
    xlabel("Resource", fontsize=15)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,3)
    iso_plot_YodInn(0.0:0.1:3.0, YodInnScalePar(R₀ = 0.8))
    ylim(0,0.5)
    xlim(0.0,3.0)
    ylabel("Consumer", fontsize=15)
    xlabel("Resource", fontsize=15)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,4)
    plot(1 ./ wn_R12_short_YodInn.xrange, wn_R12_short_YodInn.canard, color = "black", linestyle = "dashed")
    plot(1 ./ wn_R12_long_YodInn.xrange, wn_R12_long_YodInn.canard, color = "black", linestyle = "dotted")
    ylabel("Proportion with\nquasi-canard", fontsize=15)
    xlabel("1/ϵ", fontsize=15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    subplot(2,3,5)
    plot(1 ./ wn_R10_short_YodInn.xrange, wn_R10_short_YodInn.canard, color = "black", linestyle = "dashed")
    plot(1 ./ wn_R10_long_YodInn.xrange, wn_R10_long_YodInn.canard, color = "black", linestyle = "dotted")
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    xlabel("1/ϵ", fontsize=15)
    subplot(2,3,6)
    plot(1 ./ wn_R08_short_YodInn.xrange, wn_R08_short_YodInn.canard, color = "black", linestyle = "dashed")
    plot(1 ./ wn_R08_long_YodInn.xrange, wn_R08_long_YodInn.canard, color = "black", linestyle = "dotted")
    xlabel("1/ϵ", fontsize=15)
    ylim(0,1)
    xticks(fontsize=12)
    yticks(fontsize=12)
    tight_layout()
    # return SIfigure3
    savefig(joinpath(abpath(), "figs/canard_whitenoise_prop_YodInn.pdf"))
end

# Code to make resource isocline data for quasipotential figures (created in R)
resiso_data_05 = DataFrame(xrange = 0.0:0.01:3.0, resiso = [res_iso_roz_mac(R, RozMacPar(e = 0.5)) for R in 0.0:0.01:3.0])
CSV.write(joinpath(abpath(), "data/resiso_data_05.csv"), resiso_data_05)
con_iso_roz_mac(RozMacPar(e = 0.5))

resiso_data_07 = DataFrame(xrange = 0.0:0.01:3.0, resiso = [res_iso_roz_mac(R, RozMacPar(e = 0.7)) for R in 0.0:0.01:3.0])
CSV.write(joinpath(abpath(), "data/resiso_data_07.csv"), resiso_data_07)
con_iso_roz_mac(RozMacPar(e = 0.7))


text:
function trial_info!(df)
    transform!(groupby(df,[:MOUSE,:DATE, :TRIAL,:SIDE]),
        :IN => (x -> collect(1:length(x))) => :POKE_TRIAL,
        :IN => (x-> Int64.(vcat(zeros(length(x)-1), [1]))) => :LEAVE,
        :REWARD => (x -> pushfirst!(Int64.(cumsum(x)[1:end-1].+1),1)) => :BOUT_TRIAL,
        :IN => (i -> round.(i .- i[1], digits = 5)) => :IN_TRIAL,
        [:IN, :OUT] => ((i,o) -> round.(o .- i[1], digits = 5)) => :OUT_TRIAL,
    )
end

function bout_info!(df)
    transform!(groupby(df,[:MOUSE,:DATE, :TRIAL,:SIDE, :BOUT_TRIAL]),
        :IN => (i -> round.(i .- i[1], digits = 5)) => :IN_BOUT,
        [:IN, :OUT] => ((i,o) -> round.(o .- i[1], digits = 5)) => :OUT_BOUT,
        [:IN, :OUT] => ((i,o) -> round.(o .- i, digits = 5)) => :DURATION,
    )
end

function travel_info!(df)
    trav_df = combine(groupby(df,[:MOUSE,:DATE, :TRIAL]),
        :KIND => last => :TRAVEL
    )
    trav_df.TRAVEL = [x in ["short", "long"] ? x : "last" for x in trav_df.TRAVEL]
    leftjoin!(df, trav_df, on = [:MOUSE,:DATE,:TRIAL])
end

function leave_info!(df)
    transform!(groupby(df,[:MOUSE,:DATE, :TRIAL,:SIDE]),
        :IN => (x-> vcat(falses(length(x)-1), [true])) => :LEAVE
    )
end

function count_bouts(R,L,S)
    T = Bool.(S .== "travel")
    v = (R .| L) .& (.!T)
    pushfirst!(Int64.(cumsum(v)[1:end-1].+1),1)
end

function find_bouts!(df)
    df.REWARD = Bool.(df.REWARD)
    transform!(groupby(df, [:MOUSE,:DATE]),
        [:REWARD,:LEAVE, :SIDE] => ((R,L,S) -> count_bouts(R,L,S)) => :BOUT
    )
end

function cleaned_times!(df)
    transform!(groupby(df,[:MOUSE,:DATE]),
        :DURATION => (t -> round.(cumsum(t), digits = 5)) => :TIME
    )
    transform!(groupby(df,[:MOUSE,:DATE, :TRIAL]),
        :DURATION => (t -> round.(cumsum(t), digits = 5)) => :TIME_TRIAL
    )
    transform!(groupby(df,[:MOUSE,:DATE, :TRIAL, :BOUT]),
        :DURATION => (t -> round.(cumsum(t), digits = 5)) => :TIME_BOUT
    )
end

function add_time_bins!(df; binsize = 0.1)
    transform!(df,:DURATION => (ByRow(d -> Int64(round(d/binsize)))) => :BIN)
    # transform!(df,:DURATION => (ByRow(d -> Int64(ceil(d/binsize)))) => :BIN)
    filter!(r-> (r.BIN >= 1 || Bool(r.REWARD) || Bool(r.LEAVE)), df)
    replace!(df.BIN, 0 => 1)
    # filter!(r-> r.BIN >= 1, df)

end

function preprocess_pokes!(df; binsize = 0.1)
    trial_info!(df)
    bout_info!(df)
    travel_info!(df)
    find_bouts!(df)
    cleaned_times!(df)
    transform!(df,:REWARD => ByRow(Int64) => :REWARD)
    select!(df,[:MOUSE,:DATE,:SIDE,:KIND,:TRAVEL,:BOUT,:TRIAL,:BOUT_TRIAL,:POKE_TRIAL,:REWARD,:LEAVE,
        :IN, :OUT, :IN_TRIAL, :OUT_TRIAL, :IN_BOUT, :OUT_BOUT, :DURATION,:TIME, :TIME_TRIAL, :TIME_BOUT]);
    add_time_bins!(df; binsize = binsize)
    transform!(groupby(df,[:MOUSE,:DATE]),
        [:REWARD,:BIN] => ((r,b) -> sum(r)/sum(b)) => :ENV_INITIALVALUE,
        :TRIAL => (t -> Int64.(vcat([1],zeros(length(t)-1)))) => :NEWSESSION
        )
    transform!(groupby(df,[:MOUSE,:DATE,:TRIAL]),
        :REWARD => (t -> Int64.(vcat([1],zeros(length(t)-1)))) => :NEWTRIAL
        )
end

function preprocess_bouts(df)
    bouts = filter(r -> r.SIDE != "travel", df)
    transform!(groupby(bouts,[:MOUSE,:DATE,:BOUT]),
        :DURATION => (t -> round.(cumsum(t), digits = 5)) => :WAITING_TIME
    )
    filter!(r -> (Bool(r.LEAVE) .| Bool(r.REWARD)), bouts)
    bouts = bouts[:,[:DATE,:MOUSE,:SIDE,:KIND,:TRAVEL,:TRIAL, :BOUT,:REWARD,:LEAVE,:DURATION,:TIME,:WAITING_TIME,:IN,:OUT,:IN_TRIAL,:OUT_TRIAL]]
    return bouts
end
Int64.(trues(4))


text:
# Autogenerated wrapper script for BlueZ_jll for i686-linux-musl
export bccmd, bluemoon, bluetooth, bluetoothctl, bluetoothd, btattach, btmon, hex2hcd, hid2hci, l2ping, l2test, libbluetooth, mpris_proxy, obexd, rctest

using Dbus_jll
using eudev_jll
using Glib_jll
using Libical_jll
using Readline_jll
JLLWrappers.@generate_wrapper_header("BlueZ")
JLLWrappers.@declare_executable_product(bccmd)
JLLWrappers.@declare_executable_product(bluemoon)
JLLWrappers.@declare_executable_product(bluetooth)
JLLWrappers.@declare_executable_product(bluetoothctl)
JLLWrappers.@declare_executable_product(bluetoothd)
JLLWrappers.@declare_executable_product(btattach)
JLLWrappers.@declare_executable_product(btmon)
JLLWrappers.@declare_executable_product(hex2hcd)
JLLWrappers.@declare_executable_product(hid2hci)
JLLWrappers.@declare_executable_product(l2ping)
JLLWrappers.@declare_executable_product(l2test)
JLLWrappers.@declare_library_product(libbluetooth, "libbluetooth.so.3")
JLLWrappers.@declare_executable_product(mpris_proxy)
JLLWrappers.@declare_executable_product(obexd)
JLLWrappers.@declare_executable_product(rctest)
function __init__()
    JLLWrappers.@generate_init_header(Dbus_jll, eudev_jll, Glib_jll, Libical_jll, Readline_jll)
    JLLWrappers.@init_executable_product(
        bccmd,
        "bin/bccmd",
    )

    JLLWrappers.@init_executable_product(
        bluemoon,
        "bin/bluemoon",
    )

    JLLWrappers.@init_executable_product(
        bluetooth,
        "lib/cups/backend/bluetooth",
    )

    JLLWrappers.@init_executable_product(
        bluetoothctl,
        "bin/bluetoothctl",
    )

    JLLWrappers.@init_executable_product(
        bluetoothd,
        "libexec/bluetooth/bluetoothd",
    )

    JLLWrappers.@init_executable_product(
        btattach,
        "bin/btattach",
    )

    JLLWrappers.@init_executable_product(
        btmon,
        "bin/btmon",
    )

    JLLWrappers.@init_executable_product(
        hex2hcd,
        "bin/hex2hcd",
    )

    JLLWrappers.@init_executable_product(
        hid2hci,
        "lib/udev/hid2hci",
    )

    JLLWrappers.@init_executable_product(
        l2ping,
        "bin/l2ping",
    )

    JLLWrappers.@init_executable_product(
        l2test,
        "bin/l2test",
    )

    JLLWrappers.@init_library_product(
        libbluetooth,
        "lib/libbluetooth.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mpris_proxy,
        "bin/mpris-proxy",
    )

    JLLWrappers.@init_executable_product(
        obexd,
        "libexec/bluetooth/obexd",
    )

    JLLWrappers.@init_executable_product(
        rctest,
        "bin/rctest",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# it seems that currently, map operation over Nullable array is not optimzed in julia 0.4.1.

import Base: .+, .-, .*, *, ./, .\, .//, .==, .<, .!=, .<=, .%, .<<, .>>, .^, +, -, ~, &, |, $, ==, !=
import DataFrames: DataFrame

"""

`wrap_array(arr)`

Wrap an array by `AbstractArrayWrapper` if it is not `DictArray` or `labeledArray`, and not already `AbstractArrayWrapper`.

"""
wrap_array(arr::AbstractArrayWrapper) = arr
wrap_array(arr::LabeledArray) = arr
wrap_array(arr::DictArray) = arr
wrap_array(arr::AbstractArray) = AbstractArrayWrapper(arr)
wrap_array(arr::DataFrame) = arr

Base.setindex!{T,N}(arr::AbstractArrayWrapper{T,N}, v::T, arg::Int) = setindex!(arr.a, v, arg)
Base.setindex!{T,N}(arr::AbstractArrayWrapper{T,N}, v::T, args::Int...) = setindex!(arr.a, v, args...)
Base.setindex!{T,N}(arr::AbstractArrayWrapper{Nullable{T},N}, v::T, args::Int...) = setindex!(arr.a, v, args...)
Base.eltype{T,N,A}(::Type{AbstractArrayWrapper{T,N,A}}) = T
Base.linearindexing{T,N,A}(::Type{AbstractArrayWrapper{T,N,A}}) = Base.linearindexing(A)
Base.sub(arr::AbstractArrayWrapper, args::Union{Colon,Int,AbstractVector}...) = AbstractArrayWrapper(sub(arr.a, args...))
Base.slice(arr::AbstractArrayWrapper, args::Union{Colon,Int,AbstractVector}...) = AbstractArrayWrapper(slice(arr.a, args...))
Base.sub(arr::AbstractArrayWrapper, args::Tuple{Vararg{Union{Colon,Int,AbstractVector}}})= AbstractArrayWrapper(sub(arr.a, args...))
Base.slice(arr::AbstractArrayWrapper, args::Tuple{Vararg{Union{Colon,Int,AbstractVector}}}) = AbstractArrayWrapper(slice(arr.a, args...))
@delegate(AbstractArrayWrapper.a, Base.start, Base.next, Base.done, Base.size,
                           Base.ndims, Base.length, Base.setindex!, Base.find)
@delegate_and_lift(AbstractArrayWrapper.a, Base.transpose, Base.permutedims, Base.repeat,
                                   Base.repeat, Base.transpose, Base.permutedims,
                                   Base.sort, Base.sort!, Base.sortperm, Base.reverse,
                                   Base.sub, Base.slice)
Base.similar{T,N}(arr::AbstractArrayWrapper, ::Type{T}, dims::NTuple{N,Int}) = AbstractArrayWrapper(similar(arr.a, T, dims))
Base.similar{T<:AbstractFloat,U<:AbstractFloat,N,A,M}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                                                    ::Type{Nullable{U}},
                                                    dims::NTuple{M,Int}) =
  AbstractArrayWrapper(FloatNAArray(similar(arr.a.data, U, dims)))
Base.similar{U<:AbstractFloat,M}(arr::AbstractArrayWrapper,
                                 ::Type{Nullable{U}},
                                 dims::NTuple{M,Int}) =
  AbstractArrayWrapper(FloatNAArray(similar(arr.a, U, dims)))
Base.repeat(arr::AbstractArrayWrapper; kwargs...) = AbstractArrayWrapper(repeat(arr.a; kwargs...))
Base.sort(arr::AbstractArrayWrapper; kwargs...) = AbstractArrayWrapper(sort(arr.a; kwargs...))
Base.sort!(arr::AbstractArrayWrapper; kwargs...) = AbstractArrayWrapper(sort!(arr.a; kwargs...))
Base.getindex{T,N}(arr::AbstractArrayWrapper{T,N}, arg::Int) = getindex(arr.a, arg)
Base.getindex{T,N}(arr::AbstractArrayWrapper{T,N}, args::Int...) = getindex(arr.a, args...)
Base.getindex{T,N}(arr::AbstractArrayWrapper{T,N}, indices::CartesianIndex) = getindex(arr.a, indices)
Base.getindex(arr::AbstractArrayWrapper, args...) = begin
  res = getindex(arr.a, args...)
  #if is_scalar_indexing(args)
  #  res
  #else
  AbstractArrayWrapper(res)
  #end
end
Base.map(f, arr::AbstractArrayWrapper) = AbstractArrayWrapper(map(f, arr.a))
Base.map(f, arrs::AbstractArrayWrapper...) = AbstractArrayWrapper(map(f, map(x->x.a, arrs)...))

macro absarray_unary_wrapper(ops...)
  targetexpr = map(ops) do op
    quote
      $(esc(op.args[1]))(arr::DictArray) = mapvalues(x->$(esc(op.args[1]))(x), arr)
      $(esc(op.args[1]))(arr::LabeledArray) = LabeledArray($(esc(op.args[1]))(peel(arr)), pickaxis(arr))
      $(esc(op.args[1]))(x::AbstractArrayWrapper) = begin
        # AbstractArrayWrapper(map($(esc(op.args[2])), x.a))
        result = similar(x)
        $(symbol(op.args[1],naop_suffix))(result.a, x)
        result
      end
      $(symbol(op.args[1],naop_suffix))(result, x::AbstractArrayWrapper) = begin
        xa = x.a
        for i in eachindex(xa)
          @inbounds result[i] = $(esc(op.args[2]))(xa[i])
        end
      end
      $(symbol(op.args[1],naop_suffix)){T<:AbstractFloat,N,A}(result, x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = begin
        xadata = x.a.data
        resultdata = result.data
        for i in eachindex(xadata)
          @inbounds resultdata[i] = $(esc(op.args[2]))(xadata[i])
        end
      end
    end
  end
  Expr(:block, targetexpr...)
end

# Ideally, lift every possible types using some supertypes. However, a lot of annoying ambiguity warnings may occur.
# So, try to fiddle around with possible combinations that do not give any ambiguity warnings.
const LiftToNullableTypes = [Bool,
                             Integer,
                             AbstractFloat,
                             Rational,
                             Irrational{:e},
                             Real,
                             Complex,
                             AbstractString,
                             Char,
                             Symbol]

promote_nullable_types{T,U}(::Type{Nullable{T}},::Type{Nullable{U}}) = Nullable{promote_type(T,U)}
promote_nullable_types{T,U}(::Type{T},::Type{Nullable{U}}) = Nullable{promote_type(T,U)}
promote_nullable_types{T,U}(::Type{Nullable{T}},::Type{U}) = Nullable{promote_type(T,U)}
promote_nullable_types{T,U}(::Type{T},::Type{U}) = promote_type(T,U)
#promote_nullable_types(::DataType,::DataType) = Nullable{Any}

preset_nullable_type{T,U}(::Type{Nullable{T}},::Type{Nullable{U}}, tpe) = Nullable{tpe}
preset_nullable_type{T,U}(::Type{T},::Type{Nullable{U}}, tpe) = Nullable{tpe}
preset_nullable_type{T,U}(::Type{Nullable{T}},::Type{U}, tpe) = Nullable{tpe}
preset_nullable_type{T,U}(::Type{T},::Type{U}, tpe) = tpe

const naop_suffix = "!"

# some adhoc definitions to suppress ambiguity warnings clashing with irrationals.jl
.^(x::Base.Irrational{:e}, y::AbstractArrayWrapper{Real}) = AbstractArrayWrapper(.^(x, y.a))
.^{T<:Real}(x::Base.Irrational{:e}, y::AbstractArrayWrapper{T}) = AbstractArrayWrapper(.^(x, y.a))

# not used anymore.
macro absarray_binary_wrapper(ops...)
  targetexpr = map(ops) do op
    nullelem = if length(op.args) == 2
      #Expr(:curly, :Nullable, Expr(:call, :promote_type, :T, :U))
      :(promote_nullable_types(T,U))
    elseif length(op.args) == 3
      #Expr(:curly, :Nullable, op.args[3])
      :(preset_nullable_type(T,U,$(op.args[3])))
    end
    quote
      $(esc(op.args[1])){T,U}(x::AbstractArrayWrapper{T}, y::AbstractArrayWrapper{U}) = begin
        @assert(size(x) == size(y))
        #AbstractArrayWrapper(map((u,v)->$(esc(op.args[2]))(u,v), x.a, y.a))
        result = similar(x, $nullelem)
        $(symbol(op.args[1],naop_suffix))(result.a, x, y)
        result
      end
      $(symbol(op.args[1],naop_suffix)){T,U}(result, x::AbstractArrayWrapper{T}, y::AbstractArrayWrapper{U}) = begin
        xa = x.a
        ya = y.a
        for i in eachindex(xa,ya)
          @inbounds result[i] = $(esc(op.args[2]))(xa[i],ya[i])
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T,A,U,B}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                                                     y::AbstractArrayWrapper{Nullable{U},N,FloatNAArray{U,N,B}}) = begin
        xadata = x.a.data
        yadata = y.a.data
        resultdata = result.data
        for i in eachindex(xadata,yadata)
          @inbounds resultdata[i] = $(esc(op.args[2]))(xadata[i],yadata[i])
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T,A,U<:Nullable}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xadata = x.a.data
        ya= y.a
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xadata,ya)
          @inbounds yai = ya[i]
          if yai.isnull
            @inbounds resultdata[i] = na
          else
            @inbounds resultdata[i] = $(esc(op.args[2]))(xadata[i],yai.value)
          end
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T<:Nullable,U,B}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{Nullable{U},N,FloatNAArray{U,N,B}}) = begin
        xa= x.a
        yadata = y.a.data
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xa,yadata)
          @inbounds xai = xa[i]
          if xai.isnull
            @inbounds resultdata[i] = na
          else
            @inbounds resultdata[i] = $(esc(op.args[2]))(xai.value, yadata[i])
          end
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T<:Nullable,U<:Nullable}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xa= x.a
        ya = y.a
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xa,ya)
          @inbounds xai = xa[i]
          @inbounds yai = ya[i]
          if xai.isnull || yai.isnull
            @inbounds resultdata[i] = na
          else
            @inbounds resultdata[i] = $(esc(op.args[2]))(xai.value, yai.value)
          end
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T,A,U}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xadata = x.a.data
        ya= y.a
        resultdata = result.data
        for i in eachindex(xadata,ya)
          @inbounds resultdata[i] = $(esc(op.args[2]))(xadata[i],ya[i])
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T,U,B}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{Nullable{U},N,FloatNAArray{U,N,B}}) = begin
        xa= x.a
        yadata = y.a.data
        resultdata = result.data
        for i in eachindex(xa,yadata)
          @inbounds resultdata[i] = $(esc(op.args[2]))(xa[i], yadata[i])
        end
      end

      $(symbol(op.args[1],naop_suffix)){N,K,T<:Nullable,U}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xa= x.a
        ya = y.a
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xa,ya)
          @inbounds v = $(esc(op.args[2]))(xa[i], ya[i])
          @inbounds resultdata[i] = v.isnull ? na : v.value
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T,U<:Nullable}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xa= x.a
        ya = y.a
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xa,ya)
          @inbounds v = $(esc(op.args[2]))(xa[i], ya[i])
          @inbounds resultdata[i] = v.isnull ? na : v.value
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T<:Nullable,U<:Nullable}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xa= x.a
        ya = y.a
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xa,ya)
          @inbounds v = $(esc(op.args[2]))(xa[i], ya[i])
          @inbounds resultdata[i] = v.isnull ? na : v.value
        end
      end
      $(symbol(op.args[1],naop_suffix)){N,K,T,U}(result::FloatNAArray{K,N},
                                                     x::AbstractArrayWrapper{T,N},
                                                     y::AbstractArrayWrapper{U,N}) = begin
        xa= x.a
        ya = y.a
        resultdata = result.data
        na = convert(K,NaN)
        for i in eachindex(xa,ya)
          @inbounds resultdata[i] = $(esc(op.args[2]))(xa[i], ya[i])
        end
      end

      $(esc(op.args[1])){T,U<:Nullable}(x::AbstractArrayWrapper{T}, y::U) = begin
        result = similar(x, $nullelem)
        $(symbol(op.args[1],naop_suffix))(result.a, x, y)
        result
      end
      $(symbol(op.args[1],naop_suffix)){T,U}(result, x::AbstractArrayWrapper{T}, y::U) = begin
        xa = x.a
        for i in eachindex(xa)
          @inbounds result[i] = $(esc(op.args[2]))(xa[i],y)
        end
      end
      $(symbol(op.args[1],naop_suffix)){K,T,U<:Nullable,N,A}(result::FloatNAArray{K,N}, x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, y::U) = begin
        if y.isnull
          setna!(result)
        else
          xadata = x.a.data
          yvalue = y.value
          resultdata = result.data
          for i in eachindex(xadata)
            @inbounds resultdata[i] = $(esc(op.args[2]))(xadata[i],yvalue)
          end
        end
      end

      $(esc(op.args[1])){T<:Nullable,U}(x::T, y::AbstractArrayWrapper{U}) = begin
        #AbstractArrayWrapper(map(v->$(esc(op.args[2]))(x,v), y.a))
        result = similar(y, $nullelem)
        $(symbol(op.args[1],naop_suffix))(result.a, x, y)
        result
      end
      $(symbol(op.args[1],naop_suffix)){T,U}(result, x::T, y::AbstractArrayWrapper{U}) = begin
        ya = y.a
        for i in eachindex(ya)
          @inbounds result[i] = $(esc(op.args[2]))(x,ya[i])
        end
      end
      $(symbol(op.args[1],naop_suffix)){K,T<:Nullable,U,N,A}(result::FloatNAArray{K,N}, x::T, y::AbstractArrayWrapper{Nullable{U},N,FloatNAArray{U,N,A}}) = begin
        if x.isnull
          setna!(result)
        else
          xvalue = x.value
          yadata = y.a.data
          resultdata = result.data
          for i in eachindex(yadata)
            @inbounds resultdata[i] = $(esc(op.args[2]))(xvalue, yadata[i])
          end
        end
      end

      $(symbol(op.args[1],"nulltype2", naop_suffix))(result, x::AbstractArrayWrapper, y) = begin
        xa = x.a
        for i in eachindex(xa)
          @inbounds result[i] = $(esc(op.args[2]))(xa[i],y)
        end
      end
      $(symbol(op.args[1],"nulltype2", naop_suffix)){K,T,N,A}(result::FloatNAArray{K,N}, x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, y) = begin
        xadata = x.a.data
        resultdata = result.data
        for i in eachindex(xadata)
          @inbounds resultdata[i] = $(esc(op.args[2]))(xadata[i],y)
        end
      end

      $(symbol(op.args[1],"nulltype1", naop_suffix))(result, x, y::AbstractArrayWrapper) = begin
        ya = y.a
        for i in eachindex(ya)
          @inbounds result[i] = $(esc(op.args[2]))(x,ya[i])
        end
      end
      $(symbol(op.args[1],"nulltype1", naop_suffix)){K,N,T,A}(result::FloatNAArray{K,N}, x, y::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = begin
        yadata = y.a.data
        resultdata = result.data
        for i in eachindex(yadata)
          @inbounds resultdata[i] = $(esc(op.args[2]))(x,yadata[i])
        end
      end

      for nulltype in $LiftToNullableTypes
        $(esc(op.args[1]))(x::AbstractArrayWrapper{nulltype}, y::nulltype) = begin
          #AbstractArrayWrapper(map(u->$(esc(op.args[2]))(u,y), x.a))
          T = eltype(x)
          U = typeof(y)
          result = similar(x, $nullelem)
          $(symbol(op.args[1],"nulltype2",naop_suffix))(result.a, x, y)
          result
        end

        $(esc(op.args[1]))(x::nulltype, y::AbstractArrayWrapper{nulltype}) = begin
          #AbstractArrayWrapper(map(v->$(esc(op.args[2]))(x,v), y.a))
          T = typeof(x)
          U = eltype(y)
          result = similar(y, $nullelem)
          $(symbol(op.args[1],"nulltype1",naop_suffix))(result.a, x, y)
          result
        end

        $(esc(op.args[1])){T<:nulltype}(x::AbstractArrayWrapper{T}, y::nulltype) = begin
          #AbstractArrayWrapper(map(u->$(esc(op.args[2]))(u,y), x.a))
          U = typeof(y)
          result = similar(x, $nullelem)
          $(symbol(op.args[1],"nulltype2",naop_suffix))(result.a, x, y)
          result
        end

        $(esc(op.args[1])){U<:nulltype}(x::nulltype, y::AbstractArrayWrapper{U}) = begin
          #AbstractArrayWrapper(map(v->$(esc(op.args[2]))(x,v), y.a))
          T = typeof(x)
          result = similar(y, $nullelem)
          $(symbol(op.args[1],"nulltype1",naop_suffix))(result.a, x, y)
          result
        end

        $(esc(op.args[1])){T<:Nullable}(x::AbstractArrayWrapper{T}, y::nulltype) = begin
          #AbstractArrayWrapper(map(u->$(esc(op.args[2]))(u,y), x.a))
          U = typeof(y)
          result = similar(x, $nullelem)
          $(symbol(op.args[1],"nulltype2",naop_suffix))(result.a, x, y)
          result
        end

        $(esc(op.args[1])){U<:Nullable}(x::nulltype, y::AbstractArrayWrapper{U}) = begin
          #AbstractArrayWrapper(map(v->$(esc(op.args[2]))(x,v), y.a))
          T = typeof(x)
          result = similar(y, $nullelem)
          $(symbol(op.args[1],"nulltype1",naop_suffix))(result.a, x, y)
          result
        end
      end
    end
  end
  Expr(:block, targetexpr...)
end

*(x::Real, y::AbstractArrayWrapper) = x .* y
*(x::AbstractArrayWrapper, y::Real) = x .* y
*(x::Complex, y::AbstractArrayWrapper) = x .* y
*(x::AbstractArrayWrapper, y::Complex) = x .* y
*{T<:Real}(x::Nullable{T}, y::AbstractArrayWrapper) = x .* y
*{T<:Real}(x::AbstractArrayWrapper, y::Nullable{T}) = x .* y
*{T<:Complex}(x::Nullable{T}, y::AbstractArrayWrapper) = x .* y
*{T<:Complex}(x::AbstractArrayWrapper, y::Nullable{T}) = x .* y

macro nullable_unary_wrapper(ops...)
  targetexpr = map(ops) do op
    nullelem = if length(op.args) == 2
      Expr(:curly, :Nullable, :T)
    elseif length(op.args) == 3
      Expr(:curly, :Nullable, op.args[3])
    end
    quote
      $(esc(op.args[2])){T}(x::Nullable{T}) = x.isnull ? $nullelem() : Nullable($(esc(op.args[1]))(x.value))
      $(esc(op.args[2])){T}(x::T) = $(esc(op.args[1]))(x)
    end
  end
  Expr(:block, targetexpr...)
end

macro nullable_binary_wrapper(ops...)
  targetexpr = map(ops) do op
    nullelem = if length(op.args) == 2
      Expr(:curly, :Nullable, Expr(:call, :promote_type, :T, :V))
    elseif length(op.args) == 3
      Expr(:curly, :Nullable, op.args[3])
    end
    quote
      $(esc(op.args[2])){T,V}(x::Nullable{T}, y::Nullable{V}) =
        x.isnull || y.isnull ? $nullelem() : $nullelem($(esc(op.args[1]))(x.value, y.value))
      $(esc(op.args[2])){T,V}(x::Nullable{T}, y::V) =
        x.isnull ? $nullelem() : $nullelem($(esc(op.args[1]))(x.value, y))
      $(esc(op.args[2])){T,V}(x::T, y::Nullable{V}) =
        y.isnull ? $nullelem() : $nullelem($(esc(op.args[1]))(x, y.value))
      $(esc(op.args[2])){T,V}(x::T, y::V) = $(esc(op.args[1]))(x, y)
    end
  end
  Expr(:block, targetexpr...)
end

@nullable_unary_wrapper((+, naop_plus), (-, naop_minus), (~, naop_not))
@nullable_binary_wrapper((.+, naop_plus), (.-, naop_minus), (.*, naop_mul),
                         (./, naop_div),
                         # could not use .== or .!= over general types such as AbstractString.
                         # let's settle down to using == and != instead for now.
                         # at least, we don't have to provide a blanket definition (.==)(x,y) = x==y then.
                         (.\, naop_invdiv), (.//, naop_frac), (==, naop_eq, Bool), (.<, naop_lt, Bool),
                         (!=, naop_noeq, Bool), (.<=, naop_le, Bool), (.%, naop_mod), (.<<, naop_lsft),
                         (.>>, naop_rsft), (.^, naop_exp),
                         (&, naop_and), (|, naop_or), ($, naop_xor))

@absarray_unary_wrapper((+, naop_plus), (-, naop_minus), (.+, naop_plus), (.-, naop_minus), (~, naop_not))
@absarray_binary_wrapper((+, naop_plus), (-, naop_minus), (.+, naop_plus), (.-, naop_minus), (.*, naop_mul),
                         (./, naop_div),
                         (.\, naop_invdiv), (.//, naop_frac), (.==, naop_eq, Bool), (.<, naop_lt, Bool),
                         (.!=, naop_noeq, Bool), (.<=, naop_le, Bool), (.%, naop_mod), (.<<, naop_lsft),
                         (.>>, naop_rsft), (.^, naop_exp),
                         (&, naop_and), (|, naop_or), ($, naop_xor))

(==){T<:Nullable,U<:Nullable}(x::AbstractArrayWrapper{T}, y::AbstractArrayWrapper{U}) = begin
  if x === y
    return true
  else
    for (elx, ely) in zip(x.a, y.a)
      if elx.isnull && !ely.isnull
        return false
      elseif !elx.isnull && ely.isnull
        return false
      elseif !elx.isnull && !ely.isnull && elx.value != ely.value
        return false
      end
    end
    return true
  end
end

(==){T<:AbstractFloat,U<:AbstractFloat,N,A,B}(x::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, y::AbstractArrayWrapper{Nullable{U},N,FloatNAArray{U,N,B}}) = begin
  if x === y
    return true
  else
    for (elx, ely) in zip(x.a.data, y.a.data)
      if !(isnan(elx) && isnan(ely)) && elx != ely
        return false
      end
    end
    return true
  end
end

# TODO make sure this blanket definition is okay.
# remvoed in favor of using == instead of .== for naop_eq (and similarly for naop_noeq).
# (.==)(x, y) = x == y


text:
@inline function add_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:AbstractFloat}
    hi, lo = x
    hi, lo = add_2(hi, y, lo)
    return hi, lo
end

@inline function sub_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:AbstractFloat}
    hi, lo = x
    hi, lo = add_2(hi, lo, -y)
    return hi, lo
end

@inline function mul_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:AbstractFloat}
    hi, lo = x
    hihi, hilo = mul_2(y, hi)
    lohi, lolo = mul_2(y, lo)
    hi, lo = add_2(hihi, hilo, lohi, lolo)
    return hi, lo
end

@inline function dvi_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:AbstractFloat}
    yinv = inv_fp_dd(y)
    zhi, zlo = mul_dddd_dd(x, yinv)
    return zhi, zlo
end


text:
type LexicographicFitnessScheme{T} <: FitnessScheme
  randomised::Bool
  preferences::Vector{Integer}

  LexicographicFitnessScheme(r::Bool, o::Vector{Integer}) = new(r, o)
end

"""
TODO: Add description of lexicographic fitness scheme.

TODO: Could make this lazy?

**Properties:**

* `randomised::Bool`, a flag indicating whether the preference for specific
objectives should be randomised at each comparison.
* `preferences::Vector{Int}`, a list, specifying the order in which each of
objectives should be compared, given as a list of objective indices
(one-indexed).
"""
lexicographic(s::Dict{Any, Any}) =
  lexicographic(Base.get(s, "randomised", false), s["preferences"])

lexicographic(randomised::Bool, preferences::Vector{Int}) =
  LexicographicFitnessScheme(randomised, preferences)

lexicographic(preferences::Vector{Int}) =
  lexicographic(false, preferences)

uses{T}(s::LexicographicFitnessScheme{T}) = Vector{T}

assign{T}(s::LexicographicFitnessScheme{T}, v::Vector{T}) = v

function compare{T}(s::LexicographicFitnessScheme, x::Vector{T}, y::Vector{T})
  for i in (s.randomised ? shuffle(s.objectives) : s.objectives)
    if x[i] != y[i]
      return y[i] > x[i] ? 1 : -1
    end
  end
  return 0
end


text:
const FermionSite = TagType"Fermion"

function siteinds(::FermionSite,
                  N::Int; kwargs...)
  conserve_qns = get(kwargs,:conserve_qns,false)
  conserve_nf = get(kwargs,:conserve_nf,conserve_qns)
  conserve_parity = get(kwargs,:conserve_parity,conserve_qns)
  if conserve_nf
    zer = QN("Nf",0,-1) => 1
    one = QN("Nf",1,-1) => 1
    return [Index(zer,one;tags="Site,Fermion,n=$n") for n=1:N]
  elseif conserve_parity
    zer = QN("Pf",0,-2) => 1
    one = QN("Pf",1,-2) => 1
    return [Index(zer,one;tags="Site,Fermion,n=$n") for n=1:N]
  end
  return [Index(2,"Site,Fermion,n=$n") for n=1:N]
end

function state(::FermionSite,
               st::AbstractString)
  if st == "Emp" || st == "0"
    return 1
  elseif st == "Occ" || st == "1"
    return 2
  end
  throw(ArgumentError("State string \"$st\" not recognized for Fermion site"))
  return 0
end

function op(::FermionSite,
            s::Index,
            opname::AbstractString)::ITensor
  sP = prime(s)
  Emp   = s(1)
  EmpP  = sP(1)
  Occ   = s(2)
  OccP  = sP(2)

  Op = ITensor(s',dag(s))

  if opname == "N"
    Op[OccP, Occ] = 1.
  elseif opname == "C"
    Op[EmpP, Occ] = 1.
  elseif opname == "Cdag"
    Op[OccP, Emp] = 1.
  elseif opname=="F" || opname=="FermiPhase" || opname=="FP"
    Op[EmpP,Emp] =  1.
    Op[OccP,Occ] = -1.
  elseif opname == "Emp" || opname == "0"
    pEmp = ITensor(s)
    pEmp[Emp] = 1.0
    return pEmp
  elseif opname == "Occ" || opname == "1"
    pOcc = ITensor(s)
    pOcc[Occ] = 1.0
    return pOcc
  else
    throw(ArgumentError("Operator name $opname not recognized for FermionSite"))
  end
  return Op
end

function has_fermion_string(::FermionSite,
            s::Index,
            opname::AbstractString)::Bool
  if opname=="C" || opname=="Cdag"
    return true
  end
  return false
end


text:
Pkg.clone("https://github.com/dfdx/Spark.jl")
Pkg.build("Spark")
# we also need latest master of JavaCall.jl
Pkg.checkout("JavaCall")


text:
############################################################
# joAbstractParallelableLinearOperator subtypes ############
############################################################


text:
@testset "Check constructors" begin
    float_types_to_test = [Float16, Float32, Float64, ComplexF32, ComplexF64]
    int_types_to_test = [Int16, Int32, Int64]

    function test_arr(arr, f, T, s)
        @test eltype(arr) == T
        @test last(arr) ≈ sqrt(sum(abs2.((s))))
        @test arr.size == s
        @test first(arr) ≈ sqrt(sum(abs2.((ntuple(x -> 1, length(s))))))
        @test typeof(arr) <: IndexFunArray{T, length(s), F} where F
    end

    @testset "Check IndexFunArray initializer for types and boundaries" begin

        sizes = [(10, 10), (2,2,2), (3,3,3,3)]
        for s in sizes
            for T in float_types_to_test
                f(ind) = T(sqrt(sum(abs2.(ind))))
                arr = IndexFunArray(f, s)
                test_arr(arr, f, T, s)
                arr = IndexFunArray(T, f, s)
                test_arr(arr, f, T, s)
            end
        end
    end

#    @testset "Certain special cases" begin
#        @test IndexFunArray(x -> x, (1,0), (1,1), (3, 3)) == [(0, 1) (0, 2) (0, 3); (1, 1) (1, 2) (1, 3); (2, 1) (2, 2) (2, 3)]
#
#    end
end


text:
## Test
using Base.Test
using Base.Iterators
using SymEngine
using SymEngineExtras
PM = SymEngineExtras

@vars a b c x y z x_ x__ x___ y_ y__ y___ z_ z__ z___

@testset "syntactic match (non assoc, non comm)" begin
    d = Dict()
    success, d = PM.syntactic_match(a, x_, d)
    @test success

    d = Dict()
    success, d = PM.syntactic_match(sin(a), sin(x_), d)
    @test success

    d = Dict()
    success, d = PM.syntactic_match(a^b, x_^b, d)
    @test success

    d = Dict()
    success, d = PM.syntactic_match(b^a, b^x_, d)
    @test success

    d= Dict()
    success, d = PM.syntactic_match(a^a, x_^x_, d)
    @test success

    d= Dict()
    success, d = PM.syntactic_match(a^b, x_^x_, d)
    @test !success

    d= Dict(x_=>b)
    success, d = PM.syntactic_match(a^a, x_^x_, d)
    @test !success

end

@testset "match one to one" begin
    ## possibly assoc, poss comm helper
    Theta = Set((PM.aDict(),))
    theta = PM.match_one_to_one([a], x_, Theta)
    @test length(collect(theta)) == 1 && collect(theta)[1][x_] == a

    theta = PM.match_one_to_one([sin(a)], sin(x_), Theta)
    @test length(collect(theta)) == 1 && collect(theta)[1][x_] == a

    subject = [x, 2x, 3x]
    theta = PM.match_one_to_one(subject, x_, Theta, true, :Add)
    @test length(theta) == 1
    @test collect(theta)[1] == Dict(x_=>6x)

    theta = PM.match_one_to_one(subject, x_, Theta)
    @test length(theta) == 0

    theta = PM.match_one_to_one(subject, x__, Theta)
    @test length(theta) == 1 # x__ -> pattern

    theta = PM.match_one_to_one([x,x,x],x_, Theta)
    @test isempty(theta)

    theta = PM.match_one_to_one([x],x_, Theta)
    @test !isempty(theta)

end

@testset "match_sequence" begin
    ## poss assoc, not commutative
    sigma = Dict()
    subject = [a,b,c,a,b,c]
    pattern = [x__,y___,a,b,x_]
    theta = PM.match_sequence(subject, pattern, sigma, Val{false})
    d = Dict(x_=>c,x__=>[a,b,c], y___=>Basic[])
    @test any(d == d1 for d1 in theta)

    subject = [a]
    theta = PM.match_sequence(subject, [x_], sigma, Val{false})
    @test length(collect(theta)) == 1

    subject = [a,b,c]
    theta = PM.match_sequence(subject, [x_,b,z_], sigma, Val{false})
    @test length(collect(theta)) == 1

    subject = [a,b]
    theta = PM.match_sequence(subject, [x_,y_], sigma, Val{false})
    @test length(collect(theta)) == 1

    subject = [a,b,c,a,b,c]
    pattern = [x___,x_, y_, y___]
    theta = PM.match_sequence(subject, pattern, sigma, Val{false})
    d = Dict(x___=>[a], x_=>b, y_=>c,y___=>[a,b,c])
     @test any(d == d1 for d1 in theta)

    subject = Basic[1,2,3,4,3,2,1]
    pattern = [x___,x_, y_, y___]
    theta = PM.match_sequence(subject, pattern, sigma, Val{false})
    theta = Iterators.filter(m -> m[x_] < m[y_], theta)
    for m in theta
        @test m[x_] < m[y_]
    end

    subject = [a,b,c]
    pattern = [a,b,c]
    theta = PM.match_sequence(subject, pattern, sigma, Val{false})
    @test length(collect(theta)) == 1 && PM.aDict() == collect(theta)[1]  # one match, an empty substitution

    subject = Basic[a,b,c,b,a]
    pattern = [x__,c,x__] # no match as [a,b] != [b,a]
    theta = PM.match_sequence(subject, pattern, sigma, Val{false})
    @test isempty(theta)

    subject = Basic[a,b,c,b,a]
    pattern = [x__,b,y___]
    theta = PM.match_sequence(subject, pattern, sigma, Val{false})
    @test length(collect(theta)) == 2

end

@testset "comm=true" begin
    # constants
    sigma = Dict()
    out = PM.match_sequence([a,b,c], [a,b,c], sigma, Val{true})
    # out should be {\phi} (Set((Dict(),))
    @test collect(out)[1] == PM.aDict()

    # non-variable
    ss, ps = [sin(a),cos(b),c], [sin(x_),cos(y_),c]
    out = PM.match_sequence(ss, ps, sigma, Val{true})
    @test all([PM.so_far_so_good(ss, ps, sigma) for sigma in out])

    ss, ps = [sin(a),cos(b),c], [sin(x_),cos(y_),c]
    out = PM.match_sequence([sin(a),cos(b),c], [sin(x_),cos(y_),c], Dict(x_=>b), Val{true})
    @test isempty(out)  # x_->a but x_->b...

    ## regular variable
    ss, ps = [a,a,a,a,b,b,c], [x__,x__,y__,y__,z_]
    out = PM.match_sequence(ss, ps, sigma, Val{true})
    @test all([PM.so_far_so_good(ss, ps, sigma) for sigma in out])
    @test length(collect(out)) == 4 # {z->c, x=[a,b],y=a}, {z->c, x->[aa], y->c}...

    ss, ps = [sin(a),cos(b),a, a,a,b,b,c,cos(a^2)], [cos(y_), sin(x_), x_,x_,x_,y_,y_,z_, cos(a^2)]
    out = PM.match_sequence(ss, ps, sigma, Val{true})
    @test all([PM.so_far_so_good(ss, ps, sigma) for sigma in out])
    @test length(collect(out)) == 1
end

@testset "matches" begin

    ## {} isno match, {Dict()} is a match with no wildcard
    out = matches(a, b)
    @test isempty(out)

    out = matches(a,a)
    @test !isempty(out) && isempty(collect(out)[1])

    ## various matches for x
    out = matches(x_, a)
    @test  !isempty(out) && collect(out)[1] == Dict(x_ => a)

    out = matches(sin(x_), sin(a))
    @test  !isempty(out) && collect(out)[1] == Dict(x_ => a)

    out = matches([sin(x_)], [sin(a)])
    @test  !isempty(out) && collect(out)[1] == Dict(x_ => a)

    ex = x^4 * exp(x^2)
    ex1 = diff(ex, x,x,x,x,x,x,x,x,x,x)
    out = matches(ex1(x=> x_), ex1)
    @test !isempty(out) && collect(out)[1] == Dict(x_ => x)

    ## commutative/not
    out = matches([sin(x_), cos(y_)], [cos(b), sin(a)], commutative=false)
    @test isempty(out)

    out = matches([sin(x_), cos(y_)], [cos(b), sin(a)], commutative=true)
    @test !isempty(out)

    ## Associative
    out = matches([x_, c], [a,b,c], associative=false)
    @test isempty(out)

    out = matches([x_, c], [a,b,c], associative=true)
    @test !isempty(out) && length(collect(out)) == 1

    out = matches([x_, x_, y_], [a,a,a,b,b,c], associative=false)
    @test isempty(out)

    out = matches([x_, x_, y_], [a,a,a,b,b,c], associative=true)
    @test !isempty(out)

    ## associative -. x_ a plus variable
    out = matches([x_], [a,b,c], associative=false)
    @test  isempty(out)
    out = matches([x_], [a,b,c], associative=true)
    @test  !isempty(out)

    out = matches(x__ * y_,  a * b * c)
    @test  !isempty(out) && length(collect(out)) == 6 # ch(3,2) + ch(3,1) =
    out = matches(x__ + y_,  a + b + c)
    @test  !isempty(out) && length(collect(out)) == 6 # ch(3,2) + ch(3,1) =

    out = matches([x__ * y_,y_], [a, a * b * c], commutative=false)
    @test  isempty(out)
    out = matches([x__ * y_,y_], [a, a * b * c], commutative=true)
    @test  !isempty(out) && length(collect(out)) == 1 #

end


text:
using CANalyze.Decode
using CANalyze.Frames
using CANalyze.Signals
using CANalyze.Messages

signal1 = Signals.NamedSignal("ABC", nothing, nothing, Signals.Float32Signal(start=0, byte_order=:little_endian))
signal2 = Signals.NamedSignal("ABCD", nothing, nothing, Signals.Unsigned(start=40, length=17, factor=2, offset=20, byte_order=:big_endian))
signal3 = Signals.NamedSignal("ABCDE", nothing, nothing, Signals.Unsigned(start=32, length=8, factor=2, offset=20, byte_order=:little_endian))

bits1 = Signals.Bits(signal1)
println(sort(Int64[bits1.bits...]))

bits1 = Signals.Bits(signal2)
println(sort(Int64[bits1.bits...]))

bits1 = Signals.Bits(signal3)
println(sort(Int64[bits1.bits...]))

frame = Frames.CANFrame(20, [1, 2, 0xFD, 4, 5, 6, 7, 8])
m = Messages.Message(0x1FF, 8, "ABC", signal1, signal2, signal3; strict=true)

d = Decode.decode(m, frame)
println(d)


text:
using LightGraphs

function run_ttc_iteration(proposing_group::AbstractVector{Agent}, other_group::AbstractVector{Agent}, rankings::RankedPreferences)
  num_proposals = length(proposing_group)

  ## construct a  map from nodes to agents
  other_mapping = Dict()
  proposing_mapping = Dict()
  node_to_agent_map = Dict()
  for i in 1:length(other_group)
    idx = i + num_proposals
    other_mapping[other_group[i]] = idx
    node_to_agent_map[idx] = other_group[i]
  end
  for i in 1:length(proposing_group)
    proposing_mapping[proposing_group[i]] = i
    node_to_agent_map[i] = proposing_group[i]
  end

  g = SimpleDiGraph(length(proposing_group) + length(other_group))
  for i in 1:length(proposing_group)
    cur_student = proposing_group[i]
    target_school = 0
    for cur_school in rankings.ranking[cur_student]
      if haskey(other_mapping, cur_school)
        target_school = cur_school
        break
      end
    end
    add_edge!(g, i, other_mapping[target_school])
  end
  for i in 1:length(other_group)
    cur_school = other_group[i]
    target_student = 0
    for cur_student in rankings.ranking[cur_school]
      if haskey(proposing_mapping, cur_student)
        target_student = cur_student
        break
      end
    end
    cur_idx = i + num_proposals
    add_edge!(g, cur_idx, proposing_mapping[target_student])
  end

    # Find the cycles and the matches!
  matchings = Dict{Agent, Agent}()
  cycles = simplecycles(g)
  for cycle in cycles
    for i in 1:2:length(cycle)
      agent_1 = node_to_agent_map[cycle[i]]
      agent_2 = node_to_agent_map[cycle[i+1]]
      if agent_1.capacity > 0 && agent_2.capacity > 0
        agent_2.capacity = agent_2.capacity - agent_1.capacity_reduction
        agent_1.capacity = agent_1.capacity - agent_2.capacity_reduction
        if cycle[i] <= num_proposals
          matchings[agent_1] = agent_2
        else
          matchings[agent_2] = agent_1
        end
	  end
    end
  end
  new_proposing_group = Vector{Agent}()
  for proposer in proposing_group
    if proposer.capacity > 0
      push!(new_proposing_group, proposer)
    end
  end
  new_other_group = Vector{Agent}()
  for other in other_group
    if other.capacity > 0
      push!(new_other_group, other)
    end
  end
  matchings, new_proposing_group, new_other_group
end

function ttc_with_counter(proposing_group::AbstractVector{Agent}, other_group::AbstractVector{Agent}, rankings::RankedPreferences)
  matchings = Dict{Agent, Agent}()
  while length(proposing_group) != 0 && length(other_group) != 0
    new_matchings, proposing_group, other_group = run_ttc_iteration(proposing_group, other_group, rankings)
    merge!(matchings, new_matchings)
  end
  Match(matchings)
end


text:
@testset "861.score-after-flipping-matrix.jl" begin
    @test matrix_score([[0,0,1,1],[1,0,1,0],[1,1,0,0]]) == 39
end


text:
using OptimTools
using Base.Test

# test exact line search for quadratic
Q = 10*speye(40)
f(x)    = 0.5*dot(x,Q*x)
df(x)   = Q*x
d2f(x)  = Q
x0 = randn(40)

exactLS  = (f,df,fk,dfk,xk,pk)-> (-dot(pk,dfk)/dot(pk,Q*pk),1)
goldenLSf = (f,df,dk,dfk,xk,pk) -> goldenLS(f,xk,pk,0.0,5.0)

println("Test SD,NLCG,BFGS with exact line search ")
# steepest descent
x1,flag1,his1 = sd(f,df,x0,maxIter=1000,atol=1e-6,out=0,lineSearch=exactLS)
@test norm(x1)<1e-6*norm(x0)
x2,flag2,his2 = sd(f,df,x0,maxIter=1000,atol=1e-6,out=0, lineSearch=exactLS)
@test norm(x2)<1e-6
@test size(his2,1)==2 #should converge in one iteration

# nlcg
x1,flag1,his1 = nlcg(f,df,x0,maxIter=1000,atol=1e-6,out=0)
@test norm(x1)<1e-6*norm(x0)
x2,flag2,his2 = nlcg(f,df,x0,maxIter=1000,atol=1e-6,out=0, lineSearch=exactLS)
@test norm(x2)<1e-6
@test size(his2,1)==2 #should converge in one iteration

# bfgs
x1,flag1,his1 = bfgs(f,df,x0,maxIter=1000,atol=1e-6,out=0)
@test norm(x1)<1e-6*norm(x0)
x2,flag2,his2 = bfgs(f,df,x0,maxIter=1000,atol=1e-6,out=0, lineSearch=exactLS)
@test norm(x2)<1e-6
@test size(his2,1)==2 #should converge in one iteration

println("Done!")

println("Test SD,NLCG,BFGS with golden line search ")
# steepest descent
x1,flag1,his1 = sd(f,df,x0,maxIter=1000,atol=1e-6,out=0,lineSearch=goldenLSf)
@test norm(x1)<1e-6*norm(x0)
x2,flag2,his2 = sd(f,df,x0,maxIter=1000,atol=1e-6,out=0, lineSearch=goldenLSf)
@test norm(x2)<1e-6

# nlcg
x1,flag1,his1 = nlcg(f,df,x0,maxIter=1000,atol=1e-6,out=0)
@test norm(x1)<1e-6*norm(x0)
x2,flag2,his2 = nlcg(f,df,x0,maxIter=1000,atol=1e-6,out=0, lineSearch=goldenLSf)
@test norm(x2)<1e-6

# bfgs
x1,flag1,his1 = bfgs(f,df,x0,maxIter=1000,atol=1e-6,out=0)
@test norm(x1)<1e-6*norm(x0)
x2,flag2,his2 = bfgs(f,df,x0,maxIter=1000,atol=1e-6,out=0, lineSearch=goldenLSf)
@test norm(x2)<1e-6

println("Done!")