text:
include_dpo()
include(joinpath(@__DIR__, "rocket_nominal.jl"))
include(joinpath(@__DIR__, "rocket_slosh.jl"))

function fd(model::RocketNominal{Midpoint, FreeTime}, x⁺, x, u, w, h, t)
	h = u[end]
    x⁺ - (x + h * f(model, 0.5 * (x + x⁺), u, w) + w)
end

function fd(model::RocketSlosh{Midpoint, FreeTime}, x⁺, x, u, w, h, t)
	h = u[end]
    x⁺ - (x + h * f(model, 0.5 * (x + x⁺), u, w) + w)
end

# DPO
β = 1.0
δ = 1.0e-3

# Initial samples
x1_sample = resample(x1_slosh, Diagonal(ones(model_sl.n)), 1.0e-3)

prob_nom = prob_nominal.prob

# Mean problem
prob_mean = trajectory_optimization(
				model_sl,
				EmptyObjective(),
				T,
				dynamics = false,
				ul = control_bounds(model_sl, T,
					[Inf * ones(2); 0.0],
					[Inf * ones(2); 0.0])[1],
				uu = control_bounds(model_sl, T,
					[Inf * ones(2); 0.0],
					[Inf * ones(2); 0.0])[2])

# Sample problems
prob_sample = [trajectory_optimization(
				model_sl,
				EmptyObjective(),
				T,
				xl = state_bounds(model_sl, T, x1 = x1_sample[i])[1],
				xu = state_bounds(model_sl, T, x1 = x1_sample[i])[2],
				ul = ul,
				uu = uu,
				dynamics = false,
				con = con_free_time
				) for i = 1:2 * model_sl.n]

# Sample objective
Q = [(t < T ? Diagonal(100.0 * ones(model_nom.n))
	: Diagonal(1000.0 * ones(model_nom.n))) for t = 1:T]
R = [Diagonal([1.0; 1.0; 100.0]) for t = 1:T-1]

obj_sample = sample_objective(Q, R)
policy = linear_feedback(6, 2,
	idx_input = collect([1, 2, 3, 5, 6, 7]),
	idx_input_nom = (1:6),
	idx_output = (1:2))
dist = disturbances([Diagonal(δ * ones(model_sl.n)) for t = 1:T-1])
sample = sample_params(β, T)

prob_dpo = dpo_problem(
	prob_nom, prob_mean, prob_sample,
	obj_sample,
	policy,
	dist,
	sample)

# TVLQR policy
x̄_nom, ū_nom = unpack(z̄_nom, prob_nominal)
K, P = tvlqr(model_nom, x̄_nom, ū_nom, 0.0, Q, R)

# Pack
z0 = zeros(prob_dpo.num_var)
z0[prob_dpo.prob.idx.nom] = copy(z̄_nom)
z0[prob_dpo.prob.idx.mean] = copy(z̄_slosh)
for i = 1:2 * model_sl.n
	z0[prob_dpo.prob.idx.sample[i]] = copy(z̄_slosh)
end
for t = 1:T-1
	z0[prob_dpo.prob.idx.policy[prob_dpo.prob.idx.θ[t]]] = vec(copy(K[t]))
end

# Solve
if false # set to true to reoptimize
	include_snopt()
	z, info = solve(prob_dpo, copy(z0),
		nlp = :SNOPT7,
		tol = 1.0e-2, c_tol = 1.0e-2,
		time_limit = 60 * 60 * 3) # NOTE: may not reach tolerances: policy still works
	@save joinpath(@__DIR__, "sol_dpo.jld2") z
else
	println("Loading solution...")
    @load joinpath(@__DIR__, "sol_dpo.jld2") z
end

x, u = unpack(z, prob_dpo.prob.prob.nom)
@show sum([u[t][end] for t = 1:T-1]) # ~2.91


text:
using Raylib
using Raylib: RayFont, RayVector2, rayvector,
    MAROON, ORANGE, DARKGREEN, DARKBLUE, DARKPURPLE, LIME, GOLD, RED

function main()
    MAX_FONTS = 8
    screenWidth = 800
    screenHeight = 450

    Raylib.InitWindow(screenWidth, screenHeight, "raylib [text] example - raylib fonts")

    # NOTE: Textures MUST be loaded after Window initialization (OpenGL context is required)
    font_files = [
        "resources/fonts/alagard.png",
        "resources/fonts/pixelplay.png",
        "resources/fonts/mecha.png",
        "resources/fonts/setback.png",
        "resources/fonts/romulus.png",
        "resources/fonts/pixantiqua.png",
        "resources/fonts/alpha_beta.png",
        "resources/fonts/jupiter_crash.png",
    ]

    fonts = map(font_files) do file
        Raylib.LoadFont(joinpath(@__DIR__, file))
    end

    messages = [
        "ALAGARD FONT designed by Hewett Tsoi",
        "PIXELPLAY FONT designed by Aleksander Shevchuk",
        "MECHA FONT designed by Captain Falcon",
        "SETBACK FONT designed by Brian Kent (AEnigma)",
        "ROMULUS FONT designed by Hewett Tsoi",
        "PIXANTIQUA FONT designed by Gerhard Grossmann",
        "ALPHA_BETA FONT designed by Brian Kent (AEnigma)",
        "JUPITER_CRASH FONT designed by Brian Kent (AEnigma)",
    ]

    spacings = [2, 4, 8, 4, 3, 4, 4, 1]

    positions = map(1:MAX_FONTS) do i
        rayvector(
            screenWidth/2 - Raylib.MeasureTextEx(fonts[i], messages[i], 2*fonts[i].baseSize, spacings[i])[1]/2,
            60 + fonts[i].baseSize + 45(i-1)
        )
    end

    positions[4] = positions[4] .+ (0, 8)
    positions[5] = positions[5] .+ (0, 2)
    positions[8] = positions[8] .- (0, 8)

    colors = [ MAROON, ORANGE, DARKGREEN, DARKBLUE, DARKPURPLE, LIME, GOLD, RED ]

    Raylib.SetTargetFPS(60)

    while !Raylib.WindowShouldClose()
        Raylib.BeginDrawing()
        Raylib.ClearBackground(Raylib.RAYWHITE)

        Raylib.DrawText("free fonts included with raylib", 250, 20, 20, Raylib.DARKGRAY)
        Raylib.DrawLine(200, 50, 590, 50, Raylib.DARKGRAY)

        for i = 1:MAX_FONTS
            Raylib.DrawTextEx(
                fonts[i], messages[i], positions[i], 2fonts[i].baseSize, spacings[i], colors[i]
            )
        end

        Raylib.EndDrawing()
    end

    foreach(Raylib.UnloadFont, fonts)
    Raylib.CloseWindow()
end

ext: jl


text:
# See interface docs
@inline inbounds(data::Union{GridData,SimData}, I::Tuple) = inbounds(data, I...)
@inline inbounds(data::Union{GridData,SimData}, I...) =
    _inbounds(boundary(data), gridsize(data), I...)

@inline function _inbounds(boundary::BoundaryCondition, size::Tuple, I...)
    reduce(zip(size, I); init=((), true)) do (I, inbounds_acc), (s, i)
        ii, inbounds = _inbounds(boundary, s, i)
        ((I..., ii), inbounds_acc & inbounds)
    end
end
@inline _inbounds(::Remove, size::Number, i::Number) = i, _isinbounds(size, i)
@inline function _inbounds(::Wrap, size::Number, i::Number)
    if i < oneunit(i)
        size + rem(i, size), true
    elseif i > size
        rem(i, size), true
    else
        i, true
    end
end

@inline isinbounds(data::Union{SimData,GridData}, I::Tuple) = isinbounds(data, I...)
@inline isinbounds(data::Union{SimData,GridData}, I...) = _isinbounds(gridsize(data), I...)

@inline _isinbounds(size::Tuple, I...) = all(map(_isinbounds, size, I))
@inline _isinbounds(size, i) = i >= one(i) && i <= size

#= Wrap boundary where required. This optimisation allows us to ignore
bounds checks on neighborhoods and still use a wraparound grid. =#
_updateboundary!(grids::Tuple) = map(_updateboundary!, grids)
function _updateboundary!(g::GridData{Y,X,R}) where {Y,X,R}
    R < 1 && return g
    return _updateboundary!(g, boundary(g))
end
function _updateboundary!(g::GridData{Y,X,R,T}, ::Remove) where {Y,X,R,T}
    src = parent(source(g))
    # Left
    @inbounds src[1:Y, 1:R] .= Ref(padval(g))
    # Right
    @inbounds src[1:Y, X+R+1:X+2R] .= Ref(padval(g))
    # Top middle
    @inbounds src[1:R, R+1:X+R] .= Ref(padval(g))
    # Bottom middle
    @inbounds src[Y+R+1:Y+2R, R+1:X+R] .= Ref(padval(g))
    return g
end
function _updateboundary!(g::GridData{Y,X,R}, ::Wrap) where {Y,X,R}
    # TODO optimise this. Its mostly a placeholder so wrapping still works in GOL tests.
    src = parent(source(g))
    nrows, ncols = gridsize(g)
    startpadrow = startpadcol = 1:R
    endpadrow = nrows+R+1:nrows+2R
    endpadcol = ncols+R+1:ncols+2R
    startrow = startcol = 1+R:2R
    endrow = nrows+1:nrows+R
    endcol = ncols+1:ncols+R
    rows = 1+R:nrows+R
    cols = 1+R:ncols+R

    # Sides ---
    # Left
    @inbounds copyto!(src, CartesianIndices((rows, startpadcol)),
                      src, CartesianIndices((rows, endcol)))
    # Right
    @inbounds copyto!(src, CartesianIndices((rows, endpadcol)),
                      src, CartesianIndices((rows, startcol)))
    # Top
    @inbounds copyto!(src, CartesianIndices((startpadrow, cols)),
                      src, CartesianIndices((endrow, cols)))
    # Bottom
    @inbounds copyto!(src, CartesianIndices((endpadrow, cols)),
                      src, CartesianIndices((startrow, cols)))
    # Corners ---
    # Top Left
    @inbounds copyto!(src, CartesianIndices((startpadrow, startpadcol)),
                      src, CartesianIndices((endrow, endcol)))
    # Top Right
    @inbounds copyto!(src, CartesianIndices((startpadrow, endpadcol)),
                      src, CartesianIndices((endrow, startcol)))
    # Botom Left
    @inbounds copyto!(src, CartesianIndices((endpadrow, startpadcol)),
                      src, CartesianIndices((startrow, endcol)))
    # Botom Right
    @inbounds copyto!(src, CartesianIndices((endpadrow, endpadcol)),
                      src, CartesianIndices((startrow, startcol)))

    _wrapstatus!(sourcestatus(g))
    return g
end

_wrapstatus!(status::Nothing) = nothing
function _wrapstatus!(status::AbstractArray)
    # This could be further optimised.
    status[end-1, :] .|= status[1, :]
    status[:, end-1] .|= status[:, 1]
    #status[end-2, :] .|= status[1, :] .|= status[2, :]
    status[end-2, :] .= true
    status[:, end-2] .|= status[:, 1] .|= status[:, 2]
    #status[1, :] .|= status[end-2, :] .|= status[end-1, :]
    status[1, :] .= true
    status[:, 1] .|= status[:, end-2] .|= status[:, end-1]
    status .= true
end


text:
using NLS_Fit
using Test, BenchmarkTools

@testset "NLS_Fit.jl" begin

    include("misc/Misc.jl")

    include("map/Map.jl")

    include("model2fit/Model2Fit.jl")

end


text:
"""
`module JuLIP.Solve`

Contains a few geometry optimisation routines, for now see
the help for these:

* `minimise`
"""
module Solve

import Optim
import Optim: DifferentiableFunction, optimize, ConjugateGradient, LBFGS

