text:
# car plotting: 2 plots (replaces 4 old plots in paper)
# 2D sets of x1 vs x2
# show how they are never within "goal set"
# ALSO: plot of SAT trajectory! (second plot) take vals from SAT and simulate forward in time

using PGFPlots
using QHull
using JLD2
using FileIO
using LazySets
using LinearAlgebra
include("models/problems.jl")
include("OverApprox/src/overapprox_nd_relational.jl")
include("OverApprox/src/overt_parser.jl")
include("MIP/src/overt_to_mip.jl")
include("MIP/src/mip_utils.jl")
include("models/car/simple_car.jl")

controller = ARGS[1] # "smallest"

one_step_data = load("examples/jmlr/data/car_10step/car_reachability_$(controller)_controller_data_1step.jld2")
one_step_state_sets = one_step_data["symbolic_state_sets"]
one_step_meas_sets = one_step_data["symbolic_meas_sets"]

data = load("examples/jmlr/data/car_10step/car_reachability_$(controller)_controller_data.jld2")

concrete_state_sets = vcat(data["concrete_state_sets"]...)
symbolic_state_sets = vcat(data["symbolic_state_sets"]...)

# monte carlo simulate trajectories for first plot
mc_state_sets, xvec, x0, mc_meas_sets, yvec, y0 = monte_carlo_simulate(data["query"], data["input_set"], n_sim=1000000)

# clean up sets
init, reachable_state_sets = clean_up_sets(concrete_state_sets, symbolic_state_sets, data["concretization_intervals"])

#########################################
# Define Colors and Styles ##############
#########################################
define_color("concrete_color", 0x139EAB)
define_color("symbolic_color", 0x9BFF85)
define_color("mc_color", 0x38caff)
define_color("goal_color", 0xa866ff)
conc_style_solid = "solid, concrete_color, mark=none, fill=concrete_color"
conc_style_transparent = conc_style_solid*", fill opacity=0.5, opacity=0.5"
sym_style_solid = "solid, symbolic_color, mark=none, fill=symbolic_color"
sym_style_transparent =  sym_style_solid*", fill opacity=0.5, opacity=0.5"
mc_style_solid = "solid, mc_color, mark=none, fill=mc_color"
mc_style_transparent = "solid, red, mark=none, fill=red, fill opacity=0.5, opacity=0.5"#mc_style_solid*", fill opacity=0.5"

concrete_line_style = "solid, concrete_color, mark=none"
symbolic_line_style = "solid, symbolic_color, very thick, mark=none"
mc_line_style = "solid, mc_color, ultra thick, mark=none"

goal_style = "solid, goal_color, thick, mark=none, fill=goal_color"

#########################################
##### plot 1: 2D sets of x1 vs x2  ######
#########################################
dims=[1,2]
fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$y~(m)\$", xlabel="\$x~(m)\$", title="Car Position Reachable Sets")

# push!(fig, PGFPlots.Plots.Linear([9., 9.], [-4., -4.], style=mc_style_transparent, legendentry="Monte Carlo Simulations hyperrectangular hull"))

# plot init set in both concrete and hybrid colors
inputx, inputy = get_rectangle(data["input_set"], dims)
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=conc_style_solid, legendentry="Concrete Sets"))
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=sym_style_solid, legendentry="OVERT Hybrid Symbolic Sets"))

# hack for legend
push!(fig, PGFPlots.Plots.Linear([9., 9.], [-4., -4.], style=mc_style_solid, legendentry="Monte Carlo Simulations convex hull"))

# plot goal set
goal_set = Hyperrectangle(low=[-0.6, -0.2], high=[0.6, 0.2])
goal_x, goal_y = get_rectangle(goal_set, dims)
push!(fig, PGFPlots.Plots.Linear(goal_x, goal_y, style=goal_style, legendentry="Goal Set"))

### plot init set of mc simulations
points = x0[:, dims]
#plot hull
border_idx = chull(points).vertices
x = points[border_idx, 1]
y = points[border_idx, 2]
push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style_solid))

for t in 1:data["query"].ntime
    # if t == data["query"].ntime
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=conc_style_solid))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=sym_style_solid))
        # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style_transparent))
        # plot mc_sim points
        points = xvec[:, t, dims]
        #plot hull
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style_solid))
    # else
    #     push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=conc_style_transparent))
    #     push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=sym_style_transparent))
    #     # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style_transparent))
    #     # plot mc simulations hull
    #     points = xvec[:, t, dims]
    #     border_idx = chull(points).vertices
    #     x = points[border_idx, 1]
    #     y = points[border_idx, 2]
    #     push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style_solid))
    # end
end

fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/car/car_$(controller)_x12.tex", fig)
PGFPlots.save("examples/jmlr/plots/car/car_$(controller)_x12.pdf", fig)

#########################################
##### plot 2: SAT trajectory  ######
#########################################
# counter examples are "easy" here because no trajectory reaches the goal :'D
# load sat data
q1_data = load("examples/jmlr/data/car_10step/car_satisfiability_$(controller)_controller_data_q1.jld2")
q2_data = load("examples/jmlr/data/car_10step/car_satisfiability_$(controller)_controller_data_q2.jld2")
final_data = load("examples/jmlr/data/car_10step/car_satisfiability_$(controller)_controller_data_final_result.jld2")
# final data contains the SATuses "sat" or "unsat" for each timestep.
# if one of the queries has multiple timesteps that were sat, this produces multiple possible counterexample "hints"
# we take the first timestep values of any "sat" step, and simulate them forward to get a real counterexample. we could do this for every timestep that returns SAT
x0 = [q2_data["statii"][1][Symbol("$(v)_1")] for v in q2_data["query"].problem.input_vars]
xvec_mc = monte_carlo_one_simulate(q2_data["query"], x0)

fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$y~(m)\$", xlabel="\$x~(m)\$", title="Car Counter Example")

dims=[1,2]
# plot goal set
goal_set = Hyperrectangle(low=[-0.6, -0.2], high=[0.6, 0.2])
goal_x, goal_y = get_rectangle(goal_set, dims)
push!(fig, PGFPlots.Plots.Linear(goal_x, goal_y, style=goal_style, legendentry="Goal Set"))

counter_example_style = "concrete_color, solid, mark=*, mark options={concrete_color, fill=concrete_color}"
# plot the counter example
push!(fig, PGFPlots.Plots.Linear(xvec_mc[:, dims[1]], xvec_mc[:, dims[2]], style=counter_example_style, legendentry="Counter Example"))

fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/car/car_$(controller)_counter_example.tex", fig)
PGFPlots.save("examples/jmlr/plots/car/car_$(controller)_counter_example.pdf", fig)

####################
# Plot 3: 1D sets
####################

dims = [3]
conc_subsets = get_interval_subsets(one_step_state_sets, dims)
sym_subsets = get_interval_subsets(reachable_state_sets, dims)
mc_subsets = get_interval_subsets(mc_state_sets, dims)
deleteat!(mc_subsets, 1) # pop init set off start

# plot intervals using lines in PGF plots
fig = PGFPlots.Axis(style="width=8cm, height=8cm", ylabel="\$x_$(dims[1])\$", xlabel="timesteps", title="Car Reachable Sets, $controller Controller")

# plot init set
push!(fig, PGFPlots.Plots.Linear([0.,0.], [low(data["input_set"])[dims]..., high(data["input_set"])[dims]...], style=concrete_line_style, legendentry="Concrete Sets"))
# plot init set
push!(fig, PGFPlots.Plots.Linear([0.,0.], [low(data["input_set"])[dims]..., high(data["input_set"])[dims]...], style=symbolic_line_style, legendentry="OVERT Hybrid Symbolic Sets"))
# plot init set
push!(fig, PGFPlots.Plots.Linear([0.,0.], [low(data["input_set"])[dims]..., high(data["input_set"])[dims]...], style=mc_line_style, legendentry="Monte Carlo Simulations"))
# avoid set
# push!(fig, PGFPlots.Plots.Linear([-1., 26], [-0.2617, -0.2617], style="red, dashed", mark="none", legendentry="Unsafe Set"))
# plot the rest of the sets
for t in 1:data["query"].ntime
    timestep = Float64[t,t]
    push!(fig, PGFPlots.Plots.Linear(timestep, conc_subsets[t], style=concrete_line_style))
    push!(fig, PGFPlots.Plots.Linear(timestep,sym_subsets[t], style=symbolic_line_style, markSize=1))
    push!(fig, PGFPlots.Plots.Linear(timestep,mc_subsets[t], style=mc_line_style, markSize=1))
end

fig.legendStyle = "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/car/car_dim$(dims[1])_$(controller).tex", fig)
PGFPlots.save("examples/jmlr/plots/car/car_dim$(dims[1])_$(controller).pdf", fig)


text:
# Julia implementation of ROUGE score
# Ref: https://github.com/google/seq2seq/blob/master/seq2seq/metrics/rouge.py

"""
    _get_ngrams(n, text)

Calcualtes n-grams. Returns a set of n-grams.

# Arguments:
 - `n`: provide which n-grams to calculate
 - `text`: An array of tokens
"""
function _get_ngrams(n, text)
    ngrams_set = Set()
    text_length = length(text)
    max_index_ngram_start = text_length - n
    for i in 1:max_index_ngram_start+1
        push!(ngrams_set, tuple(text[i:i+n-1]))
    end
    return ngrams_set
end

"""
    _split_into_words(sentences)

Splits multiple sentences into words and flattens the result
"""
function _split_into_words(sentences)
    words = []
    for i in 1:length(sentences)
        temp = split(sentences[i], " ")
        for j in temp
            push!(words,j)
        end
    end
    return words
end

"""
    _get_word_ngrams(n, sentences)

Calculates word n-grams for multiple sentences.
"""
function _get_word_ngrams(n, sentences)
    @assert length(sentences) > 0
    @assert n>0

    words = _split_into_words(sentences)
    return _get_ngrams(n, words)
end

"""
    _lcs(x, y)

Utility function to compute the length of the longest common subsequence (lcs) between two
strings. The implementation below uses a DP programming algorithm and runs
in O(nm) time where n = len(x) and m = len(y).
"""
function _lcs(x, y)
    n, m = length(x), length(y)
    table = Dict()
    for i in 1:n+1
        for j in 1:m+1
            if i==1 || j==1
                table[i, j] = 0
            elseif x[i-1] == y[j-1]
                table[i, j] = table[i-1, j-1] + 1
            else
                table[i, j] = max(table[i-1, j], table[i, j-1])
            end
        end
    end
    return table
end

"""
    _len_lcs(x, y)

Returns the length of the Longest Common Subsequence between sequences x and y.
"""
function _len_lcs(x, y)
    table = _lcs(x, y)
    n, m = length(x), length(y)
    return table[n+1, m+1]
end

"""
    _recons_lcs(x, y)

Returns the Longest Subsequence between x and y.
"""
function _recon_lcs(x, y)
    i , j = length(x), length(y)
    table = _lcs(x, y)

    function _recon(i, j)
        if i == 1 || j == 1
            return []
        elseif x[i-1] == y[j-1]
            return push!(_recon(i - 1, j - 1) ,(x[i - 1], i))
        elseif table[i - 1, j] > table[i, j - 1]
            return _recon(i - 1, j)
        else
            return _recon(i, j - 1)
        end
    end

    recon_tuple = tuple(map(x->x[1][1], _recon(i, j))...)
end

"""
    rouge_n(evaluated_sentences, reference_sentences; n=2)

Computes ROUGE-N of two text collections of sentences. Returns f1, precision, recall for ROUGE-N.

# Arguments:
 - `evaluated_sentences`: the sentences that have been picked by the summarizer
 - `reference_sentences`: the sentences from the referene set
 - `n`: size of ngram.  Defaults to 2.

Source: (http://research.microsoft.com/en-us/um/people/cyl/download/
  papers/rouge-working-note-v1.3.1.pdf)
"""
function rouge_n(evaluated_sentences, reference_sentences; n=2)
    if length(evaluated_sentences) <= 0 || length(reference_sentences)<=0
        throw(ArgumentError())
    end

    evaluated_ngrams = _get_word_ngrams(n, evaluated_sentences)
    reference_ngrams = _get_word_ngrams(n, reference_sentences)
    reference_count = length(reference_ngrams)
    evaluated_count = length(evaluated_ngrams)

    overlapping_ngrams = intersect(evaluated_ngrams, reference_ngrams)
    overlapping_count = length(overlapping_ngrams)
    precision = 0.0
    if evaluated_count != 0
         precision = overlapping_count / evaluated_count
    end

    recall = 0.0
    if reference_count != 0
         recall = overlapping_count / reference_count
    end

    f1_score = 2.0 * ((precision * recall)) / (precision + recall + 1e-9)
    return f1_score, precision, recall
end

"""
    _f_p_r_lcs(llcs, m, n)

Computes the LCS-based F-measure score

# Arguments:
 - `llcs`: Length of LCS
 - `m`: number of words in reference summary
 - `n`: number of words in candidate summary

Source: (http://research.microsoft.com/en-us/um/people/cyl/download/papers/rouge-working-note-v1.3.1.pdf)
"""
function _f_p_r_lcs(llcs, m, n)
    r_lcs = llcs / m
    p_lcs = llcs / n
    beta = p_lcs / (r_lcs + 1e-12)
    num = (1 + (beta^2)) * r_lcs * p_lcs
    denom = r_lcs + ((beta^2) * p_lcs)
    f_lcs = num / (denom + 1e-12)
    return f_lcs, p_lcs, r_lcs
end

