text:
struct mlParams
  v_Na::Float64
  g_Na::Float64
  v_K::Float64
  g_K::Float64
  v_L::Float64
  g_L::Float64
  I_app::Float64
  gamma_na::Float64
  k_na::Float64
  beta_na::Float64
  gamma_k::Float64
  k_k::Float64
  beta_k::Float64
  M::Float64
  N::Float64
end

function ml(p)
  return mlParams(p["v_Na"] , p["g_Na"] , p["v_K"],p["g_K"] , p["v_L"] , p["g_L"] , p["I_app"] , p["gamma_na"] , p["k_na"] , p["beta_na"] , p["gamma_k"] , p["k_k"] , p["beta_k"] , p["M"] , p["N"])
end


text:
using DiscriminantAnalysis
using Base.Test

MOD = DiscriminantAnalysis

FloatingPointTypes = (Float32, Float64)
IntegerTypes = (Int32, Int64)

function sampledata{U<:Integer}(n_k::Vector{U}, p::Integer)
    p >= 3 || error("need at least 3 dimensions")
    all(n_k .> p) || error("Need p+1 obs per group")
    k = length(n_k)
    Z = vcat([vcat(vcat(eye(Float64,p), Float64[ones(p-1); 0]')/2,
                   rand(Float64, n_k[i]-(p+1), p)) .+ i for i = 1:k]...)
    refs = vcat([Int64[i for j = 1:n_k[i]] for i = 1:k]...)
    σ = sortperm(rand(sum(n_k)))
    (MOD.RefVector(refs[σ]), Z[σ,:])
end

include("test_common.jl")
include("test_lda.jl")
include("test_qda.jl")

ext: jl


text:
A = [1, 2 , 5, 6]
using Interact

@manipulate for a in slider(A, value=3, label="a")
    println(a)
end

using JLD2

@load "./test/experiments/df_point_count_save.jld2" df

df

point_counts = [50, 200, 350, 500, 650, 800, 1100, 1400]
spacings = [5, 7, 10, 14, 18, 22, 26, 30, 36, 40, 50, 60]
diag_sizes = [400, 550, 700]

sections = ["reg alg: sp lap"]
@manipulate for k in slider(1:length(spacings), value=1, label="spacing")
    spacing = spacings[k]
    fig, axs = subplots(length(diag_sizes), sharex=true, figsize=(10,10))
    for (diag_size, ax) in zip(diag_sizes, axs)
        ax.set_title("Diagonal Size: " * string(diag_size))

        local_df = df[((df.diag_size .== diag_size) .& (df.spacing .== spacing)), :]

        # get times
        timers = local_df[:, :timer]
        avg_times = map(x -> get_avg_time(get_timer(x, sections)), timers)
        println(size(avg_times))
        println(size(point_counts))
        println("--------")

        # get flow mae
        median_dicts = local_df[:, :median_results]
        median_flow_mae = map(dict -> dict["flow_mae"], median_dicts)

        # xscale("log")
        # yscale("log")
        subplots_adjust(top=1.5)
        ax.plot(point_counts, avg_times)
        # ax.plot(point_counts, median_flow_mae)
        xlabel("point count")
        ax.set_ylabel("time [s] & mae [pixel]")
    end
    gcf()
end

@manipulate for k in slider(1:length(spacings), value=1, label="spacing")
    spacing = spacings[k]
    fig, axs = subplots(3, sharex=true, figsize=(10,10), gridspec_kw=Dict("height_ratios" => [1, 1, 2]))

    axs[1].set_title("Times at spacing: " * string(spacing))
    axs[2].set_title("MAE at spacing: " * string(spacing))
    for diag_size in diag_sizes

        local_df = df[((df.diag_size .== diag_size) .& (df.spacing .== spacing)), :]

        local_point_counts = local_df[:, :point_count]

        # get times
        timers = local_df[:, :timer]
        avg_times = map(x -> get_avg_time(get_timer(x, sections)), timers)

        subplots_adjust(top=0.9)
        axs[1].plot(local_point_counts, avg_times)
        xlabel("point count")
        axs[1].set_ylabel("time [s]")

        # get flow mae
        median_dicts = local_df[:, :median_results]
        median_flow_mae = map(dict -> dict["flow_mae"], median_dicts)

        axs[2].plot(local_point_counts, median_flow_mae)
        xlabel("point count")
        axs[2].set_ylabel("mae [pixel]")

        # get points found
        points_found = local_df[:, :points_found]
        axs[3].plot(local_point_counts, points_found)
        xlabel("point count")
        axs[3].set_ylabel("points found")
        axs[3].axis("scaled")
    end
    gcf()
end

diag_sizes = [400, 550, 700];

get_avg_time(df[((df.diag_size .== 400) .& (df.spacing .== 5)), :][:, :timer][1]["reg alg: sp lap"])

##
function get_avg_time(timer)
    return (TimerOutputs.time(timer)/TimerOutputs.ncalls(timer))/10e8
end

function get_timer(timer, sections)
    private_sections = copy(sections)
    try
        timer[private_sections[1]]
    catch e
        println(e)
        return TimerOutput()
    end
    out = timer
    while private_sections != []
        out = out[popfirst!(private_sections)]
    end
    return out
end


text:
using
  PyPlot,
  JLD2,
  Printf,
  FourierFlows

using FFTW: ifft

import GeophysicalFlows.MultilayerQG
import GeophysicalFlows.MultilayerQG: energies, fluxes

# Numerical parameters and time-stepping parameters
nx = 128          # 2D resolution = nx^2
ny = nx

stepper = "FilteredRK4"   # timestepper
dt  = 5e-3     # timestep
nsteps = 40000 # total number of time-steps
nsubs  = 250   # number of time-steps for plotting
               # (nsteps must be multiple of nsubs)

# Physical parameters
 Lx  = 2π      # domain size
   μ = 5e-2    # bottom drag
beta = 5       # the y-gradient of planetary PV

nlayers = 2       # these choice of parameters give the
f0, g = 1, 1      # desired PV-streamfunction relations
  H = [0.2, 0.8]  # q1 = Δψ1 + 25*(ψ2-ψ1), and
rho = [4.0, 5.0]  # q2 = Δψ2 + 25/4*(ψ1-ψ2).
  U = zeros(nlayers)
  U[1] = 1.0
  U[2] = 0.0

gr = TwoDGrid(nx, Lx)
x, y = gridpoints(gr)
k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers

# Initialize problem
prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, Lx=Lx, f0=f0, g=g, H=H, rho=rho, U=U, dt=dt, stepper=stepper, μ=μ, beta=beta)
sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid

# Files
filepath = "."
plotpath = "./plots_2layer"
plotname = "snapshots"
filename = joinpath(filepath, "2layer.jld2")

# File management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end

# Initialize with zeros
MultilayerQG.set_q!(prob, 1e-2randn((nx, ny, nlayers)))

# Create Diagnostics
E = Diagnostic(energies, prob; nsteps=nsteps)
diags = [E]

# Create Output
get_sol(prob) = sol # extracts the Fourier-transformed solution
function get_u(prob)
  @. v.qh = sol
  streamfunctionfrompv!(v.psih, v.qh, p.invS, g)
  @. v.uh = -im*g.l *v.psih
  invtransform!(v.u, v.uh, p)
  v.u
end
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))

function plot_output(prob, fig, axs; drawcolorbar=false)

  # Plot the PV field and the evolution of energy and enstrophy.

  sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
  MultilayerQG.updatevars!(prob)

  for j in 1:nlayers
    sca(axs[j])
    pcolormesh(x, y, v.q[:, :, j])
    axis("square")
    xlim(-Lx/2, Lx/2)
    ylim(-Lx/2, Lx/2)
    title(L"$q_"*string(j)*L"$")
    if drawcolorbar==true
      colorbar()
    end

    sca(axs[j+2])
    cla()
    contourf(x, y, v.psi[:, :, j])
    contour(x, y, v.psi[:, :, j], colors="k")
    axis("square")
    xlim(-Lx/2, Lx/2)
    ylim(-Lx/2, Lx/2)
    title(L"$\psi_"*string(j)*L"$")
    if drawcolorbar==true
      colorbar()
    end
  end

  sca(axs[5])
  plot(μ*E.t[E.i], E.data[E.i][1][1], ".", color="b", label=L"$KE_1$")
  plot(μ*E.t[E.i], E.data[E.i][1][2], ".", color="r", label=L"$KE_2$")
  xlabel(L"\mu t")
  ylabel(L"KE")
  # legend()

  sca(axs[6])
  plot(μ*E.t[E.i], E.data[E.i][2][1], ".", color="k", label=L"$PE_{3/2}$")
  xlabel(L"\mu t")
  ylabel(L"PE")
  # legend()

  pause(0.001)
end

fig, axs = subplots(ncols=3, nrows=2, figsize=(15, 8))
plot_output(prob, fig, axs; drawcolorbar=true)

# Step forward
startwalltime = time()

while cl.step < nsteps
  stepforward!(prob, diags, nsubs)

  # Message
  log = @sprintf("step: %04d, t: %d, KE1: %.4f, KE2: %.4f, PE: %.4f, τ: %.2f min", cl.step, cl.t, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60)

  println(log)

  plot_output(prob, fig, axs; drawcolorbar=false)

end
println((time()-startwalltime))

plot_output(prob, fig, axs; drawcolorbar=false)

savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename, dpi=240)


text:
function getfirst(f, A)
    id = findfirst(f, A)
    return id == nothing ? nothing : A[id]
end

macro gputime(codeblock)
    quote
        local const query        = GLuint[1]
        local const elapsed_time = GLuint64[1]
        local const done         = GLint[0]
        glGenQueries(1, query)
        glBeginQuery(GL_TIME_ELAPSED, query[1])
        $(esc(codeblock))
        glEndQuery(GL_TIME_ELAPSED)

        while (done[1] != 1)
            glGetQueryObjectiv(
                query[1],
                GL_QUERY_RESULT_AVAILABLE,
                done
            )
        end
        glGetQueryObjectui64v(query[1], GL_QUERY_RESULT, elapsed_time)
        println("Time Elapsed: ", elapsed_time[1] / 1000000.0, "ms")
    end
end

function print_with_lines(out::IO, text::AbstractString)
    io = IOBuffer()
    for (i,line) in enumerate(split(text, "\n"))
        println(io, Printf.@sprintf("%-4d: %s", i, line))
    end
    write(out, take!(io))
end

print_with_lines(text::AbstractString) = print_with_lines(stdout, text)

"""
Returns the cardinality of a type. falls back to length
"""
cardinality(x) = length(x)
cardinality(x::Number) = 1
cardinality(::Nothing) = 1
cardinality(x::Type{T}) where {T <: Number} = 1

Base.length(::Type{<:Number}) = 1
Base.length(::Type{<:NTuple{N}}) where {N} = N

function glasserteltype(::Type{T}) where T
    # @assert (hasmethod(length, (T,)) || T <: DepthFormat) "Error only types with well defined lengths are allowed.\nSolution define Base.length(::Type{$T})."
end

function istexturesampler(typ::GLenum)
    return (
        typ == GL_SAMPLER_BUFFER || typ == GL_INT_SAMPLER_BUFFER || typ == GL_UNSIGNED_INT_SAMPLER_BUFFER ||
    	typ == GL_IMAGE_2D ||
        typ == GL_SAMPLER_1D || typ == GL_SAMPLER_2D || typ == GL_SAMPLER_3D ||
        typ == GL_UNSIGNED_INT_SAMPLER_1D || typ == GL_UNSIGNED_INT_SAMPLER_2D || typ == GL_UNSIGNED_INT_SAMPLER_3D ||
        typ == GL_INT_SAMPLER_1D || typ == GL_INT_SAMPLER_2D || typ == GL_INT_SAMPLER_3D ||
        typ == GL_SAMPLER_1D_ARRAY || typ == GL_SAMPLER_2D_ARRAY ||
        typ == GL_UNSIGNED_INT_SAMPLER_1D_ARRAY || typ == GL_UNSIGNED_INT_SAMPLER_2D_ARRAY ||
        typ == GL_INT_SAMPLER_1D_ARRAY || typ == GL_INT_SAMPLER_2D_ARRAY
    )
end

"""
    separate(f, A)

Separates the true and false part of `A`.
Single values get passed into `f`.
"""
function separate(f, A)
    trues = [f.(A)...]
    return A[trues], A[(!).(trues)]
end


text:
import LightGraphs:SimpleDiGraph

"""
    CompGraph

Basic graph for computation. While not strictly neccessary to compute anything,
it makes it easier to keep track of things.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> cv = (CompVertex(+, InputVertex(1), InputVertex(2)));

julia> CompGraph(inputs(cv), [CompVertex(x -> 3x, cv)])(2,3) # (2 + 3) * 3
15

julia> CompGraph(inputs(cv), [cv, CompVertex(x -> 3x, cv)])(2,3)
(5, 15)
```

"""
struct CompGraph
    inputs::AbstractVector{<:AbstractVertex}
    outputs::AbstractVector{<:AbstractVertex}
end
CompGraph(input::AbstractVertex, output::AbstractVertex) = CompGraph([input], [output])
CompGraph(input::AbstractVector{<:AbstractVertex}, output::AbstractVertex) = CompGraph(input, [output])
CompGraph(input::AbstractVertex, output::AbstractVector{<:AbstractVertex}) = CompGraph([input], output)

function (g::CompGraph)(x...) where T <:Integer
    @assert length(x) == length(g.inputs) "Must supply one input for each input vertex!"
    memo::Dict{AbstractVertex, Any} = Dict(zip(g.inputs, x))
    if length(g.outputs) == 1
        return output!(memo, g.outputs[1])
    end
    return Tuple(map(v -> output!(memo, v), g.outputs))