using JuLIP: AbstractAtoms, Preconditioner, update!, Identity, update!,
            dofs, energy, grad

export minimise!

"""
`minimise!(at::AbstractAtoms)`: geometry optimisation

`at` must have a calculator and a constraint attached.

## Keyword arguments:
* `precond = Identity()` : preconditioner
* `grtol = 1e-6`
* `ftol = 1e-32`
* `Optimiser = Optim.ConjugateGradient` (currently this is ignored)
* `verbose = 0`: 0 : no output, 1 : final, 2 : iteration
"""
function minimise!( at::AbstractAtoms;
                  precond = Identity(), gtol=1e-6, ftol=1e-32,
                  Optimiser = ConjugateGradient,
                  verbose = 1 )

   # create an objective function
   objective = DifferentiableFunction( x->energy(at, x),
                                       (x,g)->copy!(g, grad(at, x)) )
   # call Optim.jl
   # TODO: use verb flag to determine whether detailed output is wanted
   if isa(precond, Identity)
      optimiser = Optim.ConjugateGradient()
   else
      optimiser = Optim.LBFGS( P = precond,
                        precondprep! = (P, x) -> update!(P, at, x),
                        linesearch! = Optim.quadbt_linesearch! )
   end
   results = optimize( objective, dofs(at), method = optimiser,
                        f_tol = ftol, g_tol = gtol, show_trace = (verbose > 1) )
   # analyse the results
   if verbose > 0
      println(results)
   end
   return results
end

# saddle search
# include("saddlesearch.jl")

end


text:
export RedisCondition

mutable struct RedisCondition{T}
    key::String
    socket::IO
    listen::Bool

    RedisCondition{T}(channel::AbstractString) where T = RedisCondition{T}(default_connection, channel)

    RedisCondition{T}(rc::RedisConnection, key) where T = if serializeable(T)
        rc.busy && wait(rc.lock)
        rc.busy = true

        c = new{T}(key, rc.socket, false)

        finalizer(c, (x)->try
            c.listen && resp_send(rc.socket, "unsubscribe", key) |> resp_read
        finally
            rc.busy = false
            notify(rc.lock, all=false)
        end)

        c
    else
        throw(ArgumentError("RedisCondition currently not supports arbitrary element type"))
    end

    RedisCondition{T}(rcp::RedisConnectionPool, key) where T = if serializeable(T)
        socket = find_avaliable(rcp)

        c = new{T}(key, socket, false)

        finalizer(c, (x)->try
            c.listen && resp_send(socket, "unsubscribe", key) |> resp_read
        finally
            push!(rcp.queue, socket)
            notify(rcp.cond, all=false)
        end)

        c
    else
        throw(ArgumentError("RedisCondition currently not supports arbitrary element type"))
    end
end

RedisCondition(args...) = RedisCondition{String}(args...)

function wait(rc::RedisCondition{T}) where T
    if !rc.listen
        resp_send(rc.socket, "subscribe", rc.key) |> resp_read
        rc.listen = true
    end
    _, _, x = resp_read(rc.socket)
    deserialize(T, x)
end

"wake up all"
function notify(rc::RedisCondition{T}, val="") where T
    rc.listen && error("cannot notify a listening condition")
    resp_send(rc.socket, "publish", rc.key, serialize(T, val)) |> resp_read
end


text:
# Matern 1/2 ARD covariance Function

"""
    Mat12Ard <: MaternARD

ARD Matern 1/2 kernel (covariance)
```math
k(x,x') = σ² \\exp(-|x-x'|/L)
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)`` and signal standard deviation ``σ`` where
``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct Mat12Ard{T} <: MaternARD where {T<:Real}
    "Inverse squared length scale"
    iℓ2::Vector{T}
    "Signal variance"
    σ2::T
    "Priors for kernel parameters"
    priors::Array
end

"""
Matern 1/2 ARD covariance Function

    Mat12Ard(ll::Vector{T}, lσ::T)

# Arguments
  - `ll::Vector{Real}`: vector of length scales (given on log scale)
  - `lσ::Real`: signal standard deviation (given on log scale)
"""
Mat12Ard(ll::Vector{T}, lσ::T) where T = Mat12Ard{T}(exp.(-2 .* ll), exp(2 * lσ), [])

function set_params!(mat::Mat12Ard, hyp::AbstractVector)
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat12 kernel has $(num_params(mat)) parameters, received $(length(hyp))."))
    mat.iℓ2  = exp.(-2 .* @view(hyp[1:(end-1)]))
    mat.σ2 = exp(2 * hyp[end])
end

get_params(mat::Mat12Ard) = [-log.(mat.iℓ2) / 2; log(mat.σ2) / 2]
get_param_names(mat::Mat12Ard) = [get_param_names(mat.iℓ2, :ll); :lσ]
num_params(mat::Mat12Ard) = length(mat.iℓ2) + 1

cov(mat::Mat12Ard, r::Number) = mat.σ2 * exp(-r)

dk_dll(mat::Mat12Ard, r::Real, wdiffp::Real) = wdiffp / r * cov(mat,r)


text:
using BEAST
using CompScienceMeshes
using Test

# function Base.length(it::BEAST.ADIterator)
#     l = 0
#     for (m,w) in it
#         l += 1
#     end
#     l
# end
#
# tbf = BEAST.timeaxisc0d1(0,10.0,50)
#
# @test numfunctions(tbf) == 49
# @test refspace(tbf) == LagrangeRefSpace{Float64,1}()
#
# ad = assemblydata(tbf)
# @test size(ad.data) == (1,2,50)
# @test length(ad[1,1]) == 0
# @test length(ad[1,2]) == 1
# @test length(ad[50,1]) == 1
# @test length(ad[50,2]) == 0
# for c in 2:49
#     for r in 1:2
#         @test length(ad[c,r]) == 1
#     end
# end

tbf = BEAST.timebasisc0d1(1.0, 10)
timeels, timead = BEAST.assemblydata(tbf)
#ad = BEAST.temporalassemblydata(tbf)
;


text:
## (c) 2016 David A. van Leeuwen

## Unit tests for ../src/constructors.jl

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

using DataStructures
using Base.Test

print("construction, ")
## constructors
a = rand(2,3)
n1 = NamedArray(a, (OrderedDict(:a => 1, :b=>2), OrderedDict(:c=>3, :d=>2, :e=>1)), (:dim1, :dim2))
println(n1, " ", n1[1,1], " ", n1[:a, :e])
n2 = NamedArray(a, (OrderedDict(:a => 1, :b=>2), OrderedDict(:c=>3, :d=>2, :e=>1)))
println(n2)
n3 = NamedArray(a, ([:a, :b], [:c, :d, :e]), (:dim1, :dim2))
println(n3)
n4 = NamedArray(a, ([:a, :b], [:c, :d, :e]))
println(n4)
n5 = NamedArray(a, [[Char(64+i) for i in 1:d] for d in 2:3], ["a", "b"])
n6 = NamedArray(a, [[Char(64+i) for i in 1:d] for d in 2:3])
n7 = NamedArray(a)

n1 = NamedArray(Complex64, 5, 8)
n2 = NamedArray(a, (["s", "t"],[:a, :b, :c]), ("string", :symbol))

n = NamedArray(rand(2,4))
setnames!(n, ["one", "two"], 1)
setnames!(n, ["a", "b", "c", "d"], 2)

a = [1 2 3; 4 5 6]
n3 = NamedArray(a, (["a","b"],["C","D","E"]))
n4 = NamedArray(a, (@compat OrderedDict("a"=>1,"b"=>2), @compat OrderedDict("C"=>1,"D"=>2,"E"=>3)))

@test n3.array == n4.array == a
@test dimnames(n3) == dimnames(n4) == Any[:A,:B]
@test names(n3,1) == names(n4,1) == ["a","b"]
@test names(n3,2) == names(n4,2) == ["C","D","E"]

## 0-dim case #21
if VERSION ≥ v"0.4-dev"
    n0 = NamedArray(Array{Int}())
    @test size(n0) == ()
    @test n0[1] == n0.array[1]
end

## Calling constructors through convert, #38
@test convert(NamedArray, n.array) == NamedArray(n.array)
@test convert(NamedVector, [1, 2, 3]) == NamedArray([1, 2, 3])


text:
using Optimize
using NLPModels

using Plots
# pyplot()
plotly()

n_min = 2
n_max = 200

# Math Prog Base collection
include("MPBProblems.jl")
test_probs = mpb_probs

# Select solvers
using LineSearch
using Stopping
using LSDescentMethods

#stop_norm = ARCTR.stop_norm
atol = 1.0e-5
rtol = 1.0e-10

solvers = [Newton, Newton, Newton, Newton, Newton]
labels = ["Newton Armijo backtracking";"Newton HZ"; "Newton MT"; "Newton TR_SecA"; "Newton TR_Nwt"]

stp = TStopping(atol = atol, rtol = rtol, max_iter = 1000000, max_eval = 1000000, max_time = 100.0)

options = [Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => Newarmijo_wolfe)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => _hagerzhang2!)#, :τ₀ => 0.01, :mayterminate => true, :τ₁ => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => _morethuente2!)#, :τ₀ => 0.01, :τ₁ => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => TR_SecA_ls, :τ₀ => 0.000932010000000, :τ₁ => 0.206998804000000)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => TR_Nwt_ls, :τ₀ => 0.200018755000000, :τ₁ => 0.999996368310000)
           ]

include("compare_solvers.jl")
P5 = compare_solvers_with_options2(solvers, options, labels, test_probs, n_min, n_max, printskip = false)[2]


text:
using Dates, .CentralBanks

CentralBanks.ECB_FxRates("//datos02/9763-AnalisisyPlanificacion_Financiera/D. Datos/7. Tipos", Dates.Date(2005))

CentralBanks.downloadchapter(4,"c:\\Users\\Luis\\OneDrive\\Documents\\JuliaProjects\\CentralBanks")
CentralBanks.downloadchapter(2,"c:\\Users\\Luis\\OneDrive\\Documents\\JuliaProjects\\CentralBanks")
CentralBanks.downloadchapter(19,"c:\\Users\\Luis\\OneDrive\\Documents\\JuliaProjects\\CentralBanks")

CentralBanks.downloadInterestRates("c:\\Users\\Luis\\OneDrive\\Documents\\JuliaProjects\\CentralBanks")
dfs = CentralBanks.CSVtoDataFrame("c:\\Users\\Luis\\OneDrive\\Documents\\JuliaProjects\\CentralBanks\\tipos",Date(2005))
sort!(dfs["tInteres"], 1,rev = true) # Va a fallar porque no está el using DataFrame
CentralBanks.DataFrametoCSV(dfs["tInteres"], "c:\\Users\\Luis\\OneDrive\\Documents\\JuliaProjects\\CentralBanks\\tipos","tinteres.csv")

CentralBanks.combinaciones(ans)


text:
using Evolutionary
using Test
using Random
using LinearAlgebra
Random.seed!(9874984737484)

for tests in [
            "sphere.jl",
            "rosenbrock.jl",
            "schwefel.jl",
            "rastrigin.jl",
            "n-queens.jl",
            "knapsack.jl",
            "selections.jl"
]
    include(tests)
end


text:
@testset "Schwefel CMA-ES" begin

    rng = StableRNG(42)

    # Schwefel's ellipsoid test function:
    # a moderately conditioned ellipsoid with a dominating isolated eigenvalue
    function schwefel(x::AbstractVector{T}) where {T<:AbstractFloat}
        out = zero(T)
        for i = 1:length(x)
            out += sum(abs2, x[1:i])
        end
        return out
    end

    # Parameters
    N = 30
    function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::CMAES, options)
        record["σ"] = state.σ
    end

    # Testing: CMA-ES
    opts = Evolutionary.Options(rng=rng)
    λ = 14
    Random.seed!(rng, 42)
    result = Evolutionary.optimize(schwefel, ()->rand(rng, N), CMAES(lambda=λ, c_1=0.05), opts)
    println("($(λ>>1),$λ)-CMA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.converged(result)
    @test Evolutionary.minimizer(result) ≈ zeros(N) atol=1e-4
    @test minimum(result) ≈ 0.0 atol=1e-5

    bounds = Evolutionary.ConstraintBounds(fill(-1.0f0,N),fill(1.0f0,N),[],[])
    opts = Evolutionary.Options(store_trace=true, iterations=10, rng=rng)
    Random.seed!(rng, 42)
    result = Evolutionary.optimize(schwefel, bounds, CMAES(lambda=λ, weights=zeros(Float32,λ)), opts)
    println("($(λ>>1),$λ)-CMA-ES [bounds] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.iterations(result) == 10
    @test !Evolutionary.converged(result)
    @test haskey(Evolutionary.trace(result)[end].metadata, "σ")
end


text:
module Distributed

export
    MultiCPU, MultiGPU, child_architecture,
    HaloCommunication, HaloCommunicationBC,
    inject_halo_communication_boundary_conditions,
    DistributedFFTBasedPoissonSolver,
    DistributedIncompressibleModel, DistributedShallowWaterModel

using MPI

using Oceananigans.Models
using Oceananigans.Utils

include("distributed_utils.jl")
include("multi_architectures.jl")
include("halo_communication_bcs.jl")
include("halo_communication.jl")
include("distributed_fields.jl")
include("distributed_fft_based_poisson_solver.jl")
include("distributed_solve_for_pressure.jl")
include("distributed_incompressible_model.jl")
include("distributed_shallow_water_model.jl")

end # module


text:
"""
    auxiliary_variable_bounds(v::Array{JuMP.VariableRef,1})