"""
    rouge_l_sentence_level(evaluated_sentences, reference_sentences)

Computes ROUGE-L (sentence level) of two text collections of sentences.

Calculated according to:
  R_lcs = LCS(X,Y)/m,
  P_lcs = LCS(X,Y)/n,
  F_lcs = ((1 + beta^2)*R_lcs*P_lcs) / (R_lcs + (beta^2) * P_lcs)

where:
  X = reference summary
  Y = Candidate summary
  m = length of reference summary
  n = length of candidate summary

# Argumnets:
 - `evaluated_sentences`: the sentences that have been picked by the summarizer
 - `reference_sentences`: the sentences from the referene set

Source: (http://research.microsoft.com/en-us/um/people/cyl/download/papers/rouge-working-note-v1.3.1.pdf)
"""
function rouge_l_sentence_level(evaluated_sentences, reference_sentences)
    if length(evaluated_sentences) <= 0 || length(reference_sentences)<=0
        throw(ArgumentError())
    end

    reference_words = _split_into_words(reference_sentences)
    evaluated_words = _split_into_words(evaluated_sentences)
    m = length(reference_words)
    n = length(evaluated_words)
    lcs = _len_lcs(evaluated_words, reference_words)
    return _f_p_r_lcs(lcs, m, n)
end

"""
    _union_lcs(evaluated_sentences, reference_sentence)

Returns LCS_u(r_i, C) which is the LCS score of the union longest common subsequence between reference sentence ri and candidate summary C.

# Arguments:
 - `evaluated_sentences`: the sentences that have been picked by the summarizer
 - `reference_sentence`: one of the sentences in the reference summaries

For example, if r_i= w1 w2 w3 w4 w5, and C contains two sentences: c1 = w1 w2 w6 w7 w8 and c2 = w1 w3 w8 w9 w5, then the longest common subsequence of r_i and c1 is “w1 w2” and the longest common subsequence of r_i and c2 is “w1 w3 w5”. The union longest common subsequence of r_i, c1, and c2 is “w1 w2 w3 w5” and LCS_u(r_i, C) = 4/5.

"""
function _union_lcs(evaluated_sentences, reference_sentence)
    if length(evaluated_sentences) <= 0
        throw(ArgumentError())
    end

    lcs_union = Set()
    reference_words = _split_into_words([reference_sentence])
    combined_lcs_length = 0
    for eval_s in evaluated_sentences
        evaluated_words = _split_into_words([eval_s])
        lcs = Set(_recon_lcs(reference_words, evaluated_words))
        combined_lcs_length += length(lcs)
        lcs_union = union(lcs, lcs_union)
    end

    union_lcs_count = length(lcs_union)
    union_lcs_value = union_lcs_count / combined_lcs_length
    return union_lcs_value
end

"""
    rouge_l_summary_level(evaluated_sentences, reference_sentences)

Computes ROUGE-L (summary level) of two text collections of sentences.

Calculated according to:
  R_lcs = SUM(1, u)[LCS<union>(r_i,C)]/m
  P_lcs = SUM(1, u)[LCS<union>(r_i,C)]/n
  F_lcs = ((1 + beta^2)*R_lcs*P_lcs) / (R_lcs + (beta^2) * P_lcs)

where:
  SUM(i,u) = SUM from i through u
  u = number of sentences in reference summary
  C = Candidate summary made up of v sentences
  m = number of words in reference summary
  n = number of words in candidate summary

# Arguments:
 - `evaluated_sentences`: the sentences that have been picked by the summarizer
 - `reference_sentence`: the sentences in the reference summaries

Source: (http://research.microsoft.com/en-us/um/people/cyl/download/papers/rouge-working-note-v1.3.1.pdf)

"""
function rouge_l_summary_level(evaluated_sentences, reference_sentences)
    if length(evaluated_sentences) <= 0 || length(reference_sentences)<=0
        throw(ArgumentError())
    end
    m = length(_split_into_words(reference_sentences))
    n = length(_split_into_words(evaluated_sentences))
    union_lcs_sum_across_all_references = 0
    for ref_s in reference_sentences
        union_lcs_sum_across_all_references += _union_lcs(evaluated_sentences, ref_s)
    end
    return _f_p_r_lcs(union_lcs_sum_across_all_references, m, n)
end

"""
    rouge(hypotheses, references)

Calculates average rouge scores for a list of hypotheses and references.
"""
function rouge(hypotheses, references)
    rouge_1 = [rouge_n([hyp],[ref], n=1) for (hyp, ref) in zip(hypotheses, references)]
    rouge_1_f, rouge_1_p, rouge_1_r = 0.0,0.0,0.0
    for i in 1:length(rouge_1)
        rouge_1_f+=rouge_1[i][1]
        rouge_1_p+=rouge_1[i][2]
        rouge_1_r+=rouge_1[i][3]
    end
    n = length(rouge_1)
    rouge_1_f/= n
    rouge_1_p/= n
    rouge_1_r/= n

    rouge_2 = [rouge_n([hyp],[ref],n=2) for (hyp, ref) in zip(hypotheses, references)]
    rouge_2_f, rouge_2_p, rouge_2_r = 0.0,0.0,0.0
    for i in 1:length(rouge_2)
        rouge_2_f+=rouge_2[i][1]
        rouge_2_p+=rouge_2[i][2]
        rouge_2_r+=rouge_2[i][3]
    end
    n = length(rouge_1)
    rouge_2_f/= n
    rouge_2_p/= n
    rouge_2_r/= n

    rouge_l = [rouge_l_sentence_level([hyp],[ref]) for (hyp, ref) in zip(hypotheses, references)]
    rouge_l_f, rouge_l_p, rouge_l_r = 0.0,0.0,0.0
    for i in 1:length(rouge_l)
        rouge_l_f+=rouge_l[i][1]
        rouge_l_p+=rouge_l[i][2]
        rouge_l_r+=rouge_l[i][3]
    end
    n = length(rouge_l)
    rouge_l_f/= n
    rouge_l_p/= n
    rouge_l_r/= n

    return OrderedDict(
      "rouge_1 / f_score"=> rouge_1_f,
      "rouge_1 / r_score"=> rouge_1_r,
      "rouge_1 / p_score"=> rouge_1_p,
      "rouge_2 / f_score"=> rouge_2_f,
      "rouge_2 / r_score"=> rouge_2_r,
      "rouge_2 / p_score"=> rouge_2_p,
      "rouge_l / f_score"=> rouge_l_f,
      "rouge_l / r_score"=> rouge_l_r,
      "rouge_l / p_score"=> rouge_l_p,)

end

ext: jl


text:
"""
$(TYPEDEF)

A surreal number form with explicit left and right subsets.

A (valid) form represents a specific surreal number, and so instances of this
class behave as the number they represent.
"""
struct SurrealForm{L<:SurrealSet, R<:SurrealSet} <: Surreal
	l::L
	r::R

	SurrealForm(l::SurrealSet, r::SurrealSet) = new{typeof(l), typeof(r)}(l, r)
end

SurrealForm() = SurrealForm((), ())
SurrealForm(f::SurrealForm) = f
SurrealForm(l, r) = SurrealForm(SurrealSet(l), SurrealSet(r))  # Convert arguments
SurrealForm(x) = SurrealForm(left(x), right(x))

function SurrealForm(args...; check::Bool)
	f = SurrealForm(args...)
	check && check_numeric(f)
	return f
end

left(f::SurrealForm) = f.l
right(f::SurrealForm) = f.r
issurreal(f::SurrealForm) = true

Base.zero(::SurrealForm) = SurrealForm()
Base.one(::SurrealForm) = SurrealForm((0,), ())

"""
	isnumeric(f::SurrealForm)

Whether the form is *numeric* - i.e. represents a valid surreal number.
"""
isnumeric(f::SurrealForm) = all_lt(left(f), right(f))

checknumeric(f::SurrealForm) = isnumeric(f) || error("Form is not numeric")

Base.show(io::IO, mime::Union{FORM_TEXT_MIMES...}, f::SurrealForm) = show_form_outer(io, mime, f)


text:
# include("precompile_funcs.jl")

using TuringDash

# using .PrecompileFuncs
# PrecompileFuncs.precompile_turing()

app = TuringDash.make_app()

port = haskey(ENV, "PORT") ? parse(Int64, ENV["PORT"]) : 8050

TuringDash.Dash.run_server(app, "0.0.0.0", port)


text:
using LazySets.Approximations: project

for N in [Float64, Rational{Int}, Float32]
    c = N[0, 0]
    b = Ball1(c, N(1))

    # overapproximating a set of type T1 with an unsupported type T2 is the
    # identity if T1 = T2
    @test_throws MethodError overapproximate(ZeroSet{N}(2), EmptySet)
    e = EmptySet{N}()
    @test overapproximate(e, EmptySet) == e
    # the third argument is ignored
    oa = overapproximate(b, Ball1, Hyperrectangle)
    @test oa isa Ball1

    # HPolygon approximation with box directions
    p = overapproximate(b, HPolygon)
    for d in [N[1, 0], N[-1, 0]]
        @test σ(d, p)[1] ≈ σ(d, b)[1]
    end
    for d in [N[0, 1], N[0, -1]]
        @test σ(d, p)[2] ≈ σ(d, b)[2]
    end

    # Hyperrectangle approximation
    p = overapproximate(b, Hyperrectangle)
    for d in [N[1, 0], N[-1, 0]]
        @test σ(d, p)[1] ≈ σ(d, b)[1]
    end
    for d in [N[0, 1], N[0, -1]]
        @test σ(d, p)[2] ≈ σ(d, b)[2]
    end
    @test p.center ≈ c
    @test p.radius ≈ N[1, 1]

    # Interval approximation
    b = Ball1(N[0], N(1))
    p = overapproximate(b, LazySets.Interval)
    for d in [N[1], N[-1]]
        @test σ(d, p)[1] ≈ σ(d, b)[1]
    end

    # approximation with an axis-aligned hyperrectangle
    Z = Zonotope(N[-1.0, -1.0], N[-1/2 0; -1.2 -1])
    Zoa = overapproximate(Z, Hyperrectangle) # faster o.a.
    Zba = box_approximation(Z) # default o.a. implementation that uses supp function
    @test Zoa.center ≈ Zba.center && Zoa.radius ≈ Zba.radius

    # same but using cartesian product
    h1 = Hyperrectangle(N[1/2],  N[1/2])
    h2 = Hyperrectangle(N[2.5, 4.5],  N[1/2, 1/2])
    H = overapproximate(h1 × h2, Hyperrectangle) # defaults to convert method
    @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]

    # overapproximation of the lazy linear map of a hyperrectangular set
    H = Hyperrectangle(N[0, 0], N[1/2, 1])
    M = Diagonal(N[2, 2])
    OA = overapproximate(M*H, Hyperrectangle)
    @test OA isa Hyperrectangle && OA.center == N[0, 0] && OA.radius == N[1, 2]

    #overapproximation of Minkowski sum of linear maps for each block in the row block
    i1 = Interval(N[0, 1])
    h = Hyperrectangle(low=N[3, 4], high=N[5, 7])
    M = N[1 2 3; 4 5 6; 7 8 9]
    cpa = CartesianProductArray([i1, h])
    lm = M * cpa

    oa = overapproximate(lm, Hyperrectangle)
    oa_box = overapproximate(lm, Approximations.BoxDirections)
    d_oa_d_hp = overapproximate(lm, CartesianProductArray{N, Hyperrectangle{N}})
    d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections)
    oa_d_hp = overapproximate(d_oa_d_hp)
    oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections)

    @test oa == oa_d_hp
    @test oa_box == oa_d_box

    for (oax, set_type) in [(d_oa_d_hp, Hyperrectangle), (d_oa_d_box, HPolytope)]
        @test oax isa CartesianProductArray
        arr = oax.array
        @test length(arr) == 2 && dim(arr[1]) == 1 && dim(arr[2]) == 2
        @test all(X -> X isa set_type, arr)
    end

    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    i3 = Interval(N[1, 4])
    cpa = CartesianProductArray([i1, i2, i3])
    M = N[1 2 0; 0 1 0; 0 1 1]
    lm = M * cpa
    d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}})
    oa = overapproximate(lm)
    @test overapproximate(d_oa) == oa
    @test typeof(d_oa) == CartesianProductArray{N, Interval{N}}

    # interval linear map of zonotope
    M = IntervalMatrix([(N(0) ± N(1)) (N(0) ± N(0)); (N(0) ± N(0)) (N(0) ± N(1))])
    Z = Zonotope(N[1, 1], N[1 1; 1 -1])
    lm = LinearMap(M, Z)
    Zo = overapproximate(lm, Zonotope)
    @test box_approximation(Zo) == Hyperrectangle(N[0, 0], N[3, 3])

    # intersection of two polyhedra
    P = HPolyhedron([HalfSpace(N[1, 0], N(1)), HalfSpace(N[0, 1], N(1))])
    Q = HPolyhedron([HalfSpace(N[-1, 0], N(1)), HalfSpace(N[0, -1], N(1))])
    oa = overapproximate(P ∩ Q, BoxDirections)
    B = BallInf(N[0, 0], N(1))
    @test B ⊆ oa && oa ⊆ B

    # rectification
    r = Rectification(EmptySet{N}())
    @test overapproximate(r, Hyperrectangle) isa EmptySet{N}
    r = Rectification(Ball1(N[0, 0], N(1)))
    @test overapproximate(r, Hyperrectangle) ==
        Hyperrectangle(low=N[0, 0], high=N[1, 1])
end