end

"""
    output!(memo::Dict{AbstractVertex, Any}, v::AbstractVertex)

Return the output from v given any input in memo by traversing the graph.
Intermediate results from all visited vertices will be stored in memo after
function exits.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> ivs = InputVertex.(1:2);

julia> cv = CompVertex(*, ivs...);

julia> results = Dict{AbstractVertex, Any}(zip(ivs, [2,3]));

julia> output!(results, CompVertex(-, cv, ivs[1]))
4
julia> results
Dict{AbstractVertex,Any} with 4 entries:
  CompVertex(*, [InputVertex(1), InputVertex(2)], [CompVertex(-)]) => 6
  CompVertex(-, [CompVertex(*), InputVertex(1)], [])               => 4
  InputVertex(1, [CompVertex(*), CompVertex(-)])                   => 2
  InputVertex(2, [CompVertex(*)])                                  => 3
```
"""
function output!(memo::Dict{AbstractVertex, Any}, v::AbstractVertex)
    # Calculate outputs which are not already calculated
    return get!(memo, v) do
        inpt = map(iv -> output!(memo, iv), inputs(v))
        out = v(inpt...)
    end
end

"""
    SimpleDiGraph(g::CompGraph)

Return g as a SimpleDiGraph.
"""
LightGraphs.SimpleDiGraph(g::CompGraph) = SimpleDiGraph(mapfoldl(v -> flatten(v), (vs1, vs2) -> unique(vcat(vs1, vs2)), g.outputs))

"""
    SimpleDiGraph(v::AbstractVertex)

Return a SimpleDiGraph of all parents of v
"""
LightGraphs.SimpleDiGraph(v::AbstractVertex)= SimpleDiGraph(flatten(v))

"""
    SimpleDiGraph(vertices::AbstractArray{AbstractVertex,1})

Return a SimpleDiGraph of all given vertices
"""
function LightGraphs.SimpleDiGraph(vertices::AbstractArray{AbstractVertex,1})
    g = SimpleDiGraph(length(vertices))
    for (ind, v) in enumerate(vertices)
        add_edge!.([g], indexin(inputs(v), vertices), ind)
    end
    return g
end

"""
    nv(g::CompGraph)

Return the total number of vertices in the graph.
"""
LightGraphs.nv(g::CompGraph) = nv(SimpleDiGraph(g))

"""
    flatten(v::AbstractVertex)

Return an array of all input parents of v

# Examples
```julia-repl
julia> flatten(CompVertex(+, InputVertex.(1:2)...))
3-element Array{AbstractVertex,1}:
 InputVertex(1)
 InputVertex(2)
 CompVertex(+, [InputVertex(1), InputVertex(2)])
"""
function flatten(v::AbstractVertex, vertices::Vector{AbstractVertex} = Vector{AbstractVertex}(), visited::Vector{AbstractVertex} = Vector{AbstractVertex}())
    v in vertices && return vertices
    if !(v in visited)
        push!(visited, v)
        foreach(iv -> flatten(iv, vertices), inputs(v))
    end
    push!(vertices, v)
    return vertices
end

"""
    vertices(g::CompGraph)

Return an topologically sorted array of all vertices in the graph `g`.

# Examples
```julia-repl
julia> ins = InputVertex.(1:2);

julia> v1 = CompVertex(+, ins...);

julia> v2 = CompVertex(*, v1, ins[2]);

julia> graph = CompGraph(ins, v2);

julia> vertices(graph)
4-element Array{AbstractVertex,1}:
 InputVertex(1)
 InputVertex(2)
 CompVertex(+), inputs=[InputVertex(1), InputVertex(2)]
 CompVertex(*), inputs=[CompVertex(+), InputVertex(2)]
```
"""
LightGraphs.vertices(g::CompGraph) = unique(mapfoldl(flatten, vcat, g.outputs))

"""
    Base.copy(g::CompGraph, cf=clone)

Copies the given `CompGraph g` into a new instance with identical structure.

Argument `cf` may be used to alter the copy on any level, e.g to add new `MutationTraits`.

# Examples
```julia-repl
julia> ivs = InputVertex.(1:2);

julia> cv = CompVertex(+, ivs...);

julia> graph = CompGraph(ivs, [cv])
CompGraph([InputVertex(1), InputVertex(2)], [CompVertex(+)])

julia> gcopy = copy(graph)
CompGraph([InputVertex(1), InputVertex(2)], [CompVertex(+)])

julia> gcopy == graph
false
"""
function Base.copy(g::CompGraph, cf=clone)
    # Can't just have each vertex copy its inputs as a graph with multiple outputs
    # might end up with multiple copies of the same vertex

    # Instead, use the same recursion as when calculating output and store the copies
    # in the memo

    # Will contain mapping between vertex in g and its copy
    # We could initialize it with inputs, but CompGraph does not require that inputs
    # are of type InputVertex
    memo = Dict{AbstractVertex, AbstractVertex}()
    foreach(ov -> copy!(memo, ov, cf), g.outputs)
    return CompGraph(
    map(iv -> memo[iv], g.inputs),
    map(ov -> memo[ov], g.outputs)
    )
end
clone(x;cf=nothing) = deepcopy(x)

"""
    copy!(memo::Dict{AbstractVertex, AbstractVertex}, v::AbstractVertex, cf=clone)

Recursively copy the input parents of v, ensuring that each vertex gets exactly one copy.

Results will be stored in the provided dict as a mapping between original and copy.

Argument `cf` may be used to alter the copy on any level, e.g to add new `MutationTraits`.

# Examples
```julia-repl

julia> result = Dict{AbstractVertex, AbstractVertex}();

julia> NaiveNASlib.copy!(result, CompVertex(+, InputVertex.(1:2)...));

julia> result
Dict{AbstractVertex,AbstractVertex} with 3 entries:
  InputVertex(2)                                  => InputVertex(2)
  InputVertex(1)                                  => InputVertex(1)
  CompVertex(+, [InputVertex(1), InputVertex(2)]) => CompVertex(+, [InputVertex(1), InputVertex(2)])
"""
function copy!(memo::Dict{AbstractVertex, AbstractVertex}, v::AbstractVertex, cf=clone)
    return get!(memo, v) do
        # Recurse until inputs(v) is empty
        ins = map(iv -> copy!(memo, iv, cf), inputs(v))
        cf(v, ins...; cf=cf)
    end
end


text:
include("common.jl")

using .common
using Plots
using LinearAlgebra
using ProgressMeter
gr()
theme(:solarized_light)

dt = 0.00001
m1 = 0.6
m2 = 1 - m1

function step!(p::Particle)
    len_r2 = dot(p.r, p.r)
    len_r3 = len_r2 * sqrt(len_r2)
    a = - (p.r / len_r3);
    p.r += p.v .* dt
    p.v += a .* dt
end

p1 = Particle([1., 0., 0.], [0, 1., 0])

# initialize a 3D plot with 1 empty series
err_plt = plot(1, title= "Error")
plt = plot(1, xlim=(-1.2,1.2), ylim=(-1.2,1.2),
                title = "PJP <3", size = (800, 800),
                m = (:cross, 3, stroke(0)))

n = 60
prog = Progress(n,1)
# build an animated gif by pushing new points to the plot, saving every 10th frame
for i=1:n
    for j=1:100000 /10
        step!(p1)
    end
    push!(err_plt, dot(p1.r, p1.r) - 1)
    push!(plt, p1.r[1], p1.r[2])
    next!(prog)
end

plot(plt, err_plt, layout= grid(2, 1, heights=[0.7,0.3]), size = (500, 700))


text:
module STOCK

import DataFrames
import CSV

# note: excludes "dow" since the corporate entity is < 2 yrs old
# all other stocks have the full 10 years worth of data.
tickers = ["mmm", "axp", "amgn", "aapl", "ba", "cat", "cvx", "csco", "ko", "gs", "hd", "hon", "ibm", "intc", "jnj", "jpm", "mcd", "mrk", "msft", "nke", "pg", "crm", "trv", "unh", "vz", "v", "wba", "wmt", "dis"]
names = ["3M", "American Express", "Amgen", "Apple", "Boeing", "Caterpillar", "Chevron", "Cisco", "Coca-Cola", "Goldman Sachs", "Home Depot", "Honeywell", "IBM", "Intel", "Johnson & Johnson", "JP Morgan Chase", "McDonald's", "Merck", "Microsoft", "Nike", "Procter & Gamble", "Salesforce", "Travelers", "UnitedHealth", "Verizon", "Visa", "Walgreens Boots Alliance", "Walmart", "Walt Disney Company"]
sectors = ["Diversified Holding Companies", "Banking", "Biotechnology", "Computers/Consumer Electronics", "Aerospace Products/Parts", "Mobile Machinery", "Major Oil & Gas", "Networking", "Non-Alcoholic Beverages/Drinks", "Chemicals", "Securities", "Home Goods Retail", "Diversified Holding Companies", "Computer Services", "Semiconductors", "Pharmaceuticals", "Major International Banks", "Restaurants", "Pharmaceuticals", "Software", "Footwear", "Personal Care Products/Appliances", "Software", "Non-Life Insurance", "Life Insurance", "Wireless Telecommunications Services", "Consumer Finance","Drug Retail", "Mixed Retailing", "Broadcasting"]
industries = ["Business/Consumer Services", "Financial Services", "Health Care/Life Sciences", "Technology", "Industrial Goods", "Industrial Goods", "Companies on the Energy Service", "Technology", "Consumer Goods", "Financial Services", "Retail/Wholesale", "Business/Consumer Services", "Business/Consumer Services", "Technology", "Health Care/Life Sciences", "Financial Services", "Leisure/Arts/Hospitality", "Health Care/Life Sciences", "Technology", "Consumer Goods", "Consumer Goods", "Technology", "Financial Services", "Financial Services", "Telecommunication Services", "Financial Services", "Retail/Wholesale", "Retail/Wholesale", "Media/Entertainment"]
colors = Dict(
  "Technology" => "gray",
  "Financial Services" => "aqua",
  "Business/Consumer Services" => "orange",
  "Health Care/Life Sciences" => "red",
  "Industrial Goods" => "fuchsia",
  "Companies on the Energy Service" => "hotpink",
  "Consumer Goods" => "darksalmon",
  "Basic Materials/Resources" => "magenta",
  "Retail/Wholesale" => "darkgoldenrod2",
  "Leisure/Arts/Hospitality" => "yellow",
  "Telecommunication Services" => "lavenderblush",
  "Media/Entertainment" => "seagreen2",
)
frames = []

for ticker in tickers
  global frames
  # println("loading $ticker")
  df = DataFrames.DataFrame(CSV.File("data/$ticker.csv"))
  push!(frames, df)
end

backto = size(frames[1], 1)
Opens = zeros(backto, size(tickers, 1))
Closes = zeros(backto, size(tickers, 1))
Volumes = zeros(backto, size(tickers, 1))

for j in 1:size(tickers, 1)
  # println("getting open price $(tickers[j]) index $j")
  Opens[:, j] = frames[j][1:backto, 2]
  Closes[:, j] = frames[j][1:backto, 5]
  Volumes[:, j] = frames[j][1:backto, 6]
end

end


text:
# This test use audio device.
# Don't add runtests.jl

using FluidSynth

settings = Settings()
setstr(settings, "audio.driver", "pulseaudio")

synth = Synth(settings)

audio_driver = AudioDriver(settings, synth)

sfont_id = sfload(synth, "violin_sample.sf2")
program_select(synth, Int32(0), sfont_id, Int32(0), Int32(0))

noteon(synth, Int32(0), Int32(60), Int32(70))
# wait 3 second to lesten sound
sleep(3)


text:
## Product definitions

"""
Base type for all product definitions
"""
abstract type AbstractProduct end

Base.Broadcast.broadcastable(q:: AbstractProduct) = Ref(q) # treat it as a scalar in broadcasting

"""
Definition of an FX Forward product
"""
struct FXForward{T <: Real, BaseCurrency <: Currency, QuoteCurrency <: Currency} <: AbstractProduct
    notional_base_currency:: CurrencyAmount{T, BaseCurrency}
    quote_currency:: QuoteCurrency
    forward_fx_rate:: ExchangeRate{T, BaseCurrency, QuoteCurrency}
    tenor:: Tenor # relative time to maturity
end

max_issues_repo_name: lungben/MarketRiskPrototype

max_forks_repo_name: lungben/MarketRiskPrototype


text:
push!(LOAD_PATH, @__DIR__)
using Capsule, Test, TestModule

TestModule.@testmodule myTest "valueA" "valueB"

@testset "Instance registration" begin
    @test myTest.ts.fieldA == "valueA"
    @test myTest.ts.fieldB == "valueB"
end
instance = myTest.ts

@testset "function and arguments" begin
    @test myTest.noArg() == :noArg
    @test myTest.customArgumentsAreSupported("hello", "world") == (:a => "hello", :b => "world")
    @test myTest.keywordArgumentsAreSupported(; a = "hello", b = "world") == (:a => "hello", :b => "world")
    @test myTest.keywordArgumentsAreSupported("hello", "world"; c = "mundo", d = "monde") == (:a => "hello", :b => "world", :c => "mundo", :d => "monde")
end

@testset "instance auto-wireing" begin
    @test myTest.instanceArgIsAutowired() == (:ts => instance)
    @test myTest.instanceArgIsAutowired("hello", "world") == (:ts => instance, :a => "hello", :b => "world")
    @test myTest.instanceArgIsAutowired("hello"; b = "world") == (:ts => instance, :a => "hello", :b => "world")
    # @test myTest.instanceArgIsAutowired("hello") == (:ts => instance, :a => "hello")
end

@testset "scoping" begin
    @test :privateMethodIsNotAvailable ∉ names(myTest; all=true)
end


text:
using JuMP
using MathProgBase.SolverInterface
using LazySets
using Parameters
using Interpolations
using Dates
using MathOptInterface
using OVERT
using OVERT: add_overapproximate

"""
Type Definitions.
"""
mutable struct MyHyperrect
	low::Array
	high::Array
end
MyHyperrect(; low=[], high=[]) = MyHyperrect(low,high)
low(h::MyHyperrect) = h.low
high(h::MyHyperrect) = h.high

"""
----------------------------------------------
read controller
----------------------------------------------
"""

function add_controller_constraints(model, network_nnet_address, input_set, input_vars, output_vars; last_layer_activation=Id())
    network = read_nnet(network_nnet_address, last_layer_activation=last_layer_activation)
    neurons = init_neurons(model, network)
    deltas = init_deltas(model, network)
    bounds = get_bounds(network, input_set)
    encode_network!(model, network, neurons, deltas, bounds, BoundedMixedIntegerLP())
    @constraint(model, input_vars .== neurons[1])  # set inputvars
    @constraint(model, output_vars .== neurons[end])  # set outputvars
    return bounds[end]
end

function find_controller_bound(network_file, input_set, last_layer_activation)
    network = read_nnet(network_file; last_layer_activation=last_layer_activation)
    bounds = get_bounds(network, input_set)
    return bounds[end]
end

"""
----------------------------------------------
reasoning and solving concrete queries
----------------------------------------------
"""

function setup_overt_and_controller_constraints(query::OvertQuery, input_set::Hyperrectangle; t_idx::Union{Int64, Nothing}=nothing)
	"""
	This function reads overt query, computes the dictionary range, runs overt,
	  reads the controller nnet file,  and finally setup the overt
	  and controller constraints in form of a JuMP MIP model.
	inputs:
	- query: OvertQuery
	- input_set: Hyperrectangle for the initial set of variables.
	- t_idx: this is for timed dynamics, when the symbolic version is called.
	        if t_idx is an integer, a superscript will be added to all state and
	        control variable, indicating the timestep is a symbolic represenation.
	        default is nothing, which does not do the timed dynamics.
	outputs:
	- mip_model: OvertMIP file that contain the MIP formulation of all overt and
	             controller constraints.
	- oA: all OvertApproximation objects used in the dynamics,
	- oA_vars: all Overt output variables.
	"""

	# read attributes
	dynamics = query.problem.overt_dynamics
	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars

	network_file = query.network_file
	last_layer_activation = query.last_layer_activation
	N_overt = query.N_overt

	# adding subscripts if variables are timed
	if !isnothing(t_idx)
		input_vars = [Meta.parse("$(v)_$t_idx") for v in input_vars]
		control_vars = [Meta.parse("$(v)_$t_idx") for v in control_vars]
   	end

	# setup overt range dictionary
	range_dict = Dict{Symbol, Array{Float64,1}}()
	for i = 1:length(input_vars)
		range_dict[input_vars[i]] = [input_set.center[i] - input_set.radius[i],
		                           input_set.center[i] + input_set.radius[i]]
    end

    # read controller and add bounds to range dictionary
    cntr_bound = find_controller_bound(network_file, input_set, last_layer_activation)
    for i = 1:length(control_vars)
		range_dict[control_vars[i]] = [cntr_bound.center[i] - cntr_bound.radius[i],
                                     cntr_bound.center[i] + cntr_bound.radius[i]]
    end

	# call overt and setup overtMIP
	oA, oA_vars = dynamics(range_dict, N_overt, t_idx)
	mip_model = OvertMIP(oA)

	# add controller to mip
	mip_control_input_vars = [get_mip_var(v, mip_model) for v in input_vars]
	mip_control_output_vars = [get_mip_var(v, mip_model) for v in control_vars]
	controller_bound = add_controller_constraints(mip_model.model, network_file, input_set, mip_control_input_vars, mip_control_output_vars)

	mip_summary(mip_model.model)
	return mip_model, oA, oA_vars
end

function solve_for_reachability(mip_model::OvertMIP, query::OvertQuery,
	oA_vars::Array{Symbol, 1}, t_idx::Union{Int64, Nothing}; get_meas=false)
	"""
	this function sets up states in the future timestep and optimizes that.
	this will give a minimum and maximum on each of the future timesteps.
	the final answer is a hyperrectangle.
	inputs:
	- mip_model: OvertMIP object that contains all overt and controller constraints.
	- query: OvertQuery
	- oA_vars: output variables of overt.
	- t_idx: this is for timed dynamics, when the symbolic version is called.
	        if t_idx is an integer, a superscript will be added to all state and
	        control variable, indicating the timestep is a symbolic represenation.
	        default is nothing, which does not do the timed dynamics.
	outputs:
	- reacheable_set: smallest hyperrectangle that contains the reachable set.
	"""

	# inputs of reachable set hyperrectangle are initialized
	lows = Array{Float64}(undef, 0)
	highs = Array{Float64}(undef, 0)

	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars
	update_rule = query.problem.update_rule
	dt = query.dt

	if !isnothing(t_idx)
		input_vars_last = [Meta.parse("$(v)_$t_idx") for v in input_vars]
		control_vars_last = [Meta.parse("$(v)_$t_idx") for v in control_vars]
		integration_map = update_rule(input_vars_last, control_vars_last, oA_vars)
	else
		input_vars_last = input_vars
		integration_map = query.problem.update_rule(input_vars, control_vars, oA_vars)
	end

	# setup the future state and optimize.
	timestep_nplus1_vars = GenericAffExpr{Float64,VariableRef}[]
	for v in input_vars_last
	   	v_mip = mip_model.vars_dict[v]
		dv = integration_map[v]
		dv_mip = mip_model.vars_dict[dv]
		next_v_mip = v_mip + dt * dv_mip
		push!(timestep_nplus1_vars, next_v_mip)
		@objective(mip_model.model, Min, next_v_mip)
		JuMP.optimize!(mip_model.model)
		@assert termination_status(mip_model.model) == MathOptInterface.OPTIMAL
		push!(lows, objective_bound(mip_model.model))
		@debug "Objective value for MIN is: $(objective_value(mip_model.model)) and objective bound is $(objective_bound(mip_model.model))"
		@objective(mip_model.model, Max, next_v_mip)
		JuMP.optimize!(mip_model.model)
		@assert termination_status(mip_model.model) == MathOptInterface.OPTIMAL
		push!(highs, objective_bound(mip_model.model))
		@debug "Objective value for MAX is: $(objective_value(mip_model.model)) and objective bound is $(objective_bound(mip_model.model))"
	end
	# get the hyperrectangle.
	reacheable_set = Hyperrectangle(low=lows, high=highs)

	# measurement model code,
	# compute reachable set for measurements
	if query.problem.measurement_model != [] && get_meas
		@debug "Computing measurement model reachable sets."
		meas_lows = Array{Float64}(undef, 0)
		meas_highs = Array{Float64}(undef, 0)
		# deal with computing reachable set for measurements
		for measurement_matrix_row in query.problem.measurement_model
			@objective(mip_model.model, Min, sum(measurement_matrix_row.*timestep_nplus1_vars))
			JuMP.optimize!(mip_model.model)
			@assert termination_status(mip_model.model) == MathOptInterface.OPTIMAL
			push!(meas_lows, objective_bound(mip_model.model))
			@objective(mip_model.model, Max, sum(measurement_matrix_row.*timestep_nplus1_vars))
			JuMP.optimize!(mip_model.model)
			@assert termination_status(mip_model.model) == MathOptInterface.OPTIMAL
			push!(meas_highs, objective_bound(mip_model.model))
		end
		meas_reachable_set = Hyperrectangle(low=meas_lows, high=meas_highs)
	else
		meas_reachable_set = nothing
	end

	return reacheable_set, meas_reachable_set
end

function one_timestep_concretization(query::OvertQuery, input_set::Hyperrectangle; t_idx::Union{Int64, Nothing}=nothing, get_meas=false)
   """
	This function computes the next reachable set.
   inputs:
   - query: OvertQuery
	- input_set: Hyperrectangle for the initial set of variables.
   - t_idx: this is for timed dynamics, when the symbolic version is called.
            if t_idx is an integer, a superscript will be added to all state and
            control variable, indicating the timestep is a symbolic represenation.
            default is nothing, which does not do the timed dynamics.
   outputs:
   - output_set: a hyperrectangle for the output reachable set
   - oA: all OvertApproximation objects used in the dynamics,
   - oA_vars: all Overt output variables.
	"""
	# reading some query attributes
	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars
	dt = query.dt

	# setup overt and controller constraints
	mip_model, oA, oA_vars = setup_overt_and_controller_constraints(query, input_set; t_idx=t_idx)

	# solve for reacheability
	output_set, measurement_output_set = solve_for_reachability(mip_model, query, oA_vars, t_idx; get_meas=get_meas)
	return output_set, measurement_output_set, oA, oA_vars
end

function many_timestep_concretization(query::OvertQuery, input_set_0::Hyperrectangle; timed::Bool=false)
   """
	This function computes the next n reachable sets using concretization.
   inputs:
   - query: OvertQuery
	- input_set: Hyperrectangle for the initial set of variables.
   - t_idx: this is for timed dynamics, when the symbolic version is called.
            if t_idx is an integer, a superscript will be added to all state and
            control variable, indicating the timestep is a symbolic represenation.
            default is nothing, which does not do the timed dynamics.
   outputs:
   - all_sets: an array of hyperrectangle of all reachable sets, starting from the init set
   - oA: an array of OvertApproximation objects used in each timestep of concretization,
   - oA_vars: an array of Overt output variables used in each timestep of concretization,
	"""
    input_set = copy(input_set_0)
	all_sets = [input_set_0]
	all_meas_sets = []
    all_oA = Array{OverApproximation}(undef, 0)
    all_oA_vars = []
    for i = 1:query.ntime
        t1 = Dates.now()
        if timed
		   output_set, meas_output_set, oA, oA_vars = one_timestep_concretization(query, input_set; t_idx=i, get_meas=true)
        else
		   output_set, meas_output_set, oA, oA_vars = one_timestep_concretization(query, input_set; get_meas=true)
        end
        t2 = Dates.now()
        println("timestep $i computed in $((t2-t1).value/1000) seconds")
        input_set = output_set
		push!(all_sets, output_set)
		push!(all_meas_sets, meas_output_set)
        push!(all_oA, oA)
        push!(all_oA_vars, oA_vars)
    end

	return all_sets, all_meas_sets, all_oA, all_oA_vars
end

"""
----------------------------------------------
reasoning and solving symbolic queries, basic functions
----------------------------------------------
"""

function setup_mip_with_overt_constraints(query::OvertQuery, input_set::Hyperrectangle)
	"""
	this function computes the query using concrete method for the purpose
		of setting up all overt constraints.
		Ideally, all overt and controller constrains at all timesteps should be computed
		and a mip_model object is constructed for each time step. However, JuMP.model
		does not support merging different models.
		The work around this is to compute all overt constraints alltogether and
		then combine them with controller constraints.
	inputs:
    - query: OvertQuery
 	- input_set: Hyperrectangle for the initial set of variables.
    outputs:
    - mip_model: OvertMIP file that contain the MIP formulation, and overt constraints.
	             Notice that controller constraints are NOT yet added.
	- oA: all OvertApproximation objects used in the dynamics,
	- oA_vars: all Overt output variables.
	"""
	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars
	network_file = query.network_file
	ntime = query.ntime
	update_rule = query.problem.update_rule
	dt = query.dt

	# setup all overt constraints via bounds found by conc retization
	all_sets, all_meas_sets, all_oA, all_oA_vars = many_timestep_concretization(query, input_set; timed=true)
	oA_tot = add_overapproximate(all_oA)
	mip_model = OvertMIP(oA_tot)
	return mip_model, all_sets, all_meas_sets, all_oA_vars
end

function add_controllers_constraints!(mip_model, query::OvertQuery, all_sets)
	"""
	This function add controller constraints to mip_model.
	mip_model is changed in place.
	inputs:
	-mip_model: OvertMIP file that contain the MIP formulation, and overt constraints.
	             but does not yet contain controller constraints.
	- query: OvertQuery
	- all_sets: a list of all reachable sets at all timesteps computed using the concrete approach.
	"""
	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars
	network_file = query.network_file
	ntime = query.ntime

	# add controller to mip
	for i = 1:ntime
		input_set        = all_sets[i]
		input_vars_tmp   = [Meta.parse("$(v)_$i") for v in input_vars]
		control_vars_tmp = [Meta.parse("$(v)_$i") for v in control_vars]
		mip_control_input_vars  = [get_mip_var(v, mip_model) for v in input_vars_tmp]
		mip_control_output_vars = [get_mip_var(v, mip_model) for v in control_vars_tmp]
		controller_bound = add_controller_constraints(mip_model.model, network_file, input_set, mip_control_input_vars, mip_control_output_vars)
	end
end

function match_io!(mip_model, query, all_oA_vars)
	"""
	This function connects output of current timestep to the intputs of future
	timestep, via introducing new constraints in the mip formulation.
	mip_model is changed in place.
	inputs:
	-mip_model: OvertMIP file that contain the MIP formulation, and overt constraints.
	             but does not yet contain controller constraints.
	- query: OvertQuery
	- all_oA_vars:an array of Overt output variables used in each timestep
	"""
	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars
	ntime = query.ntime
	update_rule = query.problem.update_rule
	dt = query.dt

	for i = 1:ntime - 1
		oA_vars_now = all_oA_vars[i]
		input_vars_now = [Meta.parse("$(v)_$i") for v in input_vars]
		input_vars_next= [Meta.parse("$(v)_$(i+1)") for v in input_vars]
		control_vars_now = [Meta.parse("$(v)_$i") for v in control_vars]

		integration_map = update_rule(input_vars_now, control_vars_now, oA_vars_now)
		for j = 1:length(input_vars)
			v = input_vars_now[j]
			dv = integration_map[v]
			next_v = input_vars_next[j]

			v_mip = mip_model.vars_dict[v]
			dv_mip = mip_model.vars_dict[dv]
			next_v_mip = mip_model.vars_dict[next_v]

			@constraint(mip_model.model, next_v_mip == v_mip + dt * dv_mip)
     	end
   	end
end

"""
----------------------------------------------
symbolic queries, reachability.
----------------------------------------------
"""

function symbolic_reachability(query::OvertQuery, input_set::Hyperrectangle; get_meas=false)
   """
	This function computes the reachable set after n timestep symbolically.
   inputs:
   - query: OvertQuery
	- input_set: Hyperrectangle for the initial set of variables.
   outputs:
   - all_sets: an array of hyperrectangle of all reachable sets, starting from the init set
               computed with concretization
   - all_sets_symbolic: a hyperrectangle for the reachable set at t=n, computed symbolically.
	"""
	# setup all overt cosntraints
	mip_model, all_sets, all_meas_sets, all_oA_vars = setup_mip_with_overt_constraints(query, input_set)

	# read neural network and add controller constraints
	add_controllers_constraints!(mip_model, query, all_sets)

	# mip summary of constaints
    mip_summary(mip_model.model)

	# connect outputs of timestep i to inputs of timestep i-1
	match_io!(mip_model, query, all_oA_vars)

	# optimize for the output of timestep ntime.
	set_symbolic, measurement_set_symbolic = solve_for_reachability(mip_model, query, all_oA_vars[end], query.ntime, get_meas=get_meas)
	return all_sets, set_symbolic, all_meas_sets, measurement_set_symbolic
end

"""
----------------------------------------------
symbolic queries, reachability with splitting.
----------------------------------------------
"""

function split_hyperrectangle(input_set::Hyperrectangle, splits_idx)
	"""
 	This function splits the input_set into halves based on indices given in splits_idx
 	- input_set: Hyperrectangle for the initial set of variables.
	- splits_idx: indeces for splitting.
    outputs:
    - input_sets_splitted: an array of hyperrectangles which are splits of
	the original input_set.
 	"""
	if length(splits_idx) == 0
		return input_set
	else
		input_sets_splitted = []
		idx = splits_idx[1]
		center = input_set.center
		radius = input_set.radius

		left_center = copy(center)
		left_center[idx] -= radius[idx] / 2
		left_radius = copy(radius)
		left_radius[idx] /= 2
		left = Hyperrectangle(left_center, left_radius)

		right_center = copy(center)
		right_center[idx] += radius[idx] / 2
		right_radius = copy(radius)
		right_radius[idx] /= 2
		right = Hyperrectangle(right_center, right_radius)

		push!(input_sets_splitted, split_hyperrectangle(left, splits_idx[2:end]))

		push!(input_sets_splitted, split_hyperrectangle(right, splits_idx[2:end]))
	end
	input_sets_splitted = vcat(input_sets_splitted...)
	return input_sets_splitted
end

function symbolic_reachability_with_splitting(query::OvertQuery, input_sets::Array{Any, 1}, splits_idx::Array{Int, 1})
	all_concrete_sets = []
	all_symbolic_sets = []
	all_concrete_meas_sets = []
	all_symbolic_meas_sets = []
	for input_set in input_sets
		concrete_sets, symbolic_set, concrete_meas_sets, symbolic_meas_set = symbolic_reachability_with_splitting(query, input_set, splits_idx)
		push!(all_concrete_sets, concrete_sets)
		push!(all_symbolic_sets, symbolic_set)
		push!(all_concrete_meas_sets, concrete_meas_sets)
		push!(all_symbolic_meas_sets, symbolic_meas_set)
	end
	return all_concrete_sets, all_symbolic_sets, all_concrete_meas_sets, all_symbolic_meas_sets
end

function symbolic_reachability_with_splitting(query::OvertQuery, input_set::Hyperrectangle, splits_idx::Array{Int, 1})
	"""
 	This function splits the input_set into halves based on indices given in splits_idx
		and then computes the reachable set after n timestep symbolically
    inputs:
    - query: OvertQuery
 	- input_set: Hyperrectangle for the initial set of variables.
	- splits_idx: indeces for splitting.
    outputs:
    - all_sets: an array of hyperrectangle of all reachable sets, starting from the init set
                computed with concretization
    - all_sets_symbolic: a hyperrectangle for the reachable set at t=n, computed symbolically.
 	"""
	input_sets_splitted = split_hyperrectangle(input_set, splits_idx)
	all_concrete_sets = []
	all_symbolic_sets = []
	all_concrete_meas_sets = []
	all_symbolic_meas_sets = []
	for this_set in input_sets_splitted
		concrete_sets, symbolic_set, concrete_meas_sets, symbolic_meas_set = symbolic_reachability(query, this_set, get_meas=true)
		push!(all_concrete_sets, concrete_sets)
		push!(all_symbolic_sets, symbolic_set)
		push!(all_concrete_meas_sets, concrete_meas_sets)
		push!(all_symbolic_meas_sets, symbolic_meas_set)
	end
	return all_concrete_sets, all_symbolic_sets, all_concrete_meas_sets, all_symbolic_meas_sets
end

function symbolic_reachability_with_concretization(query::OvertQuery,
	input_set::Hyperrectangle, concretize_every::Union{Int, Array{Int, 1}})
   """
	This function computes the reachable set after n timestep symbolically by
		concretizing after every concretize_every timesteps.

   inputs:
   - query: OvertQuery
   - input_set: Hyperrectangle for the initial set of variables.
   - concretize_every: concretization period.
   outputs:
   - all_sets: an array of hyperrectangle of all reachable sets, starting from the init set
               computed with concretization
   - all_sets_symbolic: a hyperrectangle for the reachable set at t=n, computed symbolically.
	"""
	ntime = query.ntime # number of timesteps
	if isa(concretize_every, Int)
		@assert ntime % concretize_every == 0
		n_loops = Int(query.ntime / concretize_every)
		concretize_every = [concretize_every for i in 1:n_loops]
	elseif isa(concretize_every, Array{Int, 1})
		@assert sum(concretize_every) == query.ntime # make sure the concretization intervals add up to total time
	end

	all_concrete_sets = []
	all_symbolic_sets = []
	all_concrete_meas_sets = []
	all_symbolic_meas_sets = []
	this_set = copy(input_set)
	for n in concretize_every
		query.ntime = n
		concrete_sets, symbolic_set, concrete_meas_sets, symbolic_meas_set = symbolic_reachability(query, this_set, get_meas=true) # pass query and input set
		push!(all_concrete_sets, concrete_sets)
		push!(all_symbolic_sets, symbolic_set)
		push!(all_concrete_meas_sets, concrete_meas_sets)
		push!(all_symbolic_meas_sets, symbolic_meas_set)
		this_set = copy(symbolic_set)
	end

	query.ntime = ntime
	return all_concrete_sets, all_symbolic_sets, all_concrete_meas_sets, all_symbolic_meas_sets
end

function symbolic_reachability_with_concretization_with_splitting(query::OvertQuery,
	input_sets::Array{Any, 1},
	concretize_every::Union{Int, Array{Int, 1}},
	split_idx::Array{Int, 1})

	all_concrete_sets = []
	all_symbolic_sets = []
	all_concrete_meas_sets = []
	all_symbolic_meas_sets = []
	for input_set in input_sets
		concrete_sets, concrete_meas_sets, symbolic_set, symbolic_meas_set = symbolic_reachability_with_concretization_with_splitting(query, input_set, concretize_every_split_idx)
		push!(all_concrete_sets, concrete_sets)
		push!(all_symbolic_sets, symbolic_set)
		push!(all_concrete_meas_sets, concrete_meas_sets)
		push!(all_symbolic_meas_sets, symbolic_meas_set)
	end
	return all_concrete_sets, all_symbolic_sets, all_concrete_meas_sets, all_symbolic_meas_sets
end

function symbolic_reachability_with_concretization_with_splitting(query::OvertQuery,
	input_set::Hyperrectangle,
	concretize_every::Union{Int, Array{Int, 1}},
	split_idx::Array{Int, 1})

	ntime = query.ntime
	if isa(concretize_every, Int)
		@assert ntime % concretize_every == 0
		n_loops = Int(query.ntime / concretize_every)
		concretize_every = [concretize_every for i in 1:n_loops]
	end

	all_concrete_sets = []
	all_symbolic_sets = []
	all_concrete_meas_sets = []
	all_symbolic_meas_sets = []
	this_set = copy(input_set)
	idx = 0
	for n in concretize_every
		idx += 1
		println("concretize step: $idx")
		query.ntime = n
		concrete_sets, concrete_meas_sets, symbolic_set, symbolic_meas_set = symbolic_reachability_with_splitting(query, this_set, split_idx)
		push!(all_concrete_sets, concrete_sets)
		push!(all_symbolic_sets, symbolic_set)
		push!(all_concrete_meas_sets, concrete_meas_sets)
		push!(all_symbolic_meas_sets, symbolic_meas_set)
		this_set = copy(symbolic_set)
	end

	query.ntime = ntime
	return all_concrete_sets, all_symbolic_sets, all_concrete_meas_sets, all_symbolic_meas_sets
end

"""
----------------------------------------------
symbolic queries, satisfiability (feasibility).
----------------------------------------------
"""

function add_feasibility_constraints!(mip_model, query, oA_vars, target_constraints; apply_meas=true)
	"""
	this function adds the constraints associated with checking whether future states
	intersect with the target set or not.
	the mip_model is changed in place.
	inputs:
	- mip_model: OvertMIP object that contains all overt and controller constraints.
	- query: OvertQuery
	- oA_vars: output variables of overt.
	- target_constraints: EITHER a hyperrectangle for the target set OR a list of Constraint objects
	outputs:
	- timestep_nplus1_vars: a list of mip_variables assocaited with the states of
							next timestep this will be used for computing the
							point of feasibility (counter example)

	"""
	input_vars = query.problem.input_vars
	control_vars = query.problem.control_vars
	ntime = query.ntime
	update_rule = query.problem.update_rule
	dt = query.dt

	# setup variables current timestep
	timestep_nplus1_vars = []
	input_vars_last = [Meta.parse("$(v)_$ntime") for v in input_vars]
	control_vars_last = [Meta.parse("$(v)_$ntime") for v in control_vars]
	integration_map = update_rule(input_vars_last, control_vars_last, oA_vars)

	# setup variables at next timestep
	timestep_nplus1_vars = GenericAffExpr{Float64,VariableRef}[]
	for (i, v) in enumerate(input_vars_last)
		v_mip = mip_model.vars_dict[v]
		dv = integration_map[v]
		dv_mip = mip_model.vars_dict[dv]
		next_v_mip = v_mip + dt * dv_mip
		push!(timestep_nplus1_vars, next_v_mip)
	end

	# if a measurement model is present, check variables on measurement model
	if (query.problem.measurement_model != []) && apply_meas
		@debug "Checking feasibility of property on measurement variables, not state variables."
		output_vars = GenericAffExpr{Float64,VariableRef}[]
		for measurement_matrix_row in query.problem.measurement_model
			push!(output_vars, sum(measurement_matrix_row.*timestep_nplus1_vars))
		end
	else
		@debug "Checking feasibility of property state variables."
		output_vars = timestep_nplus1_vars
	end

	# Add constraints for feasibility.
	add_output_constraints!(target_constraints, mip_model.model, output_vars)

	return timestep_nplus1_vars
end

function add_output_constraints!(target_set::Hyperrectangle, model::JuMP.Model, vars::Array)
	# Version for hyperrectangle outputs ::Union{Hyperrectangle,MyHyperrect,InfiniteHyperrectangle}
	for (i, v) in enumerate(vars)
		v_min = low(target_set)[i]
		v_max = high(target_set)[i]
		if isfinite(v_min)
			@constraint(model, v >= v_min)
		end
		if isfinite(v_max)
			@constraint(model, v <= v_max)
		end
	end
end

function add_output_constraints!(constraint::HalfSpace, model::JuMP.Model, vars::Array)
	# ax <= b
	x = reshape(vars, (length(vars), 1))
	c = reshape(constraint.a, (1, length(constraint.a)))
	@constraint(model, (c*x)[1] <= constraint.b)
end

function add_output_constraints!(constraints::HPolyhedron, model::JuMP.Model, vars::Array)
	for constraint in constraints
		add_output_constraints!(constraint, model, vars)
	end
end

function add_output_constraints!(constraint::Constraint, model::JuMP.Model, vars::Array)
	# Constraint contains the coeffs, relation, and scalar, assuming the variables are in
	# state order
	# construct a constraint: a^T x R s
	# where a is the coefficient vector, x is the variable vector, R is the relation and s is the scalar

	# reshape
	x = reshape(vars, (length(vars), 1))
	c = reshape(constraint.coeffs, (1, length(constraint.coeffs)))
	if constraint.relation == :(<=)
		@constraint(model, (c*x)[1] <= constraint.scalar)
	elseif constraint.relation == :(>=)
		@constraint(model, (c*x)[1] >= constraint.scalar)
	elseif constraint.relation == :(==)
		@constraint(model, (c*x)[1] == constraint.scalar)
	end

end

function symbolic_satisfiability_nth(query::OvertQuery, input_set::Hyperrectangle,
	target_set, all_sets, all_oA, all_oA_vars; threads=0, apply_meas=true)
	"""
	This function computes the reachable set after n timestep symbolically.
	inputs:
	- query: OvertQuery
	- input_set: Hyperrectangle for the initial set of states.
	- target_set: Hyperrectangle for the target set of states or Constraint object specify property.
	outputs:
	- status: status of query which can be sat, unsat or error,
	- vals: if sat, returns the counter example at timestep n+1. else returns empty dictionary.
	- stats: if sat, returns the counter example at timestep 1 to n. else returns empty dictionary.
	"""
	vals, stats = Dict(), Dict()
	ntime = query.ntime
	@assert ntime == length(all_oA)

	# combine overt cosntraints
	oA_tot = add_overapproximate(all_oA)
	mip_model = OvertMIP(oA_tot, threads=threads)

	# read neural network and add controller constraints
	add_controllers_constraints!(mip_model, query, all_sets)

	# mip summary of constaints
    mip_summary(mip_model.model)

	# connect outputs of timestep i to inputs of timestep i-1
	match_io!(mip_model, query, all_oA_vars)

	# add feasibility constraints on n+1 timestep OR to measurement that is function of n+1 time states
	timestep_nplus1_vars = add_feasibility_constraints!(mip_model, query, all_oA_vars[end], target_set, apply_meas=apply_meas)

	JuMP.optimize!(mip_model.model)
	if termination_status(mip_model.model) == MathOptInterface.OPTIMAL
		# optimal
		status = "sat"
		vals = value.(timestep_nplus1_vars)
		stats = Dict()
		for i = 1:query.ntime
			input_vars_now = [Meta.parse("$(v)_$i") for v in query.problem.input_vars]
			tmp_dict = Dict((v, value(mip_model.vars_dict[v])) for v in input_vars_now)
			stats = merge(stats, tmp_dict)
		end
	elseif termination_status(mip_model.model) == MathOptInterface.INFEASIBLE
		# infeasible
		status = "unsat"
	else
  		status = "error"
	end
	return status, vals, stats
end

function symbolic_satisfiability(query::OvertQuery, input_set::Hyperrectangle, target_set; unsat_problem::Bool=false, after_n::Int=0, return_all=false, threads=0, apply_meas=true)
	"""
	Checks whether a property P is satisfied at timesteps 1 to n symbolically.
	inputs:
	- query: OvertQuery
	- input_set: Hyperrectangle for the initial set of variables.
	- target_set: Hyperrectangle for the targe set of variables.
	- unsat_problem: if true, it solves an unsatisfiability problem. i.e. the first
	                timestep at which state does not include in target set.
    - after_n: satisfiability is checked for n > after_n. n here is the time step number.
	           default value = 0. meaning that we start from begining.
	outputs:
	- status: status of query which can be sat, unsat or error,
	- vals: if sat, returns the counter example at timestep n+1. else returns empty dictionary.
	- stats: if sat, returns the counter example at timestep 1 to n. else returns empty dictionary.

	"""
	n = query.ntime
	if return_all
		SATii, valii, statii = [], [], []
	else
		SATus, vals, stats = "", Dict(), Dict() # "init" values...
	end
	problem_type = unsat_problem ? "unsat" : "sat"
	println("problem type: ", problem_type)

	input_set_tmp = copy(input_set)
	all_sets = [input_set]
	# all meas sets = [] # don't need to save measurement sets  for a satisfiability problem
	all_oA = Array{OverApproximation}(undef, 0)
	all_oA_vars = []

	for i = 1:n
		println("checking timestep ", i)
		output_set, meas_set, oA, oA_vars = one_timestep_concretization(query, input_set_tmp; t_idx=i, get_meas=false)
		input_set_tmp = output_set
		push!(all_sets, output_set)
		push!(all_oA, oA)
		push!(all_oA_vars, oA_vars)

		if i <= after_n
			println("not yet in the timestep of interest")
			continue
		end

		query.ntime = i
		SATus, vals, stats = symbolic_satisfiability_nth(query, input_set, target_set, all_sets, all_oA, all_oA_vars, threads=threads, apply_meas=apply_meas)
		println("SATus of step ", i, " is ", SATus)
		if return_all
			push!(SATii, SATus)
			push!(valii, vals)
			push!(statii, stats)
		end
		if SATus == problem_type # todo: clean up this problem type stuff
			println("Property violated at timestep $i")
			if ~return_all
				return SATus, vals, stats
			end
	 	elseif SATus == "error"
		 	throw("some error occured at timestep $i")
		end
   	end
	if return_all
		return SATii, valii, statii
	else
		println("Property holds for $n timesteps.")
		return SATus, vals, stats
	end
end

"""
----------------------------------------------
Checking Reachability Queries
----------------------------------------------
"""
function clean_up_sets(all_sets, symbolic_sets, conc_ints; dims=[4,5])
	# concretization intervals helps you parse which of all_sets are symbolic and which are one-step concrete. if the concretization intervals are [10,2] all_sets will contain:
	# init_set, c_t1, c_t2, ..., c_t10, s_t10, c_t11, c_t12
	# where c_t1 = concrete set at time t=1
	# s_t1 = symbolic set at time t=1
	# In constract, symbolic_sets contains:
	# s_t10, s_t12

	# in case all_sets is a set of sets:
	all_sets = vcat(all_sets...)
	sym_idx = cumsum(conc_ints)
	init_set = all_sets[1]
	dup_idx = findall(all_sets .∈ Ref(symbolic_sets)) # find indices of symbolic sets
	deleteat!(all_sets, dup_idx)
	deleteat!(all_sets, 1) # delete init set too. there should now be symidx[end] sets in all_sets, all concrete
	reachable_sets = deepcopy(all_sets)
	for (sym_interval, i) in enumerate(sym_idx)
		reachable_sets[i] = symbolic_sets[sym_interval]
	end
	# debug
	function debug_setcleanup()
		concrete_2d = get_subsets(all_sets, dims)
		reachable_2d = get_subsets(reachable_sets, dims)
		plot(concrete_2d, color="grey", title="debugging set cleanup")
		plot!(reachable_2d, color="red")
	end
	@debug debug_setcleanup()
	return init_set, reachable_sets
end

function clean_up_meas_sets(concrete, symbolic, conc_ints)
	sym_idx = cumsum(conc_ints)
	concrete = vcat(concrete...)
	reachable = deepcopy(concrete)
	for (i,idx) in enumerate(sym_idx)
		reachable[idx] = symbolic[i]
	end
	return reachable
end

function check_avoid_set_intersection(reachable_sets, init_set, avoid_sets)
	# if there is more than one avoid set, they are intersected separately
	# check init set, but should satisfy
	safe = true
	violations = []
	for (as_idx, as) in enumerate(avoid_sets)
		empty_intersects = [isempty(rs ∩ as) for rs in reachable_sets]
		if all(empty_intersects)
			println("Safe for avoid set $as_idx")
		else
			println("Unsafe for avoid set $as_idx. See: $empty_intersects")
			# record index of avoid set and reachable set(s) that intersected
			push!(violations, (as_idx, findall(.!empty_intersects)))
		end
		safe &= all(empty_intersects)
	end
	return safe, violations
end

function check_goal_set_achieved(reachable_sets, init_set, goal_set)
	# check if and when the goal set is check_goal_set_achieved
	# do this by checking if the reachable set is a subset of the goal set at any point in time
	achieved = [rs ⊆ goal_set for rs in reachable_sets]
	idx = findall(achieved .== true)
	if length(idx) == 0
		return false, nothing
	else
		return true, idx[1] # return first index where property is achieved!
	end
end

"""
----------------------------------------------
monte carlo simulation
----------------------------------------------
"""
function monte_carlo_one_simulate(query, x0)
	dynamics_func = query.problem.true_dynamics
	controller_nnet_address  = query.network_file
	last_layer_activation = query.last_layer_activation
	ntime = query.ntime
	dt = query.dt
	controller = read_nnet(controller_nnet_address, last_layer_activation=last_layer_activation)
	xvec = zeros(ntime+1, length(x0))
	xvec[1, :] = x0
	x = x0
	for j = 1:ntime
		u = compute_output(controller, x)
		dx = dynamics_func(x, u)
		x = x + dx*dt
		xvec[j+1, :] = x
	end
	return xvec
end

function monte_carlo_simulate(query::OvertQuery, input_set::Hyperrectangle; n_sim::Int64=1000000)
	"""
	Running monte carlo simulations for comparison.
   inputs:
   - query: OvertQuery
	- input_set: Hyperrectangle for the initial set of variables.
	- n_sim: number of mc simulations
   outputs:
   - output_sets: output sets at all times
	- xvec: simulation results.
	"""

	# unpacking query attributes
	dynamics_func = query.problem.true_dynamics
	controller_nnet_address  = query.network_file
	last_layer_activation = query.last_layer_activation
	ntime = query.ntime
	dt = query.dt
	n_states = length(input_set.center)
	n_meas = length(query.problem.measurement_model)
	min_x = [[Inf64  for n = 1:n_states] for m = 1:ntime]
	max_x = [[-Inf64 for n = 1:n_states] for m = 1:ntime]
	min_y = [[Inf64  for n = 1:n_meas] for m = 1:ntime]
	max_y = [[-Inf64  for n = 1:n_meas] for m = 1:ntime]
	controller = read_nnet(controller_nnet_address, last_layer_activation=last_layer_activation)
	xvec = zeros(n_sim, ntime, n_states)
	yvec = zeros(n_sim, ntime, n_meas)
	x0 = zeros(n_sim, n_states)
	y0 = zeros(n_sim, n_meas)
	for i = 1:n_sim
	  x  = rand(n_states)
	  x .*= input_set.radius * 2
	  x .+= input_set.center - input_set.radius
	  x0[i, :] = x
	  y0[i, :] = [LinearAlgebra.dot(x,y_i) for y_i in query.problem.measurement_model]
	  for j = 1:ntime
	      u = compute_output(controller, x)
	      dx = dynamics_func(x, u)
	      x = x + dx*dt
		  y = [LinearAlgebra.dot(x,y_i) for y_i in query.problem.measurement_model]
	      min_x[j] = min.(x, min_x[j])
	      max_x[j] = max.(x, max_x[j])
		  min_y[j] = min.(y, min_y[j])
		  max_y[j] = max.(y, max_y[j])
	      xvec[i, j, :] = x
		  yvec[i, j, :] = y
	  end
	end

	state_sets = [input_set]
	for (m1, m2) in zip(min_x, max_x)
	  println(m1, m2)
	  push!(state_sets, Hyperrectangle(low=m1, high=m2))
	end
	meas_sets = []
	for (m1, m2) in zip(min_y, max_y)
		push!(meas_sets, Hyperrectangle(low=m1, high=m2))
	end
return state_sets, xvec, x0, meas_sets, yvec, y0
end

"""
----------------------------------------------
plotting
----------------------------------------------
"""

rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

function plot_output_sets(output_sets; idx=[1,2], fig=nothing, linewidth=3,
    linecolor=:black, linestyle=:solid, fillalpha=0, fill=:red)

    fig = isnothing(fig) ? Plots.plot() : fig
    for s in output_sets
        w = s.radius[idx[1]] * 2
        h = s.radius[idx[2]] * 2
        x = s.center[idx[1]] - s.radius[idx[1]]
        y = s.center[idx[2]] - s.radius[idx[2]]
        #plot!(rectangle(w,h,x,y), fillalpha=0.0, kwargs)
        Plots.plot!(rectangle(w,h,x,y), fillalpha=fillalpha, fill=fill, legend=nothing,
                   linewidth=linewidth, linecolor=linecolor, linestyle=linestyle)
        Plots.xlabel!("x_$(idx[1])")
        Plots.ylabel!("x_$(idx[2])")
    end
    return fig
end

function plot_output_hist(data, ntime; fig=nothing, idx=[1,2], nbins=:notspecified)
    fig = isnothing(fig) ? Plots.plot() : fig
    x = data[:, ntime, idx[1]]
    y = data[:, ntime, idx[2]]
    if nbins == :notspecified
        fig = Plots.histogram2d!(x, y, legend=false)
    else
        fig = Plots.histogram2d!(x, y, nbins=nbins, legend=false)
    end
    return fig
end

function plot_mc_trajectories(data0, data; fig=nothing, idx=[1,2], n_traj = 100,
	linecolor=:red, makercolor=:orange, markersize=3)
	fig = isnothing(fig) ? Plots.plot() : fig
	x0 = data0[1:n_traj, idx[1]]
	y0 = data0[1:n_traj, idx[2]]
    x = data[1:n_traj, :, idx[1]]
    y = data[1:n_traj, :, idx[2]]
	x = hcat(x0, x)
	y = hcat(y0, y)
	for i = 1:n_traj
		plot!(x[i, :], y[i, :], marker=:o, linecolor=linecolor,
		markercolor=makercolor, markersize=markersize)
	end
	return fig
end

function plot_output_sets_pgfplot(output_sets; idx=[1,2], fig=nothing, linewidth=:thick,
    linecolor=:black, linestyle=:solid, fillalpha=0, fill=:red, labels=nothing)

	if isnothing(labels)
		labels = ["\$x_$(idx[1])\$", "\$x_$(idx[2])\$"]
	end

	if isnothing(fig)
		fig = PGFPlots.Axis(style="width=10cm, height=10cm", xlabel=labels[1], ylabel=labels[2])
	end

	line_style = "$linestyle, $linecolor, $linewidth, mark=none"
    for s in output_sets
		s1, s2 = s.center[idx[1]], s.center[idx[2]]
		r1, r2 = s.radius[idx[1]], s.radius[idx[2]]
		push!(fig, PGFPlots.Plots.Linear([s1-r1, s1+r1], [s2+r2, s2+r2], style=line_style))
		push!(fig, PGFPlots.Plots.Linear([s1-r1, s1+r1], [s2-r2, s2-r2], style=line_style))
		push!(fig, PGFPlots.Plots.Linear([s1-r1, s1-r1], [s2-r2, s2+r2], style=line_style))
		push!(fig, PGFPlots.Plots.Linear([s1+r1, s1+r1], [s2-r2, s2+r2], style=line_style))
    end

    return fig
end

function plot_output_hist_pgfplot(data, ntime; fig=nothing, idx=[1,2],
	     inner_points=false, labels=nothing)
    if isnothing(labels)
		labels = ["\$x_$(idx[1])\$", "\$x_$(idx[2])\$"]
	end

	if isnothing(fig)
		fig = PGFPlots.Axis(style="width=10cm, height=10cm", xlabel=labels[1], ylabel=labels[2])
	end

    # x = data[:, ntime, idx[1]]
    # y = data[:, ntime, idx[2]]
	# push!(fig, PGFPlots.Plots.Histogram2(x, y, density=true,
	#                                    colormap=PGFPlots.ColorMaps.Named("Jet"))

	points = data[:, ntime, idx]
	if inner_points
		dx = floor(size(points)[1] / 10000)
		dx = max(1, Int(dx))
		pp = PGFPlots.Plots.Scatter(points[1:dx:end, 1], points[1:dx:end, 2], style="mark=o, orange")
		push!(fig, pp)
	end

	border_idx = chull(points).vertices
	p = PGFPlots.Plots.Linear(points[border_idx, 1], points[border_idx, 2], style="solid, orange, line width=3pt, mark=none")
	push!(fig, p)

    return fig
end

function plot_satisfiability_pgfplot(stats, vals, query; fig=nothing, idx=[1,2], labels=nothing)
	if isnothing(labels)
		labels = ["\$x_$(idx[1])\$", "\$x_$(idx[2])\$"]
	end
	if isnothing(fig)
		fig = PGFPlots.Axis(style="width=10cm, height=10cm", xlabel=labels[1], ylabel=labels[2])
	end

	data_sat_x = Array{Float64, 1}(undef, 0)
	data_sat_y = Array{Float64, 1}(undef, 0)
    i = 1
	while true
		vx = "x$(idx[1])_$i"
		vx = Meta.parse(vx)
		if vx ∉ keys(stats)
			break
		end
		push!(data_sat_x, stats[vx])

		vy = "x$(idx[2])_$i"
		vy = Meta.parse(vy)
		push!(data_sat_y, stats[vy])
		i += 1
	end
	push!(data_sat_x, vals[idx[1]])
	push!(data_sat_y, vals[idx[2]])

	n_dim = length(query.problem.input_vars)
	x0 = [stats[Meta.parse("x$(idx[i])_1")] for i = 1:n_dim]
	data_mc = monte_carlo_one_simulate(query, x0)

	pp = PGFPlots.Plots.Linear(data_sat_x, data_sat_y, style="mark=none, orange, dashed, very thick")
	push!(fig, pp)
	pp = PGFPlots.Plots.Linear(data_mc[:, idx[1]], data_mc[:, idx[2]], style="blue")
	push!(fig, pp)
	return fig
end

# make_animation with @animate macro

function get_subsets(sets, dims; T=Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}})
	return [get_subset(s, dims)::T for s in sets]
end

function get_interval_subsets(sets, dims)
	@assert length(dims) == 1
	return [[low(h)[dims]..., high(h)[dims]...] for h::Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}} in sets]
end

# utility to extract subsets of certain dims
# for an array of hyperrectangles (input and output type)
get_subset(h::Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}, dims) = Hyperrectangle(low=low(h)[dims], high=high(h)[dims])

# for an array of HalfSpaces
# assumes that dimensions that are not selected are set to zero
function get_subset(h::HalfSpace{Float64,Array{Float64,1}}, dims)
	# if the indices NOT selected do NOT have coefficients == 0 in h.a, show @info? or @warn
	other_dims = setdiff([1:length(h.a)...], dims)
	other_dim_coeffs = h.a[other_dims]
	if !all(isapprox.(other_dim_coeffs, Ref(0.), atol=eps()))
		@warn "Halfspace representation is inaccurate. A new method must be implemented to adjust the b value of the halfspace to plot accurate projections."
	end
	return HalfSpace(h.a[dims], h.b)
end

function get_lims(sets, dims)
	lows = [low(h)[dims] for h in sets]
	highs = [high(h)[dims] for h in sets]
	lowest = minimum(vcat(lows'...), dims=1)
	highest = maximum(vcat(highs'...), dims=1)
	limits = collect(zip(lowest, highest))
	return limits[1], limits[2]
end

function expand_lims(lim_pair)
	low, high = lim_pair
	return (low - 0.1*abs(low), high + 0.1*abs(high))
end

function plot_reachable_sets(reachable_sets, target_sets, target_set_color, target_set_name, dims, plotname, dirname)
	plotly()
	rs = get_subsets(reachable_sets, dims)
	ts = get_subsets(target_sets, dims)
	xlims, ylims = get_lims(reachable_sets, dims)
	xlims = expand_lims(xlims)
	ylims = expand_lims(ylims)
	p = plot(rs, color="yellow", xlim=xlims, ylim=ylims)
	plot!(ts, color=target_set_color, label=target_set_name)
	Plots.savefig(p, dirname*plotname*".html")
end

function gen_1D_sets(sets, timesteps; width=0, exclamation=true, time_on_y=false)
	# plots line segments at each given timestep
	if time_on_y
		return [sets[i] × Interval(t-width/2,t+width/2) for (i, t) in enumerate(timesteps)]
	else
		return [Interval(t-width/2,t+width/2) × sets[i] for (i, t) in enumerate(timesteps)]
	end
end

function affine_transform(h::Hyperrectangle, scale, shift)
	low_v = low(h).*scale .+ shift
	high_v = high(h).*scale .+ shift
	return Hyperrectangle(low=low_v, high=high_v)
end

function get_rectangle(h::Hyperrectangle, dims)
	@assert length(dims) == 2
	xL, yL = low(h)[dims]
	xH, yH = high(h)[dims]

	return [xL, xH, xH, xL, xL], [yL, yL, yH, yH, yL]
end


text:
@test typeof(@inferred(NoOp())) <: NoOp <: AugmentorBase.AffineOperation

@testset "constructor" begin
    @test str_show(NoOp()) == "AugmentorBase.NoOp()"
    @test str_showconst(NoOp()) == "NoOp()"
    @test str_showcompact(NoOp()) == "No operation"
end
@testset "eager" begin
    @test_throws MethodError AugmentorBase.applyeager(NoOp(), nothing)
    @test AugmentorBase.supports_eager(NoOp) === false
    res1 = rect
    res2 = OffsetArray(rect, -2, -1)
    res3 = OffsetArray(rect, 0, 0)
    imgs = [
        (rect, res1),
        (view(rect, :, :), res1),
        (OffsetArray(rect, -2, -1), res2),
        (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res3),
    ]
    @testset "single image" begin
        for (img_in, img_out) in imgs
            res = @inferred(AugmentorBase.applyeager(NoOp(), img_in))
            @test res == img_out
            @test typeof(res) == typeof(img_out)
        end
    end
end
@testset "affine" begin
    @test_throws MethodError AugmentorBase.toaffinemap(NoOp(), nothing)
    @test AugmentorBase.supports_affine(NoOp) === true
    @test @inferred(AugmentorBase.toaffinemap(NoOp(), rect)) == AffineMap(@SMatrix([1. 0; 0 1]), @SVector([0., 0.]))
    wv = @inferred AugmentorBase.applyaffine(NoOp(), rect)
    @test wv == rect
    @test typeof(wv) <: InvWarpedView{eltype(rect),2}
end
@testset "affineview" begin
    @test AugmentorBase.supports_affineview(NoOp) === true
    wv = @inferred AugmentorBase.applyaffineview(NoOp(), rect)
    @test typeof(wv) <: SubArray{eltype(rect),2}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)) === rect
    @test wv == rect
end
@testset "lazy" begin
    @test AugmentorBase.supports_lazy(NoOp) === true
    @test @inferred(AugmentorBase.applylazy(NoOp(), rect)) === rect
    oa = OffsetArray(rect, -2, -1)
    @test @inferred(AugmentorBase.applylazy(NoOp(), oa)) === oa

    res1 = rect
    res2 = OffsetArray(rect, -2, -1)
    res3 = OffsetArray(rect, 0, 0)
    imgs = [
        (rect, res1),
        (view(rect, :, :), res1),
        (OffsetArray(rect, -2, -1), res2),
        (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res3),
    ]
    @testset "single image" begin
        for (img_in, img_out) in imgs
            res = @inferred(AugmentorBase.applylazy(NoOp(), img_in))
            @test res == img_out
        end
    end
end
@testset "view" begin
    @test AugmentorBase.supports_view(NoOp) === true
    @test @inferred(AugmentorBase.applyview(NoOp(), rect)) === view(rect, IdentityRange(1:2), IdentityRange(1:3))
end
@testset "stepview" begin
    @test AugmentorBase.supports_stepview(NoOp) === true
    @test @inferred(AugmentorBase.applystepview(NoOp(), rect)) === view(rect, 1:1:2, 1:1:3)
end
@testset "permute" begin
    @test AugmentorBase.supports_permute(NoOp) === false
end


text:
# A CKAR(2) model
# Here we demonstrate the package's utilities in the context of a simple censored and kinked AR(2) model with two regimes.
#
# Consider the model
# $$\begin{align}
# y_t^*&=\phi_1^*y_{t-1}^*+\phi_1y_{t-1}+\phi_2^*y_{t-1}^*+\phi_2y_{t-1}+\epsilon_t,\\
# y_t&=\max\{y_t^*,0\}.
# \end{align}$$
#
# This features a 'positive regime' ($y^*\geq0$) and a 'negative regime' ($y^*<0$).
##
# Noting $y_t=\mathbf{1}\{y_t^*\geq0\}y_t^*$, we can rewrite the equation for $y_t^*$ as $$y_t^*=(\phi_1^*+\phi_1\mathbf{1}\{y_{t-1}\geq0\})y_{t-1}^*+(\phi_2^*+\phi_2\mathbf{1}\{y_{t-2}\geq0\})y_{t-2}^*+\epsilon_t,$$
# where $\mathbf{1}\{y\geq0\}$ is an indicator variable.
##
using ThresholdStability

function AR2_to_TAR(ϕs, ϕ_stars)
    ϕ1, ϕ2 = ϕs
    ϕ1_star, ϕ2_star = ϕ_stars
    Σ = []
    vals = [1, -1]
    for val1 in vals
        for val2 in vals
            A = zeros(2, 2); A[2, 1] = 1
            A[1, 1] = ϕ1_star + ϕ1 * indicator(val1, 0)
            A[1, 2] = ϕ2_star + ϕ2 * indicator(val2, 0)
            push!(Σ, A)
        end
    end
    return Vector{Array{Float64, 2}}(Σ)
    # NOTE Σ is of form where Σ[1] is with yₜ₋₁,yₜ₋₂ ≥ 0, Σ[2] is with yₜ₋₁ ≥ 0 and yₜ₋₂ < 0, Σ[3] is with yₜ₋₁ < 0 and yₜ₋₂ ≥ 0
    # and Σ[4] is with yₜ₋₁,yₜ₋₂ < 0.
end
##
# This yields a set of four $2\times2$ matrices, with each matrix corresponding to a different regime pair (positive-positive, positive-negative, negative-positive and negative-negative.)
#
# For exposition, consider $\phi_1 = 0.4$, $\phi_1^* = 0.2$, $\phi_2=0.2$, and $\phi_2^*=0.1$:
##
Σ4 = AR2_to_TAR([0.4, 0.2], [0.2, 0.1])
##
# The state space constraints for this model are given by `X4`:
##
E1, E2, E3, E4 = [1 0.; 0 1.], [1 0.; 0 -1.], [-1 0.; 0 1.], [-1 0.; 0 -1.]
D1 = zeros(1,2); D2, D3, D4 = copy(D1), copy(D1), copy(D1)
X4 = [[E1, D1], [E2, D2], [E3, D3], [E4, D4]]
##
# and the automaton can be constructed using
##
G = automaton_constructor(Σ4)
##
# We then have the discrete system
s4 = discreteswitchedsystem(Σ4, G, X4)
##
# Alternatively, by substituting for $y_t$ on the first lag only we can produce a set of two $3\times3$ matrices:
##
function AR2_to_companion(ϕs, ϕ_stars)
    ϕ1, ϕ2 = ϕs
    ϕ1_star, ϕ2_star = ϕ_stars
    Σ = []
    vals = [1, -1]
    for val in vals
        A = zeros(3, 3); A[2, 1] = 1; A[1, 2] = ϕ2_star; A[1, 3] = ϕ2
        A[1, 1] = ϕ1_star + ϕ1 * indicator(val, 0)
        A[3, 1] = indicator(val, 0)
        push!(Σ, A)
    end
    return Vector{Array{Float64, 2}}(Σ)

    # NOTE Σ has form s.t. Σ[1] is with yₜ₋₁ ≥ 0 and Σ[2] is with yₜ₋₁ < 0
end

Σ2 = AR2_to_companion([0.4, 0.2], [0.2, 0.1])
Σ2st = [Σ2[1], Σ2[2], Σ2[1], Σ2[2]]
##
# `Σ2` consists of the two matrices and `Σ2st` corresponds to the same partitioning of the state space as for `Σ4`.

In this latter case, the automaton is again `G` but the state space constraints are now given by `X2`:
##
E1, E2, E3, E4 = [1 0 0.; 0 1 0.], [1 0 0.; 0 -1 0.], [-1 0 0.; 0 1 0.], [-1 0 0.; 0 -1 0.]
D1, D3 = [0 1 -1.], [0 1 -1.]; D2, D4 = [0 0 1.], [0 0 1.]
X2 = [[E1, D1], [E2, D2], [E3, D3], [E4, D4]]
##
# and we construct the discrete switched systems
##
s2 = discreteswitchedsystem(Σ2)
s2st = discreteswitchedsystem(Σ2st, G, X2)
##
# For our choice of parameter values, this model is stable. In particular, the upper bounds on the joint spectral radius (JSR) of `Σ2`, the constrained joint spectral radius (CJSR) of `(Σ4, G)`, and the state-constrained joint spectral radii (SCJSR) of `(Σ4, G, X4)` and `(Σ2st, G, X2)` all agree in value, up to 4 significant figures:
##
@show γ_jsr2 = jsr(s2)
@show γ_cjsr4 = cjsr(s4)
@show γ_scjsr4 = sosbound_γ(s4, 2)
@show γ_scjsr2 = sosbound_γ(s2st, 2)
##
# Indeed, if $\phi_i,\phi_i^*\geq0$ for $i=1,2$, then the model will be stable if and only if
# $\sum_{i=1}^2(\phi_i+\phi_i^*)<1$.
#
# As expected, the model appears stable when plotted:
##
using Plots, Distributions
using LaTeXStrings
pyplot()
function simulate_AR2(y0, Σ, T, σ)  # for companion form
    y = zeros(3, T)
    y[:, 1] = y0
    for t in 1:T-1
        ϵ_t = [rand(Normal(0, σ)), 0., 0.]
        if y[1, t] ≥ 0.
            y[:, t+1] = Σ[1]*y[:, t] + ϵ_t
        else
            y[:, t+1] = Σ[2]*y[:, t] + ϵ_t
        end
    end
    return y
end

function plot_AR2(y0, Σ, T, σ; N = 20, row=1)
    ys = []
    ens_means = zeros(T)
    for i in 1:N
        y = simulate_AR2(y0, Σ, T, σ)
        y = y[row,:]
        push!(ys, y)
        ens_means .+= y
    end
    ens_means ./= N  # ensemble means

    E_y = simulate_AR2(y0, Σ, T, 0.)  # calculating deterministic results
    E_y = E_y[row,:]

    plot(ys, color = :grey, alpha = 0.1, label = "")
    plot!(ens_means, color = :grey, linewidth = 2, label = "Ensemble mean")
    plot!(E_y, color = :blue, linewidth = 2, linestyle = :dash, label = "Deterministic")
    plot!(xlabel=L"t", legend=:topright)
end

plot_AR2(3*ones(3), Σ2, 200, 1)
plot!(ylabel=L"y^*", yguidefontrotation=-90)
##
# Introducing negative parameters, we see that the different upper bounds can diverge.
#
# Consider $\phi_1=0.5$, $\phi_1^*=0.5$, $\phi_2=-0.47$, $\phi_2^*=-0.5$. When plotted, the system appears stable:
##
Σ = AR2_to_companion([0.5, -0.47], [0.5, -0.5])
plot_AR2(3*ones(3), Σ, 200, 1, row=1)
plot!(ylabel=L"y^*", yguidefontrotation=-90)
##
# However, this is a case where the conservatism of (the upper bounds on) the JSR and CJSR would prevent us from concluding the system is stable.
##
Σ = AR2_to_TAR([0.5, -0.47], [0.5, -0.5])
E1, E2, E3, E4 = [1 0.; 0 1.], [1 0.; 0 -1.], [-1 0.; 0 1.], [-1 0.; 0 -1.]
D1 = zeros(1,2); D2, D3, D4 = copy(D1), copy(D1), copy(D1)
X = [[E1, D1], [E2, D2], [E3, D3], [E4, D4]]  # state space constraints
G = automaton_constructor(Σ)
s = discreteswitchedsystem(Σ, G, X)
@show γ_jsr = jsr(s)  # > 1
@show γ_cjsr = cjsr(s)  # > 1
##
# However, the upper bound on the SCJSR is less conservative, yielding a bound below 1. We can thus conclude the system is stable.
##
@show γ_scjsr = sosbound_γ(s, 2)  # < 1
##
# Note that (the upper bound on) the SCJSR remains conservative, even if less so than the CJSR or JSR.
#
# For example, consider $\phi_1=\phi_1^*=0.6$ and $\phi_2=\phi_2^*=-0.6$. This appears stable:
##
Σ = AR2_to_companion([0.6, -0.6], [0.6, -0.6])
plot_AR2(3*ones(3), Σ, 200, 1, row=1)
plot!(ylabel=L"y^*", yguidefontrotation=-90)
##
# Plots for other initial values $y_0$ follow a similar pattern.
#
# However, the bound on the SCJSR exceeds 1 in this case (as does the bound for the CJSR)
##
Σ = AR2_to_TAR([0.6, -0.6], [0.6, -0.6])
s = discreteswitchedsystem(Σ, G, X)
@show γ_cjsr = cjsr(s)
@show γ_scjsr = sosbound_γ(s, 2)


text:
#------------------------------------------------------------------------------
function excise(x)

  vv = vec(any(isnan(x),2))

  if any(vv)              #only keep rows with no NaNs
    x = x[!vv,:]
  end

  return x

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function exciseCopy(x)

  vv = vec(any(isnan(x),2))

  if any(vv)              #only keep rows with no NaNs
    x1 = x[!vv,:]
  else
    x1 = deepcopy(x)
  end

  return x1

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise1mPs(x)

  vv = vec(any(isnan(x),2))

  if any(vv)              #only keep rows with no NaNs
    x = x[!vv,:]
  end

  return x,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise1mCopyPs(x)

  vv = vec(any(isnan(x),2))

  if any(vv)              #only keep rows with no NaNs
    x1 = x[!vv,:]
  else
    x1 = deepcopy(x)
  end

  return x1,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise2mPs(x,y)

  vv = vec(any(isnan([x y]),2))

  if any(vv)              #only keep rows with no NaNs
    x = x[!vv,:]
    y = y[!vv,:]
  end

  return x,y,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise2mCopyPs(x,y)

  vv = vec(any(isnan([x y]),2))

  if any(vv)              #only keep rows with no NaNs
    x1 = x[!vv,:]
    y1 = y[!vv,:]
  else
    x1 = deepcopy(x)
    y1 = deepcopy(y)
  end

  return x1,y1,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise3mPs(x,y,z)

  vv = vec(any(isnan([x y z]),2))

  if any(vv)              #only keep rows with no NaNs
    x = x[!vv,:]
    y = y[!vv,:]
    z = z[!vv,:]
  end

  return x,y,z,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise3mCopyPs(x,y,z)

  vv = vec(any(isnan([x y z]),2))

  if any(vv)              #only keep rows with no NaNs
    x1 = x[!vv,:]
    y1 = y[!vv,:]
    z1 = z[!vv,:]
  else
    x1 = deepcopy(x)
    y1 = deepcopy(y)
    z1 = deepcopy(z)
  end

  return x1,y1,z1,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise4mPs(x,y,z,w)

  vv = vec(any(isnan([x y z w]),2))

  if any(vv)              #only keep rows with no NaNs
    x = x[!vv,:]
    y = y[!vv,:]
    z = z[!vv,:]
    w = w[!vv,:]
  end

  return x,y,z,w,vv

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function excise4mCopyPs(x,y,z,w)

  vv = vec(any(isnan([x y z w]),2))

  if any(vv)              #only keep rows with no NaNs
    x1 = x[!vv,:]
    y1 = y[!vv,:]
    z1 = z[!vv,:]
    w1 = w[!vv,:]
  else
    x1 = deepcopy(x)
    y1 = deepcopy(y)
    z1 = deepcopy(z)
    w1 = deepcopy(w)
  end

  return x1,y1,z1,w1,vv

end
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
module MutateIt

@windows? (const thelib = "mutateit") : (const thelib = "libmutateit")  # Name of the shared lib.

export mutateit

function mutateit(t_type, member::ASCIIString, val)
	# Mutate the member 'member' of an immutable type whose pointer is T_TYPE
	# VAL is the new value of the MEMBER field.
	# It's up to the user to guarantie that MEMBER and VAL have the same data type
	# T_TYPE can actually be either a variable of a certain type or a pointer to it.
	# In latter case, we fish the specific DataType from it.
	if (isa(t_type, Ptr))
		x_type = unsafe_load(t_type)
		p_type = t_type
	else
		x_type = t_type
		p_type = pointer([t_type])	# We need the pointer to type to later send to blind_change
	end
	dt  = typeof(x_type)			# Get the specific DataType. That's what we'll need for next inquires
	fo  = fieldoffsets(dt);		ft = dt.types
	ind = findfirst(fieldnames(dt), symbol(member))	# Find the index of the "is_continuous" member
	# This would work too
	# ind = ccall(:jl_field_index, Cint, (Any, Any, Cint), dt, symbol(member), 1) + 1
	p_val = pointer([val])
	blind_change_struct(p_type, p_val, @sprintf("%s", ft[ind]), fo[ind])
	typeof(p_type); 	typeof(p_val)		# Just to be sure that GC doesn't kill them before their due time
end

function blind_change_struct(X, what, keyword::ASCIIString, off::Integer)
	ccall((:blind_change_struct, thelib), Cint, (Ptr{Void}, Ptr{Void}, Ptr{UInt8}, Csize_t),
	       X, what, keyword, off)
end

end # module


text:
module CompleteAreaTrigger

using ..Ahorn, Maple

placements = Dict{String, Ahorn.EntityPlacement}(
    "Complete Area" => Ahorn.EntityPlacement(
        (x, y) -> Maple.Trigger("everest/getMeOutTrigger", x=x, y=y, width=16, height=16),
        "rectangle"
    )
)

end

max_stars_repo_name: thefoxbot/MoreEntities

max_issues_repo_name: thefoxbot/MoreEntities

max_forks_repo_name: thefoxbot/MoreEntities


text:
@testset "Custom Fourier Types" begin
    N = 5
    x = randn((N, N))
    fs = FourierTools.FourierSplit(x, 2, 2, 4)
    @test FourierTools.parenttype(fs) == typeof(x)

    fj = FourierTools.FourierJoin(x, 2, 2, 4)

    @test FourierTools.parenttype(fj) == typeof(x)
    @test FourierTools.parenttype(typeof(fj)) == typeof(x)

    @test FourierTools.IndexStyle(typeof(fj)) == IndexStyle(typeof(fj))
end


text:
module MPIExtensions

using Adapt: Adapt
using CUDA: CUDA
using MPI: MPI

const mpi_is_cuda_aware = Ref(false)

# Data Movement between CPU and GPU. `Lux` & `Flux` have much better support for this but this is all we need here
# If not an abstract array don't do anything
# Maybe we want to have a central repository containing these device transfer utilities
cpu(x) = x
gpu(x) = x
cpu(x::AbstractArray) = Adapt.adapt(Array, x)
gpu(x::AbstractArray) = Adapt.adapt(CUDA.CuArray, x)

function __init__()
    # If using `mpi_is_cuda_aware` anywhere use the @static macro since we anyways need to recompile the code when MPI implementation changes
    disable_cuda_aware_support = parse(Bool, get(ENV, "FLUXMPI_DISABLE_CUDAMPI_SUPPORT", "false"))
    mpi_is_cuda_aware[] = !disable_cuda_aware_support && MPI.has_cuda()
    if disable_cuda_aware_support && MPI.has_cuda()
        @info "CUDA-aware MPI support disabled using `FLUXMPI_DISABLE_CUDAMPI_SUPPORT=true`" maxlog = 1
    elseif !mpi_is_cuda_aware[]
        @warn "MPI Implementation is not CUDA Aware" maxlog = 1
    end
end

# Certain Non-blocking collective communication primitives
# TODO: Should be moved to MPI.jl?
"""
    Iallreduce!(sendbuf, recvbuf, op, comm)
    Iallreduce!(sendrecvbuf, op, comm)

Performs non-blocking elementwise reduction using the operator `op` on the buffer `sendbuf`. `sendbuf` can also be a scalar, in which case `recvbuf` will be a value of the same type.

`recvbuf` and an MPI_Request object are returned. The value in `recvbuf` is only valid after the request has been completed. (`MPI.Wait!`)

!!! warning
    OpenMPI doesn't support Iallreduce! with CUDA. See https://github.com/open-mpi/ompi/issues/9845
"""
function Iallreduce!(rbuf::MPI.RBuffer, op::Union{MPI.Op,MPI.MPI_Op}, comm::MPI.Comm)
    req = MPI.Request()
    # int MPI_Iallreduce(const void *sendbuf, void *recvbuf, int count,
    #                    MPI_Datatype datatype, MPI_Op op, MPI_Comm comm,
    #                    MPI_Request * request)
    MPI.@mpichk ccall(
        (:MPI_Iallreduce, MPI.libmpi),
        Cint,
        (MPI.MPIPtr, MPI.MPIPtr, Cint, MPI.MPI_Datatype, MPI.MPI_Op, MPI.MPI_Comm, Ptr{MPI.MPI_Request}),
        rbuf.senddata,
        rbuf.recvdata,
        rbuf.count,
        rbuf.datatype,
        op,
        comm,
        req,
    )
    req.buffer = rbuf
    finalizer(MPI.free, req)
    return rbuf.recvdata, req
end

Iallreduce!(rbuf::MPI.RBuffer, op, comm::MPI.Comm) = Iallreduce!(rbuf, MPI.Op(op, eltype(rbuf)), comm)

Iallreduce!(sendbuf, recvbuf, op, comm::MPI.Comm) = Iallreduce!(MPI.RBuffer(sendbuf, recvbuf), op, comm)

Iallreduce!(buf, op, comm::MPI.Comm) = Iallreduce!(MPI.IN_PLACE, buf, op, comm)

"""
    Ibcast!(buf, root, comm)

Non-blocking broadcast of the buffer `buf` to all processes in `comm`.

`recvbuf` and an MPI_Request object are returned. The value in `recvbuf` is only valid after the request has been completed. (`MPI.Wait!`)
"""
function Ibcast!(buf::MPI.Buffer, root::Integer, comm::MPI.Comm)
    req = MPI.Request()
    # int MPI_Ibcast(void *buffer, int count, MPI_Datatype datatype, int root,
    #                MPI_Comm comm, MPI_Request *request)
    MPI.@mpichk ccall(
        (:MPI_Ibcast, MPI.libmpi),
        Cint,
        (MPI.MPIPtr, Cint, MPI.MPI_Datatype, Cint, MPI.MPI_Comm, Ptr{MPI.MPI_Request}),
        buf.data,
        buf.count,
        buf.datatype,
        root,
        comm,
        req,
    )
    req.buffer = buf
    finalizer(MPI.free, req)
    return buf.data, req
end

Ibcast!(buf, root::Integer, comm::MPI.Comm) = Ibcast!(MPI.Buffer(buf), root, comm)

"""
    allreduce!(v, op, comm)

Perform `MPI.Allreduce!` ensuring that CuArrays are safely transfered to CPU if CUDA-aware MPI is unavailable/disabled.
"""
function allreduce!(v::CUDA.CuArray, op::Any, comm::MPI.Comm)
    @static if !mpi_is_cuda_aware[]
        v = cpu(v)
    end
    MPI.Allreduce!(v, op, comm)
    @static if !mpi_is_cuda_aware[]
        v = gpu(v)
    end
    return v
end

function allreduce!(v::AbstractArray, op::Any, comm::MPI.Comm)
    MPI.Allreduce!(v, op, comm)
    return v
end

"""
    bcast!(v, op, comm)

Perform `MPI.Bcast!` ensuring that CuArrays are safely transfered to CPU if CUDA-aware MPI is unavailable/disabled.
"""
function bcast!(v::CUDA.CuArray, root::Integer, comm::MPI.Comm)
    @static if !mpi_is_cuda_aware[]
        v = cpu(v)
    end
    MPI.Bcast!(v, root, comm)
    @static if !mpi_is_cuda_aware[]
        v = gpu(v)
    end
    return v
end

function bcast!(v::AbstractArray, root::Integer, comm::MPI.Comm)
    MPI.Bcast!(v, root, comm)
    return v
end

"""
    reduce!(v, op, comm)

Perform `MPI.Reduce!` ensuring that CuArrays are safely transfered to CPU if CUDA-aware MPI is unavailable/disabled.
"""
function reduce!(v::CUDA.CuArray, op::Any, root::Integer, comm::MPI.Comm)
    @static if !mpi_is_cuda_aware[]
        v = cpu(v)
    end
    MPI.Reduce!(v, op, root, comm)
    @static if !mpi_is_cuda_aware[]
        v = gpu(v)
    end
    return v
end

function reduce!(v::AbstractArray, op::Any, root::Integer, comm::MPI.Comm)
    MPI.Reduce!(v, op, root, comm)
    return v
end

end


text:
include("help.jl")
function mark_field!(r::Robot)
    num_vert = moves!(r, Sud)
    num_hor = moves!(r, West)
    #УТВ: Робот - в Юго-Западном углу

    putmarkers!(r, Nord)
    #УТВ: Все клетки поля промаркированы

    moves!(r, Sud)
    moves!(r, West)
    moves!(r, Nord, num_vert)
    moves!(r, Ost, num_hor)
    #УТВ: Робот - в исходном положении
end

function move_if_possible!(r::Robot, direct_side::HorizonSide)::Bool
    orthogonal_side = next(direct_side)
    reverse_side = inverse(orthogonal_side)
    num_steps=0
    while isborder(r, direct_side) == true
        if !isborder(r, orthogonal_side)
            move!(r, orthogonal_side)
            num_steps += 1
        else
            break
        end
    end
    #УТВ: Робот или уперся в угол внешней рамки поля, или готов сделать шаг (или несколько) в направлении
    # direct_side
    if isborder(r,direct_side) == false
        move!(r,direct_side)
        if num_steps > 0
            while isborder(r,reverse_side) == true
                move!(r,direct_side)
            end
        end
        result = true
    else
        result = false
    end

    for _ in 1:num_steps
        move!(r, reverse_side)
    end

    return result
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while move_if_possible!(r,side)
        num_steps+=1
    end
    return num_steps
end

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move_if_possible!(r,side)
    end
end

function putmarkers!(r::Robot, side::HorizonSide)
    putmarker!(r)
    while move_if_possible!(r,side)
        putmarker!(r)
    end

    if !isborder(r,Ost)
        move!(r,Ost)
        return putmarkers!(r, inverse(side))
    end
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
function run_instance(w,ξ,η; epsconv=1.0e-4,callback=(x...)->nothing,ρ0=0.25,
    λ0=1.0,δρ=0.0,η0=1.0,δη=0.0,damp=0.99,maxiter=50000,EPInputState=nothing)
    println("ρ0=$ρ0, δρ=$δρ, η0=$η0, δη=$δη")
    ret,(Γ,Λ),runtime,MSEdB,Overlap=EP_Regress_ZeroTemp(ξ,w,
        η,λ=λ0,ρ0=ρ0,λ0=λ0,δρ=δρ,δλ=0.0,η0=η0,δη=δη,damp=damp,maxiter=maxiter,
        epsconv=epsconv,callback=callback,EPInputState=EPInputState)
    println("Teacher-student overlap: $Overlap")
    println("Teacher-student mean squared error (dB): $MSEdB")
    return ret,(Γ,Λ),runtime,MSEdB,Overlap
end

function run_instance_gaussian_iid(w,ρ,α,η; callback=(x...)->nothing,
        ρ0=ρ,λ0=1.0,δρ=0.0,η0=η,δη=0.0,damp=0.99,epsconv=1e-4,maxiter=50000)
    prior=Lambda_prior(η)
    N=length(w)
    ξ=generate_gaussian_iid_patterns(N,α)
    M=size(ξ)[1];
    initial_EP_state=initialize(M,N,η)
    ret,(Γ,Λ),runtime,MSEdB,Overlap = run_instance(w,ξ,η,epsconv=epsconv,
        callback=callback,ρ0=ρ,λ0=λ0,δρ=δρ,η0=η0,δη=δη,damp=damp,maxiter=maxiter,EPInputState=initial_EP_state)
    return ret,(Γ,Λ),runtime,ρ0,MSEdB,Overlap
end

function run_instance_gaussian_cor(w,ρ,α,η; callback=(x...)->nothing,
        ρ0=ρ,λ0=1.0,δρ=0.0,η0=η,δη=0.0,damp=0.99,epsconv=1e-4,maxiter=50000)
    prior=Lambda_prior(η)
    N=length(w)
    ξ=generate_gaussian_cor_patterns(N,α)
    M=size(ξ)[1];
    initial_EP_state=initialize(M,N,η)
    ret,(Γ,Λ),runtime,MSEdB,Overlap = run_instance(w,ξ,η,epsconv=epsconv,
        callback=callback,ρ0=ρ,λ0=λ0,δρ=δρ,η0=η0,δη=δη,damp=damp,maxiter=maxiter,EPInputState=initial_EP_state)
    return ret,(Γ,Λ),runtime,ρ0,MSEdB,Overlap
end

max_stars_repo_name: mpieropan/EPDilutedClassifier

max_issues_repo_name: mpieropan/EPDilutedClassifier

max_forks_repo_name: mpieropan/EPDilutedClassifier


text:
module ApproxFunOrthogonalPolynomials
using Base, LinearAlgebra, Reexport, BandedMatrices, BlockBandedMatrices, AbstractFFTs, FFTW, BlockArrays, FillArrays, FastTransforms, IntervalSets,
            DomainSets, Statistics, SpecialFunctions, FastGaussQuadrature

@reexport using ApproxFunBase

import AbstractFFTs: Plan, fft, ifft
import FFTW: plan_r2r!, fftwNumber, REDFT10, REDFT01, REDFT00, RODFT00, R2HC, HC2R,
                r2r!, r2r,  plan_fft, plan_ifft, plan_ifft!, plan_fft!

import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSpace, VFun, RowVector,
                    UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, WeightSpace, NoSpace, Space,
                    HeavisideSpace, PointSpace,
                    IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix,
                    AnyDomain, ZeroSpace, ArraySpace, TrivialInterlacer, BlockInterlacer,
                    AbstractTransformPlan, TransformPlan, ITransformPlan,
                    ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral,
                    ConcreteVolterra, Volterra, VolterraWrapper,
                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation, EvaluationWrapper,
                    Conversion, defaultConversion, defaultcoefficients, default_Fun, Multiplication, Derivative, Integral, bandwidths,
                    ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                    DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral, IntegralWrapper,
                    ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven, Dirichlet, ConcreteDirichlet,
                    TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                    hasfasttransform, canonicalspace, domain, setdomain, prectype, domainscompatible,
                    plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform,
                    CanonicalTransformPlan, ICanonicalTransformPlan,
                    Integral,
                    domainspace, rangespace, boundary,
                    union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points,
                    rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp,
                    linesum, differentiate, integrate, linebilinearform, bilinearform,
                    UnsetNumber, coefficienttimes, subspace_coefficients, sumspacecoefficients, specialfunctionnormalizationpoint,
                    Segment, IntervalOrSegmentDomain, PiecewiseSegment, isambiguous, Vec, eps, isperiodic,
                    arclength, complexlength,
                    invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
                    reverseorientation, checkpoints, evaluate, extrapolate, mul_coefficients, coefficients, isconvertible,
                    clenshaw, ClenshawPlan, sineshaw,
                    toeplitz_getindex, toeplitz_axpy!, sym_toeplitz_axpy!, hankel_axpy!, ToeplitzOperator, SymToeplitzOperator, hankel_getindex,
                    SpaceOperator, ZeroOperator, InterlaceOperator,
                    interlace!, reverseeven!, negateeven!, cfstype, pad!, alternatesign!, mobius,
                    extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots,splitatroots,
                    chebmult_getindex, intpow, alternatingsum,
                    domaintype, diagindshift, rangetype, weight, isapproxinteger, default_Dirichlet, scal!, dotu,
                    components, promoterangespace, promotedomainspace,
                    block, blockstart, blockstop, blocklengths, isblockbanded, pointscompatible, affine_setdiff, complexroots,
                    ℓ⁰, recα, recβ, recγ, ∞

import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
            Interval, ChebyshevInterval, boundary, ∂, rightendpoint, leftendpoint,
            dimension, WrappedDomain

import BandedMatrices: bandrange, bandshift,
                inbands_getindex, inbands_setindex!, bandwidth, AbstractBandedMatrix,
                colstart, colstop, colrange, rowstart, rowstop, rowrange,
                bandwidths, _BandedMatrix, BandedMatrix

import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype, iterate,
                >=, /, ^, \, ∪, transpose, size, tail, broadcast, broadcast!, copyto!, copy, to_index, (:),
                similar, map, vcat, hcat, hvcat, show, summary, stride, sum, cumsum, sign, imag, conj, inv,
                complex, reverse, exp, sqrt, abs, abs2, sign, issubset, values, in, first, last, rand, intersect, setdiff,
                isless, union, angle, join, isnan, isapprox, isempty, sort, merge, promote_rule,
                minimum, maximum, extrema, argmax, argmin, findmax, findmin, isfinite,
                zeros, zero, one, promote_rule, repeat, length, resize!, isinf,
                getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
                split

import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, dot, cross,
                qr, qr!, rank, isdiag, istril, istriu, issymmetric, ishermitian,
                Tridiagonal, diagm, diagm_container, factorize, nullspace,
                Hermitian, Symmetric, adjoint, transpose, char_uplo

import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_chebyshevtransform,
                        plan_chebyshevtransform!, plan_ichebyshevtransform, plan_ichebyshevtransform!,
                        pochhammer, lgamma, chebyshevtransform!, ichebyshevtransform!

import BlockBandedMatrices: blockbandwidths, subblockbandwidths

# we need to import all special functions to use Calculus.symbolic_derivatives_1arg
# we can't do importall Base as we replace some Base definitions
import SpecialFunctions: sinpi, cospi, airy, besselh,
                    asinh, acosh,atanh, erfcx, dawson, erf, erfi,
                    sin, cos, sinh, cosh, airyai, airybi, airyaiprime, airybiprime,
                    hankelh1, hankelh2, besselj, besselj0, bessely, besseli, besselk,
                    besselkx, hankelh1x, hankelh2x, exp2, exp10, log2, log10,
                    tan, tanh, csc, asin, acsc, sec, acos, asec,
                    cot, atan, acot, sinh, csch, asinh, acsch,
                    sech, acosh, asech, tanh, coth, atanh, acoth,
                    expm1, log1p, lfact, sinc, cosc, erfinv, erfcinv, beta, lbeta,
                    eta, zeta, gamma, polygamma, invdigamma, digamma, trigamma,
                    abs, sign, log, expm1, tan, abs2, sqrt, angle, max, min, cbrt, log,
                    atan, acos, asin, erfc, inv

points(d::IntervalOrSegmentDomain{T},n::Integer; kind::Int=1) where {T} =
    fromcanonical.(Ref(d), chebyshevpoints(float(real(eltype(T))), n; kind=kind))  # eltype to handle point
bary(v::AbstractVector{Float64},d::IntervalOrSegmentDomain,x::Float64) = bary(v,tocanonical(d,x))

include("bary.jl")

include("ultraspherical.jl")
include("Domains/Domains.jl")
include("Spaces/Spaces.jl")
include("roots.jl")
include("specialfunctions.jl")
include("fastops.jl")

end