Given a vector of JuMP variables (maximum 4 variables), this function returns the worst-case
bounds, the product of these input variables can admit.
"""
function auxiliary_variable_bounds(v::Array{JuMP.VariableRef,1})

    v_l = [JuMP.lower_bound(v[1]), JuMP.upper_bound(v[1])]
    v_u = [JuMP.lower_bound(v[2]), JuMP.upper_bound(v[2])]
    M = v_l * v_u'

    if length(v) == 2
       # Bilinear
        return minimum(M), maximum(M)

    elseif (length(v) == 3) || (length(v) == 4)
       M1 = zeros(Float64, (2, 2, 2))
       M1[:,:,1] = M * JuMP.lower_bound(v[3])
       M1[:,:,2] = M * JuMP.upper_bound(v[3])
       if length(v) == 3
          # Trilinear
          return minimum(M1), maximum(M1)
       elseif length(v) == 4
          # Quadrilinear
          M2 = zeros(Float64, (2, 2, 4))
          M2[:,:,1:2] = M1 * JuMP.lower_bound(v[4])
          M2[:,:,3:4] = M1 * JuMP.upper_bound(v[4])
          return minimum(M2), maximum(M2)
       end
    end
end

"""
    gate_element_bounds(M::Array{Float64,3})

Given a set of elementary gates, `{G1, G2, ... ,Gn}`, this function evaluates
the range of every co-ordinate of the superimposed gates, over all possible gates.
"""
function gate_element_bounds(M::Array{Float64,3})

    M_l = zeros(size(M)[1], size(M)[2])
    M_u = zeros(size(M)[1], size(M)[2])

    for i = 1:size(M)[1]
        for j = 1:size(M)[2]

            M_l[i,j] = minimum(M[i,j,:])
            M_u[i,j] = maximum(M[i,j,:])

            if M_l[i,j] > M_u[i,j]
                Memento.error(_LOGGER, "Lower and upper bound conflict in the elements of input elementary gates")
            elseif (M_l[i,j] in [-Inf, Inf]) || (M_u[i,j] in [-Inf, Inf])
                Memento.error(_LOGGER, "Unbounded entry detected in the input elementary gates")
            end
        end
    end

    return M_l, M_u
end

"""
    get_commutative_gate_pairs(M::Dict{String,Any}; identity_in_pairs = true)

Given a dictionary of elementary quantum gates, this function returns all pairs of commuting
gates. Optional argument, `identity_pairs` can be set to `false` if identity matrix need not be part of the commuting pairs.
"""
function get_commutative_gate_pairs(M::Dict{String,Any}; identity_in_pairs = true)

    num_gates = length(keys(M))
    commute_pairs = Array{Tuple{Int64,Int64},1}()
    commute_pairs_prodIdentity = Array{Tuple{Int64,Int64},1}()

    for i = 1:(num_gates-1)
        for j = (i+1):num_gates
            M_i = M["$i"]["matrix"]
            M_j = M["$j"]["matrix"]

            if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
                continue
            end

            M_ij = M_i*M_j
            M_ji = M_j*M_i
            Id = Matrix(LA.I, size(M_ij)[1], size(M_ij)[2])

            # Commuting pairs == Identity
            if isapprox(M_ij, Id, atol=1E-6)
                push!(commute_pairs_prodIdentity, (i,j))

            # Commuting pairs != Identity
            elseif isapprox(M_ij, M_ji, atol = 1E-4)
                push!(commute_pairs, (i, j))

            end

        end
    end

    if identity_in_pairs
        # Commuting pairs involving Identity
        identity_idx = []
        for i in keys(M)
            if "Identity" in M[i]["type"]
                push!(identity_idx, parse(Int64, i))
            end
        end

        if length(identity_idx) > 0
            for i = 1:length(identity_idx)
                for j = 1:num_gates
                    if j != identity_idx[i]
                        push!(commute_pairs, (j, identity_idx[i]))
                    end
                end
            end
        end

    end

    return commute_pairs, commute_pairs_prodIdentity
end

"""
    get_redundant_gate_product_pairs(M::Dict{String,Any})

Given a dictionary of elementary quantum gates, this function returns all pairs of gates whose product is
one of the input elementary gates. For example, let `G_basis = {G1, G2, G3}` be the elementary gates. If `G1*G2 ∈ G_basis`,
then `(1,2)` is considered as a redundant pair.
"""
function get_redundant_gate_product_pairs(M::Dict{String,Any})
    num_gates = length(keys(M))
    redundant_pairs_idx = Array{Tuple{Int64,Int64},1}()

    # Non-Identity redundant pairs
    for i = 1:(num_gates-1)
        for j = (i+1):num_gates
            M_i = M["$i"]["matrix"]
            M_j = M["$j"]["matrix"]

            if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
                continue
            end

            for k = 1:num_gates
                if (k != i) && (k != j) && !("Identity" in M["$k"]["type"])
                    M_k = M["$k"]["matrix"]

                    if isapprox(M_i*M_j, M_k, atol = 1E-4)
                        push!(redundant_pairs_idx, (i, j))
                    end

                    if isapprox(M_j*M_i, M_k, atol = 1E-4)
                        push!(redundant_pairs_idx, (j, i))
                    end
                end
            end

        end
    end

    return redundant_pairs_idx
end

"""
    get_idempotent_gates(M::Dict{String,Any})

Given the dictionary of complex quantum gates, this function returns the indices of matrices which are self-idempotent
or idempotent with other set of input gates, excluding the Identity gate.
"""
function get_idempotent_gates(M::Dict{String,Any})
    num_gates = length(keys(M))
    idempotent_gates_idx = Vector{Int64}()

    # Excluding Identity gate in input
    for i = 1:num_gates
        M_i = M["$i"]["matrix"]
        for j = 1:num_gates
            M_j = M["$j"]["matrix"]
            if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
                continue
            end

            if isapprox(M_i^2, M_j, atol=1E-4)
                push!(idempotent_gates_idx, i)
            end

        end
    end

    return idempotent_gates_idx
end

"""
    get_involutory_gates(M::Dict{String,Any})

Given the dictionary of complex gates `G_1, G_2, ..., G_n`, this function returns the indices of these gates
which are involutory, i.e, `G_i^2 = Identity`, excluding the Identity gate.
"""
function get_involutory_gates(M::Dict{String,Any})
    num_gates = length(keys(M))
    involutory_gates_idx = Vector{Int64}()

    if num_gates > 0
        n_r = size(M["1"]["matrix"])[1]
        n_c = size(M["1"]["matrix"])[2]
    end
    Id = Matrix{ComplexF64}(LA.I, n_r, n_c)

    # Excluding Identity gate in input
    for i=1:num_gates
        if !("Identity" in M["$i"]["type"]) && (isapprox((M["$i"]["matrix"])^2, Id, atol=1E-5))
            push!(involutory_gates_idx, i)
        end
    end

    return involutory_gates_idx
end

"""
    complex_to_real_gate(M::Array{Complex{Float64},2})

Given a complex-valued two-dimensional quantum gate of size NxN, this function returns a real-valued gate
of dimensions 2Nx2N.
"""
function complex_to_real_gate(M::Array{Complex{Float64},2})

    n = size(M)[1]
    M_real = zeros(2*n, 2*n)

    ii = 1; jj = 1;
    for i = collect(1:2:2*n)
        for j = collect(1:2:2*n)

            if isapprox(real(M[ii,jj]), 0, atol=1E-6)
                M_real[i,j] = 0
                M_real[i+1,j+1] = 0
            else
                M_real[i,j] = real(M[ii,jj])
                M_real[i+1,j+1] = real(M[ii,jj])
            end

            if isapprox(imag(M[ii,jj]), 0, atol=1E-6)
                M_real[i,j+1] = 0
                M_real[i+1,j] = 0
            else
                M_real[i,j+1] = imag(M[ii,jj])
                M_real[i+1,j] = -imag(M[ii,jj])
            end

            jj += 1
        end
        jj = 1
        ii += 1
    end

    return M_real
end

"""
    real_to_complex_gate(M::Array{Complex{Float64},2})

Given a real-valued two-dimensional quantum gate of size 2Nx2N, this function returns a complex-valued gate
of size NxN, if the input gate is in a valid complex form.
"""
function real_to_complex_gate(M::Array{Float64,2})

    n = size(M)[1]

    if !iseven(n)
        Memento.error(_LOGGER, "Specified gate can admit only even numbered columns and rows")
    end

    M_complex = zeros(Complex{Float64}, (Int(n/2), Int(n/2)))

    ii = 1; jj = 1;
    for i = collect(1:2:n)
        for j = collect(1:2:n)

            if !isapprox(M[i,j], M[i+1, j+1], atol = 1E-5) || !isapprox(M[i+1,j], -M[i,j+1], atol = 1E-5)
                Memento.error(_LOGGER, "Specified real form of the complex gate is invalid")
            end

            M_re = M[i,j]
            M_im = M[i,j+1]
            (isapprox(M_re, 0, atol=1E-6)) && (M_re = 0)
            (isapprox(M_im, 0, atol=1E-6)) && (M_im = 0)

            M_complex[ii,jj] = complex(M_re, M_im)
            jj += 1
        end
        jj = 1
        ii += 1
    end

    return M_complex
end

"""
    round_complex_values(M::Array{Complex{Float64},2})

Given a complex-valued gate, this function returns a complex-valued gate which
rounds the values closest to 0 and 1. This is useful to avoid numerical issues.
"""
function round_complex_values(M::Array{Complex{Float64},2})
    if size(M)[1] == 0
        Memento.error(_LOGGER, "Input cannot be a scalar")
    end

    if length(size(M)) == 2
        n_r = size(M)[1]
        n_c = size(M)[2]

        M_round = Array{Complex{Float64},2}(zeros(n_r,n_c))

        for i=1:n_r
            for j=1:n_c
                M_round[i,j] = complex(QCO.round_real_value(real(M[i,j])), QCO.round_real_value(imag(M[i,j])))
            end
        end
        return M_round
    end

end

"""
    round_real_value(x::T) where T <: Number