# tests that do not work with Rational{Int}
for N in [Float64, Float32]
    # useful for benchmarking overapproximate and LinearMap's support vector
    # (see #290)
    function overapproximate_lmap(n)
        B = BallInf(ones(N, n), N(2))
        π = sparse(N[1, 2], N[1, 2], ones(N, 2), 2, n)
        return Approximations.overapproximate(π*B)
    end
    o = overapproximate_lmap(50)
    @test o.center == N[1, 1] && o.radius == N[2, 2]

    # Approximation of a 2D centered unit ball in norm 1
    # All vertices v should be like this:
    # ‖v‖ >= 1 and ‖v‖ <= 1+ε
    # Where ε is the given error bound
    b = Ball1(N[0, 0], N(1))
    ε = N(0.01)
    p = tovrep(overapproximate(b, ε))
    for v in vertices_list(p)
    @test norm(v) >= N(1)
    @test norm(v) <= N(1 + ε)
    end

    # Check that there are no redundant constraints for a ballinf
    b = BallInf(N[0.5, 0.5], N(0.1))
    lcl = overapproximate(b, N(0.001)).constraints
    @test length(lcl) == 4
    @test lcl[1].a == N[1, 0]
    @test lcl[1].b == N(0.6)
    @test lcl[2].a == N[0, 1]
    @test lcl[2].b == N(0.6)
    @test lcl[3].a == N[-1, 0]
    @test lcl[3].b == N(-0.4)
    @test lcl[4].a == N[0, -1]
    @test lcl[4].b == N(-0.4)

    # Check that there are no redundant constraints for a HPolygon (octagon)
    p = HPolygon{N}()
    addconstraint!(p, LinearConstraint(N[1, 0], N(1)))
    addconstraint!(p, LinearConstraint(N[0, 1], N(1)))
    addconstraint!(p, LinearConstraint(N[-1, 0], N(1)))
    addconstraint!(p, LinearConstraint(N[0, -1], N(1)))
    addconstraint!(p, LinearConstraint(N[sqrt(2)/2, sqrt(2)/2], N(1)))
    addconstraint!(p, LinearConstraint(N[-sqrt(2)/2, sqrt(2)/2], N(1)))
    addconstraint!(p, LinearConstraint(N[sqrt(2)/2, -sqrt(2)/2], N(1)))
    addconstraint!(p, LinearConstraint(N[-sqrt(2)/2, -sqrt(2)/2], N(1)))
    lcl = overapproximate(p, N(.001)).constraints
    @test length(lcl) == 8
    @test lcl[1].a ≈ N[1, 0]
    @test lcl[1].b ≈ N(1)
    @test lcl[2].a ≈ N[sqrt(2)/2, sqrt(2)/2]
    @test lcl[2].b ≈ N(1)
    @test lcl[3].a ≈ N[0, 1]
    @test lcl[3].b ≈ N(1)
    @test lcl[4].a ≈ N[-sqrt(2)/2, sqrt(2)/2]
    @test lcl[4].b ≈ N(1)
    @test lcl[5].a ≈ N[-1, 0]
    @test lcl[5].b ≈ N(1)
    @test lcl[6].a ≈ N[-sqrt(2)/2, -sqrt(2)/2]
    @test lcl[6].b ≈ N(1)
    @test lcl[7].a ≈ N[0, -1]
    @test lcl[7].b ≈ N(1)
    @test lcl[8].a ≈ N[sqrt(2)/2, -sqrt(2)/2]
    @test lcl[8].b ≈ N(1)

    # Zonotope approximation of convex hull of two zonotopes

    # same order (mean method only)
    Z1 = Zonotope(zeros(N, 2), hcat(N[1, 0]))
    Z2 = Zonotope(zeros(N, 2), hcat(N[0, 1]))
    Zch = overapproximate(ConvexHull(Z1, Z2), Zonotope; algorithm="mean")
    # the result is a diamond (tight)
    @test Zch == Zonotope(N[0, 0], N[0.5 0.5; 0.5 -0.5])

    # different order (mean method only)
    Z1 = Zonotope(zeros(N, 2), hcat(N[1, 0]))
    Z2 = Zonotope(N[0, 0], N[1 0; 0 1])
    # the result is a box (tight)
    Zch = overapproximate(ConvexHull(Z1, Z2), Zonotope; algorithm="mean")
    @test Zch == Zonotope(N[0, 0], N[1 0; 0 1])

    # comparison of different methods
    Z1 = Zonotope(ones(N, 2), [N[1, 0], N[0, 1], N[1, 1]])
    Z2 = Zonotope(-ones(N, 2), [N[0.5, 1], N[-0.1, 0.9], N[1, 4]])
    Y = ConvexHull(Z1, Z2)
    # overapproximate with a polygon
    Y_polygon = overapproximate(Y, N(1e-3))
    # overapproximate with a zonotope (different methods)
    overapproximate(Y, Zonotope) # default
    Y_zonotope_1 = overapproximate(Y, Zonotope; algorithm="mean")
    Y_zonotope_2 = overapproximate(Y, Zonotope; algorithm="join")
    @test Y_polygon ⊆ Y_zonotope_1
    @test Y_polygon ⊆ Y_zonotope_2
    # check property of 'join' algorithm: box approximation equality
    X1 = box_approximation(Y)
    X2 = box_approximation(Y_zonotope_2)
    @test center(X1) ≈ center(X2) &&
          radius_hyperrectangle(X1) ≈ radius_hyperrectangle(X2)

    # example from paper
    Z1 = Zonotope(N[3, 0], N[1 2; 1 1])
    Z2 = Zonotope(N[1, 0], N[-2 1; 1 1])
    @test overapproximate(ConvexHull(Z1, Z2), Zonotope; algorithm="join") ==
          Zonotope(N[2, 0], N[0 1 3; 1 1 0])

    # different orders
    Z1 = Zonotope(N[3, 0], N[1 2 1; 1 1 2])
    Z2 = Zonotope(N[1, 0], N[-2 1; 1 1])
    Y = ConvexHull(Z1, Z2)
    Y_polygon = overapproximate(Y, N(1e-3))
    Y_zonotope_1 = overapproximate(Y, Zonotope; algorithm="mean")
    Y_zonotope_2 = overapproximate(Y, Zonotope; algorithm="join")
    @test Y_polygon ⊆ Y_zonotope_1
    @test Y_polygon ⊆ Y_zonotope_2

    # ResetMap and CPA
    X = CartesianProductArray([Ball1(N[0, 0], N(1)), Ball1(N[0, 0], N(2)),
                               Ball1(N[0, 0], N(3))])
    resets = Dict(3 => N(0), 6 => N(0))
    rm = ResetMap(X, resets)
    Y = overapproximate(rm, CartesianProductArray, Hyperrectangle)
    @test array(Y) == [Hyperrectangle(N[0, 0], N[1, 1]),
        Hyperrectangle(N[0, 0], N[0, 2]), Hyperrectangle(N[0, 0], N[3, 0])]
end

for N in [Float64]
    # decomposed linear map approximation
    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    M = N[1 2; 0 1]
    cpa = CartesianProductArray([i1, i2])
    lm = M * cpa
    d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}})
    oa = overapproximate(lm, OctDirections)
    @test oa ⊆ d_oa

    # decomposed intersection between Cartesian product array and polyhedron
    h1 = Hyperrectangle(low=N[3, 4], high=N[5, 7])
    h2 = Hyperrectangle(low=N[5, 5], high=N[6, 8])
    cpa1 = CartesianProductArray([i1, i2, h1])
    o_cpa1 = overapproximate(cpa1)
    cpa2 = CartesianProductArray([i1, i2, h2])
    o_cpa2 = overapproximate(cpa2)
    G = HalfSpace(N[1, 0, 0, 0], N(1))
    G_3 = HalfSpace(N[0, 0, 1, 0], N(1))
    G_3_neg = HalfSpace(N[0, 0, -1, 0], N(0))
    G_comb = HalfSpace(N[1, 1, 0, 0], N(2.5))

    d_int_g = Intersection(cpa1, G)
    d_int_g_3 = Intersection(cpa1, G_3)
    d_int_g_3_neg = Intersection(cpa1, G_3_neg)
    d_int_cpa = intersection(cpa1, cpa2)
    l_int_cpa = intersection(o_cpa1, o_cpa2)
    o_d_int_g = overapproximate(d_int_g, CartesianProductArray, Hyperrectangle)
    o_d_int_g_3 = overapproximate(d_int_g_3, CartesianProductArray, Hyperrectangle)
    o_d_int_g_3_neg = overapproximate(d_int_g_3_neg, CartesianProductArray, Hyperrectangle)

    @test overapproximate(o_d_int_g) ≈ overapproximate(d_int_g)
    @test isempty(d_int_g_3) && o_d_int_g_3 == EmptySet{N}()
    @test overapproximate(o_d_int_g_3_neg) ≈ overapproximate(d_int_g_3_neg)
    @test overapproximate(d_int_cpa, Hyperrectangle) == l_int_cpa
    @test all([X isa CartesianProductArray for X in [d_int_cpa, o_d_int_g, o_d_int_g_3_neg]])

    int_g_comb = Intersection(cpa1, G_comb)
    o_d_int_g_comb = overapproximate(int_g_comb, CartesianProductArray, Hyperrectangle)
    o_bd_int_g_comb = overapproximate(int_g_comb, BoxDirections)
    @test overapproximate(o_d_int_g_comb) ≈ overapproximate(o_bd_int_g_comb)
    projection = project(o_bd_int_g_comb, [1, 2], BoxDirections)
    @test vertices_list(projection) ≈ [[0.5, 2.5], [0., 2.5], [0., 2], [0.5, 2.]]
    projection = project(overapproximate(int_g_comb, OctDirections), [1, 2], OctDirections)
    @test vertices_list(projection) ≈ [[0., 2.5], [0., 2.], [0.5, 2.]]

    # Zonotope approximation
    Z1 = Zonotope(ones(N, 2), [N[1, 0], N[0, 1], N[1, 1]])
    Z2 = Zonotope(-ones(N, 2), [N[0.5, 1], N[-0.1, 0.9], N[1, 4]])
    Y = ConvexHull(Z1, Z2)
    Y_polygon = overapproximate(Y, N(1e-3)) # overapproximate with a polygon
    Y_zonotope = overapproximate(Y, Zonotope) # overapproximate with a zonotope
    @test Y_polygon ⊆ Y_zonotope

    #decomposed linear map approximation
    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    M = N[1 2; 0 1]
    cpa = CartesianProductArray([i1, i2])
    lm = M * cpa
    d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}})
    oa = overapproximate(lm, OctDirections)
    @test oa ⊆ d_oa

    # =======================================
    # Zonotope overapprox. of a Taylor model
    # =======================================
    x₁, x₂, x₃ = set_variables(N, ["x₁", "x₂", "x₃"], order=5)
    Dx₁ = IA.Interval(N(1.0), N(3.0))
    Dx₂ = IA.Interval(N(-1.0), N(1.0))
    Dx₃ = IA.Interval(N(-1.0), N(0.0))
    D = Dx₁ × Dx₂ × Dx₃   # domain
    x0 = IntervalBox(IA.mid.(D)...)
    I = IA.Interval(N(0.0), N(0.0)) # interval remainder
    p₁ = 1 + x₁ - x₂
    p₂ = x₃ - x₁
    vTM = [TaylorModels.TaylorModelN(pi, I, x0, D) for pi in [p₁, p₂]]
    Z1 = overapproximate(vTM, Zonotope)
    @test center(Z1) == N[3, -2.5]
    @test Matrix(genmat(Z1)) == N[1 -1 0; -1 0 0.5]
end


text:
@testset "5_rotation" begin
  @testset "rotation_ckks" begin
    @testset "EncryptionParameters" begin
      @test_nowarn EncryptionParameters(SchemeType.ckks)
    end
    enc_parms = EncryptionParameters(SchemeType.ckks)

    @testset "polynomial modulus degree" begin
      @test_nowarn set_poly_modulus_degree!(enc_parms, 8192)
      @test poly_modulus_degree(enc_parms) == 8192
    end

    @testset "coefficient modulus" begin
      @test_nowarn coeff_modulus_create(8192, [40, 40, 40, 40, 40])
      @test_nowarn set_coeff_modulus!(enc_parms, coeff_modulus_create(8192, [40, 40, 40, 40, 40]))
    end

    @testset "SEALContext" begin
      @test_nowarn SEALContext(enc_parms)
    end
    context = SEALContext(enc_parms)

    @testset "KeyGenerator" begin
      @test_nowarn KeyGenerator(context)
    end
    keygen = KeyGenerator(context)

    @testset "PublicKey" begin
      @test_nowarn PublicKey()
    end
    public_key_ = PublicKey()

    @testset "create_public_key" begin
      @test_nowarn create_public_key!(public_key_, keygen)
    end

    @testset "SecretKey" begin
      @test_nowarn secret_key(keygen)
    end
    secret_key_ = secret_key(keygen)

    @testset "RelinKeys" begin
      @test_nowarn RelinKeys()
    end
    relin_keys_ = RelinKeys()

    @testset "create_relin_keys" begin
      @test_nowarn create_relin_keys!(relin_keys_, keygen)
    end

    @testset "GaloisKeys" begin
      @test_nowarn GaloisKeys()
    end
    galois_keys_ = GaloisKeys()

    @testset "create_galois_keys" begin
      @test_nowarn create_galois_keys!(galois_keys_, keygen)
    end

    @testset "Encryptor" begin
      @test_nowarn Encryptor(context, public_key_)
    end
    encryptor = Encryptor(context, public_key_)

    @testset "Evaluator" begin
      @test_nowarn Evaluator(context)
    end
    evaluator = Evaluator(context)

    @testset "Decryptor" begin
      @test_nowarn Decryptor(context, secret_key_)
    end
    decryptor = Decryptor(context, secret_key_)

    @testset "CKKSEncoder" begin
      @test_nowarn CKKSEncoder(context)
    end
    encoder = CKKSEncoder(context)

    slot_count_ = 4096
    @testset "slot_count" begin
      @test slot_count(encoder) == slot_count_
    end

    input = collect(range(0.0, 1.0, length=slot_count_))
    plain = Plaintext()
    encrypted = Ciphertext()
    initial_scale = 2.0^50
    @testset "encode and encrypt" begin
      @test_nowarn encode!(plain, input, initial_scale, encoder)
      @test_nowarn encrypt!(encrypted, plain, encryptor)
    end

    rotated = Ciphertext()
    @testset "rotate 2 steps left" begin
      @test_nowarn rotate_vector!(rotated, encrypted, 2, galois_keys_, evaluator)
    end

    result = similar(input)
    @testset "decrypt and decode" begin
      @test_nowarn decrypt!(plain, rotated, decryptor)
      @test_nowarn decode!(result, plain, encoder)
    end

    @testset "compare results" begin
      true_result = circshift(input, -2)
      @test isapprox(result, true_result)
    end
  end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: ec2_instance_connect
using AWS.Compat
using AWS.UUIDs

"""
    send_serial_console_sshpublic_key(instance_id, sshpublic_key)
    send_serial_console_sshpublic_key(instance_id, sshpublic_key, params::Dict{String,<:Any})

Pushes an SSH public key to the specified EC2 instance. The key remains for 60 seconds,
which gives you 60 seconds to establish a serial console connection to the instance using
SSH. For more information, see EC2 Serial Console in the Amazon EC2 User Guide.

# Arguments
- `instance_id`: The ID of the EC2 instance.
- `sshpublic_key`: The public key material. To use the public key, you must have the
  matching private key. For information about the supported key formats and lengths, see
  Requirements for key pairs in the Amazon EC2 User Guide.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"SerialPort"`: The serial port of the EC2 instance. Currently only port 0 is supported.
  Default: 0
"""
send_serial_console_sshpublic_key(InstanceId, SSHPublicKey; aws_config::AbstractAWSConfig=global_aws_config()) = ec2_instance_connect("SendSerialConsoleSSHPublicKey", Dict{String, Any}("InstanceId"=>InstanceId, "SSHPublicKey"=>SSHPublicKey); aws_config=aws_config)
send_serial_console_sshpublic_key(InstanceId, SSHPublicKey, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ec2_instance_connect("SendSerialConsoleSSHPublicKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("InstanceId"=>InstanceId, "SSHPublicKey"=>SSHPublicKey), params)); aws_config=aws_config)

"""
    send_sshpublic_key(availability_zone, instance_id, instance_osuser, sshpublic_key)
    send_sshpublic_key(availability_zone, instance_id, instance_osuser, sshpublic_key, params::Dict{String,<:Any})

Pushes an SSH public key to the specified EC2 instance for use by the specified user. The
key remains for 60 seconds. For more information, see Connect to your Linux instance using
EC2 Instance Connect in the Amazon EC2 User Guide.

# Arguments
- `availability_zone`: The Availability Zone in which the EC2 instance was launched.
- `instance_id`: The ID of the EC2 instance.
- `instance_osuser`: The OS user on the EC2 instance for whom the key can be used to
  authenticate.
- `sshpublic_key`: The public key material. To use the public key, you must have the
  matching private key.

"""
send_sshpublic_key(AvailabilityZone, InstanceId, InstanceOSUser, SSHPublicKey; aws_config::AbstractAWSConfig=global_aws_config()) = ec2_instance_connect("SendSSHPublicKey", Dict{String, Any}("AvailabilityZone"=>AvailabilityZone, "InstanceId"=>InstanceId, "InstanceOSUser"=>InstanceOSUser, "SSHPublicKey"=>SSHPublicKey); aws_config=aws_config)
send_sshpublic_key(AvailabilityZone, InstanceId, InstanceOSUser, SSHPublicKey, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ec2_instance_connect("SendSSHPublicKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AvailabilityZone"=>AvailabilityZone, "InstanceId"=>InstanceId, "InstanceOSUser"=>InstanceOSUser, "SSHPublicKey"=>SSHPublicKey), params)); aws_config=aws_config)


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã ``2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Tree-based grammar guided genetic programming that uses a context-free
grammar in BNF.
"""
module GP

export GPESParams, GPESResult, gp_search, exprsearch, SearchParams, SearchResult, get_derivtree
export GPPopulation, GPIndividual
export crossover, mutate, max_depth

using Reexport
using ExprSearch, MinDepths
using RLESUtils, GitUtils, CPUTimeUtils, RandUtils, Observers, LogSystems, TreeIterators, TreeUtils
using GrammaticalEvolution
@reexport using DerivationTrees  #for pretty strings
using CPUTime
using DerivTreeVis
using RLESTypes.SymbolTable

import Compat.view
import DerivationTrees: get_children, max_depth
import ..ExprSearch: exprsearch, get_fitness, get_derivtree, get_expr
import Base: copy!, sort!, isless, resize!, empty!, push!, pop!, length, start, next, done,
    rand, getindex

immutable RuleNotFoundException <: Exception end
immutable DepthExceededException <: Exception end

const DEFAULT_EXPR = :()

include("logdefs.jl")

immutable GPESParams <: SearchParams
    pop_size::Int64
    maxdepth::Int64
    iterations::Int64
    tournament_size::Int64
    top_keep::Float64
    crossover_frac::Float64
    mutate_frac::Float64
    rand_frac::Float64
    default_expr
    logsys::LogSystem
    userargs::SymbolTable #passed to get_fitness if specified
end
GPESParams(pop_size::Int64, maxdepth::Int64, iterations::Int64, tournament_size::Int64,
    top_keep::Float64, crossover_frac::Float64, mutate_frac::Float64, rand_frac::Float64,
    default_expr, logsys::LogSystem=logsystem(); userargs::SymbolTable=SymbolTable()) =
        GPESParams(pop_size, maxdepth, iterations, tournament_size, top_keep,
        crossover_frac, mutate_frac, rand_frac, default_expr, logsys, userargs)

type GPESResult <: SearchResult
    tree::DerivationTree
    fitness::Float64
    expr
    best_at_eval::Int64
    totalevals::Int64
end
GPESResult(grammar::Grammar) = GPESResult(DerivationTree(DerivTreeParams(grammar)), realmax(Float64),
    DEFAULT_EXPR, 0, 0)

exprsearch(p::GPESParams, problem::ExprProblem) = gp_search(p, problem::ExprProblem)

get_derivtree(result::GPESResult) = result.tree
get_expr(result::GPESResult) = result.expr
get_fitness(result::GPESResult) = result.fitness

type GPIndividual
    derivtree::DerivationTree
    expr
    fitness::Nullable{Float64}
end
GPIndividual(grammar::Grammar) = GPIndividual(DerivationTree(DerivTreeParams(grammar)),
    DEFAULT_EXPR, Nullable{Float64}())

type GPPopulation
    individuals::Vector{GPIndividual}
end
GPPopulation() = GPPopulation(Array(GPIndividual, 0))

function gp_search(p::GPESParams, problem::ExprProblem)
    @notify_observer(p.logsys.observer, "verbose1", ["Starting GP search"])
    @notify_observer(p.logsys.observer, "computeinfo", ["starttime", string(now())])

    initialize!(problem)
    grammar = get_grammar(problem)
    mdr = min_depth_rule(grammar)
    mda = min_depth_actions(mdr, grammar)
    result = GPESResult(grammar)

    tstart = CPUtime_start()

    #first iter is ramped init
    iter = 1
    pop = GPPopulation()
    ramped_initialize!(pop, grammar, mdr, mda, p.pop_size, p.maxdepth)
    parallel_evaluate(p, pop, result, problem, p.default_expr)
    sort!(pop)
    fitness = get(pop[1].fitness)
    code = string(pop[1].expr)
    nevals = iter * p.pop_size
    @notify_observer(p.logsys.observer, "elapsed_cpu_s", [nevals, CPUtime_elapsed_s(tstart)])
    @notify_observer(p.logsys.observer, "fitness", Any[iter, fitness])
    @notify_observer(p.logsys.observer, "code", Any[iter, code])
    @notify_observer(p.logsys.observer, "population", Any[iter, pop])
    @notify_observer(p.logsys.observer, "current_best", [nevals, fitness, code])
    iter += 1

    while iter <= p.iterations
        pop = generate(p, grammar, mda, pop, result, problem)
        ind = best_ind(pop)
        fitness = get(ind.fitness)
        code = string(ind.expr)
        nevals = iter * p.pop_size
        @notify_observer(p.logsys.observer, "elapsed_cpu_s", [nevals, CPUtime_elapsed_s(tstart)])
        @notify_observer(p.logsys.observer, "fitness", Any[iter, fitness])
        @notify_observer(p.logsys.observer, "code", Any[iter, code])
        @notify_observer(p.logsys.observer, "population", Any[iter, pop])
        @notify_observer(p.logsys.observer, "current_best", [nevals, fitness, code])
        iter += 1
    end
    @assert result.fitness <= get(best_ind(pop).fitness)

    #dealloc pop
    for ind in pop
        rm_tree!(ind.derivtree)
    end

    @notify_observer(p.logsys.observer, "result", [result.fitness, string(result.expr),
        result.best_at_eval, result.totalevals])

    #meta info
    @notify_observer(p.logsys.observer, "computeinfo", ["endtime",  string(now())])
    @notify_observer(p.logsys.observer, "computeinfo", ["hostname", gethostname()])
    @notify_observer(p.logsys.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
    @notify_observer(p.logsys.observer, "computeinfo", ["cpu_time", CPUtime_elapsed_s(tstart)])
    @notify_observer(p.logsys.observer, "parameters", ["pop_size", p.pop_size])
    @notify_observer(p.logsys.observer, "parameters", ["maxdepth", p.maxdepth])
    @notify_observer(p.logsys.observer, "parameters", ["iterations", p.iterations])
    @notify_observer(p.logsys.observer, "parameters", ["tournament_size", p.tournament_size])
    @notify_observer(p.logsys.observer, "parameters", ["top_keep", p.top_keep])
    @notify_observer(p.logsys.observer, "parameters", ["crossover_frac", p.crossover_frac])
    @notify_observer(p.logsys.observer, "parameters", ["mutate_frac", p.mutate_frac])
    @notify_observer(p.logsys.observer, "parameters", ["rand_frac", p.rand_frac])
    @notify_observer(p.logsys.observer, "parameters", ["default_expr", string(p.default_expr)])

    result
end

start(pop::GPPopulation) = start(pop.individuals)
next(pop::GPPopulation, s) = next(pop.individuals, s)
done(pop::GPPopulation, s) = done(pop.individuals, s)

function isless(x::GPIndividual, y::GPIndividual)
   get(x.fitness, realmax(eltype(x.fitness))) < get(y.fitness, realmax(eltype(y.fitness)))
end

function best_ind(pop::GPPopulation)
    @assert !isnull(pop.individuals[1].fitness)
    @assert pop.individuals[1] <= pop.individuals[end]
    pop.individuals[1]
end

function resize!(pop::GPPopulation, N::Int64)
    resize!(pop.individuals, N)
end

pop!(pop::GPPopulation) = pop!(pop.individuals)
push!(pop::GPPopulation, ind::GPIndividual) = push!(pop.individuals, ind)
empty!(pop::GPPopulation) = empty!(pop.individuals)
length(pop::GPPopulation) = length(pop.individuals)
getindex(pop::GPPopulation, index) = getindex(pop.individuals, index)

function copy!(dst::GPIndividual, src::GPIndividual)
   copy!(dst.derivtree, src.derivtree)
   dst.expr = src.expr
   dst.fitness = src.fitness #isbits, so can just copy
end

max_depth(ind::GPIndividual) = max_depth(ind.derivtree)

"""
Initialize population.  Ramped initialization.
"""
function ramped_initialize!(pop::GPPopulation, grammar::Grammar, mdr::MinDepthByRule,
    mda::MinDepthByAction, pop_size::Int64, maxdepth::Int64)

    empty!(pop)
    startmindepth = mdr[:start]
    for d in cycle(startmindepth:maxdepth)
        try
            ind = rand(grammar, mda, d)
            push!(pop, ind)
        catch e
            if !isa(e, IncompleteException)
                rethrow(e)
            end
        end
        length(pop) < pop_size || break
    end
    pop
end

function resample_subtree!(tree::DerivationTree, node::DerivTreeNode, mda::MinDepthByAction,
    target_depth::Int64)
    rm_node(tree, node.children) #return child nodes
    opennodes = DerivTreeNode[node]
    tree.nopen = 1
    rand_subtree!(tree, node, opennodes, mda, target_depth)
end

"""
Random individual to a target_depth
"""
function rand(grammar::Grammar, mda::MinDepthByAction, target_depth::Int64)
    opennodes = DerivTreeNode[]
    ind = GPIndividual(grammar)
    tree = ind.derivtree
    node = initialize!(tree)
    push!(opennodes, node)
    rand_subtree!(tree, node, opennodes, mda, target_depth)
    ind
end

function rand_subtree!(tree::DerivationTree, node::DerivTreeNode, opennodes::Vector{DerivTreeNode},
    mda::MinDepthByAction, target_depth::Int64)
    while !isempty(opennodes)
        node = pop!(opennodes)
        if is_decision(node)
            depths = mda[Symbol(node.rule.name)]
            actions = find(x->x <= target_depth-node.depth, depths)
            if isempty(actions)
                warn("init(): no valid actions")
                throw(IncompleteException())
            end
            a = rand_element(actions)
            nodes = expand_node!(tree, node, a)
        else
            nodes = expand_node!(tree, node)
        end
        append!(opennodes, nodes)
    end
end
"""
Evaluate fitness function
"""
function seq_evaluate(p::GPESParams, pop::GPPopulation, result::GPESResult, problem::ExprProblem,
    default_expr)
    for ind in pop
        if isnull(ind.fitness)
            try
                ind.expr = get_expr(ind.derivtree)
                fitness = get_fitness(problem, ind.derivtree, p.userargs)
                ind.fitness = Nullable{Float64}(fitness)
                result.totalevals += 1
                if fitness < result.fitness
                    result.fitness = fitness
                    copy!(result.tree, ind.derivtree)
                    result.best_at_eval = result.totalevals
                    result.expr = ind.expr
                end
            catch e
                if !isa(e, IncompleteException)
                    rethrow(e)
                end
                ind.expr = default_expr
                ind.fitness = Nullable{Float64}(realmax(Float64))
            end
        end
    end
end

"""
Threaded evaluation of fitnesses
"""
function parallel_evaluate(p::GPESParams, pop::GPPopulation, result::GPESResult,
    problem::ExprProblem, default_expr)
    #evaluate fitness in parallel
    for i = 1:length(pop) #use this to disable threads
    #Threads.@threads for i = 1:length(pop)
        ind = pop[i]
        if isnull(ind.fitness)
            try
                fitness = get_fitness(problem, ind.derivtree, p.userargs)
                ind.fitness = Nullable{Float64}(fitness)
                ind.expr = get_expr(ind.derivtree)
            catch e
                if !isa(e, IncompleteException)
                    println("Exception caught! ", e)
                    rethrow(e)
                end
                ind.fitness = Nullable{Float64}(realmax(Float64))
            end
        end
    end
    #update global result
    for i = 1:length(pop)
        ind = pop[i]
        fitness = get(ind.fitness)
        result.totalevals += 1
        if fitness < result.fitness
            result.fitness = fitness
            copy!(result.tree, ind.derivtree)
            result.best_at_eval = result.totalevals
            result.expr = get_expr(ind.derivtree)
        end
    end
end

"""
Tournament selection
Best fitness wins deterministically
"""
function selection(pop::GPPopulation, tournament_size::Int64)
    #randomly choose N inds and compare their fitness
    ids = randperm(length(pop))
    #sorted, just return the one with the lowest index
    id = minimum(view(ids, 1:tournament_size))
    pop[id]
end

TreeUtils.get_children(node::DerivTreeNode) = DerivationTrees.get_children(node)

function findbyrule(ind::GPIndividual, name::AbstractString)
    candidates = DerivTreeNode[]
    for node in tree_iter(ind.derivtree.root) #TODO: replace with AbstractTrees.jl interface
        if node.rule.name == name
            push!(candidates, node)
        end
    end
    candidates
end

"""
Crossover
"""
function crossover(ind1::GPIndividual, ind2::GPIndividual, grammar::Grammar, maxdepth::Int64)
    ind3 = GPIndividual(grammar)
    ind4 = GPIndividual(grammar)
    copy!(ind3.derivtree, ind1.derivtree)
    copy!(ind4.derivtree, ind2.derivtree)

    node3 = rand_node(n->!is_terminal(n), ind3.derivtree.root) #choose crossover point on ind3
    candidates = findbyrule(ind4, node3.rule.name) #find corresponding on ind4
    isempty(candidates) && throw(RuleNotFoundException())
    node4 = rand_element(candidates)
    swap_children!(node3, node4)
    if max_depth(ind3) > maxdepth || max_depth(ind4) > maxdepth
        throw(DepthExceededException())
    end
    (ind3, ind4)
end #module

"""
Mutation
"""
function mutate(ind1::GPIndividual, grammar::Grammar, mda::MinDepthByAction, target_depth::Int64)
    ind2 = GPIndividual(grammar)
    copy!(ind2.derivtree, ind1.derivtree)

    #choose random node
    tree = ind2.derivtree
    node = rand_node(tree.root)
    resample_subtree!(tree, node, mda, target_depth)
    ind2
end

"""
Adjust to target pop_size
"""
function adjust_to_size!(newpop::GPPopulation, pop::GPPopulation, grammar::Grammar,
    pop_size::Int64, tournament_size::Int64)
    #reproduce if undersize
    while length(newpop) < pop_size
        ind1 = selection(pop, tournament_size)
        ind2 = GPIndividual(grammar)
        copy!(ind2, ind1)
        push!(newpop, ind2)
    end

    #trim to size
    while length(newpop) > pop_size
        ind = pop!(newpop)
        rm_tree!(ind.derivtree)
    end
    newpop
end

function sort!(pop::GPPopulation) #lowest (best) fitness first
    sort!(pop.individuals)
end

"""
Create next generation
"""
function generate(p::GPESParams, grammar::Grammar, mda::MinDepthByAction,
    pop::GPPopulation, result::GPESResult, problem::ExprProblem)

    #println("enter generate")
    n_keep = floor(Int64, p.top_keep*p.pop_size)
    n_crossover = floor(Int64, p.crossover_frac*p.pop_size)
    n_mutate = floor(Int64, p.mutate_frac*p.pop_size)
    n_rand = floor(Int64, p.rand_frac*p.pop_size)

    parallel_evaluate(p, pop, result, problem, p.default_expr)
    sort!(pop)
    newpop = GPPopulation()

    #pass through the top ones
    for i = 1:n_keep
        ind1 = GPIndividual(grammar)
        copy!(ind1, pop[i])
        push!(newpop, ind1)
    end

    #crossover
    n = 0
    while n < n_crossover
        #println("crossover: $n")
        ind1 = selection(pop, p.tournament_size)
        ind2 = selection(pop, p.tournament_size)
        try
            (ind3, ind4) = crossover(ind1, ind2, grammar, p.maxdepth)
            push!(newpop, ind3)
            push!(newpop, ind4)
            n += 2
        catch e
            if !isa(e, Union{RuleNotFoundException,DepthExceededException})
                rethrow(e)
            end
        end
    end

    #mutate
    n = 0
    while n < n_mutate
        #println("mutate: $n")
        ind1 = selection(pop, p.tournament_size)
        try
            ind2 = mutate(ind1, grammar, mda, p.maxdepth)
            push!(newpop, ind2)
            n += 1
        catch e
            if !isa(e, IncompleteException)
                rethrow(e)
            end
        end
    end

    #random
    n = 0
    while n < n_rand
        #println("random: $n")
        try
            ind1 = rand(grammar, mda, p.maxdepth)
            push!(newpop, ind1)
            n += 1
        catch e
            if !isa(e, IncompleteException)
                rethrow(e)
            end
        end
    end

    #adjust size: trim or reproduce to size
    adjust_to_size!(newpop, pop, grammar, p.pop_size, p.tournament_size)

    #dealloc old pop
    for ind in pop
        rm_tree!(ind.derivtree)
    end

    parallel_evaluate(p, newpop, result, problem, p.default_expr)
    sort!(newpop)

    newpop
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
module LambdaTAME

using LinearAlgebra
using NPZ
using JSON
using Distributed
using SparseArrays
using Random
using MatrixNetworks
using DataStructures
using NearestNeighbors
using Parameters

using Distributions
using Metis

#  https://github.com/charlescolley/DistributedTensorConstruction.jl.git
using DistributedTensorConstruction
# imports SymTensorUnweighted type, tensors_from_graph, contraction_divide_out!
import DistributedTensorConstruction: contraction

using NetworkAlign

import Statistics: mean
import Arpack: svds
import Combinatorics: permutations, combinations
import StatsBase: sample

#store path to project repo here
PROJECT_PATH = "."

abstract type returnType end

# https://github.com/eikmeier/TGPA.git
#TGPA_PATH = "Users/ccolley/Code/TGPA"
#include(TGPA_PATH*"/TGPA_generator.jl")

#
# https://github.com/nassarhuda/NetworkAlignment.jl.git
#NETWORK_ALIGNMENT_PATH = "/Users/charlie/Documents/Code/NetworkAlignment.jl/src/"
#NETWORK_ALIGNMENT_PATH = "/homes/ccolley/Documents/Software/NetworkAlignment.jl/src/" #nilpotent
#include(NETWORK_ALIGNMENT_PATH*"NetworkAlignment.jl")

#LOWRANK_EIGENALIGN_PATH = "/Users/charlie/Documents/Code/lowrank_spectral_v1_julia/all_code/include_all.jl" #local
#LOWRANK_EIGENALIGN_PATH = "/homes/ccolley/Documents/Software/lowrank_spectral_v1_julia/all_code/include_all.jl"#nilpotent
#include(LOWRANK_EIGENALIGN_PATH)
include("../include/lowrank_spectral/align_networks_eigenalign.jl")
#TODO: fix 'WARNING: Method definition' messages

#TODO: Check to see if this can be run straight out of the box.
#local graph repos
MULTIMAGNA=PROJECT_PATH*"/data/sparse_tensors/"

#no checks are run, will lead to undefined behavior if variables aren't as specified
struct ThirdOrderSymTensor
    n::Int                  #assuming > 0
    indices::Array{Int,2}   #assuming no repeated permutations, i != k , k != j , i != j
    values::Array{Float64,1}
end

struct UnweightedThirdOrderSymTensor #hyperedge weights are 1.0
    n::Int                  #assuming > 0
    indices::Vector{Vector{Tuple{Int,Int}}}   #assuming no repeated permutations, i != k , k != j , i != j
end

#TODO:
#   Note that to use SparseSymmetricTensors.jl, must update experiment drivers to
#   take a flag that use the local load function

include("Contraction.jl")
include("Matchings.jl")
include("PostProcessing.jl")
include("AlignmentDrivers.jl")
include("RandomGraphs.jl")
include("fileio.jl")
include("Experiments.jl")
#include("TAME_Implementations.jl")
include("LambdaTAME_implementation.jl")
include("LowRankTAME_implementation.jl")
include("TAME_implementation.jl")
#include("Experimental_code.jl")
include("AdditionalAlignmentMethods.jl")

#include("SparseSymmetricTensorCode.jl")

#  --  From Experiments.jl --  #

#all alignment method flag types
export ΛTAME_M, ΛTAME_MultiMotif_M, LowRankTAME_M, LowRankTAME_MultiMotif_M, TAME_M, TAME_MultiMotif_M, EigenAlign_M, LowRankEigenAlign_M, LowRankEigenAlignOnlyEdges_M, Degree_M ,Random_M
export ΛTAME_GramMatching, ΛTAME_rankOneMatching
#random graph flag types
export ErdosRenyi, RandomGeometric, HyperKron
#noise model flag types
export ErdosRenyiNoise, DuplicationNoise

export distributed_pairwise_alignment
export distributed_random_trials
export align_tensors # this is the file based version
export align_matrices
export random_graph_exp

export load_ThirdOrderSymTensor, load_UnweightedThirdOrderSymTensor
export graph_to_ThirdOrderTensor

export KlauAlgo, noPostProcessing, SuccessiveKlauAlgo, LocalSearch

#  --  From Matchings.jl --  #
#flag types
export returnTimings,noTimings

export bipartite_matching_primal_dual

#  --  From TAME_Implementations.jl --  #
export align_tensors_profiled
export ΛTAME
export LowRankTAME
export LowRankTAME_profiled
export TAME
export TAME_profiled

end #module

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
module tiddlywinks

import tiddlywinks

include("functions.jl")
export tiddly_greet
# Write your package code here.

end


text:
module PDE

using GPUBenchmarks, BenchmarkTools, Primes

description = """
Kuramoto-Sivashinsky algorithm benchmark ([original benchmark](https://github.com/johnfgibson/julia-pde-benchmark/blob/master/1-Kuramoto-Sivashinksy-benchmark.ipynb)).

This benchmark is dominated by the cost of the FFT, leading to worse results for OpenCL with
CLFFT compared to the faster CUFFT.
Similarly the multithreaded backend doesn't improve much over base with the same FFT implementation.
Result of the benchmarked PDE:
![](https://github.com/JuliaGPU/GPUBenchmarks.jl/blob/master/results/plots/pde_result.png?raw=true)
"""

sq3(x) = Complex64(x * x * x) # cudanative doesn't like Complex64 directly broadcasted.

# source: http://nbviewer.jupyter.org/url/homepages.warwick.ac.uk/staff/C.Ortner/julia/PlaneWaves.ipynb
# an optimised implementation of CH_vectorised!

function inner_loop(
        nruns,
        planv, planw, planwi, u0, u, v, w, Â,
        c1, c2, c3, c4
    )
    u .= identity.(u0)
    w .= tocomplex64.(u0)
    for n = 1:nruns
        v .= sq3.(u)
        planv * v
        planw * w
        @. w = ((1f0+c1*Â)*w - (c2*Â)*v) / ((1f0+c3*Â+c4)*Â)
        planwi * w
        u .= real.(w)
    end
    synchronize(u)
    return u
end

function makeresult(bench, N, device, hardware, mdiff)
    BenchResult(
        "PDE",
        bench,
        N,
        Float32,
        string(device),
        hardware,
        @__FILE__,
        mdiff
    )
end

function execute2(device)
    hardware, array_type = init(device)
    T = Float32
    results = BenchResult[]
    for x in (4,)
        N = 2 ^ x
        # initialisations
        h = 2*π/N; epsn = h * 3; C = 2/epsn; tau = epsn * h
        k = [0:N/2; -N/2+1:-1]
        acpu = kron(k.^2, ones(1,N)) + kron(ones(N), k'.^2)
        ucpu = 2*(rand(N, N)-0.5)

        Â = array_type(T.(acpu))
        u = array_type(T.(ucpu))

        w = tocomplex64.(u)
        v = copy(w)
        c1 = (C*tau + tau/epsn)
        c2 = (tau / epsn)
        c3 = (epsn * tau)
        c4 = (C*tau)
        planv = plan_fft!(v)
        planw = plan_fft!(w)
        planwi = plan_ifft!(w)
        nruns = 5
        u0 = copy(u)
        bench = @benchmark $(inner_loop)(
                $nruns,
                $planv, $planw, $planwi,$u0, $u, $v, $w, $Â,
                $c1, $c2, $c3, $c4
        )
        wcpu = tocomplex64.(ucpu)
        vcpu = copy(wcpu)
        planvcpu = plan_fft!(vcpu)
        planwcpu = plan_fft!(wcpu)
        planwicpu = plan_ifft!(wcpu)
        inner_loop(
                nruns,
                planvcpu, planwcpu, planwicpu, copy(ucpu), ucpu, vcpu, wcpu, acpu,
                c1, c2, c3, c4
        )
        println(Array(ucpu))
        # mdiff = meandifference(ucpu, u)
        # @show mdiff
        # push!(results, makeresult(bench, N, device, hardware, mdiff))
        # free(Â); free(u0); free(u); free(v); free(w)
    end
    results
end

# execute(Float32, 512, :opencl)

function CH_memory_af!(nruns, N, AT)
    # initialisations
    h = Float32(2*π/N); epsn = Float32(h * 3); C = Float32(2/epsn); tau = Float32(epsn * h)
    k = [0:N/2; -N/2+1:-1]
    Â = AT((Float32.(kron(k.^2, ones(1,N)) + kron(ones(N), k'.^2))))
    u = AT(Float32.((2*(rand(N, N)-0.5))))

    # ============= ACTUAL CODE THAT IS BEING TESTED ======================
    # allocate arrays and define constants
    w = AT{Complex64}(u)
    v = copy(w)
    c1 = (C*tau + tau/epsn)
    c2 = (tau / epsn)
    c3 = (epsn * tau)
    c4 = (C*tau)
    tic()
    for n = 1:nruns
        v .= complex.(u .* u .* u)
        v = fft(v)
        w = fft(w)
        w .= (( (1f0 + c1) .* Â) .* w .- (c2 .* Â) .* v) ./ ((1f0 + c3 .* Â .+ c4) .* Â)
        w = ifft(w)
        u .= real.(w)
    end
    GPUArrays.synchronize(u)
    toc()
    # ======================================================================
    u
end

# source: https://github.com/johnfgibson/julia-pde-benchmark/blob/master/1-Kuramoto-Sivashinksy-benchmark.ipynb
function inner_ks(IFFT!, FFT!, Nt, Nn, Nn1, u, G, A_inv, B, dt2, dt32)
    for n = 1:Nt
        Nn1 .= Nn       # shift nonlinear term in time
        Nn .= u         # put u into Nn in prep for comp of nonlinear term

        IFFT! * Nn
            # transform Nn to gridpt values, in place
        Nn .= Nn .* Nn   # collocation calculation of u^2
        FFT!*Nn        # transform Nn back to spectral coeffs, in place

        Nn .= G .* Nn    # compute Nn == -1/2 d/dx (u^2) = -u u_x

        # loop fusion! Julia translates the folling line of code to a single for loop.
        u .= A_inv .* (B .* u .+ dt32 .* Nn .- dt2 .* Nn1)
    end
    synchronize(u)
end

function makeresult(bench, N, device, hardware, mdiff)
    BenchResult(
        "PDE",
        bench,
        N,
        Complex64,
        string(device),
        hardware,
        @__FILE__,
        mdiff
    )
end
function execute(device)
    hardware, AT = init(device)
    results = BenchResult[]
    is_gpuarrays(device) || device == :julia_base || return results
    T = Float32
    for i = 1:7
        N = 10^i
        Lx = T(64*pi)
        Nx = T(N)
        dt = T(1/16)
        Nt = 50

        x = Lx*(0:Nx-1)/Nx
        u = T.(cos.(x) + 0.1*sin.(x/8) + 0.01*cos.((2*pi/Lx)*x))

        u = AT((T(1)+T(0)im)*u)             # force u to be complex
        Nx = length(u)                      # number of gridpoints
        kx = T.(vcat(0:Nx/2-1, 0:0, -Nx/2+1:-1))# integer wavenumbers: exp(2*pi*kx*x/L)
        alpha = T(2)*pi*kx/Lx                  # real wavenumbers:    exp(alpha*x)

        D = T(1)im*alpha                       # spectral D = d/dx operator

        L = alpha.^2 .- alpha.^4            # spectral L = -D^2 - D^4 operator

        G = AT(T(-0.5) .* D)               # spectral -1/2 D operator, to eval -u u_x = 1/2 d/dx u^2

        # convenience variables
        dt2  = T(dt/2)
        dt32 = T(3*dt/2)
        A_inv = AT((ones(T, Nx) - dt2*L).^(-1))
        B = AT(ones(T, Nx) + dt2*L)

        # compute in-place FFTW plans
        FFT! = plan_fft!(u)
        IFFT! = plan_ifft!(u)

        # compute nonlinear term Nn == -u u_x
        powed = u .* u
        Nn = G .* fft(powed);    # Nn == -1/2 d/dx (u^2) = -u u_x
        Nn1 = copy(Nn);        # Nn1 = Nn at first time step
        FFT! * u;

        # timestepping loop
        bench = @benchmark $inner_ks($(IFFT!), $(FFT!), $Nt, $Nn, $Nn1, $u, $G, $A_inv, $B, $dt2, $dt32)
        push!(results, makeresult(bench, N, device, hardware, 0.0))
    end
    results
end

function ksintegrateNaive(u, Lx, dt, Nt, nsave)
    Nx = length(u)                  # number of gridpoints
    x = collect(0:(Nx-1)/Nx)*Lx
    kx = vcat(0:Nx/2-1, 0, -Nx/2+1:-1)  # integer wavenumbers: exp(2*pi*kx*x/L)
    alpha = 2*pi*kx/Lx              # real wavenumbers:    exp(alpha*x)
    D = 1im*alpha;                  # D = d/dx operator in Fourier space
    L = alpha.^2 - alpha.^4         # linear operator -D^2 - D^4 in Fourier space
    G = -0.5*D                      # -1/2 D operator in Fourier space

    Nsave = div(Nt, nsave)+1        # number of saved time steps, including t=0
    t = (0:Nsave)*(dt*nsave)        # t timesteps
    U = zeros(Nsave, Nx)            # matrix of u(xⱼ, tᵢ) values
    U[1,:] = u                      # assign initial condition to U
    s = 2                           # counter for saved data

    dt2  = dt/2
    dt32 = 3*dt/2;
    A_inv = (ones(Nx) - dt2*L).^(-1)
    B     =  ones(Nx) + dt2*L

    Nn  = G.*fft(u.*u) # -u u_x (spectral), notation Nn = N^n     = N(u(n dt))
    Nn1 = copy(Nn)     #                   notation Nn1 = N^{n-1} = N(u((n-1) dt))
    u  = fft(u)        # transform u to spectral

    # timestepping loop
    for n = 1:Nt
        Nn1 = copy(Nn)                 # shift nonlinear term in time: N^{n-1} <- N^n
        Nn  = G.*fft(real(ifft(u)).^2) # compute Nn = -u u_x

        u = A_inv .* (B .* u + dt32*Nn - dt2*Nn1)

        if mod(n, nsave) == 0
            U[s,:] = real(ifft(u))
            s += 1
        end
    end
    t,U
end
#
# using FileIO, Interpolations, Colors, GPUBenchmarks, GPUArrays, ColorVectorSpace, FixedPointNumbers
# Lx = 64*pi
# Nx = 1024
# dt = 1/16
# nsave = 8
# Nt = 3200
#
# x = Lx*(0:Nx-1)/Nx
# u = cos.(x) + 0.1*sin.(x/8) + 0.01*cos.((2*pi/Lx)*x);
# t,U = ksintegrateNaive(u, Lx, dt, Nt, nsave)
# cn = 100
# cmap = interpolate(colormap("Oranges", cn), BSpline(Linear()), OnCell());
# mini, maxi = extrema(U)
# img_color = map(U) do val
#     val = (val - mini) / (maxi - mini)
#     val = 1 - clamp(val, 0f0, 1f0);
#     idx = (val * (cn - 1)) + 1.0
#     RGB{N0f8}(cmap[idx])
# end
# save(GPUBenchmarks.dir("results", "plots", "pde_result.png"), img_color)

end


text:
@testset "ObsView constructor" begin
    @test_throws DimensionMismatch ObsView((rand(2,10), rand(9)))
    @test_throws DimensionMismatch ObsView((rand(2,10), rand(9)), 1:2)
    @test_throws DimensionMismatch ObsView((rand(2,10), rand(4,9,10), rand(9)))

    @testset "bounds check" begin
        for var in (vars..., tuples..., CustomType())
            @test_throws BoundsError ObsView(var, -1:100)
            @test_throws BoundsError ObsView(var, 1:16)
            @test_throws BoundsError ObsView(var, [1, 10, 0, 3])
            @test_throws BoundsError ObsView(var, [1, 10, -10, 3])
            @test_throws BoundsError ObsView(var, [1, 10, 18, 3])
        end
    end

    @testset "Tuples" begin
        @test typeof(@inferred(ObsView((X,X)))) <: ObsView
        @test typeof(@inferred(ObsView((X,X), 1:15))) <: ObsView
        d = ObsView((X, y), 5:10)
        @test @inferred(getobs(d , 2)) == getobs((X, y), 6)
        @test numobs(d) == 6
    end

    @testset "Array, SubArray, SparseArray" begin
        for var in (Xs, ys, vars...)
            subset = @inferred(ObsView(var))
            @test subset.data === var
            @test subset.indices === 1:15
            @test typeof(subset) <: ObsView
            @test @inferred(numobs(subset)) === numobs(var)
            @test @inferred(getobs(subset)) == getobs(var)
            @test @inferred(ObsView(subset)) === subset
            @test @inferred(ObsView(subset, 1:15)) === subset
            @test subset[end] == obsview(var, 15)
            @test @inferred(subset[15]) == obsview(var, 15)
            @test @inferred(subset[2:5]) == obsview(var, 2:5)
            for idx in (1:10, [1,10,15,3], [2])
                @test ObsView(var)[idx] == obsview(var, idx)
                @test ObsView(var)[idx] == obsview(var, collect(idx))
                subset = @inferred(ObsView(var, idx))
                @test typeof(subset) <: ObsView{typeof(var), typeof(idx)}
                @test subset.data === var
                @test subset.indices === idx
                @test @inferred(numobs(subset)) === length(idx)
                @test @inferred(getobs(subset)) == getobs(var, idx)
                @test @inferred(ObsView(subset)) === subset
                @test @inferred(subset[1]) == obsview(var, idx[1])
                @test numobs(subset[1:1]) == numobs(ObsView(var, obsview(idx, 1:1)))
            end
        end
    end

    @testset "custom types" begin
        @test_throws MethodError ObsView(EmptyType())
        @test_throws MethodError ObsView(EmptyType(), 1:10)
        @test_throws BoundsError getobs(ObsView(CustomType(), 11:20), 11)
        @test typeof(@inferred(ObsView(CustomType()))) <: ObsView
        @test numobs(ObsView(CustomType())) === 15
        @test numobs(ObsView(CustomType(), 1:10)) === 10
        @test getobs(ObsView(CustomType())) == collect(1:15)
        @test getobs(ObsView(CustomType(),1:10),10) == 10
        @test getobs(ObsView(CustomType(),1:10),[3,5]) == [3,5]
    end
end

@testset "ObsView getindex and getobs" begin
    @testset "Matrix and SubArray{T,2}" begin
        for var in (X, Xv)
            subset = @inferred(ObsView(var, 5:12))
            @test typeof(@inferred(getobs(subset))) <: Array{Float64,2}
            @test @inferred(numobs(subset)) == length(subset) == 8
            @test @inferred(subset[2:5]) == obsview(X, 6:9)
            @test @inferred(subset[3:6]) != obsview(X, 6:9)
            @test @inferred(getobs(subset, 2:5)) == X[:, 6:9]
            @test @inferred(getobs(subset, [3,1,4])) == X[:, [7,5,8]]
            @test typeof(subset[2:5]) <: SubArray
            @test @inferred(subset[collect(2:5)]) == obsview(X, collect(6:9))
            @test typeof(subset[collect(2:5)]) <: SubArray
            @test @inferred(getobs(subset)) == getobs(subset[1:end]) == X[:, 5:12]
        end
    end

    @testset "Vector and SubArray{T,1}" begin
        for var in (y, yv)
            subset = @inferred(ObsView(var, 6:10))
            @test typeof(getobs(subset)) <: Array{String,1}
            @test @inferred(numobs(subset)) == length(subset) == 5
            @test @inferred(subset[2:3]) == obsview(y, 7:8)
            @test @inferred(getobs(subset, 2:3)) == y[7:8]
            @test @inferred(getobs(subset, [2,1,4])) == y[[7,6,9]]
            @test typeof(subset[2:3]) <: SubArray
            @test @inferred(subset[collect(2:3)]) == obsview(y, collect(7:8))
            @test typeof(subset[collect(2:3)]) <: SubArray
            @test @inferred(getobs(subset)) == getobs(subset[1:end]) == y[6:10]
        end
    end
end

@testset "getobs!" begin
    @test getobs!(nothing, ObsView(y, 1)) == "setosa"

    @testset "ObsView" begin
        xbuf1 = zeros(4,8)
        s1 = ObsView(X, 2:9)
        @test @inferred(getobs!(xbuf1,s1)) === xbuf1
        @test xbuf1 == getobs(s1)
        xbuf1 = zeros(4,5)
        s1 = ObsView(X, 5:12)
        @test @inferred(getobs!(xbuf1,s1,2:6)) === xbuf1
        @test xbuf1 == getobs(s1,2:6) == getobs(X,6:10)

        s3 = ObsView(Xs, 6:10)
        @test @inferred(getobs!(nothing,s3)) == getobs(Xs,6:10)

        s4 = ObsView(CustomType(), 6:10)
        @test @inferred(getobs!(nothing,s4)) == getobs(s4)
        s5 = ObsView(CustomType(), 3:10)
        @test @inferred(getobs!(nothing,s5,2:6)) == getobs(s5,2:6)
    end

    @testset "Tuple with ObsView" begin
        xbuf = zeros(4,2)
        ybuf = ["foo", "bar"]
        s1 = ObsView(Xs, 5:9)
        s2 = ObsView(X, 5:9)
        @test getobs!((nothing,xbuf),(s1,s2), 2:3) == (getobs(Xs,6:7),xbuf)
        @test xbuf == getobs(X,6:7)
        @test getobs!((nothing,xbuf),(s1,s2), 2:3) == (getobs(Xs,6:7),xbuf)
        @test getobs!((nothing,xbuf),(s1,s2), 2:3) == (getobs(Xs,6:7),xbuf)
    end
end

@testset "ObsView other" begin

    @testset "constructor" begin
        # @test_throws DimensionMismatch ObsView((rand(2,10),rand(9)))
        # @test_throws DimensionMismatch ObsView((rand(2,10),rand(4,9,10),rand(9)))
        @test_throws MethodError ObsView(EmptyType())
        @test_throws MethodError ObsView((EmptyType(),EmptyType()))
        @test_throws MethodError ObsView((EmptyType(),EmptyType()))
        for var in (vars..., Xs, ys)
            A = @inferred(ObsView(var))
        end
        for var in (vars..., tuples..., Xs, ys)
            A = ObsView(var)
            @test @inferred(parent(A)) === var
            @test @inferred(ObsView(A)) == A
            @test @inferred(ObsView(var)) == A
        end
    end

    @testset "typestability" begin
        for var in (vars..., tuples..., Xs, ys)
            @test typeof(@inferred(ObsView(var))) <: ObsView
            @test typeof(@inferred(ObsView(var))) <: ObsView
        end
        for tup in tuples
            @test typeof(@inferred(ObsView(tup))) <: ObsView
        end
        @test typeof(@inferred(ObsView(CustomType()))) <: ObsView
    end

    @testset "AbstractArray interface" begin
        for var in (vars..., tuples..., Xs, ys)
            # for var in (vars[1], )
            A = ObsView(var)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[16]
            @test @inferred(numobs(A)) == 15
            @test @inferred(length(A)) == 15
            @test @inferred(size(A)) == (15,)
            @test @inferred(A[2:3]) == obsview(var, 2:3)
            @test @inferred(A[[1,3]]) == obsview(var, [1,3])
            @test @inferred(A[1]) == obsview(var, 1)
            @test @inferred(A[11]) == obsview(var, 11)
            @test @inferred(A[15]) == obsview(var, 15)
            @test A[end] == A[15]
            @test @inferred(getobs(A,1)) == getobs(var, 1)
            @test @inferred(getobs(A,11)) == getobs(var, 11)
            @test @inferred(getobs(A,15)) == getobs(var, 15)
            @test typeof(@inferred(collect(A))) <: Vector
        end
        for var in vars
            A = ObsView(var)
            @test getobs(A) == var
            @test getobs.(A) == [getobs(var,i) for i in 1:numobs(var)]
        end
        for var in tuples
            A = ObsView(var)
            @test getobs(A) == var
            @test getobs.(A) == [map(x->getobs(x,i), var) for i in 1:numobs(var)]
        end
    end

    @testset "subsetting" begin
        for var_raw in (vars..., tuples..., Xs, ys)
            for var in (var_raw, ObsView(var_raw))
                A = ObsView(var)
                @test getobs(@inferred(ObsView(A))) == @inferred(getobs(A))

                S = @inferred(ObsView(A, 1:5))
                @test typeof(S) <: ObsView
                @test @inferred(length(S)) == 5
                @test @inferred(size(S)) == (5,)
                @test @inferred(A[1:5]) == S[:]
                @test @inferred(getobs(A,1:5)) == getobs(S)
                @test @inferred(getobs(S)) == getobs(ObsView(ObsView(var,1:5)))

                S = @inferred(ObsView(A, 1:5))
                @test typeof(S) <: ObsView
            end
        end
        A = ObsView(X)
        @test typeof(A.data) <: Array
        S = @inferred(ObsView(A))
        @test typeof(S) <: ObsView
        @test @inferred(length(S)) == 15
        @test typeof(S.data) <: Array
    end

    @testset "iteration" begin
        count = 0
        for (i,x) in enumerate(ObsView(X1))
            @test all(i .== x)
            count += 1
        end
        @test count == 15
    end
end


text:
using Dates: now
using NCDatasets

using Oceananigans.Fields
using Oceananigans.Fields: cpudata
using Oceananigans.Diagnostics: Average
using Oceananigans.Utils: validate_intervals, versioninfo_with_gpu, oceananigans_versioninfo
using Oceananigans.Grids: topology, interior_x_indices, interior_y_indices, interior_z_indices,
                          all_x_indices, all_y_indices, all_z_indices

# Possibly should
collect_face_nodes(topo, ξF) = collect(ξF)[1:end-1]
collect_face_nodes(::Bounded, ξF) = collect(ξF)

const default_output_attributes = Dict(
    "u" => Dict("longname" => "Velocity in the x-direction", "units" => "m/s"),
    "v" => Dict("longname" => "Velocity in the y-direction", "units" => "m/s"),
    "w" => Dict("longname" => "Velocity in the z-direction", "units" => "m/s"),
    "b" => Dict("longname" => "Buoyancy",                    "units" => "m/s²"),
    "T" => Dict("longname" => "Conservative temperature",    "units" => "°C"),
    "S" => Dict("longname" => "Absolute salinity",           "units" => "g/kg")
)

"""
    NetCDFOutputWriter{D, O, I, T, S} <: AbstractOutputWriter

An output writer for writing to NetCDF files.
"""
mutable struct NetCDFOutputWriter{D, O, I, T, S} <: AbstractOutputWriter
              filename :: String
               dataset :: D
               outputs :: O
    iteration_interval :: I
         time_interval :: T
                  mode :: String
                slices :: S
              previous :: Float64
               verbose :: Bool
end

"""
    NetCDFOutputWriter(model, outputs; filename, iteration_interval=nothing, time_interval=nothing,
                       global_attributes=Dict(), output_attributes=Dict(), dimensions=Dict(),
                       mode="c", compression=0, with_halos=false, verbose=false, slice_kwargs...)

Construct a `NetCDFOutputWriter` that writes `(label, output)` pairs in `outputs` (which should
be a `Dict`) to a NetCDF file, where `label` is a string that labels the output and `output` is
either a field from the model (e.g. `model.velocities.u`) or a function `f(model)` that returns
something to be written to disk. Custom output requires the spatial `dimensions` to be manually
specified.

Keyword arguments
=================
- `iteration_interval`: Save output every `n` model iterations.

- `time_interval`: Save output every `t` units of model clock time.

- `filename`: Filepath to save output to.

- `global_attributes`: Dict of model properties to save with every file (deafult: `Dict()`)

- `output_attributes`: Dict of attributes to be saved with each field variable (reasonable
  defaults are provided for velocities, buoyancy, temperature, and salinity).

- `dimensions`: A `Dict` of dimension tuples to apply to outputs (useful for function outputs
  as field dimensions can be inferred).

- `with_halos`: Include the halo regions in the grid coordinates and output fields
  (default: `false`).

- `mode`: "a" (for append) and "c" (for clobber or create). Default: "c". See NCDatasets.jl
  documentation for more information on the `mode` option.

- `compression`: Determines the compression level of data (0-9, default 0)

- `slice_kwargs`: `dimname = Union{OrdinalRange, Integer}` will slice the dimension `dimname`.
  All other keywords are ignored. E.g. `xC = 3:10` will only produce output along the dimension
  `xC` between indices 3 and 10 for all fields with `xC` as one of their dimensions. `xC = 1`
  is treated like `xC = 1:1`. Multiple dimensions can be sliced in one call. Not providing slices
  writes output over the entire domain (including halo regions if `with_halos=true`).

Examples
========
Saving the u velocity field and temperature fields, the full 3D fields and surface 2D slices
to separate NetCDF files:

```jldoctest netcdf1
using Oceananigans, Oceananigans.OutputWriters

grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1));

model = IncompressibleModel(grid=grid);

simulation = Simulation(model, Δt=12, stop_time=3600);

fields = Dict("u" => model.velocities.u, "T" => model.tracers.T);

simulation.output_writers[:field_writer] =
    NetCDFOutputWriter(model, fields, filename="fields.nc", time_interval=60)

# output
NetCDFOutputWriter (time_interval=60): fields.nc
├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0)
└── 2 outputs: ["T", "u"]
```

```jldoctest netcdf1
simulation.output_writers[:surface_slice_writer] =
    NetCDFOutputWriter(model, fields, filename="surface_xy_slice.nc",
                       time_interval=60, zC=grid.Nz, zF=grid.Nz+1)

# output
NetCDFOutputWriter (time_interval=60): surface_xy_slice.nc
├── dimensions: zC(1), zF(1), xC(16), yF(16), xF(16), yC(16), time(0)
└── 2 outputs: ["T", "u"]
```

Writing a scalar, profile, and slice to NetCDF:

```jldoctest
using Oceananigans, Oceananigans.OutputWriters

grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 2, 3));

model = IncompressibleModel(grid=grid);

simulation = Simulation(model, Δt=1.25, stop_iteration=3);

f(model) = model.clock.time^2; # scalar output

g(model) = model.clock.time .* exp.(znodes(Cell, grid)); # vector/profile output

h(model) = model.clock.time .* (   sin.(xnodes(Cell, grid, reshape=true)[:, :, 1])
                            .*     cos.(ynodes(Face, grid, reshape=true)[:, :, 1])); # xy slice output

outputs = Dict("scalar" => f, "profile" => g, "slice" => h);

dims = Dict("scalar" => (), "profile" => ("zC",), "slice" => ("xC", "yC"));

output_attributes = Dict(
    "scalar"  => Dict("longname" => "Some scalar", "units" => "bananas"),
    "profile" => Dict("longname" => "Some vertical profile", "units" => "watermelons"),
    "slice"   => Dict("longname" => "Some slice", "units" => "mushrooms")
);

global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7);

simulation.output_writers[:things] =
    NetCDFOutputWriter(model, outputs,
                       iteration_interval=1, filename="things.nc", dimensions=dims, verbose=true,
                       global_attributes=global_attributes, output_attributes=output_attributes)

# output
NetCDFOutputWriter (iteration_interval=1): things.nc
├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0)
└── 3 outputs: ["profile", "slice", "scalar"]
```
"""
function NetCDFOutputWriter(model, outputs; filename,
    iteration_interval = nothing,
         time_interval = nothing,
     global_attributes = Dict(),
     output_attributes = Dict(),
            dimensions = Dict(),
                  mode = "c",
           compression = 0,
            with_halos = false,
               verbose = false,
    xC = with_halos ? all_x_indices(Cell, model.grid) : interior_x_indices(Cell, model.grid),
    xF = with_halos ? all_x_indices(Face, model.grid) : interior_x_indices(Face, model.grid),
    yC = with_halos ? all_y_indices(Cell, model.grid) : interior_y_indices(Cell, model.grid),
    yF = with_halos ? all_y_indices(Face, model.grid) : interior_y_indices(Face, model.grid),
    zC = with_halos ? all_z_indices(Cell, model.grid) : interior_z_indices(Cell, model.grid),
    zF = with_halos ? all_z_indices(Face, model.grid) : interior_z_indices(Face, model.grid)
    )

    validate_intervals(iteration_interval, time_interval)

    # Ensure we can add metadata to the global attributes later by converting to pairs of type {Any, Any}.
    global_attributes = Dict{Any, Any}(k => v for (k, v) in global_attributes)

    # Generates a dictionary with keys "xC", "xF", etc, whose values give the slices to be saved.
    slice_keywords = Dict(name => a for (name, a) in zip(("xC", "yC", "zC", "xF", "yF", "zF"),
                                                         ( xC,   yC,   zC,   xF,   yF,   zF )))

    # Allow values to be of Any type as OffsetArrays may get modified below.
    grid = model.grid
    dims = Dict{String,Any}(
        "xC" => with_halos ? grid.xC : collect(xnodes(Cell, grid)),
        "xF" => with_halos ? grid.xF : collect(xnodes(Face, grid)),
        "yC" => with_halos ? grid.yC : collect(ynodes(Cell, grid)),
        "yF" => with_halos ? grid.yF : collect(ynodes(Face, grid)),
        "zC" => with_halos ? grid.zC : collect(znodes(Cell, grid)),
        "zF" => with_halos ? grid.zF : collect(znodes(Face, grid))
    )

    dim_attribs = Dict(
        "xC" => Dict("longname" => "Locations of the cell centers in the x-direction.", "units" => "m"),
        "xF" => Dict("longname" => "Locations of the cell faces in the x-direction.",   "units" => "m"),
        "yC" => Dict("longname" => "Locations of the cell centers in the y-direction.", "units" => "m"),
        "yF" => Dict("longname" => "Locations of the cell faces in the y-direction.",   "units" => "m"),
        "zC" => Dict("longname" => "Locations of the cell centers in the z-direction.", "units" => "m"),
        "zF" => Dict("longname" => "Locations of the cell faces in the z-direction.",   "units" => "m")
    )

    # Add useful metadata as global attributes
    global_attributes["date"] = "This file was generated on $(now())."
    global_attributes["Julia"] = "This file was generated using " * versioninfo_with_gpu()
    global_attributes["Oceananigans"] = "This file was generated using " * oceananigans_versioninfo()

    # Slice coordinate arrays stored in the dims dict
    for (dim, indices) in slice_keywords
        dim = string(dim) # convert symbol to string
        dims[dim] = dims[dim][get_slice(indices)] # overwrite entries in dims Dict
    end

    # Open the NetCDF dataset file
    dataset = Dataset(filename, mode, attrib=global_attributes)

    # Define variables for each dimension and attributes if this is a new file.
    if mode == "c"
        for (dim_name, dim_array) in dims
            defVar(dataset, dim_name, dim_array, (dim_name,),
                   compression=compression, attrib=dim_attribs[dim_name])
        end

        # Creates an unliimited dimension "time"
        defDim(dataset, "time", Inf)
        defVar(dataset, "time", Float64, ("time",))

        # Ensure we have an attribute for every output. Use reasonable defaults if
        # none were specified by the user.
        for c in keys(outputs)
            if !haskey(output_attributes, c)
                output_attributes[c] = default_output_attributes[c]
            end
        end

        # Initiates empty variables for fields from the user-supplied variable outputs
        for (name, output) in outputs
            if output isa Field
                FT = eltype(output.grid)
                defVar(dataset, name, FT, (netcdf_spatial_dimensions(output)..., "time"),
                       compression=compression, attrib=output_attributes[name])
            else
                defVar(dataset, name, Float64, (dimensions[name]..., "time"),
                       compression=compression, attrib=output_attributes[name])
            end
        end

        sync(dataset)
    end

    # extract outputs whose values are Fields
    field_outputs = filter(o -> o.second isa Field, outputs)

    # Store a slice specification for each field.
    slices = Dict(name => slice_indices(field; xC=xC, yC=yC, zC=zC, xF=xF, yF=yF, zF=zF)
                  for (name, field) in field_outputs)

    return NetCDFOutputWriter(filename, dataset, outputs, iteration_interval, time_interval,
                              mode, slices, 0.0, verbose)
end

Base.open(ow::NetCDFOutputWriter) = Dataset(ow.filename, "a")
Base.close(ow::NetCDFOutputWriter) = close(ow.dataset)

"""
    netcdf_spatial_dimensions(::field)

Returns the NetCDF dimensions associated with a field.

Examples
========
julia> netcdf_spatial_dimensions(model.velocities.u)
("xF", "yC", "zC")

julia> netcdf_spatial_dimensions(model.tracers.T)
("xC", "yC", "zC")
"""
netcdf_spatial_dimensions(::Field{LX, LY, LZ}) where {LX, LY, LZ} = xdim(LX), ydim(LY), zdim(LZ)

xdim(::Type{Face}) = "xF"
ydim(::Type{Face}) = "yF"
zdim(::Type{Face}) = "zF"

xdim(::Type{Cell}) = "xC"
ydim(::Type{Cell}) = "yC"
zdim(::Type{Cell}) = "zC"

# This function allows users to specify slices with integers; eg xC=3.
# Note: size(a[3:3, :, :]) = (1, Ny, Nz) versus size(a[3, :, :]) = (Ny, Nz)
get_slice(n::Integer) = n:n
get_slice(n::UnitRange) = n

"""
    slice_indices(field; slice_specs...)

Returns an array of indices that specify a view over a field's data.
"""
slice_indices(field; slice_specs...) =
    [get_slice(slice_specs[Symbol(dim)]) for dim in netcdf_spatial_dimensions(field)]

function save_output_to_netcdf!(nc, model, f::Field, name, time_index)
    data = cpudata(f) # Transfer data to CPU if parent(output) is a CuArray
    nc.dataset[name][:, :, :, time_index] = view(data, nc.slices[name]...)
end

function save_output_to_netcdf!(nc, model, avg::Average, name, time_index)
    data = avg(model)
    data = dropdims(data, dims=avg.dims)
    colons = Tuple(Colon() for _ in 1:ndims(data))
    nc.dataset[name][colons..., time_index] = data
end

function save_output_to_netcdf!(nc, model, output, name, time_index)
    data = output(model)
    colons = Tuple(Colon() for _ in 1:ndims(data))
    nc.dataset[name][colons..., time_index] = data
end

"""
    write_output!(output_writer, model)

Writes output to netcdf file `output_writer.filename` at specified intervals. Increments the `time` dimension
every time an output is written to the file.
"""
function write_output!(ow::NetCDFOutputWriter, model)
    ds, verbose, filepath = ow.dataset, ow.verbose, ow.filename

    time_index = length(ds["time"]) + 1
    ds["time"][time_index] = model.clock.time

    if verbose
        @info "Writing to NetCDF: $filepath..."
        @info "Computing NetCDF outputs for time index $(time_index): $(keys(ow.outputs))..."

        # Time and file size before computing any outputs.
        t0, sz0 = time_ns(), filesize(filepath)
    end

    for (name, output) in ow.outputs
        # Time before computing this output.
        verbose && (t0′ = time_ns())

        save_output_to_netcdf!(ow, model, output, name, time_index)

        if verbose
            # Time after computing this output.
            t1′ = time_ns()
            @info "Computing $name done: time=$(prettytime((t1′-t0′)/1e9))"
        end
    end

    sync(ow.dataset)

    if verbose
        # Time and file size after computing and writing all outputs.
        t1, sz1 = time_ns(), filesize(filepath)
        verbose && @info begin
            @sprintf("Writing done: time=%s, size=%s, Δsize=%s",
                    prettytime((t1-t0)/1e9), pretty_filesize(sz1), pretty_filesize(sz1-sz0))
        end
    end

    return nothing
end

function Base.show(io::IO, ow::NetCDFOutputWriter)
    freq_int = isnothing(ow.iteration_interval) && isnothing(ow.time_interval) ? "" :
               isnothing(ow.iteration_interval) ? "(time_interval=$(ow.time_interval))" :
               isnothing(ow.time_interval)      ? "(iteration_interval=$(ow.iteration_interval))" :
               "(iteration_interval=$(ow.iteration_interval), time_interval=$(ow.time_interval))"

    dims = join([dim * "(" * string(length(ow.dataset[dim])) * "), "
                 for dim in keys(ow.dataset.dim)])[1:end-2]

    print(io, "NetCDFOutputWriter $freq_int: $(ow.filename)\n",
        "├── dimensions: $dims\n",
        "└── $(length(ow.outputs)) outputs: $(keys(ow.outputs))")
end


text:
using DataStructures

using Soss: Let, Follows, Return, LineNumber
using MacroTools
using MacroTools: @q, striplines
using Parameters

export tocube
function tocube(m::Model)

    m = canonical(m) |> dropLines

    dims = Dict([k => length(x) for (k,x) in pairs(rand(m))])
    dimTotal = sum(values(dims))

    function proc(m, st::Let; ctx=ctx)
        x = st.x
        val = st.rhs
        :($x = @. $val)
    end

    function proc(m, st::Follows; ctx=ctx)
        x = st.x
        val = st.rhs
        u = ctx[:u]
        lo = ctx[:j][1] + 1
        hi = push!(ctx[:j],1,dims[x])

        res = @q begin
            $u[$lo:$hi,:] .= cdf.($val, $x)
        end
        res
    end

    proc(m, st::Return; ctx=ctx)  = :(return $(st.rhs))
    proc(m, st::LineNumber; ctx=ctx) = nothing

    ctx = Dict{Symbol, Any}()
    ctx[:j] = counter(Int)
    ctx[:m] = m
    @gensym u
    ctx[:u] = u

    body = Soss.buildSource(m, proc; ctx=ctx) |> striplines

    f = gensym(:tocube)

    argsExpr = Expr(:tuple,variables(m)...)

    MacroTools.flatten(@q (
        function $f(; kwargs...)
            $u = zeros($dimTotal,10)
            @unpack $argsExpr = kwargs
            $body
            return $u
        end
    ))
end


text:
export POP!, PUSH!, COPYPUSH!, COPYPOP!, GLOBAL_STACK

const GLOBAL_STACK = []

############# global stack operations ##########
@inline function PUSH!(x)
    push!(GLOBAL_STACK, x)
    _zero(x)
end

@inline function POP!(x::T) where T
    NiLangCore.deanc(x, _zero(x))
    loaddata(T, pop!(GLOBAL_STACK))
end

UNSAFE_PUSH!(args...) = PUSH!(args...)
@inline function UNSAFE_POP!(x::T) where T
    loaddata(T, pop!(GLOBAL_STACK))
end

############# local stack operations ##########
@inline function PUSH!(stack, x::T) where T
    push!(stack, x)
    stack, _zero(T)
end

@inline function POP!(stack, x::T) where T
    NiLangCore.deanc(x, _zero(T))
    stack, loaddata(T, pop!(stack))
end

@inline function UNSAFE_POP!(stack, x::T) where T
    stack, loaddata(T, pop!(stack))
end

loaddata(::Type{T}, x::T) where T = x

@dual POP! PUSH!
@dual UNSAFE_POP! UNSAFE_PUSH!

############# copied push/pop stack operations ##########
@inline function COPYPUSH!(stack, x)
    push!(stack, copy(x))
    stack, x
end

@inline function COPYPOP!(stack, x)
    y = pop!(stack)
    NiLangCore.deanc(x, y)
    stack, y
end

UNSAFE_COPYPUSH!(args...) = COPYPUSH!(args...)
@inline function UNSAFE_COPYPOP!(stack, x::XT) where XT
    y = pop!(stack)
    stack, convert(XT, y)
end

@inline function COPYPUSH!(x)
    push!(GLOBAL_STACK, x)
    x
end

@inline function COPYPOP!(x::XT) where XT
    y = pop!(GLOBAL_STACK)
    NiLangCore.deanc(x, y)
    convert(XT, y)
end

@inline function UNSAFE_COPYPOP!(x)
    y = pop!(GLOBAL_STACK)
    y
end

@dual COPYPOP! COPYPUSH!
@dual UNSAFE_COPYPOP! UNSAFE_COPYPUSH!

loaddata(::Type{T}, x::TX) where {T<:NullType, TX} = T(x)


text:
#
# Do work here ...
#

v = sin(1.234567) + cos(0.122122)
println("Julia side: I will return a Float64: ", v)

return v

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
module SfObjects

using Unitful, ..SfGeometry, ..SfMatter

import ..SfMatter: mass, density, bouyancy
import ..SfGeometry: radius, volume, area, shape
import PhysicalConstants.CODATA2018: g_n

export Object
export mass, bouyancy, area, volume, radius, sphere_of, cube_of

struct Object
	shape::Shape
	material::Material
end

shape(object::Object) = object.shape

radius(object::Object) = radius(object.shape)
area(object::Object) = area(object.shape)
volume(object::Object) = volume(object.shape)

density(object::Object) = density(object.material)

mass(object::Object) = mass(object.material, object.shape)

"""
    mass(x, s::Shape)

Calculate the mass of an object with shape `s` made of `x`, which must have [`density`](@ref)
"""
mass(x, s::Shape) = density(x) * volume(s) |> u"kg"

"""
    bouyancy(s::Shape, fluid; g = g_n)

Bouyancy force felt by an object of shape `s` in a medium made of `fluid` (which must have [`density`](@ref)),
undergoing acceleration `g` (which defaults to Earth gravity).
"""
bouyancy(s::Shape, fluid; g = g_n) = bouyancy(volume(s), density(fluid), g = g)

"""
    bouyancy(o, fluid; g = g_n)

Net force felt by an object `o` (which must have [`shape`](@ref) and [`mass`](@ref)) in a medium made of
`fluid` (which must have [`density`](@ref)), undergoing acceleration `g` (which defaults to Earth gravity).

Negative return values imply sinking.
"""
bouyancy(o, fluid; g = g_n) = bouyancy(shape(o), fluid, g = g) - (mass(o) * g)

"""
	sphere_of(material::Material, m::Unitful.Mass)

Create a sphere of `material` that has mass `m`.
"""
sphere_of(material::Material, m::Unitful.Mass) = Object(Sphere(upreferred(m / density(material))), material)

"""
	cube_of(material::Material, m::Unitful.Mass)

Create a cube of `material` that has mass `m`.
"""
cube_of(material::Material, m::Unitful.Mass) = Object(Cube(upreferred(m / density(material))), material)

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics