text:
using Test
using IncrementalInference

@testset "Testing basic marginalization" begin
# linear hex
N=6
fg = generateCanonicalFG_lineStep(N;
                                  graphinit=false,
                                  poseEvery=1,
                                  landmarkEvery=N+1,
                                  posePriorsAt=[0],
                                  landmarkPriorsAt=[],
                                  sightDistance=N+1)

deleteFactor!.(fg, [Symbol("x$(i)lm0f1") for i=1:(N-1)])

# tree = buildTreeReset!(fg, show=true, drawpdf=true)
tree, smtasks, hists = solveTree!(fg)

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

defaultFixedLagOnTree!(fg, 6)
# Back up data from these two poses so we can compare them once we solve again.
lm0 = deepcopy(getVal(fg, :lm0))
X0 = deepcopy(getVal(fg, :x0))
X1 = deepcopy(getVal(fg, :x1))

fifoFreeze!(fg)

tree, smtasks, hists = solveTree!(fg; recordcliqs=ls(fg));
for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

# although variables are marginalized, no cliques are
@test calcCliquesRecycled(tree) == (6,0,0,0)

lm0cmp = deepcopy(getVal(fg, :lm0))
X0cmp = deepcopy(getVal(fg, :x0))
X1cmp = deepcopy(getVal(fg, :x1))

@test X0 == X0cmp #Frozen
@test lm0 == lm0cmp #Frozen
@test X1 != X1cmp #Recalculated

deleteVariable!(fg, :x0)

addVariable!.(fg, [Symbol("x$i") for i = 7:9], ContinuousScalar)
addFactor!(fg, [:x6,:x7], LinearRelative(Normal(1.0, 0.1)))
addFactor!(fg, [:x7,:x8], LinearRelative(Normal(1.0, 0.1)))
addFactor!(fg, [:x8,:x9], LinearRelative(Normal(1.0, 0.1)))
addFactor!(fg, [:lm0, :x9], LinearRelative(Normal(9,0.1)))

tree, smtasks, hists = solveTree!(fg; recordcliqs=ls(fg));

#clique 7 should be marginalized and therefor not do up or downsolve
@test calcCliquesRecycled(tree) == (7,1,0,0)
@test !(IIF.solveDown_StateMachine in getindex.(hists[7], 3))
@test !(IIF.solveUp_StateMachine in getindex.(hists[7], 3))

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

@test lm0 == getVal(fg, :lm0) #Still Frozen
@test X1cmp == getVal(fg, :x1) #also now Frozen

unfreezeVariablesAll!(fg)

# freeze again
defaultFixedLagOnTree!(fg, 9)

tree, smtasks, hists = solveTree!(fg)

@test lm0 == getVal(fg, :lm0) #Still Frozen
@test X1cmp != getVal(fg, :x1) #not frozen

# freeze 2,4,6 to all marginalize clique 2
setfreeze!(fg, [:x2, :x4, :x6])
tree, smtasks, hists = solveTree!(fg; recordcliqs=ls(fg));

#clique 2 should be marginalized and therefor not do up or downsolve
@test calcCliquesRecycled(tree) == (7,1,0,0)
@test !(IIF.solveDown_StateMachine in getindex.(hists[2], 3))
@test !(IIF.solveUp_StateMachine in getindex.(hists[2], 3))
@test areCliqVariablesAllMarginalized(fg, tree.cliques[2])

tree, smtasks, hists = solveTree!(fg, tree; recordcliqs=ls(fg));
for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.25)
end

@test calcCliquesRecycled(tree) == (7,1,6,0)

X1 = deepcopy(getVal(fg, :x1))

setfreeze!(fg, [:x3, :x5])
tree, smtasks, hists = solveTree!(fg, tree; recordcliqs=ls(fg));
# csmAnimate(tree, hists, frames=1)

@test calcCliquesRecycled(tree) == (7,3,4,0)

@test lm0 == getVal(fg, :lm0) #Still Frozen
@test X1 != getVal(fg, :x1) #not frozen

for i = [2,4,6]
    @test areCliqVariablesAllMarginalized(fg, tree.cliques[i])
end

end

@testset "Testing basic incremental recycle" begin

fg = generateCanonicalFG_lineStep(3;
                                   poseEvery=1,
                                   landmarkEvery=3,
                                   posePriorsAt=[],
                                   landmarkPriorsAt=[0],
                                   sightDistance=2,
                                   solverParams=SolverParams(algorithms=[:default, :parametric]))

getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
# getSolverParams(fg).dbg = true

# tree = buildTreeReset!(fg, drawpdf=true, show=true)

smtasks = Task[]
tree, smt, hists = solveTree!(fg; smtasks=smtasks, recordcliqs=ls(fg));

addFactor!(fg, [:lm3], Prior(Normal(3, 0.1)), graphinit=false)
_, _, hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

@test !(IIF.solveUp_StateMachine in getindex.(hists[4], 3))

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

end

@testset "Testing incremental hex" begin

N=6
sfg = generateCanonicalFG_lineStep(N;
                                   graphinit=false,
                                   poseEvery=1,
                                   landmarkEvery=N+1,
                                   posePriorsAt=[0],
                                   landmarkPriorsAt=[],
                                   sightDistance=N+1,
                                   solverParams=SolverParams(algorithms=[:default, :parametric]))

deleteFactor!.(sfg, [Symbol("x$(i)lm0f1") for i=1:(N-1)])

vsyms = sortDFG(ls(sfg))
fsyms = sortDFG(lsf(sfg))

fg = deepcopyGraph(LightDFG, sfg, vsyms[1:3], fsyms[1:3])

getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
getSolverParams(fg).dbg = true
getSolverParams(fg).useMsgLikelihoods = true

tree = buildTreeReset!(fg)#, drawpdf=true, show=true)

smtasks = Task[]
_,_,hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    # println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

# add some more
deepcopyGraph!(fg, sfg, vsyms[4:6], fsyms[4:6])

_,_,hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    # println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

# add some more
deepcopyGraph!(fg, sfg, vsyms[7:8], fsyms[7:8])

_,_,hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    # println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

# add some more,close the loop
deepcopyGraph!(fg, sfg, Symbol[], [fsyms[9]])

_,_,hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

# force a reshuffle
addFactor!(fg, [:x4], Prior(Normal(4.1,0.1)))
_,_,hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    # println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

# and another reuse
addFactor!(fg, [:x4], Prior(Normal(3.9,0.1)))

_,_,hists = solveTree!(fg, tree; smtasks=smtasks, recordcliqs=ls(fg));

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    # println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.25)
end

# printCSMHistoryLogical(hists)
# all except clique 1 should go for UPRECYCLED directly
#     | 1   x4             | 2   x2             | 3   x5             | 4   x3             | 5   lm0            | 6   x1
# ----+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------
# 1   | 1   setClique NULL | 3   setClique NULL | 4   setClique NULL | 5   setClique NULL | 6   setClique NULL | 7   setClique NULL
# 2   | 2   buildCliq NULL | 8   buildCliq UPRE | 9   buildCliq UPRE | 10  buildCliq UPRE | 11  buildCliq UPRE | 12  buildCliq UPRE
# 3   | 13  waitForUp NULL | 14  waitForUp UPRE | 15  waitForUp UPRE | 17  waitForUp UPRE | 19  waitForUp UPRE | 21  waitForUp UPRE
# 4   | 34  preUpSolv NULL | 31  preUpSolv UPRE | 16  preUpSolv UPRE | 18  preUpSolv UPRE | 20  preUpSolv UPRE | 22  preUpSolv UPRE
# 5   | 35  solveUp_S NULL | 32  postUpSol UPRE | 23  postUpSol UPRE | 24  postUpSol UPRE | 25  postUpSol UPRE | 26  postUpSol UPRE
# 6   | 36  postUpSol UPSO | 33  waitForDo UPRE | 27  waitForDo UPRE | 28  waitForDo UPRE | 29  waitForDo UPRE | 30  waitForDo UPRE
# 7   | 37  waitForDo UPSO | 39  solveDown UPRE | 40  solveDown UPRE | 43  solveDown UPRE | 41  solveDown UPRE | 44  solveDown UPRE
# 8   | 38  solveDown DOWN | 47  updateFro DOWN | 46  updateFro DOWN | 52  updateFro DOWN | 50  updateFro DOWN | 54  updateFro DOWN
# 9   | 42  updateFro DOWN | 49  exitState DOWN | 48  exitState DOWN | 53  exitState DOWN | 51  exitState DOWN | 55  exitState DOWN
# 10  | 45  exitState DOWN |

end


text:
module KhepriLibrary
using KhepriBase

include("TestCells.jl")

end

ext: jl


text:
# Confidence intervals should be tested component-wise against
# expected results; splatting them into arrays and checking for
# approximate array equality is incorrect when either limit is
# infinite.
macro test_ci_approx(x::Expr, y::Expr)
    return quote
        Base.Test.@test typeof($(esc(x))) <: Tuple{Real,Real}
        Base.Test.@test typeof($(esc(x))) == typeof($(esc(y)))
        Base.Test.@test all(map(isapprox, $(esc(x)), $(esc(y))))
    end
end

include("binomial.jl")
include("circular.jl")
include("fisher.jl")
include("kolmogorov_smirnov.jl")
include("kruskal_wallis.jl")
include("mann_whitney.jl")
include("t.jl")
include("z.jl")
include("wilcoxon.jl")
include("power_divergence.jl")
include("anderson_darling.jl")
include("common.jl")


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

# Compute coboundary_1 in 2D via product FV * EV^t with fixing of redundancies

FV = [[1,2,3,4,5,17,16,12],
[1,2,3,4,6,7,8,9,10,11,12,13,14,15],
[4,5,9,11,12,13,14,15,16,17],
[2,3,6,7], [8,9,10,11]]

EV = [[1,2],[2,3],[3,4],[4,5],[1,12],[2,6],[3,7],[4,9],[5,17],[6,7],[8,9],
[8,10],[9,11],[10,11],[11,15],[12,13],[12,16],[13,14],[14,15],[16,17]]

V = Float64[ 0   2   5   7  10   2   5   3   7  3  7  0  3  3  7  0  10;
            16  16  16  16  16  13  13  11  11  8  8  5  5  2  2  0   0]

# non-convex cells (convex=false), including outer cell (exterior=true)
copFE = Lar.coboundary_1( V::Lar.Points, FV::Lar.Cells, EV::Lar.Cells, false,true );
Matrix(copFE)

VV = [[k] for k in 1:size(V,2)];
using Plasm
Plasm.view( Plasm.numbering(3)((V,[VV, EV, FV])) )

# final solid visualization
copEV = Lar.coboundary_0(EV)
W = convert(Lar.Points, V')
triangulated_faces = Lar.triangulate2D(W, [copEV, copFE])
#triangulated_faces = Lar.triangulate2d(V, EV)
FVs = convert(Array{Lar.Cells}, triangulated_faces)
Plasm.viewcolor(V::Lar.Points, FVs::Array{Lar.Cells})

#TODO: debug ...


text:
# Compares the time of simulation and
include("../../src/ttv.jl")
include("/Users/ericagol/Computer/Julia/regress.jl")

# Turn off fastkicks:
nopair = false

elements = readdlm("elements_rebound.txt",',')

#n = 4
rstar = -1e12
for n=2:11
# Now, try calling with kicks between planets rather than -drift+Kepler:
pair_input = ones(Bool,n,n)
# We want Keplerian between star & planets, and impulses between
# planets.  Impulse is indicated with 'true', -drift+Kepler with 'false':
for i=2:n
  pair_input[1,i] = false
  # We don't need to define this, but let's anyways:
  pair_input[i,1] = false
end
# Now, only include Kepler solver for adjacent planets:
for i=2:n-1
  pair_input[i,i+1] = false
  pair_input[i+1,i] = false
end

# Loop over the number of planets in the simulation from 2 to 8:
t0 = 0.0
h  = 0.05
# Rebound integration is for 5000 code units;  inner planet has
# a period of 2pi, while in our case 1-day, so differs by 5000/(2pi):
tmax = 2500/pi

ntt = zeros(Int64,n)
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+2
end
#println("ntt: ",ntt)
tt1 = zeros(n,maximum(ntt))
tt = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count1 = zeros(Int64,n)
count = zeros(Int64,n)
# Call the ttv function, without gradient:
rstar = -1e12
if nopair
  ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
else
  ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar;pair=pair_input)
end

tic()
if nopair
  ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
else
  ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar;pair=pair_input)
end
var_0 = toq()

# Now call with gradient:

dtdq0 = zeros(n,maximum(ntt),7,n)
dtdq0 = zeros(n,maximum(ntt),7,n)
#dtdelements0 = zeros(n,maximum(ntt),7,n)
#@time dtdelements0 = ttv_elements!(n,t0,h,tmax,elements,tt2,count1,dtdq0,rstar)
#dtdelements0 = ttv_elements!(n,t0,h,tmax,elements,tt2,count1,dtdq0,rstar)
#println("n ",n," t0 ",t0," h ",h," tmax ",tmax," elements ",elements)
if nopair
  dtdelements = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar)
else
  dtdelements = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar;pair=pair_input)
end
tic()
if nopair
  dtdelements = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar)
else
  dtdelements = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar;pair=pair_input)
end
var_1 = toq()

println("Nplanet: ",n-1," No gradient: ",var_0, " gradient: ", var_1, " ratio: ", var_1/var_0)
end

max_forks_repo_name: langfzac/nbg-papers


text:
include("dae_initialization.jl")
include("../../components/ThirdOrderEq.jl")

include("system.jl")

powergrid = PowerGrid(node_list,line_list)
operationpoint = find_operationpoint(powergrid)

timespan = (0.0, 30)

# I choose du0 = 0 because I assume that we start from a stable fixpoint for now
# but this should definitely change later on....
# Seems like a valid DAEProblem is found. The function is of type nd_ODE_ODE{...}
DAEProblem(dae_rhs(powergrid),zeros(systemsize(powergrid)),operationpoint.vec,timespan)

# The function looks the same as in the DAEProblem
ODEProblem{iipfunc}(rhs(powergrid),operationpoint.vec,timespan)

solve(powergrid, operationpoint,timespan)

# Fails because there is no method for nd_ODE_ODE{...}.
# But exactly this type seems to work in the case of ODE solvers...
dae_solve(powergrid, operationpoint, timespan)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
# This implements the L0 EM algorithm for regularized regression as described in:
#
#  Zhenqiu Liu and Gang Li. "Efficient Regularized Regression for Variable
#  Selection with L0 Penalty" Submitted to arXiv on July 29th 2014.
#
# theta = l0EM(X, y, lambda; kws...)
#
# Find coefficients 'theta' for regularized regression with L0 penalty, i.e.
# theta minimizes 0.5*norm(y - X*theta) + lambda/2*sum(theta .> 0.0).
#
# Outputs:
#  'theta' is an M*1 vector with the (sparse) coefficients.
#
# Inputs:
#  'X' is an N by M matrix.
#
#  'y' is an N by 1 array.
#
#  'lambda' is the regularization parameter. By default it is set via AIC, i.e. to 2.
#    Other, typically valid values are log(N) (BIC) and 2*log(M) (RIC).
#
# Keyword args:
#  'epsilon' is the minimum coefficient value that is allowed, below which the
#    coefficient is set to zero. Defaults to 1e-4.
#
#  'delta_treshold' is the treshold value below which the iteration is deemed to have
#    converged. Defaults to 1e-5.
#
#  'nonnegative' is true iff only non-negative coefficients are allowed. Defaults to false.
#
#  'maxIterations' is max number of iterations (in case a lambda given so it doesn't converge).
#
function l0EM(X, y, lambda = 2;
  epsilon = 1e-4, deltaTreshold = 1e-5,
  nonnegative = false, maxIterations = 10000)

  N, M = size(X)
  #lambda_eye = UniformScaling(lambda) # Same as but often faster than: lambda * eye(N)
  lambda_eye = lambda * eye(N)
  xt = X'
  theta = xt * ((X * xt .+ lambda_eye) \ y) # Same as but often faster / more stable than: theta = xt * inv(X * xt .+ lambda_eye) * y

  if nonnegative
    set_negative_to_zero!(theta)
  end

  iterations = 0

  while iterations < maxIterations
    iterations += 1

    # E-step:
    eta = theta

    # M-step:
    eta_squared = eta .^ 2
    xt_eta = broadcast(*, eta_squared, xt)
    theta = xt_eta * ((X * xt_eta .+ lambda_eye) \ y) # Same as but often faster / more stable than: theta = xt_eta * inv(X * xt_eta .+ lambda_eye) * y

    # Ensure non-negativity if required
    if nonnegative
      set_negative_to_zero!(theta)
    end

    # We have converged if change is too small => break.
    if norm(theta - eta) < deltaTreshold
      break
    end
  end

  theta[abs(theta) .< epsilon] = 0.0

  theta
end

function set_negative_to_zero!(v)
  v[v .< 0.0] = 0.0
  v
end

# Given mild constraints on the co-linearity of features there is
# a max lambda value we can calculate.
function find_max_lambda(X, y)
  lambda_from_col_j(j) = (sum(X[:,j] .* y)^2) / (4 * sum(X[:,j].^2))
  N, M = size(X)
  maximum(map(lambda_from_col_j, 1:M))
end

num_selected(v) = sum(v .> 0.0)

# Instead of checking a fixed number of lambda values in a range from 1e-8
# to max_lambda we do a binary search to identify lambda values where there
# is a transition in the number of selected variables.
#
#  cs, ns = adaptive_lambda_regularized_regression(X, y, regressor)
#
# Outputs:
#  'cs' is a dict mapping the tried lambda values to the regressed coefficients for
#    each lambda.
#
#  'ns' is a dict mapping an integer number of selected variables to the lambda value
#    for which that number of selected variables was first encountered.
#
# Inputs:
#  'X' is an N by M matrix corresponding to N cases each with M measured features.
#
#  'y' is an N by 1 array corresponding to the dependent variable for the N cases.
#
#  'regressor' is a Julia function to call to perform regularized regression.
#    It needs to have an interface like regressor(X, y, lambda; kws...).
#    Defaults to using the L0 EM algorithm for regularized regression.
#
# Optional keyword arguments:
#  'minLambda' is the minimum lambda value to use. Defaults to 1e-7 and might not
#    converge at all if a smaller value is selected.
#
#  'numLambdas' is an indication (not exact) of the number of calls to the regressor
#    that we are allowed to do. The actual number of calls can deviate somewhat but not a lot.
#
#  'stepDivisor' determines how small lambda increments are considered when searching in
#    between the standard log-spaced values of LASSO and similar algorithms. Defaults to
#    2^12 which tends to give a sensible number of alternative models to choose from.
#    Increase it to find a large number of selected number of variables (returned in 'ns').
#
#  'logSpace' is true iff the search should be in log space i.e. if taking the logarithm
#    into account when selecting mid points in the binary search.
#
#  'includeAICetal' is true iff we ensure AIC, BIC and RIC are included among the lambdas
#    to be tested. The stepDivisor should typically be set lower (2^6) is this is set to true.
#    Defaults to false.
#
#  'kws' are any other keyword arguments that will be passed on to the regressor.
function adaptive_lambda_regularized_regression(X, y, regressor = l0EM;
  minLambda = 1e-7, numLambdas = 100,
  stepDivisor = 2^12,
  logSpace = true,
  includeAICetal = false,
  kws...)

  max_lambda = find_max_lambda(X, y)
  log_increment = (log(max_lambda + 1) - log(minLambda + 1)) / (numLambdas-1)

  num_vars_to_lambda = Dict{Int64, Float64}()
  coeffs = Dict{Float64, Matrix{Float64}}()

  num_calls_to_regularized_regressor = 0

  update_if_not_there!(l) = begin
    if haskey(coeffs, l)
      cs = coeffs[l]
    else
      coeffs[l] = cs = regressor(X, y, l; kws...)
      num_calls_to_regularized_regressor += 1
    end
    nvars = num_selected(cs)
    if !haskey(num_vars_to_lambda, nvars)
      num_vars_to_lambda[nvars] = l
    end
    return (cs, nvars)
  end

  midpoint(minL, maxL) = begin
    if logSpace
      log_minL = log(minL)
      exp( log_minL + (log(maxL) - log_minL) / 2 )
    else
      minL + (maxL - minL) / 2
    end
  end

  # Search the lambda values from minLambda to maxLambda
  binary_index_search(minLambda, maxLambda, minDelta) = begin
    #println("bis($minLambda, $maxLambda, $minDelta)")
    if num_calls_to_regularized_regressor <= numLambdas && (maxLambda - minLambda > minDelta)
      csmin, nvarmin = update_if_not_there!(minLambda)
      #println("  nvarmin = $nvarmin ($minLambda)")
      csmax, nvarmax = update_if_not_there!(maxLambda)
      #println("  nvarmax = $nvarmax ($maxLambda)")
      if (nvarmin - nvarmax > 1)
        mid = midpoint(minLambda, maxLambda)
        csmid, nvarmid = update_if_not_there!(mid)
        #println("  nvarmid = $nvarmid ($mid)")
        if nvarmin - nvarmid > 1
          binary_index_search(minLambda, mid, minDelta)
        end
        if nvarmid - nvarmax > 1
          binary_index_search(mid, maxLambda, minDelta)
        end
      end
    end
  end

  if includeAICetal

    # Always include the AIC, BIC and RIC endpoints and search them first to ensure
    # they are included.
    a, b = 2.0 < log(N) ? (2.0, log(N)) : (log(N), 2.0)
    binary_index_search(a, b, log_increment/stepDivisor)
    b, c = b < (2*log(M)) ? (b, 2*log(M)) : (2*log(M), b)
    binary_index_search(b, c, log_increment/stepDivisor)

    # And then from minLambda up to AIC and from RIC up to max lambda.
    binary_index_search(minLambda, a, log_increment/stepDivisor)
    binary_index_search(c, max_lambda, log_increment/stepDivisor)

  else
    binary_index_search(minLambda, max_lambda, log_increment/stepDivisor)
  end

  return (coeffs, num_vars_to_lambda)
end

mse(theta, X, y) = mean( (X * theta - y).^2 )

# Basic test:
function basic_l0EM_test(N, M, NumReps = 30)
  for rep in 1:NumReps
    k = rand(1:iceil(2*log(M)))
    X = randn(N, M)
    actual_theta = [linspace(1,k,k), zeros(M-k)]
    y = X * actual_theta + 0.10 * randn(N, 1)

    cs, ns = adaptive_lambda_regularized_regression(X, y, l0EM;
      stepDivisor = 2^4, includeAICetal = true);

    # Find lambda with lowest MSE on test set
    Xtest = randn(N, M)
    ytest = Xtest * actual_theta + 0.10 * randn(N, 1)
    best_mse = Inf
    local best_lambda, best_nvars
    for lambda in values(ns)
      newmse = mse(cs[lambda], Xtest, ytest)
      if newmse < best_mse
        best_nvars = num_selected(cs[lambda])
        best_lambda = lambda
        best_mse = newmse
      end
    end
    println("k = $k, Lambda = $best_lambda, num vars = $best_nvars, MSE = $best_mse")
  end
end

basic_l0EM_test(50, 100, 10)

# Using log space binary search seems to speed the search up about 10-30% depending on problem size.
# Potentially more for big data sets and with cross validation.
#@time cs2, ns2 = adaptive_lambda_regularized_regression(X, y, l0EM; stepDivisor = 2^2, includeAICetal = true);
#@time cs4, ns4 = adaptive_lambda_regularized_regression(X, y, l0EM; stepDivisor = 2^4, includeAICetal = true);
#@time cs6, ns6 = adaptive_lambda_regularized_regression(X, y, l0EM; stepDivisor = 2^6, includeAICetal = true);
#@time cs4b, ns4b = adaptive_lambda_regularized_regression(X, y, l0EM; stepDivisor = 2^4, includeAICetal = false);
#@time csf, nsf = adaptive_lambda_regularized_regression(X, y, l0EM; logSpace = false);
#sort(collect(keys(ns)))
#length(cs)

# If not using the adaptive (binary) search the standard method is to do 100 lambda values
# in a log scale:
#function log_split_lambdas(X, y, numLambdas = 10)
#  exp(linspace(0.0, log(find_max_lambda(X, y) + 1), numLambdas)) - 1.0 + 1e-6
#end
#ts100 = map((l) -> l0EM(X, y, l), log_split_lambdas(X, y, 100))
#num_vars = map((t) -> sum(t .> 0.0), ts100)


text:
 using Distributed
@everywhere using DifferentialEquations
@everywhere using ProgressMeter

include("PSS_function_parallel.jl")
include("Yfind.jl")
# include("is_it.jl")

@everywhere using PyCall
pygui(:qt)
@everywhere using PyPlot
pygui(true)
#PSS_function(Q,P, Energy, t_end, Max_hit)
@everywhere max_hit=10
@everywhere Energy=.25
@everywhere H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
@everywhere Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
@everywhere ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
@everywhere ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
@everywhere function Eq_of_M(du,u,p,t)
     du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
     du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
     du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
     du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
     du[5]=0
 return
end
@everywhere condition(u,t,integrator)= u[5]>max_hit
@everywhere function condition2(u,t,integrator) # Event when event_f(u,t) == 0
 u[1]
end
@everywhere affect!(integrator) = terminate!(integrator)
@everywhere function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
 end
@everywhere cb2 = ContinuousCallback(condition2,nothing,affect2!,  save_positions = (false,true))
@everywhere cb1 = DiscreteCallback(condition,affect!)

@everywhere cb=CallbackSet(cb2,cb1)
@everywhere H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))

n_iter_P=200
n_iter_Q=200

N=n_iter_P*n_iter_Q

T=@showprogress pmap(PSS_function_parallel,1:N)
COLOR=["b", "g", "c", "m", "y", "k", "r"]
#
# outfileQ = "PSSoutfileQ_"*string(Energy)*"nP"*string(n_iter_P)*"nQ"*string(n_iter_Q)*".dat"
# outfileP = "PSSoutfileP_"*string(Energy)*"nP"*string(n_iter_P)*"nQ"*string(n_iter_Q)*".dat"
# outfilePAIR = "PSSoutfileDuple_"*string(Energy)*"nP"*string(n_iter_P)*"nQ"*string(n_iter_Q)*".dat"
#
# a = open(outfileQ, "w")
# b = open(outfileP, "w")
# c = open(outfilePAIR, "w")
# Q_tabs=zeros(0)
# P_tabs=zeros(0)
# QQ=zeros(max_hit+2,N)
# PP=zeros(max_hit+2,N)
figure()
axis([ -2.5, 2.5,-1.5, 1.5])
# axis([ -.02,.02,-.2, .2])

@showprogress 1 "Computing..." for k=1:N
    if T[k]!=nothing
        current_color=COLOR[mod(k,length(COLOR))+1]
        plot(T[k][2],T[k][1], color=current_color,".",markersize=1, markeredgewidth=.1)
        plot(T[k][2],-T[k][1], color=current_color,".",markersize=1, markeredgewidth=.1)
        plot(-T[k][2],T[k][1], color=current_color,".",markersize=1, markeredgewidth=.1)
        plot(-T[k][2],-T[k][1], color=current_color,".",markersize=1, markeredgewidth=.1)
        # for j=1:length(T[k][1])
        #     # push!(Q_tabs,T[k][1][j])
        #     # push!(P_tabs,T[k][2][j])
        #     QQ[j,k]=T[k][1][j]
        #     PP[j,k]=T[k][2][j]
        # end
    end
end
# Q=vcat(Q_tabs,-Q_tabs,Q_tabs,-Q_tabs)
# P=vcat(P_tabs,P_tabs,-P_tabs,-P_tabs)
# # PAIR=hcat(Q,P)
# plot(P,Q ,color="k",".",markersize=1, markeredgewidth=.1)
# writedlm(outfileQ, QQ)
# writedlm(outfileP, PP)
#
# # writedlm(outfilePAIR, PAIR)
#
# close(a)
# close(b)


text:
# Configuration for nearly neutral simuation
const nn_simtype = 2    # nearly neutral infinite sites
const type_str = "N2"
const mu_list_flag = true
@everywhere const N_list     = [25, 50, 100,200,400,800,1600, 3200, 6400]     # popsize
const mu_list= [0.025,0.01,0.005,0.0025]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 500000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
dfe = dfe_neutral
dfe_str = "dfe_neutral"

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
import ChemistryFeaturization: encode, encodable_elements, decode, output_shape, features
using ChemistryFeaturization.Data
using ChemistryFeaturization.ElementFeature

using DataFrames

const default_nbins = 10

"""
    GraphNodeFeaturization(atom_features, codecs)
    GraphNodeFeaturization(feature_names, lookup_table; nbins, logspaced, categorical)

A featurization for AtomGraph objects that encodes features associated with each node. Contains a collection of `AbstractAtomFeatureDescriptor` objects, and can be initialized by passing those, or by passing parameters for constructing them.

## Required Arguments
- `feature_names::Vector{String}`: Names of each feature
- `lookup_table::DataFrame` (optional): source of data for features, if not included in built-in `element_data_df`

## Keyword Arguments
- `nbins::Union{Vector{Integer},Integer}`: Number of bins to use for one-hot encoding of continuous-valued features. Will be ignored for categorical features.
- `logspaced::Union{Vector{Bool},Bool}`: Whether to logarithmically space the bins
- `categorical::Union{Vector{Bool},Bool}`: Whether each feature is categorical or continous-valued.
"""
struct GraphNodeFeaturization <: AbstractFeaturization
    features::Vector{<:AbstractAtomFeatureDescriptor}
    codecs::Vector{<:AbstractCodec}
end

GraphNodeFeaturization(features::Vector{<:AbstractAtomFeatureDescriptor}) =
    GraphNodeFeaturization(features, default_codec.(features))

# NB: this constructor only works if every feature is an ElementFeature
function GraphNodeFeaturization(
    feature_names::Vector{String},
    lookup_table::Union{DataFrame,Nothing} = nothing;
    nbins::Union{Vector{<:Integer},Integer,Nothing} = nothing,
    logspaced::Union{Vector{Bool},Bool,Nothing} = nothing,
    categorical::Union{Vector{Bool},Bool,Nothing} = nothing,
)
    num_features = length(feature_names)
    local lookup_table_here, logspaced_here, categorical_here, nbins_here
    if isnothing(lookup_table)
        lookup_table_here = element_data_df
    else
        # need to merge them in case some data is in one place and some the other
        lookup_table_here = outerjoin(element_data_df, lookup_table, on = :Symbol)
    end

    if isnothing(logspaced)
        logspaced_here = map(fn -> default_log(fn, lookup_table_here), feature_names)
    else
        logspaced_here = get_param_vec(logspaced, num_features)
    end

    if isnothing(categorical)
        categorical_here =
            map(fn -> default_categorical(fn, lookup_table_here), feature_names)
    else
        categorical_here = get_param_vec(categorical, num_features)
    end

    if isnothing(nbins)
        nbins_here = [default_nbins for i = 1:num_features]
    else
        nbins_here = get_param_vec(nbins, num_features, pad_val = default_nbins)
    end

    bins = map(zip(feature_names, nbins_here, logspaced_here, categorical_here)) do args
        get_bins(
            args[1],
            lookup_table_here,
            nbins = args[2],
            logspaced = args[3],
            categorical = args[4],
        )
    end

    codecs = OneHotOneCold.(categorical_here, bins)
    efs = ElementFeatureDescriptor.(feature_names, Ref(lookup_table_here))

    GraphNodeFeaturization(efs, codecs)
end

# pretty printing, short
function Base.show(io::IO, fzn::GraphNodeFeaturization)
    st = "GraphNodeFeaturization encoding $(length(fzn.features)) features"
    print(io, st)
end

# pretty printing, long
function Base.show(io::IO, ::MIME"text/plain", fzn::GraphNodeFeaturization)
    st = "GraphNodeFeaturization encoding $(length(fzn.features)) features:\n"
    for feature in fzn.features
        st = string(st, "\t", feature, "\n")
    end
    st = st[1:end-1]
    print(io, st)
end

output_shape(fzn::GraphNodeFeaturization) = sum(output_shape.(fzn.codecs))

features(fzn::GraphNodeFeaturization) = fzn.features

"""
    chunk_vec(vec, nbins)

Helper function that divides up an already-constructed feature vector into "chunks" (one for each feature) of lengths specified by the vector nbins.

Sum of nbins should be equal to the length of vec.

# Examples
```jldoctest
julia> chunk_vec([1,0,0,1,0], [3,2])
2-element Array{Array{Bool,1},1}:
 [1, 0, 0]
 [1, 0]
 ```
"""
function chunk_vec(vec::Vector{<:Real}, nbins::Vector{<:Integer})
    chunks = fill(Bool[], size(nbins, 1))
    @assert length(vec) == sum(nbins) "Total number of bins doesn't match length of feature vector."
    for i = 1:size(nbins, 1)
        if i == 1
            start_ind = 1
        else
            start_ind = sum(nbins[1:i-1]) + 1
        end
        chunks[i] = vec[start_ind:start_ind+nbins[i]-1]
    end
    return chunks
end

function encode(ag::AtomGraph, fzn::GraphNodeFeaturization)
    encoded = reduce(hcat, map(x -> encode(ag, x[1], x[2]), zip(fzn.features, fzn.codecs)))
    return encoded
end

function decode(encoded::Matrix{<:Real}, fzn::GraphNodeFeaturization)
    num_atoms = size(encoded, 1)
    nbins = [output_shape(c) for c in fzn.codecs]
    local decoded = Dict{Integer,Dict{String,Any}}()
    for i = 1:num_atoms
        #println("atom $i")
        chunks = chunk_vec(encoded[i, :], nbins)
        decoded[i] = Dict{String,Any}()
        for (chunk, f, c) in zip(chunks, fzn.features, fzn.codecs)
            #println("    $(f.name): $(decode(f, chunk))")
            decoded[i][f.name] = decode(chunk, c)
        end
    end
    return decoded
end


text:
module ControlSystemIdentification

using DSP,
    LinearAlgebra,
    Statistics,
    StatsBase,
    Random,
    ComponentArrays,
    Optim,
    ControlSystems,
    FillArrays,
    Parameters,
    TotalLeastSquares,
    RecipesBase,
    FFTW,
    Roots,
    MonteCarloMeasurements,
    LowLevelParticleFilters,
    BandedMatrices
import StatsBase: predict
import MatrixEquations
import Optim: minimizer, Options
import ControlSystems: ninputs, noutputs, nstates
import StatsBase.residuals

export iddata,
    noutputs,
    ninputs,
    nstates,
    input,
    output,
    sampletime,
    hasinput,
    apply_fun,
    resample,
    timevec
export StateSpaceNoise,
    pem, simulation_errors, prediction_errors, predict, simulate, noise_model
export n4sid, era, okid
export getARXregressor,
    getARregressor,
    find_na,
    arx,
    ar,
    arxar,
    residuals,
    arma,
    arma_ssa,
    armax,
    impulseest,
    bodeconfidence,
    tls,
    wtls_estimator,
    plr,
    estimate_residuals
export FRD, tfest, coherence, coherenceplot, simplot, simplot!, predplot, predplot!, Hz, rad

export model_spectrum

export KalmanFilter

export weighted_estimator, Bandpass, Bandstop, Lowpass, Highpass

export kautz, laguerre, laguerre_oo, adhocbasis, sum_basis, basis_responses, filter_bank, basislength, ωζ2complex, add_poles, minimum_phase

include("utils.jl")
include("types.jl")
include("pem.jl")
include("frd.jl")
include("arx.jl")
include("subspace.jl")
include("spectrogram.jl")
include("frequency_weights.jl")
include("basis_functions.jl")

predict(sys, d::AbstractIdData, args...) =
    hasinput(sys) ? predict(sys, output(d), input(d), args...) :
    predict(sys, output(d), args...)

function predict(sys, y, u, x0 = zeros(sys.nx))
    model = SysFilter(sys, copy(x0))
    yh = [model(yt, ut) for (yt, ut) in observations(y, u)]
    oftype(y, yh)
end
predict(sys::ControlSystems.TransferFunction, args...) = predict(ss(sys), args...)

"""
	yh = predict(ar::TransferFunction, y)

Predict AR model
"""
function predict(G::ControlSystems.TransferFunction, y)
    _, a, _, _ = params(G)
    yr, A = getARregressor(output(y), length(a))
    yh = A * a
    oftype(output(y), yh)
end

function simulate(sys, u, x0 = zeros(sys.nx))
    u = input(u)
    model = SysFilter(sys, copy(x0))
    yh = map(observations(u, u)) do (ut, _)
        model(ut)
    end
    oftype(u, yh)
end
simulate(sys::ControlSystems.TransferFunction, args...) = simulate(ss(sys), args...)

@userplot Simplot
"""
	simplot(sys, data, x0=zeros(sys.nx); ploty=true)

Plot system simulation and measured output to compare them.
`ploty` determines whether or not to plot the measured signal
"""
simplot
@recipe function simplot(p::Simplot; ploty = true)
    sys, d = p.args[1:2]
    y = oftype(randn(2, 2), output(d))
    u = oftype(randn(2, 2), input(d))
    x0 = length(p.args) > 3 ? p.args[4] : zeros(sys.nx)
    yh = simulate(sys, u, x0)
    xguide --> "Time [s]"
    yguide --> "Output"
    t = timevec(d)
    err = nrmse(y, yh)
    ploty && @series begin
        label --> "y"
        t, y'
    end
    @series begin
        label --> "sim fit :$(round.(err, digits=2))%"
        t, yh'
    end
    nothing
end

@userplot Predplot
"""
	predplot(sys, data, x0=zeros(sys.nx); ploty=true)

Plot system simulation and measured output to compare them.
`ploty` determines whether or not to plot the measured signal
"""
predplot
@recipe function predplot(p::Predplot; ploty = true)
    sys, d = p.args[1:2]
    y = oftype(randn(2, 2), output(d))
    u = oftype(randn(2, 2), input(d))
    x0 = length(p.args) > 3 ? p.args[4] : zeros(sys.nx)
    yh = predict(sys, y, u, x0)
    xguide --> "Time [s]"
    yguide --> "Output"
    t = timevec(d)
    err = nrmse(y, yh)
    ploty && @series begin
        label --> "y"
        t, y'
    end
    @series begin
        label --> "pred fit :$(round.(err, digits=2))%"
        t, yh'
    end
    nothing
end

function ControlSystems.lsim(sys::StateSpaceNoise, u; x0 = zeros(sys.nx))
    simulate(sys, input(u), x0)
end

ControlSystems.innovation_form(sys::Union{StateSpaceNoise,N4SIDStateSpace}) =
    ss(sys.A, sys.K, sys.C, Matrix(Eye(sys.ny)), sys.Ts) # innovation model

end # module


text:
using Pkg
Pkg.update()

try
    using UnicodePlots
    using PyCall
    using NSGAII
    using OrderedCollections
    #using Plotly
catch
    #import Pkg
    #Pkg.add("Plotly")
    import Pkg; Pkg.add("Plots")
    import Pkg; Pkg.add("UnicodePlots")
    Pkg.add("OrderedCollections")

    #Pkg.add("PackageCompiler.jl")
    #Pkg.clone("https://github.com/JuliaCN/Py2Jl.jl")
    Pkg.clone("https://github.com/gsoleilhac/NSGAII.jl")
    Pkg.add("ProgressMeter")
    Pkg.add("UnicodePlots")
    Pkg.add("Plots")
    Pkg.add("OrderedCollections")
    Pkg.add("PyCall")
    Pkg.build("PyCall")
    #Pkg.clone("https://github.com/gsoleilhac/NSGAIII.jl")
    #Pkg.clone("https://github.com/JuliaCN/Py2Jl.jl")

end
Pkg.update()
using PyCall
using UnicodePlots

include("src/SpikingNeuralNetworks.jl")
include("src/units.jl")
include("src/plot.jl")
using UnicodePlots
using PyCall
using NSGAII
SNN = SpikingNeuralNetworks.SNN
#=
py"""
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
matplotlib.get_backend()
import sys
import os
sys.path.append(os.getcwd())
"""
py"""
from neuronunit import tests
from neuronunit.tests import fi
from examples.simple_with_injection import SimpleModel
from neo import AnalogSignal
import quantities as pq
try:
   from julia import Main
except:
   import julia
   julia.install()
   from julia import Main

"""
#=
simple = py"SimpleModel()"

py"""
import pickle
from neuronunit.tests.fi import RheobaseTestP

cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
for test in cell_tests.values():
    if "Rheobase test" in test.keys():
        temp_test = {k:v for k,v in test.items()}
        break
rt = temp_test["Rheobase test"]

rtest = RheobaseTestP(observation=rt.observation)
#                        name='RheobaseTest')
JHH = {
'Vr': -68.9346,
'Cm': 0.0002,
'gl': 1.0 * 1e-5,
'El': -65.0,
'EK': -90.0,
'ENa': 50.0,
'gNa': 0.02,
'gK': 0.006,
'Vt': -63.0
}
JHH1 = {k:(v-0.01*v) for k,v in JHH.items()}
JHH2 = {k:(v+0.01*v) for k,v in JHH.items()}
ranges = {k:[v-0.01*v,v+0.01*v] for k,v in JHH.items()}

"""
#
using OrderedCollections
ranges = OrderedDict(py"ranges")
H1=[values(ranges)]

current_params = py"rt.params"a
print(current_params)
simple.attrs = py"JHH"
=#
=#

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
function load_image(
    filename::AbstractString;
    dims::Tuple{Int, Int}, grayscale=true
)
    image = grayscale ? Gray.(load(filename)) : load(filename)
    image[1:dims[1], 1:dims[2]]
end

function load_images(
    image_dir::AbstractString, v::UnitRange{Int};
    dims::Tuple{Int, Int}, grayscale=true
)
    files = readdir(image_dir, join=true, sort=false)[v]
    [
        (
            load_image(file; dims=dims, grayscale=grayscale),
            match(r"\/([a-zA-Z0-9_]+)\.", file).captures[1]
        )
        for file in files
    ]
end

function load_images(
    image_dir::AbstractString;
    dims::Tuple{Int, Int}, grayscale=true
)
    v = 1:length(readdir(image_dir, join=true, sort=false))
    load_images(image_dir, v; dims=dims, grayscale=grayscale)
end

function images_to_image_matrix(images)
    reduce(hcat, map(x->reshape(x, :, 1), convert.(Array{Float64}, images)))
end

function image_matrix_to_images(image_matrix, dims)
    images = [
        reshape(image_matrix[:, i], dims[1], dims[2])
        for i in 1:size(image_matrix, 2)
    ]
    [Gray.(image) for image in images]
end

function tile_images(images, cols)
    rows = size(images, 1) ÷ cols
    reduce(
        vcat,
        (reduce(hcat, (images[i*cols+j] for j = 1:cols)) for i = 0:rows-1)
    )
end

function rescale(data)
    min = minimum(data)
    max = maximum(data)
    scale(x) = (x-min)/(max-min)
    scale.(data)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/evaluator/passes.jl
# Functions used to compute forward pass of nonlinear functions which include:
# set_value_post, overwrite_or_intersect, forward_pass_kernel, associated blocks
#############################################################################

const FORWARD_DEBUG = false

"""
$(FUNCTIONNAME)

Post process set_value operator. By default, performs the affine interval cut on
a MC structure.
"""
function set_value_post(x_values::Vector{Float64}, val::MC{N,T}, lower_variable_bounds::Vector{Float64},
                        upper_variable_bounds::Vector{Float64}, sparsity::Vector{Int},
                        subgrad_tol::Float64) where {N, T<:RelaxTag}

    lower = val.cv
    upper = val.cc
    lower_refinement = true
    upper_refinement = true

    for i = 1:N

        cv_val = @inbounds val.cv_grad[i]
        cc_val = @inbounds val.cc_grad[i]

        i_sol = @inbounds sparsity[i]
        x_val = @inbounds x_values[i_sol]
        lower_bound = @inbounds lower_variable_bounds[i_sol]
        upper_bound = @inbounds upper_variable_bounds[i_sol]

        if lower_refinement
            if cv_val > 0.0
                if isinf(lower_bound)
                    !upper_refinement && break
                    lower_refinement = false
                else
                    lower += cv_val*(lower_bound - x_val)
                    #delX = sub_round(lower_bound, x_val, RoundDown)
                    #lower = add_round(lower, mul_round(cv_val, delX, RoundDown), RoundDown)
                end
            else
                if isinf(upper_bound)
                    !upper_refinement && break
                    lower_refinement = false
                else
                    lower += cv_val*(upper_bound - x_val)
                    #delX = sub_round(upper_bound, x_val, RoundUp)
                    #lower = add_round(lower, mul_round(cv_val, delX, RoundDown), RoundDown)
                end
            end
        end

        if upper_refinement
            if cc_val > 0.0
                if isinf(lower_bound)
                    !lower_refinement && break
                    upper_refinement = false
                else
                    upper += cc_val*(upper_bound - x_val)
                    #delX = sub_round(upper_bound, x_val, RoundUp)
                    #upper = add_round(upper, mul_round(cc_val, delX, RoundUp), RoundUp)
                end
            else
                if isinf(upper_bound)
                    !lower_refinement && break
                    upper_refinement = false
                else
                    upper += cc_val*(lower_bound - x_val)
                    #delX = sub_round(lower_bound, x_val, RoundDown)
                    #upper = add_round(upper, mul_round(cc_val, delX, RoundUp), RoundUp)
                end
            end
        end
    end

    if lower_refinement && (val.Intv.lo + subgrad_tol > lower)
        lower = val.Intv.lo
    elseif !lower_refinement
        lower = val.Intv.lo
    else
        lower -= subgrad_tol #sub_round(lower, subgrad_tol, RoundDown)
    end

    if upper_refinement && (val.Intv.hi - subgrad_tol < upper)
        upper = val.Intv.hi
    elseif !upper_refinement
        upper = val.Intv.hi
    else
        upper += subgrad_tol #add_round(upper, subgrad_tol, RoundUp)
    end

    return MC{N,T}(val.cv, val.cc, Interval{Float64}(lower, upper), val.cv_grad, val.cc_grad, val.cnst)
end

"""
$(FUNCTIONNAME)

Intersects the new set valued operator with the prior and performs affine bound tightening

- First forward pass: `is_post` should be set by user option, `is_intersect` should be false
  so that the tape overwrites existing values, and the `interval_intersect` flag could be set
  to either value.
- Forward CP pass (assumes same reference point): `is_post` should be set by user option,
  `is_intersect` should be true so that the tape intersects with  existing values, and the
  `interval_intersect` flag should be false.
- Forward CP pass (assumes same reference point): `is_post` should be set by user option,
  `is_intersect` should be true so that the tape intersects with existing values, and the
  `interval_intersect` flag should be false.
- Subsequent forward passes at new points: is_post` should be set by user option,
  `is_intersect` should be true so that the tape intersects with existing values, and the
  `interval_intersect` flag should be `true` as predetermined interval bounds are valid but
   the prior values may correspond to different points of evaluation.
"""
function overwrite_or_intersect(xMC::MC{N,T}, past_xMC::MC{N,T}, x::Vector{Float64}, lbd::Vector{Float64},
                                ubd::Vector{Float64}, subgrad_tol::Float64, sparsity::Vector{Int}, is_post::Bool,
                                is_intersect::Bool,
                                interval_intersect::Bool) where {N,T<:RelaxTag}

    if is_post && is_intersect && interval_intersect
        return set_value_post(x, xMC ∩ past_xMC.Intv, lbd, ubd, sparsity, subgrad_tol)

    elseif is_post && is_intersect && !interval_intersect
        return set_value_post(x, xMC ∩ past_xMC, lbd, ubd, sparsity, subgrad_tol)

    elseif is_post && !is_intersect
        return set_value_post(x, xMC, lbd, ubd, sparsity, subgrad_tol)

    elseif !is_post && is_intersect && interval_intersect
        return xMC ∩ past_xMC.Intv

    elseif !is_post && is_intersect && !interval_intersect
        return xMC ∩ past_xMC

    end
    return xMC
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution of node representing `n = x + y`.
"""
function forward_plus_binary!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                              numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                              x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                              sparsity::Vector{Int},
                              is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool) where {N,T<:RelaxTag}

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index =  children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index = children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    output_is_number = arg1_is_number && arg2_is_number

    # a + b
    if output_is_number
         numberstorage[k] = num1 + num2

    # x + b
    elseif !arg1_is_number && arg2_is_number
        outset = set1 + num2
        # is_first_eval ? (set1 + num2) : plus_kernel(set1, num2, setstorage[k].Intv)

    # a + y
    elseif arg1_is_number && !arg2_is_number
        outset = num1 + set2
        # is_first_eval ? (num1 + set2) : plus_kernel(num1, set2, setstorage[k].Intv)

    # x + y
    else
        outset = set1 + set2
        # is_first_eval ? (set1 + set2) : plus_kernel(set1, set2, setstorage[k].Intv)

    end

     numvalued[k] = output_is_number
    if !output_is_number
         setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, false,
                                                         is_intersect, interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution of node representing `n = +(x, y, z,...)`.
"""
function forward_plus_narity!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                              numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                              x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64},
                              subgrad_tol::Float64, sparsity::Vector{Int},
                              is_post::Bool, is_intersect::Bool, interval_intersect::Bool) where {N,T<:RelaxTag}

    # get row indices
    idx = first(children_idx)

    # extract values for argument 1
    arg_index =  children_arr[idx]
    output_is_number = numvalued[arg_index]
    if output_is_number
       tmp_set = zero(MC{N,T})
       tmp_num = numberstorage[arg_index]
    else
       tmp_num = 0.0
       tmp_set = setstorage[arg_index]
    end
    output_is_number = true

    for idx = 2:length(children_idx)
        cidx = children_idx[idx]
        arg_index =  children_arr[cidx]
        arg_is_number = numvalued[arg_index]
        if arg_is_number
            tmp_num += numberstorage[arg_index]
        else
            tmp_set += setstorage[arg_index]
        end
        output_is_number &= arg_is_number
    end

     numvalued[k] = output_is_number
    if output_is_number
         numberstorage[k] = tmp_num
    else
        tmp_set += tmp_num
        setstorage[k] = overwrite_or_intersect(tmp_set, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, false, is_intersect,
                                               interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = x*y`.
"""
function forward_multiply_binary!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                                  numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                                  x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                                  sparsity::Vector{Int},
                                  is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool) where {N,T<:RelaxTag}
    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = one(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 1.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index =  children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = one(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 1.0
    end

    output_is_number = arg1_is_number && arg2_is_number

    # a * b
    if output_is_number
         numberstorage[k] = num1 * num2

    # x * b
    elseif !arg1_is_number && arg2_is_number
        outset = set1*num2 #is_first_eval ? (set1 * num2) : mult_kernel(set1, num2, setstorage[k].Intv)

    # a * y
    elseif arg1_is_number && !arg2_is_number
        outset = num1*set2 #is_first_eval ? (num1 * set2) : mult_kernel(set2, num1, setstorage[k].Intv)

    # x * y
    else
        outset = set1*set2 #is_first_eval ? (set1 * set2) : mult_kernel(set1, set2, setstorage[k].Intv)

    end

    numvalued[k] = output_is_number
    if !output_is_number
         setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                                         interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution of node representing `n = *(x, y, z,...)`.
"""
function forward_multiply_narity!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                                  numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                                  x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                                  sparsity::Vector{Int},
                                  is_post::Bool, is_intersect::Bool, interval_intersect::Bool) where {N,T<:RelaxTag}
    # get row indices
    idx = first(children_idx)

    # extract values for argument 1
    arg_index =  children_arr[idx]
    output_is_number =  numvalued[arg_index]
    if output_is_number
        tmp_set = 1.0#one(MC{N,T})
        tmp_num = numberstorage[arg_index]
    else
        tmp_num = 1.0
        tmp_set = setstorage[arg_index]
    end

    for idx = 2:length(children_idx)
        cidx =  children_idx[idx]
        arg_index_t =  children_arr[cidx]
        arg_is_number_t = numvalued[arg_index_t]
        if arg_is_number_t
            tmp_num = tmp_num*numberstorage[arg_index_t]
        else
            tmp_set = tmp_set*setstorage[arg_index_t]
        end
        output_is_number &= arg_is_number_t
    end

     numvalued[k] = output_is_number
    if output_is_number
         numberstorage[k] = tmp_num
    else
       tmp_set *= tmp_num
       setstorage[k] = overwrite_or_intersect(tmp_set, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                              interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = x-y`.
"""
function forward_minus!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                        numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                        x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                        sparsity::Vector{Int},
                        is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool) where {N,T<:RelaxTag}

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index = children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 =  setstorage[arg2_index]
        num2 = 0.0
    end

    output_is_number = arg1_is_number && arg2_is_number

    # a - b
    if output_is_number
         numberstorage[k] = num1 - num2

    # x - b
    elseif !arg1_is_number && arg2_is_number
        outset = set1 - num2 #is_first_eval ? (set1 - num2) : minus_kernel(set1, num2, setstorage[k].Intv)

    # a - y
    elseif arg1_is_number && !arg2_is_number
        outset = num1 - set2 #is_first_eval ? (num1 - set2) : minus_kernel(num1, set2, setstorage[k].Intv)

    # x - y
    else
        outset = set1 - set2 #is_first_eval ? (set1 - set2) : minus_kernel(set1, set2, setstorage[k].Intv)

    end

    numvalued[k] = output_is_number
    if !output_is_number
         setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                                         interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = x^y`.
"""
function forward_power!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                        numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                        x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64},
                        subgrad_tol::Float64, sparsity::Vector{Int},
                        is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool,
                        ctx::GuardCtx) where {N,T<:RelaxTag}

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 =  numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 =  setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index =  children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    # is output a number (by closure of the reals)
    output_is_number = arg1_is_number && arg2_is_number
    numvalued[k] = output_is_number

    # x^1 = x
    if arg2_is_number && (num2 == 1.0)
        if arg1_is_number
             numberstorage[k] = num1
        else
             setstorage[k] = set1
        end
        return nothing

    # x^0 = 1
    elseif arg2_is_number && (num2 == 0.0)
        if arg1_is_number
             numberstorage[k] = 1.0
        else
             setstorage[k] = zero(MC{N,T})
        end
        return nothing

    else
        # a^b
        if arg1_is_number && arg2_is_number
             numberstorage[k] = num1^num2

        # x^b
        elseif !arg1_is_number && arg2_is_number
            outset = set1^num2
            #is_first_eval ? pow(set1, num2) : ^(set1, num2, setstorage[k].Intv)

        # a^y
        elseif arg1_is_number  && !arg2_is_number
            guard_on = ctx.metadata.guard_on
            outset = guard_on ? overdub(ctx, ^, num1, set2) : num1^set2 # overdub(ctx, pow, num1, set2)
            #is_first_eval ? overdub(ctx, pow, num1, set2) : overdub(ctx, ^, num1, set2, setstorage[k].Intv)

        # x^y
        elseif !arg1_is_number && !arg2_is_number
            guard_on = ctx.metadata.guard_on
            outset = guard_on ? overdub(ctx, ^, set1, set2) : set1^set2
            #is_first_eval ? overdub(ctx, pow, set1, set2) : overdub(ctx, ^, set1, set2, setstorage[k].Intv)

        end
    end

    if !output_is_number
        setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                               interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = x/y`.
"""
function forward_divide!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                         numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                         x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                         sparsity::Vector{Int},
                         is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool,
                         ctx::GuardCtx) where {N,T<:RelaxTag}

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index = children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    # is output a number (by closure of the reals)?
    output_is_number = arg1_is_number && arg2_is_number
    numvalued[k] = output_is_number

    # a/b
    if output_is_number
         numberstorage[k] = num1/num2

    # x/b
    elseif !arg1_is_number && arg2_is_number

        guard_on = ctx.metadata.guard_on
        if guard_on
            outset = Cassette.overdub(ctx, /, set1, num2)
        else
            outset = set1/num2
        end
        # is_first_eval ? set1/num2 : div_kernel(set1, num2, setstorage[k].Intv)

    # a/y
    elseif arg1_is_number && !arg2_is_number

        guard_on = ctx.metadata.guard_on
        if guard_on
            outset = Cassette.overdub(ctx, /, num1, set2)
        else
            outset = num1/set2
        end
        # is_first_eval ? num1/set2 : div_kernel(num1, set2, setstorage[k].Intv)

    # x/y
    else

        guard_on = ctx.metadata.guard_on
        if guard_on
            outset = Cassette.overdub(ctx, /, set1, set2)
        else
            outset = set1/set2
        end
        # is_first_eval ? set1/set2 : div_kernel(set1, set2, setstorage[k].Intv)

    end

    if !output_is_number
        setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post,
                                               is_intersect, interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = user_f(x, y...)`.
"""
function forward_user_multivariate!(k::Int64, op::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                                    numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                                    x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                                    sparsity::Vector{Int},
                                    is_post::Bool, is_intersect::Bool, interval_intersect::Bool, ctx::GuardCtx,
                                    user_operators::JuMP._Derivatives.UserOperatorRegistry,
                                    num_mv_buffer::Vector{Float64}) where {N, T<:RelaxTag}

    n = length(children_idx)
    evaluator = user_operators.multivariate_operator_evaluator[op - JuMP._Derivatives.USER_OPERATOR_ID_START + 1]
    set_input = zeros(MC{N,T}, n)
    num_input = view(num_mv_buffer, 1:n)
    fill!(num_input, -Inf)

    buffer_count = 1
    output_is_number = true
    for c_idx in children_idx
        arg_index = children_arr[c_idx]
        arg_is_number = numvalued[arg_index]
        if arg_is_number
            num_input[buffer_count] = numberstorage[arg_index]
        else
            set_input[buffer_count] = setstorage[arg_index]
        end
        output_is_number &= arg_is_number
        buffer_count += 1
    end

    if output_is_number
        numberstorage[k] = MOI.eval_objective(evaluator, num_input)
    else
        for i = 1:(buffer_count - 1)
            if !isinf(num_input[i])
                 set_input[buffer_count] = MC{N,T}(num_input[buffer_count])
            end
        end
        guard_on = ctx.metadata.guard_on
        if guard_on
            outset = Cassette.overdub(ctx, MOI.eval_objective, evaluator, set_input)
        else
            outset = MOI.eval_objective(evaluator, set_input)
        end
        setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                               interval_intersect)
    end
    numvalued[k] = output_is_number

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = f(c)` where f is standard function
and `c` is a number.
"""
function forward_univariate_number!(k::Int64, op::Int64, arg_idx::Int, numvalued::Vector{Bool}, numberstorage::Vector{Float64})

    tmp_num = numberstorage[arg_idx]
    outnum = eval_univariate_set(op, tmp_num)

    numberstorage[k] = outnum
    numvalued[k] = true

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = f(x)` where f is standard function
that requires a single tiepoint calculation per convex/concave relaxation (e.g. tan).
"""
function forward_univariate_tiepnt_1!(k::Int64, op::Int64, child_idx::Int64, setstorage::Vector{V},
                                      x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64,
                                      sparsity::Vector{Int},
                                      tpdict::Dict{Int64, Tuple{Int64,Int64,Int64,Int64}},
                                      tp1storage::Vector{Float64}, tp2storage::Vector{Float64},
                                      is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool, ctx::GuardCtx) where V

    tmp_set = setstorage[child_idx]

    tidx1, tidx2 = tpdict[k]
    tp1 =  tp1storage[tidx1]
    tp2 =  tp2storage[tidx1]
    new_tie_points = tp1 === Inf

    guard_on = ctx.metadata.guard_on
    if guard_on
        outset, tp1, tp2 = Cassette.overdub(ctx, single_tp_set, op, tmp_set, setstorage[k], tp1, tp2, is_first_eval)
    else
        outset, tp1, tp2 = single_tp_set(op, tmp_set, setstorage[k], tp1, tp2, is_first_eval)
    end

    if new_tie_points
         tp1storage[tidx1] = tp1
         tp2storage[tidx1] = tp2
    end

    setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                           interval_intersect)
    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = f(x)` where f is standard function
that requires a two tiepoint calculations per convex/concave relaxation (e.g. sin).
"""
function forward_univariate_tiepnt_2!(k::Int64, op::Int64, child_idx::Int64, setstorage::Vector{V},
                                      x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64},
                                      subgrad_tol::Float64, sparsity::Vector{Int},
                                      tpdict::Dict{Int64, Tuple{Int64,Int64,Int64,Int64}},
                                      tp1storage::Vector{Float64}, tp2storage::Vector{Float64},
                                      tp3storage::Vector{Float64}, tp4storage::Vector{Float64},
                                      is_post::Bool, is_intersect::Bool, is_first_eval::Bool,
                                      interval_intersect::Bool, ctx::GuardCtx) where V

    tmp_set = setstorage[child_idx]

    # retreive previously calculated tie-points
    # These are re-initialize to Inf for each box
    tidx1, tidx2 = tpdict[k]
    tp1 = tp1storage[tidx1]
    tp2 = tp2storage[tidx1]
    tp3 = tp3storage[tidx2]
    tp4 = tp4storage[tidx2]

    new_tie_points = tp1 === Inf

    # Perform an evaluation of the univariate function overdubbed with Cassette.jl
    guard_on = ctx.metadata.guard_on
    if guard_on
        outset, tp1, tp2, tp3, tp4 = Cassette.overdub(ctx, double_tp_set, op, tmp_set, setstorage[k],
                                                      tp1, tp2, tp3, tp4, is_first_eval)
    else
        outset, tp1, tp2, tp3, tp4 = double_tp_set(op, tmp_set, setstorage[k], tp1, tp2, tp3, tp4, is_first_eval)
    end
    #Cassette.overdub(ctx, double_tp_set, op, tmp_set, setstorage[k], tp1, tp2, tp3, tp4, is_first_eval)

    # Store new tiepoints if new evaluation
    if new_tie_points
         tp1storage[tidx1] = tp1
         tp2storage[tidx1] = tp2
         tp3storage[tidx2] = tp3
         tp4storage[tidx2] = tp4
    end

    setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                           interval_intersect)
    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = user_f(x)`.
"""
function forward_univariate_user!(k::Int64, op::Int64, child_idx::Int64, arg_is_number::Bool,
                                  setstorage::Vector{V}, x::Vector{Float64}, lbd::Vector{Float64},
                                  ubd::Vector{Float64}, subgrad_tol::Float64, sparsity::Vector{Int},
                                  is_post::Bool, is_intersect::Bool, is_first_eval::Bool,
                                  interval_intersect::Bool, ctx::GuardCtx, user_operators) where V

    userop = op - JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START + 1
    f = user_operators.univariate_operator_f[userop]

    if arg_is_number
        tmp_num = setstorage[child_idx]
        outnum = f(tmp_num)
        numberstorage[k] = outnum

    else
        tmp_set = setstorage[child_idx]

        guard_on = ctx.metadata.guard_on
        if guard_on
            outset = Cassette.overdub(ctx, f, tmp_set)
        else
            outset = f(tmp_set)
        end

        setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity,
                                               is_post, is_intersect, interval_intersect)
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Updates storage tapes with forward evalution for node representing `n = f(x)` where `f` is a standard function
that does not require a tiepoint evaluation (e.g. exp).
"""
function forward_univariate_other!(k::Int64, op::Int64, child_idx::Int64, setstorage::Vector{V},
                                   x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64},
                                   subgrad_tol::Float64, sparsity::Vector{Int}, is_post::Bool,
                                   is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool,
                                   ctx::GuardCtx) where V

    tmp_set = setstorage[child_idx]

    guard_on = ctx.metadata.guard_on
    if guard_on
        outset = Cassette.overdub(ctx, eval_univariate_set, op, tmp_set)
    else
        outset = eval_univariate_set(op, tmp_set)
    end
    setstorage[k] = overwrite_or_intersect(outset, setstorage[k], x, lbd, ubd, subgrad_tol, sparsity,
                                           is_post, is_intersect, interval_intersect)

    return nothing
end

"""
$(FUNCTIONNAME)

"""
function expand_set(::Type{MC{N2,T}}, x::MC{N1,T}, fsparse::Vector{Int64},
                    subsparse::Vector{Int64}, cv_buffer::Vector{Float64},
                    cc_buffer::Vector{Float64}) where {N1 ,N2, T<:RelaxTag}

    # TODO: Perform this in a manner that does not allocate via an ntuple
    # operator or via generated code
    cvg = x.cv_grad
    ccg = x.cc_grad
    xcount = 1
    xcurrent = subsparse[1]
    for i = 1:N2
        if fsparse[i] === xcurrent
            cv_buffer[i] = cvg[xcount]
            cc_buffer[i] = ccg[xcount]
            xcount += 1
            if xcount <= N1
                xcurrent = subsparse[xcount]
            else
                break
            end
        else
            cv_buffer[i] = 0.0
        end
    end
    cv_grad = SVector{N2,Float64}(cv_buffer)
    cc_grad = SVector{N2,Float64}(cc_buffer)
    return MC{N2,T}(x.cv, x.cc, x.Intv, cv_grad, cc_grad, x.cnst)
end

"""
$(FUNCTIONNAME)

Unpacks the `MC{N1,T}` in the subexpression to a `MC{N2,T}` where `N1` is the sparsity of the
subexpression and `N2` is the sparsity of the function tape. Note that the sparsity of the
subexpression is less than the sparsity of the function itself. This limits the storage
required by tapes but prevents reverse mode subgradient propagation through expressions
[This may be subject to change in the future once the reverse mode propagation becomes
more robust].
"""
function forward_get_subexpression!(k::Int64, active_subexpr::NonlinearExpression{MC{N,T}},
                                    setstorage::Vector{MC{Q,T}}, numberstorage::Vector{Float64},
                                    numvalued::Vector{Bool}, fsparsity::Vector{Int64},
                                    cv_buffer::Vector{Float64}, cc_buffer::Vector{Float64}) where {N,Q,T<:RelaxTag}

    is_number = active_subexpr.isnumber[1]
    if is_number
        numberstorage[k] = active_subexpr.setstorage[1]
    else
        set_arg = active_subexpr.setstorage[1]
        setstorage[k] = expand_set(MC{Q,T}, set_arg, fsparsity, active_subexpr.grad_sparsity, cv_buffer, cc_buffer)
    end
    numvalued[k] = is_number

    return nothing
end

const id_to_operator = Dict(value => key for (key, value) in JuMP.univariate_operator_to_id)

"""
$(FUNCTIONNAME)

Performs a forward pass using the tape information passed as arguments. Each variety of node calls an associated
forward_xxx function where xxx is a descriptor.
"""
function forward_pass_kernel!(nd::Vector{JuMP.NodeData}, adj::SparseMatrixCSC{Bool,Int64}, x::Vector{Float64},
                              lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                              setstorage::Vector{MC{N,T}}, numberstorage::Vector{Float64}, numvalued::Vector{Bool},
                              tpdict::Dict{Int64,Tuple{Int64,Int64,Int64,Int64}}, tp1storage::Vector{Float64},
                              tp2storage::Vector{Float64}, tp3storage::Vector{Float64}, tp4storage::Vector{Float64},
                              user_operators::JuMP._Derivatives.UserOperatorRegistry, subexpressions::Vector{NonlinearExpression},
                              func_sparsity::Vector{Int64}, reverse_sparsity::Vector{Int64},
                              num_mv_buffer::Vector{Float64}, ctx::GuardCtx,
                              is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool,
                              cv_buffer::Vector{Float64}, cc_buffer::Vector{Float64},
                              treat_x_as_number::Vector{Bool}, subgrad_tol::Float64) where {N, T<:RelaxTag}

    children_arr = rowvals(adj)

    FORWARD_DEBUG && println(" ")
    for k = length(nd):-1:1

        oldset = setstorage[k]
        nod = nd[k]
        op = nod.index

        if nod.nodetype == JuMP._Derivatives.VALUE
            numvalued[k] = true
            FORWARD_DEBUG && println("value[$op]    at k = $k -> $(numberstorage[k])")

        elseif nod.nodetype == JuMP._Derivatives.PARAMETER
            numvalued[k] = true
            FORWARD_DEBUG && println("parameter[$op] at k = $k -> $(numberstorage[k])")

        elseif nod.nodetype == JuMP._Derivatives.VARIABLE
            isa_number =  treat_x_as_number[op]
            numvalued[k] = isa_number
            xval = x[op]
            if isa_number
                 numberstorage[k] = xval
            else
                seed_index = reverse_sparsity[op]
                seed_grad = seed_gradient(seed_index, Val{N}())
                lower_bnd_val = lbd[op]
                upper_bnd_val = ubd[op]
                xcv_eps = xval - (xval - lower_bnd_val)*1E-8
                xcc_eps = xval + (upper_bnd_val - xval)*1E-8
                xMC = MC{N,T}(xcv_eps, xcc_eps, Interval{Float64}(lower_bnd_val, upper_bnd_val),
                              seed_grad, seed_grad, false)
                setstorage[k] = is_first_eval ? xMC : (xMC ∩ setstorage[k].Intv)
            end
            FORWARD_DEBUG && println("variable[$op] at k = $k -> $(setstorage[k])")
        elseif nod.nodetype == JuMP._Derivatives.SUBEXPRESSION
            active_subexpr = subexpressions[op]
            forward_get_subexpression!(k, active_subexpr, setstorage, numberstorage, numvalued, func_sparsity,
                                       cv_buffer, cc_buffer)

        elseif nod.nodetype == JuMP._Derivatives.CALL

            children_idx = nzrange(adj, k)
            n_children = length(children_idx)

            # :+ with arity two or greater
            if op === 1
                n = length(children_idx)
                if n === 2
                    forward_plus_binary!(k, children_arr, children_idx, numvalued, numberstorage,
                                         setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect, is_first_eval,
                                         interval_intersect)
                else
                    forward_plus_narity!(k, children_arr, children_idx, numvalued, numberstorage,
                                         setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect, interval_intersect)
                end
                FORWARD_DEBUG && println("plus[$n]     at k = $k -> $(setstorage[k])")

            # :- with arity two
            elseif op === 2
                forward_minus!(k, children_arr, children_idx, numvalued, numberstorage,
                               setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect, is_first_eval,
                               interval_intersect)
                FORWARD_DEBUG && println("minus        at k = $k -> $(setstorage[k])")

            # :* with arity two or greater
            elseif op === 3
                n = length(children_idx)
                if n === 2
                    forward_multiply_binary!(k, children_arr, children_idx, numvalued,
                                             numberstorage, setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post,
                                             is_intersect, is_first_eval, interval_intersect)
                else
                    forward_multiply_narity!(k, children_arr, children_idx, numvalued,
                                             numberstorage, setstorage, x, lbd, ubd, subgrad_tol,
                                             sparsity, is_post, is_intersect, interval_intersect)
                end
                FORWARD_DEBUG && println("mult[$n]     at k = $k -> $(setstorage[k])")

            # :^
            elseif op === 4
                forward_power!(k, children_arr, children_idx, numvalued, numberstorage,
                               setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect, is_first_eval,
                               interval_intersect, ctx)

                FORWARD_DEBUG && println("power       at k = $k -> $(setstorage[k])")

            # :/
            elseif op === 5
                forward_divide!(k, children_arr, children_idx, numvalued, numberstorage,
                                setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect, is_first_eval,
                                interval_intersect, ctx)

                FORWARD_DEBUG && println("divide      at k = $k -> $(setstorage[k])")

            # user multivariate function
            elseif op >= JuMP._Derivatives.USER_OPERATOR_ID_START
                forward_user_multivariate!(k, op, children_arr, children_idx, numvalued, numberstorage,
                                           setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post, is_intersect,
                                           interval_intersect, ctx, user_operators, num_mv_buffer)
                FORWARD_DEBUG && println("user_mult   at k = $k -> $(setstorage[k])")

            else
               error("Unsupported operation $(operators[op])")
            end

        elseif nod.nodetype == JuMP._Derivatives.CALLUNIVAR

            # checks to see if operator is a number
            child_idx = first(nzrange(adj, k))
            arg_idx = children_arr[adj.colptr[k]]
            arg_is_number =  numvalued[arg_idx]
            numvalued[k] = arg_is_number

            # performs univariate operators on number valued inputs
            if op >= JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START
                forward_univariate_user!(k, op, arg_idx, arg_is_number, setstorage, x, lbd, ubd, subgrad_tol,
                                         sparsity, is_post, is_intersect, is_first_eval, interval_intersect, ctx,
                                         user_operators)

            elseif arg_is_number
                forward_univariate_number!(k, op, arg_idx, numvalued, numberstorage)

            # performs set valued operators that require a single tiepoint calculation
            elseif single_tp(op)
                forward_univariate_tiepnt_1!(k, op, arg_idx, setstorage, x, lbd, ubd,  subgrad_tol, sparsity, tpdict,
                                             tp1storage, tp2storage, is_post, is_intersect,
                                             is_first_eval, interval_intersect, ctx)

            # performs set valued operators that require two tiepoint calculations
            elseif double_tp(op)
                forward_univariate_tiepnt_2!(k, op, arg_idx, setstorage, x, lbd, ubd,  subgrad_tol, sparsity, tpdict,
                                             tp1storage, tp2storage, tp3storage, tp4storage,
                                             is_post, is_intersect, is_first_eval,
                                             interval_intersect, ctx)

            # performs set valued operator on other functions in base library
            else
                forward_univariate_other!(k, op, arg_idx, setstorage, x, lbd, ubd, subgrad_tol, sparsity, is_post,
                                          is_intersect, is_first_eval, interval_intersect, ctx)

            end
            FORWARD_DEBUG && println("fop[$op]   at k = $k -> $(setstorage[k])")
        else
            error("Unrecognized node type $(nod.nodetype).")

        end

    end

    return nothing
end

###
### Define forward evaluation pass
###
function forward_pass!(evaluator::Evaluator, d::NonlinearExpression{V}) where V
    # check that prior subexpressions have been evaluated
    # i.e. box_id is same and reference point is the same
    for i = 1:d.dependent_subexpression_count
        if !prior_eval(evaluator, i)
            subexpr = evaluator.subexpressions[i]
            forward_pass!(evaluator, subexpr)
        end
    end

    forward_pass_kernel!(d.nd, d.adj, evaluator.x, evaluator.lower_variable_bounds,
                         evaluator.upper_variable_bounds, d.grad_sparsity,
                         d.setstorage,
                         d.numberstorage, d.isnumber, d.tpdict,
                         d.tp1storage, d.tp2storage, d.tp3storage, d.tp4storage,
                         evaluator.user_operators, evaluator.subexpressions,
                         d.grad_sparsity, d.reverse_sparsity,
                         evaluator.num_mv_buffer, evaluator.ctx,
                         evaluator.is_post, evaluator.is_intersect,
                         evaluator.is_first_eval, evaluator.interval_intersect,
                         d.cv_grad_buffer, d.cc_grad_buffer,
                         evaluator.treat_x_as_number, evaluator.subgrad_tol)
    return nothing
end

function forward_pass!(evaluator::Evaluator, d::BufferedNonlinearFunction{V}) where V

    forward_pass!(evaluator, d.expr)
    d.has_value = true
    d.last_past_reverse = false

    return nothing
end


text:
"""
Add numbers.
"""
function add_nums end

"""
    add_nums(a,b)

Add two numbers together.
"""
function add_nums(a,b)
  a+b
end

"""
    add_nums(a,b,c)

Add three numbers together.
"""
function add_nums(a,b,c)
  a+b+c
end


text:
#Default board size if no argument given
#This also makes it convenient to run inside an IDE
N = 5

if length(ARGS) >= 1
  N = parse(Int64,ARGS[1])
end

function displaysquarebits(board)
  board_size = size(board,1)
  for i = 1:board_size
    for j = 1:board_size
      if board[i,j]
        print("1")
      else
        print("0")
      end
    end
    println()
  end
end

# Create all the square sizes
function create_square_indexes(n)
  if n < 1
    error("The board must have at least one square")
  end
  if n == 1
    return Array{NTuple{8,Int64},1}()
  end

  square_indexes = Array{NTuple{8,Int64},1}()
  # Generate "flat" square bitmasks
  for square_size::UInt64=1:n-1
    for y=1:n-square_size
      for x=1:n-square_size
        push!(square_indexes,(x,y,
                              x+square_size,y,
                              x,y+square_size,
                              x+square_size,y+square_size))
      end
    end
  end

  #Generate "diamond" square bitmasks
  for square_size::UInt64=1:n÷2
    for y=1:n-square_size*2
      for x=square_size+1:n-square_size
        push!(square_indexes,(x,y,
                              x+square_size,y+square_size,
                              x-square_size,y+square_size,
                              x,y+2*square_size))
      end
    end
  end

  return square_indexes
end

function has_square(board, square_indexes)
  for square_index in square_indexes
    first_corner::Bool = board[square_index[1],square_index[2]]
    if board[square_index[3],square_index[4]] == first_corner &&
       board[square_index[5],square_index[6]] == first_corner &&
       board[square_index[7],square_index[8]] == first_corner
       return true
     end
  end
  return false
end

println("START")
println()

#TODO: implement this
#iterate_valid_boards(N)

function has_square_20_million_times(board, square_indexes)
  x = false
  for i=1:20000000
    x = has_square(board, square_indexes)
  end
  return x
end

function time_has_square(board, n)
  square_indexes = create_square_indexes(n)
  @timev has_square(board, square_indexes)
  @timev has_square(board, square_indexes)
  @timev has_square(board, square_indexes)
  @timev has_square(board, square_indexes)
  @timev has_square(board, square_indexes)
  @timev has_square_20_million_times(board, square_indexes)
  #@code_native has_square(board, square_indexes)
end

solution_5 = [false false false false false;
              false true  false true  true;
              false true  true  true  false;
              true  true  false false true;
              true  false true  true  true]

solution_7 = [true  true  true  true  false true  true;
              false true  false false false false false;
              false false false true  true  false true;
              true  true  false false true  true  false;
              false false true  true  true  false false;
              true  false true  false true  false true;
              false true  true  false false true  true]

time_has_square(solution_5,5)
time_has_square(solution_7,7)

#Verify correctness
#square_indexes_5 = create_square_indexes(5)
#println(has_square(solution_5, square_indexes_5))
#square_indexes_7 = create_square_indexes(7)
#println(has_square(solution_7, square_indexes_7))


text:
import Pkg;
Pkg.add("FastGaussQuadrature")
Pkg.add("Plots")
Pkg.add("LaTeXStrings")
Pkg.add("PyPlot")
Pkg.add("WriteVTK")


text:
"""
	gmtwhich(cmd0::String; kwargs...)

Find full path to specified files

Full option list at [`gmtwhich`]($(GMTdoc)gmtwhich.html)

Parameters
----------

- **A** | **with_permissions** :: [Type => Bool]

	Only consider files that the user has permission to read [Default consider all files found].
    ($(GMTdoc)gmtwhich.html#a)
- **C** | **confirm** :: [Type => Bool]

	Instead of reporting the paths, print the confirmation Y if the file is found and N if it is not.
    ($(GMTdoc)gmtwhich.html#c)
- **D** | **report_dir** :: [Type => Bool]

	Instead of reporting the paths, print the directories that contains the files.
    ($(GMTdoc)gmtwhich.html#d)
- **G** | **download** :: [Type => Str | []]      ``Arg = [c|l|u]``

	If a file argument is a downloadable file (either a full URL, a @file for downloading from
	the GMT Site Cache, or @earth_relief_*.grd) we will try to download the file if it is not
	found in your local data or cache dirs.
    ($(GMTdoc)gmtwhich.html#g)
- $(GMT.opt_V)
"""
function gmtwhich(cmd0::String; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmtwhich", cmd0)

	d = KW(kwargs)
	cmd = parse_V_params("", d)
    cmd = parse_these_opts(cmd, d, [[:A :with_permissions], [:C :confirm], [:D :report_dir], [:G :download]])

	common_grd(d, cmd0, cmd, "gmtwhich ", nothing)		# Finish build cmd and run it
end


text:
import TriangleMesh, UnicodePlots, Plots
using Printf, LoopVectorization, LinearAlgebra, SparseArrays

include("CreateMeshFCFV.jl")
include("VisuFCFV.jl")
include("DiscretisationFCFV.jl")
include("DiscretisationFCFV_Poisson_o2.jl")
include("SolversFCFV_Poisson.jl")

#--------------------------------------------------------------------#

function SetUpProblem!( mesh, T, Tdir, Tneu, se, a, b, c, d, alp, bet )
    # Evaluate T analytic on cell faces
    @avx for in=1:mesh.nf
        x        = mesh.xf[in]
        y        = mesh.yf[in]
        Tdir[in] = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        dTdx     = Tdir[in] * (a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y))
        dTdy     = Tdir[in] * (alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y))
        Tneu[in] = -dTdy # nx*dTdx + nt*dTdy on SOUTH face
    end
    # Evaluate T analytic on barycentres
    @avx for iel=1:mesh.nel
        x       = mesh.xc[iel]
        y       = mesh.yc[iel]
        T       = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        T[iel]  = T
        se[iel] = T*(-a*alp*cos(a*x + c*y) + b*bet*sin(b*x + d*y))*(a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y)) + T*(a^2*alp*sin(a*x + c*y) + b^2*bet*cos(b*x + d*y)) + T*(-alp*c*cos(a*x + c*y) + bet*d*sin(b*x + d*y))*(alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y)) + T*(alp*c^2*sin(a*x + c*y) + bet*d^2*cos(b*x + d*y))
    end
    return
end

#--------------------------------------------------------------------#

function ComputeError( mesh, Te, qx, qy, a, b, c, d, alp, bet )
    eT  = zeros(mesh.nel)
    eqx = zeros(mesh.nel)
    eqy = zeros(mesh.nel)
    Ta  = zeros(mesh.nel)
    qxa = zeros(mesh.nel)
    qya = zeros(mesh.nel)
    @avx for iel=1:mesh.nel
        x        = mesh.xc[iel]
        y        = mesh.yc[iel]
        Ta[iel]  = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        qxa[iel] = -Ta[iel] * (a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y))
        qya[iel] = -Ta[iel] * (alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y))
        eT[iel]  = Te[iel] - Ta[iel]
        eqx[iel] = qx[iel] - qxa[iel]
        eqy[iel] = qy[iel] - qya[iel]
    end
    errT  = norm(eT)/norm(Ta)
    errqx = norm(eqx)/norm(qxa)
    errqy = norm(eqy)/norm(qya)
    return errT, errqx, errqy
end

#--------------------------------------------------------------------#

function StabParam(tau, dA, Vol, mesh_type)
    if mesh_type=="Quadrangles";        taui = tau;    end
    # if mesh_type=="UnstructTriangles";  taui = tau*dA; end
    if mesh_type=="UnstructTriangles";  taui = tau end
    return taui
end

#--------------------------------------------------------------------#

@views function main( n )

    println("\n******** FCFV POISSON ********")

    # Create sides of mesh
    xmin, xmax = 0, 1
    ymin, ymax = 0, 1
    nx, ny     = n*8, n*8
    R          = 0.5
    inclusion  = 0
    solver     = 0
    o2         = 1
    BC         = [1; 1; 1; 1] # S E N W --- 1: Dirichlet / 2: Neumann
    # mesh_type  = "Quadrangles"
    mesh_type  = "UnstructTriangles"

    # Generate mesh
    if mesh_type=="Quadrangles"
        if o2==0 tau  = 1e0 end
        if o2==1 tau  = 1e4 end
        mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC )
    elseif mesh_type=="UnstructTriangles"
        if o2==0 tau  = 1e2 end
        if o2==1 tau  = 1e6 end
        mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC )
    end
    println("Number of elements: ", mesh.nel)

    # Source term and BCs etc...
    Tanal  = zeros(mesh.nel)
    se     = zeros(mesh.nel)
    Tdir   = zeros(mesh.nf)
    Tneu   = zeros(mesh.nf)
    alp = 0.1; bet = 0.3; a = 5.1; b = 4.3; c = -6.2; d = 3.4;
    println("Model configuration :")
    @time SetUpProblem!(mesh , Tanal, Tdir, Tneu, se, a, b, c, d, alp, bet)

    # Compute some mesh vectors
    println("Compute FCFV vectors:")
    @time ae, be, be_o2, ze, pe, mei, pe, rj  = ComputeFCFV_o2(mesh, se, Tdir, tau, o2)

    # Assemble element matrices and RHS
    println("Compute element matrices:")
    @time Kv, fv = ElementAssemblyLoop_o2(mesh, ae, be, be_o2, ze, mei, pe, rj, Tdir, Tneu, tau, o2)

    # Assemble triplets and sparse
    println("Assemble triplets and sparse:")
    @time K, f = CreateTripletsSparse(mesh, Kv, fv)

    # Solve
    println("Solve:")
    @time Th = SolvePoisson(mesh, K, f, solver)

    # Reconstruct element values
    println("Compute element values:")
    @time Te, qx, qy = ComputeElementValues_o2(mesh, Th, ae, be, be_o2, ze, rj, mei, Tdir, tau, o2)

    # Compute discretisation errors
    err_T, err_qx, err_qy = ComputeError( mesh, Te, qx, qy, a, b, c, d, alp, bet )
    println("Error in T:  ", err_T )
    println("Error in qx: ", err_qx)
    println("Error in qy: ", err_qy)

    # # Visualise
    # println("Visualisation:")
    # @time PlotMakie( mesh, Te, xmin, xmax, ymin, ymax, :viridis )

end

n = 4
main( n )
# main( n )
# main( n )


text:
"""
`PyPinYin` package to interface with Python's `pypinyin` through `PyCall`.
On loading, three functions are provided as generic functions:
    + pinyin
    + lazypinyin
    + hanzi2pinyin
To find documentation for those functions, one should go into
Julia's `?` REPL mode. README.md also provides many examples.
"""
module PyPinYin

include("./styles.jl")

export pinyin
export lazypinyin
export printpinyin
export hanzi2pinyin

using PyCall

function __init__()
    py"""
    from pypinyin import lazy_pinyin as lazypy
    from pypinyin import pinyin as pinyins
    """
end

"""
    lazypinyin(hanzi::AbstractString;    # 汉字串
               style::Int=NORMAL,        # 声调风格
               errors::String="default", # 如何处理没有拼音的字符(default/ignore/replace)
               strict::Bool=true,        # 是否严格遵照《汉语拼音方案》处理声母和韵母
               v2u::Bool=false,          # 是否用 `ü` 代替 `v`
               tone5::Bool=false,        # 是否用 5 标识轻声
               sandhi::Bool=false)       # 是否变调，nǐ hǎo -> ní hǎo
"""
function lazypinyin(hanzi::AbstractString;
                    style::Int=NORMAL,
                    errors::String="default",
                    strict::Bool=true,
                    v2u::Bool=false,
                    tone5::Bool=false,
                    sandhi::Bool=false)

    return py"lazypy"(hanzi,
                      style=style,
                      errors=errors,
                      strict=strict,
                      v_to_u=v2u,
                      neutral_tone_with_five=tone5,
                      tone_sandhi=sandhi)
end

"""
    pinyin(hanzi::AbstractString;   # 汉字串
           style::Int=TONE,         # 声调风格
           heteronym::Bool=false,   # 多音字
           errors::String="default" # 如何处理没有拼音的字符(default/ignore/replace)
           strict::Bool=true,       # 是否严格遵照《汉语拼音方案》处理声母和韵母
           v2u::Bool=false,         # 是否用 `ü` 代替 `v`
           tone5::Bool=false)       # 是否用 5 标识轻声
"""
function pinyin(hanzi::AbstractString;
                style::Int=TONE,
                heteronym::Bool=false,
                errors::String="default",
                strict::Bool=true,
                v2u::Bool=false,
                tone5::Bool=false)

    return py"pinyins"(hanzi,
                       style=style,
                       heteronym=heteronym,
                       errors=errors,
                       strict=strict,
                       v_to_u=v2u,
                       neutral_tone_with_five=tone5)
end

"""
    printpinyin(pinyins)
将汉字转拼音的结果输出到一行并以空格分开
"""
function printpinyin(pinyins)
    n = length(pinyins)
    for i = 1:n-1
        print("$(pinyins[i]) ")
    end
    print(pinyins[n])
end

"""
    hanzi2pinyin(hanzi::AbstractString) = lazypinyin(hanzi, style=TONEONTOP, v2u=true, sandhi=true)
常见转拼音需求的快捷函数
# Example
julia> hanzi2pinyin("你好啊Julia，战略！") |> printpinyin
ní hǎo a Julia， zhàn lüè ！
"""
function hanzi2pinyin(hanzi::AbstractString)
    return lazypinyin(hanzi, style=TONEONTOP, v2u=true, sandhi=true)
end

end # module


text:
"""Pareto smoothed importance sampling (PSIS)

This module implements Pareto smoothed importance sampling (PSIS) and PSIS
leave-one-out cross-validation for Python (Numpy).

Included functions
------------------
psisloo
    Pareto smoothed importance sampling leave-one-out log predictive densities.

psislw
    Pareto smoothed importance sampling.

gpdfitnew
    Estimate the paramaters for the Generalized Pareto Distribution (GPD).

gpinv
    Inverse Generalised Pareto distribution function.

sumlogs
    Sum of vector where numbers are represented by their logarithms.

References
----------
Aki Vehtari, Andrew Gelman and Jonah Gabry (2015). Efficient implementation
of leave-one-out cross-validation and WAIC for evaluating fitted Bayesian
models. arXiv preprint arXiv:1507.04544.

Aki Vehtari and Andrew Gelman (2015). Pareto smoothed importance sampling.
arXiv preprint arXiv:1507.02646.
"""

module ParetoSmoothedImportanceSampling

using StatsFuns, Statistics

psis_path = @__DIR__

include("psisloo.jl")
include("psislw.jl")
include("gpdfitnew.jl")
include("gpinv.jl")
include("waic.jl")
include("pk_utilities.jl")
include("dic.jl")

export
    psis_path

end