Given a real-valued number, this function returns a real-value which rounds the values closest to 0 and 1.
"""
function round_real_value(x::T) where T <: Number
    if isapprox(abs(x), 0, atol=1E-6)
        x = 0
    elseif isapprox(x, 1, atol=1E-6)
        x = 1
    end

    return x
end

"""
    unique_idx(x::AbstractArray{T})

This function returns the indices of unique elements in a given array of scalar or vector inputs. Overall,
this function computes faster than Julia's built-in `findfirst` command.
"""
function unique_idx(x::AbstractArray{T}) where T
    uniqueset = Set{T}()
    ex = eachindex(x)
    idxs = Vector{eltype(ex)}()
    for i in ex
        xi = x[i]
        if !(xi in uniqueset)
            push!(idxs, i)
            push!(uniqueset, xi)
        end
    end
    idxs
end

"""
    unique_matrices(M::Array{Float64, 3})

This function returns the unique set of matrices and the corresponding indices
of unique matrices from the given set of matrices.
"""
function unique_matrices(M::Array{Float64, 3})
    M[isapprox.(M, 0, atol=1E-6)] .= 0

    M_reshape = [];
    for i=1:size(M)[3]
        push!(M_reshape, round.(reshape(M[:,:,i], size(M)[1]*size(M)[2]), digits=5))
    end

    idx = QCO.unique_idx(M_reshape)

    return M[:,:,idx], idx
end

"""
    kron_single_qubit_gate(num_qubits::Int64, M::Array{Complex{Float64},2}, qubit_loc::String)

Given number of qubits of the circuit, the complex-valued one-qubit gate and the qubit location ("q1","q2',"q3",...),
this function returns a full-sized gate after applying appropriate kronecker products. This function supports any number
integer-valued qubits.
"""
function kron_single_qubit_gate(num_qubits::Int64, M::Array{Complex{Float64},2}, qubit_loc::String)

    if size(M)[1] != 2
        Memento.error(_LOGGER, "Input should be an one-qubit gate")
    end

    qubit = parse(Int, qubit_loc[2:end])

    if !(qubit in 1:num_qubits)
        Memento.error(_LOGGER, "Specified qubit location, $qubit, has to be ∈ [q1,...,q$num_qubits]")
    end

    I = QCO.IGate(1)
    M_kron = 1

    for i = 1:num_qubits
        M_iter = I
        if i == qubit
            M_iter = M
        end
        M_kron = kron(M_kron, M_iter)
    end

    if size(M_kron)[1] == 2^(num_qubits)
        return M_kron
    else
        Memento.error(_LOGGER, "Dimensions mismatch in evaluation of Kronecker of single-qubit gate")
    end

end

"""
    kron_two_qubit_gate(num_qubits::Int64, M::Array{Complex{Float64},2}, c_qubit_loc::String, t_qubit_loc::String)

Given number of qubits of the circuit, the complex-valued two-qubit gate and the control and target qubit locations ("q1","q2',"q3",...),
this function returns a full-sized gate after applying appropriate kronecker products. This function supports any number integer-valued qubits.
"""
function kron_two_qubit_gate(num_qubits::Int64, M::Array{Complex{Float64},2}, c_qubit_loc::String, t_qubit_loc::String)

    if size(M)[1] != 4
        Memento.error(_LOGGER, "Input should be a two-qubit gate")
    end

    c_qubit = parse(Int, c_qubit_loc[2:end])
    t_qubit = parse(Int, t_qubit_loc[2:end])

    if !(c_qubit in 1:num_qubits) || !(t_qubit in 1:num_qubits)
        Memento.error(_LOGGER, "Specified control and target qubit locations have to be ∈ [q1,...,q$num_qubits]")
    elseif isapprox(c_qubit, t_qubit, atol = 1E-6)
        Memento.error(_LOGGER, "Control and target qubits cannot be identical for a multi-qubit elementary gate")
    end

    I = QCO.IGate(1)
    Swap = QCO.SwapGate()

    # If on adjacent qubits
    if abs(c_qubit - t_qubit) == 1

        M_kron = 1
        for i = 1:(num_qubits-1)
            M_iter = I
            if (i in [c_qubit, t_qubit]) && (i+1 in [c_qubit, t_qubit])
                M_iter = M
            end
            M_kron = kron(M_kron, M_iter)
        end

    # If on non-adjacent qubits
    # Assuming the qubit numbering starts at 1, and not 0
    elseif abs(c_qubit - t_qubit) >= 2
        ct = abs(c_qubit - t_qubit)

        M_sub_depth = 2*ct - 1
        M_sub_loc = ct
        M_sub_swap_id = ct
        M_sub = Matrix{Complex{Float64}}(LA.I, 2^(ct+1),2^(ct+1))

        # Idea is to represent gate_14 = swap_34 * swap_23 * gate_12 * swap_23 * swap_34
        for d=1:M_sub_depth
            M_sub_kron = 1
            swap_qubits = true

            for i=1:ct
                M_sub_iter = I

                if (d == M_sub_loc) && (i == 1)
                    M_sub_iter = M
                elseif !(d == M_sub_loc) && (i == M_sub_swap_id) && swap_qubits
                    M_sub_iter = Swap
                    if (d < M_sub_loc) && (M_sub_swap_id > 2)
                        M_sub_swap_id -= 1
                    elseif (d > M_sub_loc)
                        M_sub_swap_id += 1
                    end
                    swap_qubits = false
                end

                M_sub_kron = kron(M_sub_kron, M_sub_iter)
            end

            M_sub *= M_sub_kron
        end

        M_kron = 1
        i = 1
        while i <= num_qubits
            M_iter = I
            if (i in [c_qubit, t_qubit])
                M_iter = M_sub
                i += (ct+1)
            else
                i += 1
            end
            M_kron = kron(M_kron, M_iter)
        end

    end

    if size(M_kron)[1] == 2^(num_qubits)
        return M_kron
    else
        Memento.error(_LOGGER, "Dimensions mismatch in evaluation of Kronecker of two-qubit gate")
    end

end

"""
    _parse_gates_with_kron_symbol(s::String)

Given a string with gates separated by kronecker symbols (x), this function parses and returns the vector of gates. For
example, if the input string is `H_1xCNot_2_3xT_4`, the output will be `Vector{String}(["H_1", "CNot_2_3", "T_4"])`.
"""
function _parse_gates_with_kron_symbol(s::String)

    gates = Vector{String}()
    gate_id = string()

    for i = 1:length(s)
       if s[i] != kron_symbol
          gate_id = gate_id * s[i]
       else
          push!(gates, gate_id)
          (i != length(s)) && (gate_id = string())
       end

       if i == length(s)
          push!(gates, gate_id)
       end
    end

    return gates
 end

"""
    _parse_gate_string(s::String)

Given a string representing a single gate with qubit numbers separated by symbol `_`, this function parses and returns the vector of qubits on
which the input gate is located. For example, if the input string is `CRX_2_3`, the output will be `Vector{Int64}([2,3])`.
"""
function _parse_gate_string(s::String; type=false, qubits=false)

    if occursin(kron_symbol, s)
        Memento.error(_LOGGER, "Kron symbol is not supported for parsing qubit numbers in $s")
    end

    gates = Vector{String}()
    gate_id = string()

    for i = 1:length(s)
       if s[i] != qubit_separator
          gate_id = gate_id * s[i]
       else
          push!(gates, gate_id)
          (i != length(s)) && (gate_id = string())
       end

       if (i == length(s)) && (s[i] != qubit_separator)
          push!(gates, gate_id)
       end
    end

    if type && qubits
        return gates[1], parse.(Int, gates[2:end]) # Assuming 1st element is the gate type/name
    elseif type
        return gates[1]
    elseif qubits
        return parse.(Int, gates[2:end])
    end

 end

"""
    is_gate_real(M::Array{Complex{Float64},2})

Given a complex-valued quantum gate, M, this function returns if M has purely real parts or not as it's elements.
"""
function is_gate_real(M::Array{Complex{Float64},2})
    M_imag = imag(M)
    n_r = size(M_imag)[1]
    n_c = size(M_imag)[2]

    if sum(isapprox.(M_imag, zeros(n_r, n_c), atol=1E-6)) == n_r*n_c
        return true
    else
        return false
    end
 end

"""
    _get_constraint_slope_intercept(vertex1::Vector{<:Number}, vertex2::Vector{<:Number})

Given co-ordinates of two points in a plane, this function returns the slope (m) and intercept (c) of the
line joining these two points.
"""
function _get_constraint_slope_intercept(vertex1::Tuple{<:Number, <:Number}, vertex2::Tuple{<:Number, <:Number})

    if isapprox.(vertex1, vertex2, atol=1E-6) == [true, true]
        Memento.warn(_LOGGER, "Invalid slope and intercept for two identical vertices")
        return
    end

    if isapprox(vertex1[1], vertex2[1], atol = 1E-6)
        return Inf, Inf
    else
        m = QCO.round_real_value((vertex2[2] - vertex1[2]) / (vertex2[1] - vertex1[1]))
        c = QCO.round_real_value(vertex1[2] - (m * vertex1[1]))

        return m,c
    end

 end

"""
    is_multi_qubit_gate(gate::String)

Given the input gate string, this function returns a boolean if the input gate is a multi qubit gate or not.
For example, for a 2-qubit gate `CRZ_1_2`, output is `true`.
"""
 function is_multi_qubit_gate(gate::String)

    if occursin(kron_symbol, gate)
        return true
    end

    qubit_loc = QCO._parse_gate_string(gate, qubits = true)

    if length(qubit_loc) > 1
        return true
    elseif length(qubit_loc) == 1
        return false
    else
        Memento.error(_LOGGER, "Atleast one qubit has to be specified for an input gate")
    end

end

function _verify_θ_bounds(angle::Number)
    if !(-2*π <= angle <= 2*π)
        Memento.error(_LOGGER, "θ angle is not within valid bounds")
    end
end

function _verify_ϕ_bounds(angle::Number)
    if !(-2*π <= angle <= 2*π)
        Memento.error(_LOGGER, "ϕ angle is not within valid bounds")
    end
end

function _verify_λ_bounds(angle::Number)
    if !(-2*π <= angle <= 2*π)
        Memento.error(_LOGGER, "λ angle is not within valid bounds")
    end
end

"""
    _determinant_test_for_infeasibility(data::Dict{String,Any})

Given the processed data dictionary, this function performs a few simple tests based on the determinant values of
elementary and target gates to detect MIP infeasibility.
"""
function _determinant_test_for_infeasibility(data::Dict{String,Any})

    if data["are_gates_real"]
        det_target = LA.det(data["target_gate"])
    else
        det_target = LA.det(QCO.real_to_complex_gate(data["target_gate"]))
    end

    if isapprox(imag(det_target), 0, atol = 1E-6)
        if isapprox(det_target, -1, atol=1E-6)
            sum_det = 0
            for k = 1:length(keys(data["gates_dict"]))
                det_val = LA.det(data["gates_dict"]["$k"]["matrix"])
                if isapprox(imag(det_val), 0, atol = 1E-6)
                    sum_det += det_val
                end
            end

            if (isapprox(sum_det, length(keys(data["gates_dict"])), atol = 1E-6)) && (data["decomposition_type"] == "exact")
                Memento.error(_LOGGER, "Infeasible decomposition: det.(elementary_gates) = 1, while det(target_gate) = -1")
            end
        end
    else
        det_gates_real = true
        for k = 1:length(keys(data["gates_dict"]))
            if !(isapprox(imag(LA.det(data["gates_dict"]["$k"]["matrix"])), 0, atol=1E-6))
                det_gates_real = false
                continue
            end
        end

        if det_gates_real && (data["decomposition_type"] == "exact")
            Memento.error(_LOGGER, "Infeasible decomposition: det.(elementary_gates) = real, while det(target_gate) = complex")
        end

    end

end


text:
using UnicodePlots: heatmap
part1(input) = length(fold(input, part = 1))

function part2(input)
    coords = fold(input, part = 2)
    m = zeros(Int, 6, 39)
    for (x, y) in coords
        m[6 - y, x + 1] = 1
    end
    return heatmap(m)
end

function fold(input; part)
    coords = Tuple{Int, Int}[]
    for line in eachline(input)
        ',' in line && push!(coords, Tuple(parse.(Int, split(line, ","))))
        startswith(line, "fold along x") && fold_x!(coords, parse(Int, last(split(line, "="))))
        startswith(line, "fold along y") && fold_y!(coords, parse(Int, last(split(line, "="))))
        part == 1 && startswith(line, "fold") && break
    end
    return coords
end

fold_x!(coords, pos) = unique!(map!(c -> (mirror(c[1], pos), c[2]), coords, coords))
fold_y!(coords, pos) = unique!(map!(c -> (c[1], mirror(c[2], pos)), coords, coords))
mirror(x, pos) = x < pos ? x : 2 * pos - x


text:
# Included by symbolic_solver.jl

# In this module are methods to compute rows of the equilibrium matrix.
# Each row method has it's equation documented in LaTeX. The more expanded
# notation is designed to better represent the way the code is written. In this
# form, the left-most sum index in each term on the left-hand-side corresponds
# to one column of the solution vector, i.e.:
# 1. free chemical potentials for pure elements
# 2. free potentials
# 3. free phase amounts
# Note that the right-most term in each sum is the term in the solution vector,
# i.e. μ_B, ΔPot, Δℵ

# TODO: all the row methods could use some simplification in the arguments to
# not depend on the fixed/free amounts that are not relevant. It would also be
# nice to make sure they work regardless of whether the inputs are symbolic or
# numeric.

# TODO: The residual terms for N, N_A, and X(B) subtract the current value from the
# prescribed value. It's not clear to me mathematically why this is, but it is easy to
# verify for an ideal system that is not mass balanced that this should be the case.

@doc raw"""
    c_iG

# Equation

```math
c^\alpha_{iG} = -\sum_j e^\alpha_{ij} \frac{\partial G^\alpha_M}{\partial y^\alpha_j}
```
"""
function c_iG(phase_record, i)
    # c^\alpha_{iG} = -\sum_j e^\alpha_{ij} \frac{\partial G^\alpha_M}{\partial y^\alpha_j}
    total = 0.0
    state_vars_offset = length(phase_record.state_variables)
    for j in 1:length(phase_record.site_fractions)
        total += phase_record.inv_phase_matrix[i,j]*phase_record.grad[state_vars_offset+j]
    end
    return -total
end

@doc raw"""
    c_iPot

# Equation

```math
c^\alpha_{i\mathrm{Pot}} = -\sum_j e^\alpha_{ij} \frac{\partial^2 G^\alpha_M}{\partial \mathrm{Pot} \partial y^\alpha_j}
```
"""
function c_iPot(phase_record, i, pot_idx)
    total = 0.0
    state_vars_offset = length(phase_record.state_variables)
    for j in 1:length(phase_record.site_fractions)
        total += phase_record.inv_phase_matrix[i,j]*phase_record.hess[pot_idx,state_vars_offset+j]
    end
    return -total
end

@doc raw"""
    c_iA

# Equation

```math
c^\alpha_{iA} = \sum_j e^\alpha_{ij} \frac{\partial M^\alpha_A}{\partial y^\alpha_j}
```
"""
function c_iA(phase_record, i, A)
    total = 0.0
    state_vars_offset = length(phase_record.state_variables)
    for j in 1:length(phase_record.site_fractions)
        total += phase_record.inv_phase_matrix[i,j]*phase_record.mass_jac[A, state_vars_offset+j]
    end
    return total
end

@doc raw"""
    get_stable_phase_row_rhs

# Examples
```
srow, srhs = get_stable_phase_row_rhs([prx], 1, [], [1, 2], [P, T], [], [], [1])
```

# Equation

```math
\sum_{B_{\mathrm{free}}} M_B^\alpha \mu_B + \sum_\mathrm{Pot} -\frac{\partial G_M^\alpha}{\partial \mathrm{Pot}} \Delta \mathrm{Pot} + \sum_\beta 0 = G_M^\alpha + \sum_{B_{\mathrm{fixed}}} -M_B^\alpha \mu_B
```
"""
function get_stable_phase_row_rhs(phase_records, phase_idx,
                                  fixed_chempot_symbols, free_chempot_idxs,
                                  fixed_pot_symbols, free_pot_idxs,
                                  fixed_phase_symbols, free_phase_idxs,
                                  )
    phase_rec = phase_records[phase_idx]
    # Construct the row in the equilibrium matrix
    soln_size = (length(free_chempot_idxs) + length(free_pot_idxs) + length(free_phase_idxs))
    row = Array{Num}(undef, soln_size)
    for B in 1:length(free_chempot_idxs)
        row[B] = phase_rec.mass[B]
    end
    # ΔPotential columns
    col_offset = length(free_chempot_idxs)
    for pot in 1:length(free_pot_idxs)
        row[col_offset+pot] = -phase_rec.grad[free_pot_idxs[pot]]
    end
    # Δℵ columns
    col_offset += length(free_pot_idxs)
    for β in 1:length(free_phase_idxs)
        row[col_offset+β] = 0.0
    end

    # Construct right-hand-side (RHS)
    rhs = 0.0
    # Gibbs energy
    rhs += phase_rec.obj
    # Fixed chemical potentials
    for B in 1:length(fixed_chempot_symbols)
        rhs -= phase_rec.mass[B] * fixed_chempot_symbols[B]
    end

    return (row, rhs)
end

@doc raw"""
    get_N_A_row_rhs

# Examples
```
using Symbolics
@variables T P N_A
@variables Y_BETA_A Y_BETA_B
G_BETA_A = 8000.0-10.0*T;
G_BETA_B = 12000.0-10.0*T;
G_BETA = Y_BETA_A*G_BETA_A + Y_BETA_B*G_BETA_B + R*T*(Y_BETA_A*log(Y_BETA_A) + Y_BETA_B*log(Y_BETA_B));
mass_BETA = [Y_BETA_A, Y_BETA_B];
state_variables = [P, T];
site_fractions = [Y_BETA_A, Y_BETA_B];
prx = PhaseRecord("BETA", G_BETA, mass_BETA, state_variables, site_fractions);

r, rrhs = get_N_A_row_rhs([prx], 1, N_A, [], [1, 2], [P, T], [], [], [1])
```

# Equation
```math
\sum_{B_{\mathrm{free}}} \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha}  c_{iB} \mu_B + \sum_\mathrm{Pot} \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha} c_{i\mathrm{Pot}} \Delta \mathrm{Pot} + \sum_\beta M_A^\beta \Delta \aleph^\beta  \\
= - \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha} c_{iG} - \sum_{B_{\mathrm{fixed}}} \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha}  c_{iB} \mu_B + \left( \tilde{N}_A - \sum_\alpha \aleph^\alpha M^\alpha_A \right)
```

"""
function get_N_A_row_rhs(phase_records, el_idx, N_A_prescribed,
                         fixed_chempot_symbols, free_chempot_idxs,
                         fixed_pot_symbols, free_pot_idxs,
                         fixed_phase_symbols, free_phase_idxs,
                         )
    # Construct the row in the equilibrium matrix
    soln_size = (length(free_chempot_idxs) + length(free_pot_idxs) + length(free_phase_idxs))
    row = Array{Num}(undef, soln_size)
    # Chemical potential columns
    for B in 1:length(free_chempot_idxs)
        total = 0.0
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.site_fractions)
                total += prx.ℵ * prx.mass_jac[el_idx,statevar_offset+i] * c_iA(prx,i,B)
            end
        end
        row[B] = total
    end
    # ΔPotential columns
    col_offset = length(free_chempot_idxs)
    for pot in 1:length(free_pot_idxs)
        total = 0.0
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.state_variables)
                total += prx.ℵ * prx.mass_jac[el_idx,statevar_offset+i] * c_iPot(prx,i,free_pot_idxs[pot])
            end
        end
        row[col_offset+pot] = total
    end
    # Δℵ columns
    col_offset += length(free_pot_idxs)
    for β in 1:length(free_phase_idxs)
        row[col_offset+β] = phase_records[β].mass[el_idx]
    end

    # construct the right-hand-side term
    rhs = 0.0
    # C_iG terms
    for α in 1:length(phase_records)
        prx = phase_records[α]
        statevar_offset = length(prx.state_variables)
        for i in 1:length(prx.state_variables)
            rhs -= prx.ℵ * prx.mass_jac[el_idx,statevar_offset+i] * c_iG(prx, i)
        end
    end
    # Fixed chemical potentials
    for B in 1:length(fixed_chempot_symbols)
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.state_variables)
                rhs -= prx.ℵ * prx.mass_jac[el_idx,statevar_offset+i] * fixed_chempot_symbols[B]
            end
        end
    end
    # Prescribed N_A
    rhs += N_A_prescribed
    # System current N_A
    for α in 1:length(phase_records)
        prx = phase_records[α]
        rhs -= prx.mass[el_idx]
    end

    return (row, rhs)
end

@doc raw"""
    get_x_A_row_rhs

See the function `get_N_A_row_rhs`. This is the corresponding mole fraction condition.

# Equation

```math
\sum_{B_{\mathrm{free}}} \sum_\alpha \sum_i \frac{\aleph^\alpha c_{iB}}{N} \left( \frac{\partial M_A^\alpha}{\partial y_i^\alpha} - x_A \sum_C \frac{\partial M_C^\alpha}{\partial y_i^\alpha} \right)  \mu_B \\
+ \sum_\mathrm{Pot} \sum_\alpha \sum_i \frac{\aleph^\alpha c_{i\mathrm{Pot}}}{N} \left( \frac{\partial M_A^\alpha}{\partial y_i^\alpha} - x_A \sum_C \frac{\partial M_C^\alpha}{\partial y_i^\alpha} \right)  \Delta \mathrm{Pot} \\
+ \sum_{\beta_\mathrm{free}} \frac{1}{N} \left( M_A^\beta - x_A \sum_C M_C^\beta \right) \Delta \aleph^\beta \\
= - \sum_\alpha \sum_i \frac{\aleph^\alpha c_{iG}}{N} \left( \frac{\partial M_A^\alpha}{\partial y_i^\alpha} - x_A \sum_C \frac{\partial M_C^\alpha}{\partial y_i^\alpha} \right) \\
- \sum_{B_{\mathrm{fixed}}} \sum_\alpha \sum_i \frac{\aleph^\alpha c_{iB}}{N} \left( \frac{\partial M_A^\alpha}{\partial y_i^\alpha} - x_A \sum_C \frac{\partial M_C^\alpha}{\partial y_i^\alpha} \right)  \mu_B \\
+ \left( \tilde{x}_A - x_A \right)
```
"""
function get_x_A_row_rhs(phase_records, el_idx, x_A_prescribed,
                         fixed_chempot_symbols, free_chempot_idxs,
                         fixed_pot_symbols, free_pot_idxs,
                         fixed_phase_symbols, free_phase_idxs,
                         )
    num_elements = length(phase_records[1].mass)  # assumes number of elements

    # Constuct N and x_A, which will be used in almost all the expresions below
    N_A = 0.0  # Intermediate, just used for computing x_A
    for α in 1:length(phase_records)
        N_A += phase_records[α].ℵ * phase_records[α].mass[el_idx]
    end
    N = 0.0
    for B in 1:num_elements
        for α in 1:length(phase_records)
            N += phase_records[α].ℵ * phase_records[α].mass[B]
        end
    end
    x_A = N_A / N

    # TODO: conditions other than chemical potential
    # Construct the row in the equilibrium matrix
    soln_size = (length(free_chempot_idxs) + length(free_pot_idxs) + length(free_phase_idxs))
    row = Array{Num}(undef, soln_size)
    # Chemical potential columns
    for B in 1:length(free_chempot_idxs)
        total = 0.0
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.site_fractions)
                factor = prx.ℵ * c_iA(prx,i,B) / N
                inner_total = prx.mass_jac[el_idx,statevar_offset+i]
                for C in 1:num_elements
                    inner_total -= x_A * prx.mass_jac[C,statevar_offset+i]
                end
                total += factor * inner_total
            end
        end
        row[B] = total
    end
    # ΔPotential columns
    col_offset = length(free_chempot_idxs)
    for pot in 1:length(free_pot_idxs)
        total = 0.0
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.state_variables)
                factor = prx.ℵ * c_iPot(prx,i,free_pot_idxs[pot]) / N
                inner_total = prx.mass_jac[el_idx,statevar_offset+i]
                for C in 1:num_elements
                    inner_total -= x_A * prx.mass_jac[C,statevar_offset+i]
                end
                total += factor * inner_total
            end
        end
        row[col_offset+pot] = total
    end
    # Δℵ columns
    col_offset += length(free_pot_idxs)
    for β in 1:length(free_phase_idxs)
        prx = phase_records[β]
        factor = 1 / N
        inner_total = prx.mass[el_idx]
        for C in 1:num_elements
            inner_total -= x_A * prx.mass[C]
        end
        row[col_offset+β] = factor * inner_total
    end

    # construct the right-hand-side term
    rhs = 0.0
    # c_iG terms
    for α in 1:length(phase_records)
        prx = phase_records[α]
        statevar_offset = length(prx.state_variables)
        for i in 1:length(prx.state_variables)
            factor = prx.ℵ * c_iG(prx, i) / N
            inner_total = prx.mass_jac[el_idx,statevar_offset+i]
            for C in 1:num_elements
                inner_total -= x_A * prx.mass_jac[C,statevar_offset+i]
            end
            rhs -= factor * inner_total
        end
    end
    # Fixed chemical potentials
    for B in 1:length(fixed_chempot_symbols)
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.state_variables)
                factor = prx.ℵ * c_iA(prx, i, B) / N
                inner_total = prx.mass_jac[el_idx,statevar_offset+i]
                for C in 1:num_elements
                    inner_total -= x_A * prx.mass_jac[C,statevar_offset+i]
                end
                rhs -= factor * inner_total
            end
        end
    end
    # Prescribed amount
    rhs += x_A_prescribed
    rhs -= x_A

    return (row, rhs)
end

@doc raw"""
    get_N_row_rhs

See the function `get_N_A_row_rhs`. This is conceptually the same as an N_A
condition with an inner loop over all the elements in each column.

# Equation
```math
\sum_{B_{\mathrm{free}}} \sum_A \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha}  c_{iB} \mu_B
+ \sum_\mathrm{Pot} \sum_A \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha} c_{i\mathrm{Pot}} \Delta \mathrm{Pot}
+ \sum_\beta \sum_A M_A^\beta \Delta \aleph^\beta  \\
= \sum_A \left(- \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha} c_{iG} - \sum_{B_{\mathrm{fixed}}} \sum_\alpha \sum_i \aleph^\alpha \frac{\partial M_A^\alpha}{\partial y_i^\alpha}  c_{iB} \mu_B + \right) + (\tilde{N} - N)
```

"""
function get_N_row_rhs(phase_records, N_prescribed,
                       fixed_chempot_symbols, free_chempot_idxs,
                       fixed_pot_symbols, free_pot_idxs,
                       fixed_phase_symbols, free_phase_idxs,
                       )
    # Construct the row in the equilibrium matrix
    soln_size = (length(free_chempot_idxs) + length(free_pot_idxs) + length(free_phase_idxs))
    num_elements = length(phase_records[1].mass)  # assume mass length is the same for every phase record
    row = Array{Num}(undef, soln_size)
    # Chemical potential columns
    for B in 1:length(free_chempot_idxs)
        total = 0.0
        for A in 1:num_elements
            for α in 1:length(phase_records)
                prx = phase_records[α]
                statevar_offset = length(prx.state_variables)
                for i in 1:length(prx.site_fractions)
                    total += prx.ℵ * prx.mass_jac[A,statevar_offset+i] * c_iA(prx,i,B)
                end
            end
        end
        row[B] = total
    end
    # ΔPotential columns
    col_offset = length(free_chempot_idxs)
    for pot in 1:length(free_pot_idxs)
        total = 0.0
        for A in 1:num_elements
            for α in 1:length(phase_records)
                prx = phase_records[α]
                statevar_offset = length(prx.state_variables)
                for i in 1:length(prx.state_variables)
                    total += prx.ℵ * prx.mass_jac[A,statevar_offset+i] * c_iPot(prx,i,free_pot_idxs[pot])
                end
            end
        end
        row[col_offset+pot] = total
    end
    # Δℵ columns
    col_offset += length(free_pot_idxs)
    for β in 1:length(free_phase_idxs)
        total = 0.0
        for A in 1:num_elements
            total += phase_records[β].mass[A]
        end
        row[col_offset+β] = total
    end

    # construct the right-hand-side term
    rhs = 0.0
    for A in 1:num_elements
        # C_iG terms
        for α in 1:length(phase_records)
            prx = phase_records[α]
            statevar_offset = length(prx.state_variables)
            for i in 1:length(prx.state_variables)
                rhs -= prx.ℵ * prx.mass_jac[A,statevar_offset+i] * c_iG(prx, i)
            end
        end
        # Fixed chemical potential terms
        for B in 1:length(fixed_chempot_symbols)
            for α in 1:length(phase_records)
                prx = phase_records[α]
                statevar_offset = length(prx.state_variables)
                for i in 1:length(prx.state_variables)
                    rhs -= prx.ℵ * prx.mass_jac[A,statevar_offset+i] * fixed_chempot_symbols[B]
                end
            end
        end
    end
    # Prescribed N
    rhs += N_prescribed
    # Current system total N
    for A in 1:num_elements
        for α in 1:length(phase_records)
            prx = phase_records[α]
            rhs -= prx.ℵ * prx.mass[A]
        end
    end

    return (row, rhs)
end

# This is the only "public" function of this file.
"""
    cond_row_rhs

This function takes a condition symbol and dispatches on the correct method that
condition that returns the elements in the row and the right hand side.

The number of columns are determined by the `fixed_free_terms` - there's one
column for each free chemical potential, free potential, and free phase amount.
"""
function cond_row_rhs(cond, elements, phases, phase_records, fixed_free_terms)
    # assumes elements, phases, phase_records are sorted
    str_cond = string(cond)
    if str_cond == "N"
        row, rhs = get_N_row_rhs(phase_records, cond, fixed_free_terms...)
    elseif startswith(str_cond, "N_")
        el = str_cond[3:end]
        el_idx = findfirst(x -> x == el, elements)
        row, rhs = get_N_A_row_rhs(phase_records, el_idx, cond, fixed_free_terms...)
    elseif startswith(str_cond, "X_")
        el = str_cond[3:end]
        el_idx = findfirst(x -> x == el, elements)
        row, rhs = get_x_A_row_rhs(phase_records, el_idx, cond, fixed_free_terms...)
    else
        throw("Condition for $str_cond is not yet implemented")
    end

    return row, rhs
end


text:
# TODO: NonSeekableDataStream
# TODO: Remove coltypes(ds::AbstractDataStream)
# TODO: Remove colnames(ds::AbstractDataStream)
# TODO: Implement MatrixStream, which reads CSV's into a Matrix of Float64's

export openstream

abstract AbstractDataStream

type SeekableDataStream <: AbstractDataStream
    io::IO
    p::ParsedCSV
    o::ParseOptions
    nrows::Int

    function SeekableDataStream(io::IO, p::ParsedCSV, o::ParseOptions, nrows::Integer)
        r = new(io, p, o, int(nrows))
        finalizer(r, r -> close(r.io))
        return r
    end
end

function openstream(pathname::String;
                    nrows::Int = 1,
                    header::Bool = true,
                    separator::Char = ',',
                    allowquotes::Bool = true,
                    quotemark::Char = '"',
                    decimal::Char = '.',
                    nastrings::Vector = ASCIIString["", "NA"],
                    truestrings::Vector = ASCIIString["T", "t", "TRUE", "true"],
                    falsestrings::Vector = ASCIIString["F", "f", "FALSE", "false"],
                    makefactors::Bool = false,
                    colnames::Vector = UTF8String[],
                    cleannames::Bool = false,
                    coltypes::Vector{Any} = Any[],
                    allowcomments::Bool = false,
                    commentmark::Char = '#',
                    ignorepadding::Bool = true,
                    skipstart::Int = 0,
                    skiprows::Vector{Int} = Int[],
                    skipblanks::Bool = true,
                    encoding::Symbol = :utf8)
    io = open(pathname, "r")
    nbytes = 2^20
    p = ParsedCSV(Array(Uint8, nbytes), Array(Int, 1),
                  Array(Int, 1), BitArray(1))
    o = ParseOptions(header, separator, allowquotes, quotemark,
                     decimal, nastrings, truestrings, falsestrings,
                     makefactors, colnames, cleannames, coltypes,
                     allowcomments, commentmark, ignorepadding,
                     skipstart, skiprows, skipblanks, encoding)
    return SeekableDataStream(io, p, o, nrows)
end

function Base.show(io::IO, ds::SeekableDataStream)
    @printf io "SeekableDataStream"
    @printf io "Minibatch Size: %d\n" ds.nrows
    return
end

function Base.start(s::SeekableDataStream)
    seek(s.io, 0)
    return DataFrame()
end

function Base.next(s::SeekableDataStream, df::DataFrame)
    if position(s.io) == 0
        df = readtable!(s.p, s.io, s.nrows, s.o)
        return df, df
    else
        bytes, fields, rows = DataFrames.readnrows!(s.p, s.io, s.nrows, s.o)
        cols = fld(fields, rows)
        DataFrames.filldf!(df, rows, cols, bytes, fields, s.p, s.o)
        return df, df
    end
end

Base.done(s::SeekableDataStream, df::DataFrame) = eof(s.io)

##############################################################################
#
# Streaming data functions
#
##############################################################################

function Base.sum(ds::AbstractDataStream, dim::Integer)
    if dim != 2
        error("Only column sums available for AbstractDataStream")
    end

    df1 = start(ds)

    df1, df1 = next(ds, df1)

    cnames = colnames(df1)
    p = length(cnames)
    sums = zeros(p)
    counts = zeros(Int, p)

    n = nrow(df1)
    for j in 1:p
        c = df1[j]
        if eltype(c) <: Real
            for i in 1:n
                if !isna(c[i])
                    sums[j] += c[i]
                    counts[j] += 1
                end
            end
        end
    end

    while !done(ds, df1)
        df1, df1 = next(ds, df1)
        n = nrow(df1)
        for j in 1:p
            c = df1[j]
            if eltype(c) <: Real
                for i in 1:n
                    if !isna(c[i])
                        sums[j] += c[i]
                        counts[j] += 1
                    end
                end
            end
        end
    end

    res = DataFrame({Float64 for i in 1:p}, cnames, 1)

    for j in 1:p
        if counts[j] != 0
            res[1, j] = sums[j]
        end
    end

    return res
end

function Base.prod(ds::AbstractDataStream, dim::Integer)
    if dim != 2
        error("Only column sums available for AbstractDataStream")
    end

    df1 = start(ds)

    df1, df1 = next(ds, df1)

    cnames = colnames(df1)
    p = length(cnames)
    prods = ones(p)
    counts = zeros(Int, p)

    n = nrow(df1)
    for j in 1:p
        c = df1[j]
        if eltype(c) <: Real
            for i in 1:n
                if !isna(c[i])
                    prods[j] *= c[i]
                    counts[j] += 1
                end
            end
        end
    end

    while !done(ds, df1)
        df1, df1 = next(ds, df1)
        n = nrow(df1)
        for j in 1:p
            c = df1[j]
            if eltype(c) <: Real
                for i in 1:n
                    if !isna(c[i])
                        prods[j] *= c[i]
                        counts[j] += 1
                    end
                end
            end
        end
    end

    res = DataFrame({Float64 for i in 1:p}, cnames, 1)

    for j in 1:p
        if counts[j] != 0
            res[1, j] = prods[j]
        end
    end

    return res
end

function Base.mean(ds::AbstractDataStream, dim::Integer)
    if dim != 2
        error("Only column sums available for AbstractDataStream")
    end

    df1 = start(ds)

    df1, df1 = next(ds, df1)

    cnames = colnames(df1)
    p = length(cnames)
    sums = zeros(p)
    counts = zeros(Int, p)

    n = nrow(df1)
    for j in 1:p
        c = df1[j]
        if eltype(c) <: Real
            for i in 1:n
                if !isna(c[i])
                    sums[j] += c[i]
                    counts[j] += 1
                end
            end
        end
    end

    while !done(ds, df1)
        df1, df1 = next(ds, df1)
        n = nrow(df1)
        for j in 1:p
            c = df1[j]
            if eltype(c) <: Real
                for i in 1:n
                    if !isna(c[i])
                        sums[j] += c[i]
                        counts[j] += 1
                    end
                end
            end
        end
    end

    res = DataFrame({Float64 for i in 1:p}, cnames, 1)

    for j in 1:p
        if counts[j] != 0
            res[1, j] = sums[j] / counts[j]
        end
    end

    return res
end

# function colvars(ds::AbstractDataStream)
#   p = length(colnames(ds))
#   means = zeros(p)
#   deltas = zeros(p)
#   m2s = zeros(p)
#   vars = zeros(p)
#   ns = zeros(Int, p)

#   for minibatch in ds
#     for row_index in 1:nrow(minibatch)
#       for column_index in 1:p
#         if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
#           ns[column_index] += 1
#           deltas[column_index] = minibatch[row_index, column_index] - means[column_index]
#           means[column_index] += deltas[column_index] / ns[column_index]
#           m2s[column_index] = m2s[column_index] + deltas[column_index] * (minibatch[row_index, column_index] - means[column_index])
#           vars[column_index] = m2s[column_index] / (ns[column_index] - 1)
#         end
#       end
#     end
#   end

#   result_types = {Float64 for i in 1:p}
#   results = DataFrame(result_types, colnames(ds), 1)

#   for column_index in 1:p
#     if ns[column_index] != 0
#       results[1, column_index] = vars[column_index]
#     end
#   end

#   return results
# end

# function colstds(ds::AbstractDataStream)
#   vars = colvars(ds)
#   stds = deepcopy(vars)
#   column_types = coltypes(vars)
#   for j in 1:length(column_types)
#     if column_types[j] <: Real
#       stds[1, j] = sqrt(vars[1, j])
#     end
#   end
#   return stds
# end

# function colmins(ds::AbstractDataStream)
#   p = length(colnames(ds))
#   mins = [Inf for i in 1:p]
#   ns = zeros(Int, p)

#   for minibatch in ds
#     for row_index in 1:nrow(minibatch)
#       for column_index in 1:p
#         if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
#           if minibatch[row_index, column_index] < mins[column_index]
#             mins[column_index] = minibatch[row_index, column_index]
#             ns[column_index] += 1
#           end
#         end
#       end
#     end
#   end

#   result_types = {Float64 for i in 1:p}
#   df = DataFrame(result_types, colnames(ds), 1)

#   for column_index in 1:p
#     if ns[column_index] != 0
#       df[1, column_index] = mins[column_index]
#     end
#   end

#   return df
# end

# function colmaxs(ds::AbstractDataStream)
#   p = length(colnames(ds))
#   maxs = [-Inf for i in 1:p]
#   ns = zeros(Int, p)

#   for minibatch in ds
#     for row_index in 1:nrow(minibatch)
#       for column_index in 1:p
#         if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
#           if minibatch[row_index, column_index] > maxs[column_index]
#             maxs[column_index] = minibatch[row_index, column_index]
#             ns[column_index] += 1
#           end
#         end
#       end
#     end
#   end

#   result_types = {Float64 for i in 1:p}
#   df = DataFrame(result_types, colnames(ds), 1)

#   for column_index in 1:p
#     if ns[column_index] != 0
#       df[1, column_index] = maxs[column_index]
#     end
#   end

#   return df
# end

# function colranges(ds::AbstractDataStream)
#   p = length(colnames(ds))
#   mins = [Inf for i in 1:p]
#   maxs = [-Inf for i in 1:p]
#   ns = zeros(Int, p)

#   for minibatch in ds
#     for row_index in 1:nrow(minibatch)
#       for column_index in 1:p
#         if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
#           ns[column_index] += 1
#           if minibatch[row_index, column_index] < mins[column_index]
#             mins[column_index] = minibatch[row_index, column_index]
#           end
#           if minibatch[row_index, column_index] > maxs[column_index]
#             maxs[column_index] = minibatch[row_index, column_index]
#           end
#         end
#       end
#     end
#   end

#   result_types = {Float64 for i in 1:p}
#   df_mins = DataFrame(result_types, colnames(ds), 1)
#   df_maxs = DataFrame(result_types, colnames(ds), 1)

#   for column_index in 1:p
#     if ns[column_index] != 0
#       df_mins[1, column_index] = mins[column_index]
#       df_maxs[1, column_index] = maxs[column_index]
#     end
#   end

#   return (df_mins, df_maxs)
# end

# # Two-pass algorithm for covariance and correlation
# function cov(ds::AbstractDataStream)
#   p = length(colnames(ds))

#   # Make one pass to compute means
#   means = colmeans(ds)

#   # Now compute covariances during second pass
#   ns = zeros(Int, p, p)
#   covariances = datazeros(p, p)

#   for minibatch in ds
#     for row_index in 1:nrow(minibatch)
#       for column_index in 1:p
#         for alt_column_index in 1:p
#           if coltypes(minibatch)[column_index] <: Real &&
#                 !isna(minibatch[row_index, column_index]) &&
#                 coltypes(minibatch)[alt_column_index] <: Real &&
#                 !isna(minibatch[row_index, alt_column_index])
#             ns[column_index, alt_column_index] += 1
#             n = ns[column_index, alt_column_index]
#             a = minibatch[row_index, column_index] - means[1, column_index]
#             b = minibatch[row_index, alt_column_index] - means[1, alt_column_index]
#             covariances[column_index, alt_column_index] = ((n - 1) / n) * covariances[column_index, alt_column_index] + (a * b) / n
#           end
#         end
#       end
#     end
#   end

#   # Scale estimates by (n / (n - 1))
#   for i in 1:p
#     for j in 1:p
#       if ns[i, j] <= 2
#         covariances[i, j] = NA
#       else
#         n = ns[i, j]
#         covariances[i, j] *= (n / (n - 1))
#       end
#     end
#   end

#   return covariances
# end

# function cor(ds::AbstractDataStream)
#   covariances = cov(ds)
#   correlations = deepcopy(covariances)
#   p = nrow(correlations)
#   for i in 1:p
#     for j in 1:p
#       correlations[i, j] = covariances[i, j] / sqrt(covariances[i, i] * covariances[j, j])
#     end
#   end
#   return correlations
# end

function Base.select(ds::AbstractDataStream, query::Integer)
    i = 0
    for df in ds
        u = nrow(df)
        if i + u > query
            return df[query - i, :]
        end
        i += u
    end
    error("Did not find requested row")
end

# # TODO: Stop returning empty DataFrame at the end of a stream
# #       (NOTE: Probably not possible because we don't know nrows.)
# # TODO: Implement
# #        * colentropys
# #        * colcardinalities
# #        * colmedians
# #        * colffts
# #        * colnorms


text:
using SWIP
using Test

@testset "SWIP.jl" begin
    # Write your own tests here.
end


text:
function mass_matrix_tg_entries!(
    mass_matrix,
    tg::TG,
    global_index::Base.ImmutableDict{Symbol, Int64},
) where {TG <: PSY.TurbineGov}
    @debug "Using default mass matrix entries $TG"
end

function mass_matrix_avr_entries!(
    mass_matrix,
    tg::PSY.HydroTurbineGov,
    global_index::Base.ImmutableDict{Symbol, Int64},
)
    mass_matrix[global_index[:x_g1], global_index[:x_g1]] = PSY.get_Tf(avr)
    mass_matrix[global_index[:x_g3], global_index[:x_g3]] = PSY.get_Tg(avr)
end

function mdl_tg_ode!(
    ::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.TGFixed, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Update inner vars
    P_ref = get_P_ref(device)
    inner_vars[τm_var] = P_ref * PSY.get_efficiency(PSY.get_prime_mover(device))

    return
end

function mdl_tg_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.TGTypeI, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain references
    ω_ref = get_ω_ref(device)
    P_ref = get_P_ref(device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, PSY.TGTypeI)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    x_g1 = internal_states[1]
    x_g2 = internal_states[2]
    x_g3 = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, PSY.TGTypeI)
    ω = @view device_states[external_ix]

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    inv_R = PSY.get_R(tg) < eps() ? 0.0 : (1.0 / PSY.get_R(tg))
    Ts = PSY.get_Ts(tg)
    Tc = PSY.get_Tc(tg)
    T3 = PSY.get_T3(tg)
    T4 = PSY.get_T4(tg)
    T5 = PSY.get_T5(tg)
    V_min, V_max = PSY.get_valve_position_limits(tg)

    #Compute auxiliary parameters
    P_in_sat = clamp(P_ref + inv_R * (ω_ref - ω[1]), V_min, V_max)

    #Compute block derivatives
    _, dxg1_dt = low_pass(P_in_sat, x_g1, 1.0, Ts)
    y_ll, dxg2_dt = lead_lag(x_g1, x_g2, 1.0, T3, Tc)
    τ_m, dxg3_dt = lead_lag(y_ll, x_g3, 1.0, T4, T5)

    #Compute 3 States TG ODE:
    output_ode[local_ix[1]] = dxg1_dt
    output_ode[local_ix[2]] = dxg2_dt
    output_ode[local_ix[3]] = dxg3_dt

    #Update mechanical torque
    inner_vars[τm_var] = τ_m

    return
end

function mdl_tg_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.TGTypeII, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain references
    ω_ref = get_ω_ref(device)
    P_ref = get_P_ref(device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, PSY.TGTypeII)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    xg = internal_states[1]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, PSY.TGTypeII)
    ω = @view device_states[external_ix]

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    inv_R = PSY.get_R(tg) < eps() ? 0.0 : (1.0 / PSY.get_R(tg))
    T1 = PSY.get_T1(tg)
    T2 = PSY.get_T2(tg)

    #Compute block derivatives
    y_ll1, dxg_dt = lead_lag(inv_R * (ω_ref - ω[1]), xg, 1.0, T1, T2)
    τ_m = y_ll1 + P_ref / 1.0

    #Compute 1 State TG ODE:
    output_ode[local_ix[1]] = dxg_dt

    #Update mechanical torque
    inner_vars[τm_var] = τ_m

    return
end

function mdl_tg_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.SteamTurbineGov1, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain TG
    tg = PSY.get_prime_mover(device)
    #Obtain references
    P_ref = get_P_ref(device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, typeof(tg))

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    x_g1 = internal_states[1]
    x_g2 = internal_states[2]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, typeof(tg))
    ω = @view device_states[external_ix]

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    inv_R = PSY.get_R(tg) < eps() ? 0.0 : (1.0 / PSY.get_R(tg))
    T1 = PSY.get_T1(tg)
    T2 = PSY.get_T2(tg)
    V_min, V_max = PSY.get_valve_position_limits(tg)
    T3 = PSY.get_T3(tg)
    D_T = PSY.get_D_T(tg)

    #Compute auxiliary parameters
    ref_in = inv_R * (P_ref - (ω[1] - 1.0))

    #Compute block derivatives
    x_g1_sat, dxg1_dt = low_pass_nonwindup(ref_in, x_g1, 1.0, T1, V_min, V_max)
    y_ll, dxg2_dt = lead_lag(x_g1_sat, x_g2, 1.0, T2, T3)
    τ_m = y_ll - D_T * (ω[1] - 1.0)

    #Compute 2 State TG ODE:
    output_ode[local_ix[1]] = dxg1_dt
    output_ode[local_ix[2]] = dxg2_dt

    #Update mechanical torque
    inner_vars[τm_var] = τ_m

    return
end

function mdl_tg_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.GasTG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain references
    P_ref = get_P_ref(device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, PSY.GasTG)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    x_g1 = internal_states[1]
    x_g2 = internal_states[2]
    x_g3 = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, PSY.GasTG)
    ω = @view device_states[external_ix]

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    inv_R = PSY.get_R(tg) < eps() ? 0.0 : (1.0 / PSY.get_R(tg))
    T1 = PSY.get_T1(tg)
    T2 = PSY.get_T2(tg)
    T3 = PSY.get_T3(tg)
    D_turb = PSY.get_D_turb(tg)
    AT = PSY.get_AT(tg)
    KT = PSY.get_Kt(tg)
    V_min, V_max = PSY.get_V_lim(tg)

    #Compute auxiliary parameters
    x_in = min((P_ref - inv_R * (ω[1] - 1.0)), (AT + KT * (AT - x_g3)))

    #Compute block derivatives
    x_g1_sat, dxg1_dt = low_pass_nonwindup(x_in, x_g1, 1.0, T1, V_min, V_max)
    _, dxg2_dt = low_pass(x_g1_sat, x_g2, 1.0, T2)
    _, dxg3_dt = low_pass(x_g2, x_g3, 1.0, T3)

    #Compute output torque
    τ_m = x_g2 - D_turb * (ω[1] - 1.0)

    #Compute 1 State TG ODE:
    output_ode[local_ix[1]] = dxg1_dt
    output_ode[local_ix[2]] = dxg2_dt
    output_ode[local_ix[3]] = dxg3_dt

    #Update mechanical torque
    inner_vars[τm_var] = τ_m

    return
end

function mdl_tg_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.HydroTurbineGov, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}

    #Obtain references
    P_ref = get_P_ref(device)
    ω_ref = get_ω_ref(device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(device, PSY.HydroTurbineGov)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    x_g1 = internal_states[1] # Filter
    x_g2 = internal_states[2] # PI State
    x_g3 = internal_states[3] # Gate
    x_g4 = internal_states[4] # Flow

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(device, PSY.HydroTurbineGov)
    ω = device_states[external_ix][1]
    Δω = ω - ω_ref

    #Get Parameters
    tg = PSY.get_prime_mover(device)
    R = PSY.get_R(tg)
    r = PSY.get_r(tg)
    Tr = PSY.get_Tr(tg)
    Tf = PSY.get_Tf(tg)
    Tg = PSY.get_Tg(tg)
    #Gate velocity limits not implemented
    #VELM = PSY.get_VELM(tg)
    G_min, G_max = PSY.get_gate_position_limits(tg)
    Tw = PSY.get_Tw(tg)
    At = PSY.get_At(tg)
    D_T = PSY.get_D_T(tg)
    q_nl = PSY.get_q_nl(tg)

    #Compute block derivatives
    c, dxg2_dt = pi_block_nonwindup(x_g1, x_g2, 1.0 / r, 1.0 / (r * Tr), G_min, G_max)
    P_in = P_ref - Δω - R * c
    _, dxg1_dt = low_pass_mass_matrix(P_in, x_g1, 1.0, Tf)
    _, dxg3_dt = low_pass_mass_matrix(c, x_g3, 1.0, Tg)
    h = (x_g4 / x_g3)^2

    #Compute 4 States TG ODE:
    output_ode[local_ix[1]] = dxg1_dt
    output_ode[local_ix[2]] = dxg2_dt
    output_ode[local_ix[3]] = dxg3_dt
    output_ode[local_ix[4]] = (1.0 - h) / Tw

    #Update mechanical torque
    inner_vars[τm_var] = (x_g4 - q_nl) * h * At - D_T * Δω * x_g3
    return
end


text:
# Code to implement common AbstractVector operations, such as concatenation and resizing

function Base.empty(t::Table)
    return Table(map(empty, columns(t)))
end

function Base.empty!(t::Table)
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(empty!, columns(t))
    return t
end

function Base.resize!(t::Table, n::Integer)
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(col -> resize!(col, n), columns(t))
    return t
end

function Base.vcat(t1::Table{names}, t2::Table{names}) where {names}
    c1 = columns(t1)
    c2 = columns(t2)
    c3 = map(vcat, c1, c2)
    return Table(c3)
end

function Base.push!(t::Table{names}, row::NamedTuple{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(push!, columns(t), row)
    return t
end

function Base.pop!(t::Table)
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return map(pop!, columns(t), row)
end

function Base.pushfirst!(t::Table{names}, row::NamedTuple{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(pushfirst!, columns(t), row)
    return t
end

function Base.popfirst!(t::Table{names})
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return map(popfirst!, columns(t), row)
end

function Base.append!(t1::Table{names}, t2::Table{names}) where {names}
    if accelerators(t1) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(append!, columns(t1), columns(t2))
    return t1
end

function Base.prepend!(t1::Table{names}, t2::Table{names}) where {names}
    if accelerators(t1) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(prepend!, columns(t1), columns(t2))
    return t1
end

function Base.deleteat!(t::Table, i)
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map(col -> deleteat!(col, i), columns(t))
    return t
end

function Base.insert!(t::Table{names}, i, row::NamedTuple{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    map((col, value) -> insert!(col, i, value), columns(t), row)
    return t
end

function Base.splice!(t::Table, index::Integer)
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return map(col -> splice!(col, index), columns(t))
end

function Base.splice!(t::Table{names}, index::Integer, row::NamedTuple{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return map((col, value) -> splice!(col, index, value), columns(t), row)
end

function Base.splice!(t::Table{names}, index::Integer, t2::Table{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return map((col, values) -> splice!(col, index, values), columns(t), columns(t2))
end

function Base.splice!(t::Table, inds::UnitRange)
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return Table(map(col -> splice!(col, inds), columns(t)))
end

function Base.splice!(t::Table{names}, inds::UnitRange, row::NamedTuple{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return Table(map((col, value) -> splice!(col, inds, value), columns(t), row))
end

function Base.splice!(t::Table{names}, inds::UnitRange, t2::Table{names}) where {names}
    if accelerators(t) !== ()
        error("Mutating tables with acceleration indexes currently unsupported")
    end
    return Table(map((col, values) -> splice!(col, inds, values), columns(t), columns(t2)))
end


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

module PostProcess

export PostProcessing, StandardPostProc, postprocess

import Compat.ASCIIString

using ..TrajSaveReplay
using ..AddSupplementary
using ..SaveHelpers
using ..Log_To_CSV
using ..Log_To_Scripted
using ..Log_To_Waypoints
using ..Summarize
using ..Label270_To_Text
import ...RLESCAS

abstract PostProcessing

type StandardPostProc <: PostProcessing
  formats::Vector{ASCIIString}
  filters::Vector{ASCIIString}
end
StandardPostProc() = StandardPostProc(ASCIIString[], ASCIIString[])

function postprocess(filename::AbstractString, opts::StandardPostProc)
  formats = opts.formats
  filters = opts.filters

  sort!(formats)
  #sorting prevents pdf from appearing after tex.  This works around tex being deleted as
  #an intermediate file during the pdf process

  #fill and add supplementary to all files
  fill_replay(filename, overwrite=true)
  add_supplementary(filename)
  #filters
  for f in filters
    if f == "nmacs_only" && !nmacs_only(filename)
      return #if it fails any of the filters, we're done
    end
  end

  for f in formats
    if f == "pdf"
      RLESCAS.include_visualize()
      RLESCAS.trajPlot(filename, format=:PDF)
    elseif f == "tex"
      RLESCAS.include_visualize()
      RLESCAS.trajPlot(filename, format=:TEX)
    elseif f == "scripted"
      log_to_scripted(filename)
    elseif f == "waypoints"
      log_to_waypoints(filename)
    elseif f == "csv"
      log_to_csv(filename)
    elseif f == "label270_text"
      label270_to_text(filename)
    elseif f == "summary"
      summarize(filename)
    else
      warn("config: unrecognized output")
    end
  end
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
#######################################################################
#
# An example showing all 48 default chart styles available in Excel 2007
# using Python and XlsxWriter. Note, these styles are not the same as
# the styles available in Excel 2013.

# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_styles.xlsx")

	# sythesis Python's title string function
	title(words) = join([uppercasefirst(t) for t in split(words, ' ')], " ")

	# Show the styles for all of these chart types.
	for chart_type in ["column" "area" "line" "pie"]
		# Add a worksheet for each chart type.
		ws = add_worksheet!(wb, title(chart_type))
		set_zoom!(ws, 30)
		style_number = 1

		# Create 48 charts, each with a different style.
		for row_num in 0:16:90, col_num in 0:9:64
			chart = add_chart!(wb, Dict("type"=> chart_type))
			add_series!(chart, Dict("values"=> "=Data!\$A\$1:\$A\$6"))
			set_title!(chart, Dict("name"=> "Style $style_number"))
			set_legend!(chart, Dict("none"=> true))
			set_style!(chart, style_number)
			insert_chart!(ws, row_num, col_num , chart)
			style_number += 1
		end
	end
	# Create a worksheet with data for the charts.
	ws = add_worksheet!(wb, "Data")
	data = [10, 40, 50, 20, 10, 50]
	write_column!(ws, "A1", data)
	hide!(ws)

	close(wb)
	isfile("chart_styles.xlsx")
end

test()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.57782063716216 <= q <= 17.16108472444359)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.6503895050792288 + 0.5000448231799557*1/(1 + exp(-(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])))) + 0.8223864208279696*1/(1 + exp(-(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])))) + -0.39198040344556784*1/(1 + exp(-(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])))) + 0.4922116220556738*1/(1 + exp(-(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))))) + 1/(1 + exp(-(0.1539410450835783 + -0.07568393712605781*1/(1 + exp(-(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])))) + -0.902375002810281*1/(1 + exp(-(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])))) + 0.7665443117707627*1/(1 + exp(-(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])))) + -0.49047645888895497*1/(1 + exp(-(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))))) + 1/(1 + exp(-(0.6062617117178619 + -0.6503507748090844*1/(1 + exp(-(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])))) + -0.9378359170676176*1/(1 + exp(-(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])))) + -0.23896524724444168*1/(1 + exp(-(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])))) + -0.3347498056547815*1/(1 + exp(-(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))))) + 1/(1 + exp(-(-0.5413645906370745 + 0.4055343512605676*1/(1 + exp(-(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])))) + -0.13388786735965885*1/(1 + exp(-(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])))) + -0.9731450409901217*1/(1 + exp(-(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])))) + 0.17849747613123057*1/(1 + exp(-(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions