text:
# Exception returned by bad tabular indexing
struct IndexError <: Exception
    str::String
end

# Dims trait for dimensionality of a Tabular
struct Dims{n}; end
@pure Dims(n::Int) = Dims{n}()
Dims(::Tabular{n}) where {n} = Dims(n)
get(::Dims{n}) where {n} = n

@pure Base.:+(::Dims{n}, m::Int) where {n} = Dims{n+m}()
@pure Base.:+(n::Int, ::Dims{m}) where {m} = Dims{n+m}()
@pure Base.:+(::Dims{n}, ::Dims{m}) where {n,m} = Dims{n+m}()

@pure Base.:-(::Dims{n}, m::Int) where {n} = Dims{n-m}()
@pure Base.:-(n::Int, ::Dims{m}) where {m} = Dims{n-m}()
@pure Base.:-(::Dims{n}, ::Dims{m}) where {n,m} = Dims{n-m}()

@pure Base.:(==)(::Dims{n}, m::Int) where {n} = n == m
@pure Base.:(==)(n::Int, ::Dims{m}) where {m} = n == m
@pure Base.:(==)(::Dims{n}, ::Dims{m}) where {n,m} = n == m

@pure Base.:(<)(::Dims{n}, m::Int) where {n} = n < m
@pure Base.:(<)(n::Int, ::Dims{m}) where {m} = n < m
@pure Base.:(<)(::Dims{n}, ::Dims{m}) where {n,m} = n < m

@pure Base.:(<=)(::Dims{n}, m::Int) where {n} = n <= m
@pure Base.:(<=)(n::Int, ::Dims{m}) where {m} = n <= m
@pure Base.:(<=)(::Dims{n}, ::Dims{m}) where {n,m} = n <= m

@pure Base.:(>=)(::Dims{n}, m::Int) where {n} = n >= m
@pure Base.:(>=)(n::Int, ::Dims{m}) where {m} = n >= m
@pure Base.:(>=)(::Dims{n}, ::Dims{m}) where {n,m} = n >= m

@pure Base.:(>)(::Dims{n}, m::Int) where {n} = n > m
@pure Base.:(>)(n::Int, ::Dims{m}) where {m} = n > m
@pure Base.:(>)(::Dims{n}, ::Dims{m}) where {n,m} = n > m

struct ScalarIndex; end
struct NonscalarIndex; end

"""
    index_shape(tabular, inds...)

Returns a tuple indicating which `inds` are a `ScalarIndex()`, and which are a
`NonScalarIndex()`. The output is useful to determine dimensionality of the output of
indexing and other operations.
"""
@inline function index_shape(x::Series, i)
    shape = _index_shape(Dims(1), typeof(get(x)), typeof(i))
    if shape === (nothing,)
        # TODO: An explicit path using dat from elements?
        throw(IndexError("Tried to index a $(typeof(x)) with a $(typeof(i))"))
    else
        return shape # fast path for fully-typed, inferred case
    end
end

@inline function index_shape(x::Table, i1, i2)
    shape = _index_shape(Dims(2), typeof(get(x)), typeof(i1), typeof(i2))
    if shape[1] === nothing || shape[2] === nothing
        # TODO: An explicit path using data from elements?
        throw(IndexError("Tried to index a $(typeof(x)) with $(typeof(i1)), $(typeof(i2))"))
    else
        return shape # fast path for fully-typed, inferred case
    end
end

@pure function Dims(x::Tuple{Vararg{Union{ScalarIndex, NonscalarIndex}}})
    Dims{count(i -> (i isa NonscalarIndex), x)}()
end

# ============================
#  index_shape implementation
# ============================
_index_shape(::Dims{N}, container_type, index_types...) where {N} = ntuple(i->nothing, Val{N})

# Series

# associatives
_index_shape(::Dims{1}, ::Type{<:Associative{K}}, ::Type{K}) where {K} = (ScalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:Associative{K}}, ::Type{<:AbstractVector{K}}) where {K} = (NonscalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:Associative}, ::Type{Colon}) = (NonscalarIndex(),)

# vectors
_index_shape(::Dims{1}, ::Type{<:AbstractVector}, ::Type{<:Integer}) = (ScalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:AbstractVector}, ::Type{<:AbstractVector{<:Integer}}) = (NonscalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:AbstractVector}, ::Type{Colon}) = (NonscalarIndex(),)

# tuples of pairs
_index_shape(::Dims{1}, ::Type{<:Tuple{Vararg{Pair}}}, ::Type) = (ScalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:Tuple{Vararg{Pair}}}, ::Type{<:Label}) = (ScalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:Tuple{Vararg{Pair}}}, ::Type{<:Tuple}) = (NonscalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:Tuple{Vararg{Pair}}}, ::Type{<:Tuple{Vararg{Label}}}) = (NonscalarIndex(),)
_index_shape(::Dims{1}, ::Type{<:Tuple{Vararg{Pair}}}, ::Type{Colon}) = (NonscalarIndex(),)

# everything else - structs
_index_shape(::Dims{1}, ::Type, ::Type{<:Label}) = (ScalarIndex(),)
_index_shape(::Dims{1}, ::Type, ::Type{<:Tuple{Vararg{Label}}}) = (NonscalarIndex(),)
_index_shape(::Dims{1}, ::Type, ::Type{Colon}) = (NonscalarIndex(),)

# Table

# associatives
function _index_shape(::Dims{2}, ::Type{<:Associative{K,V}}, ::Type{I}, ::Type{K}) where {K,V,I}
    (_index_shape(Dims(1), V, I)..., ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:Associative{K,V}}, ::Type{I}, ::Type{<:AbstractVector{K}}) where {K,V,I}
    (_index_shape(Dims(1), V, I)..., NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:Associative{K,V}}, ::Type{I}, ::Type{Colon}) where {K,V,I}
    (_index_shape(Dims(1), V, I)..., NonscalarIndex())
end

# vectors
function _index_shape(::Dims{2}, ::Type{<:AbstractVector{V}}, ::Type{I}, ::Type{<:Integer}) where {V,I}
    (_index_shape(Dims(1), V, I)..., ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractVector{V}}, ::Type{I}, ::Type{<:AbstractVector{<:Integer}}) where {V,I}
    (_index_shape(Dims(1), V, I)..., NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractVector{V}}, ::Type{I}, ::Type{Colon}) where {V,I}
    (_index_shape(Dims(1), V, I)..., NonscalarIndex())
end

# matrices
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{<:Integer}, ::Type{<:Integer})
    (ScalarIndex(), ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{<:Integer}, ::Type{<:AbstractVector{<:Integer}})
    (ScalarIndex(), NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{<:Integer}, ::Type{Colon})
    (ScalarIndex(), NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{<:AbstractVector{<:Integer}}, ::Type{<:Integer})
    (NonscalarIndex(), ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{Colon}, ::Type{<:Integer})
    (NonscalarIndex(), ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{<:AbstractVector{<:Integer}}, ::Type{<:AbstractVector{<:Integer}})
    (NonscalarIndex(), NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{<:AbstractVector{<:Integer}}, ::Type{Colon})
    (NonscalarIndex(), NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{Colon}, ::Type{<:AbstractVector{<:Integer}})
    (NonscalarIndex(), NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{<:AbstractMatrix}, ::Type{Colon}, ::Type{Colon})
    (NonscalarIndex(), NonscalarIndex())
end

# tuples of pairs
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type) where {I, T <: Tuple{Vararg{Pair}}}
    (_index_shape(Dims(1), first_val_type(T), I)..., ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{<:Label}) where {I, T <: Tuple{Vararg{Pair}}}
    (_index_shape(Dims(1), first_val_type(T), I)..., ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{<:Tuple}) where {I, T <: Tuple{Vararg{Pair}}}
    (_index_shape(Dims(1), first_val_type(T), I)..., NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{<:Tuple{Vararg{Label}}}) where {I, T <: Tuple{Vararg{Pair}}}
    (_index_shape(Dims(1), first_val_type(T), I)..., NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{Colon}) where {I, T <: Tuple{Vararg{Pair}}}
    (_index_shape(Dims(1), first_val_type(T), I)..., NonscalarIndex())
end
@pure function first_val_type(T::Type{<:Tuple{Vararg{Pair}}})
    T.parameters[1].parameters[2]
end

# everything else - structs
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{<:Label}) where {I, T}
    (_index_shape(Dims(1), first_field_type(T), I)..., ScalarIndex())
end
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{<:Tuple{Vararg{Label}}}) where {I, T}
    (_index_shape(Dims(1), first_field_type(T), I)..., NonscalarIndex())
end
function _index_shape(::Dims{2}, ::Type{T}, ::Type{I}, ::Type{Colon}) where {I, T}
    (_index_shape(Dims(1), first_field_type(T), I)..., NonscalarIndex())
end
@pure function first_field_type(T::Type)
    T.types[1]
end

# =========
#  indices
# =========
@inline indices(t::Tabular{N}) where {N} = _indices(Dims(N), get(t))

# Series

@inline _indices(::Dims{1}, d::Associative) = ((keys(d)),)
# @inline _indices(::Dims{1}, v::AbstractVector{<:Pair}) = ((_map(first, v)),)

# Flat arrays...
@inline _indices(::Dims{1}, a::AbstractVector) = indices(a)
@inline _indices(::Dims{2}, a::AbstractMatrix) = indices(a)

# tuples of pairs
@inline _indices(::Dims{1}, t::Tuple{Vararg{Pair}}) = (_map(first, t),)

# everything else - structs
@inline _indices(::Dims{1}, t) = (field_indices(typeof(t)),)
@pure field_indices(::Type{T}) where {T} = (map(name -> Label{name}(), fieldnames(T))...)

# Tables

@inline function _indices(::Dims{2}, d::Associative)
    (_indices(Dims(1), first(values(d)))..., keys(d))
end

# @inline function _indices(::Dims{2}, x::AbstractVector{<:Pair}) where {N}
#     f = first(x).second
#     (_indices(Dims(1), f)..., _map(first, x))
# end

@inline function _indices(::Dims{2}, v::AbstractVector)
    (_indices(Dims(1), first(v))..., indices(v)...)
end

@inline function _indices(::Dims{2}, t::Tuple{Vararg{Pair}})
    (_indices(Dims(1), first(t).second)..., _map(first, t))
end

@inline function _indices(::Dims{2}, t)
    (_indices(Dims(1), getfield(t,1))..., field_indices(typeof(t)))
end

# ==========
#  Indexing
# ==========

@propagate_inbounds function getindex(t::Tabular{N}, i::Vararg{Any, N}) where N
    dims = Dims(index_shape(t, i...))
    data = _getindex(Dims(N), get(t), i...)
    if dims == 0
        return data
    else
        return Tabular{get(dims)}(data)
    end
end

@propagate_inbounds function setindex!(t::Tabular{N}, val, i::Vararg{Any, N}) where N
    if index_shape(t, i...) isa Tuple{Vararg{ScalarIndex}}
        _setindex!(Dims(N), t.data, val, i...)
    else
        error("Non-scalar setindex! has not yet been implemented.")
    end
    return val
end

# Series

# associatives
@propagate_inbounds function _getindex(::Dims{1}, dict::Associative{K}, i::K) where {K}
    dict[i]
end
@propagate_inbounds function _getindex(::Dims{1}, dict::Associative{K}, inds::AbstractVector{K}) where {K}
    return Dict(map(k -> Pair(k, dict[k]), inds))
end
@propagate_inbounds function _getindex(::Dims{1}, d::Associative, ::Colon)
    return copy(d)
end

@propagate_inbounds function _setindex!(::Dims{1}, dict::Associative{K}, val, i::K) where {K}
    dict[i] = val
end

# vectors
@propagate_inbounds function _getindex(::Dims{1}, v::AbstractVector, i::Integer)
    v[i]
end
@propagate_inbounds function _getindex(::Dims{1}, v::AbstractVector{T}, i::AbstractVector{I}) where {I <: Integer,T}
    # Unlike AbstractArray, must preserve index values (not just element values)
    # TODO: some sort-based dictionary would be better for this kind of data,
    #       and might be very efficient for ranges
    d = Dict{I, T}()
    for ind ∈ i
        d[ind] = v[ind]
    end
    return d
end
@propagate_inbounds function _getindex(::Dims{1}, v::AbstractVector, ::Colon)
    v[:]
end

@propagate_inbounds function _setindex!(::Dims{1}, v::AbstractVector, val, i::Integer)
    v[i] = val
end

# tuples of pairs
@propagate_inbounds function _getindex(::Dims{1}, t::Tuple{Vararg{Pair}}, i)
    get_tuple_value(t, i)
end
@propagate_inbounds function _getindex(::Dims{1}, t::Tuple{Vararg{Pair}}, i::Label)
    get_tuple_value(t, i)
end
@propagate_inbounds function _getindex(::Dims{1}, t::Tuple{Vararg{Pair}}, inds::Tuple)
    _map(i -> get_tuple_pair(t, i), inds)
end
@propagate_inbounds function _getindex(::Dims{1}, t::Tuple{Vararg{Pair}}, inds::Tuple{Vararg{Label}})
    _map(i -> get_tuple_pair(t, i), inds)
end
@propagate_inbounds function _getindex(::Dims{1}, t::Tuple{Vararg{Pair}}, ::Colon)
    t
end

# everything else - structs
@inline function _getindex(::Dims{1}, x, ::Label{L}) where L
    getfield(x, L)
end
@inline function _getindex(::Dims{1}, x, inds::Tuple{Vararg{Label}})
    _map(l -> l => _getfield(x, l), inds)
end
@inline function _getindex(::Dims{1}, x, ::Colon)
    # Try return the same struct type when possible,
    # and make a copy since this isn't a view (x may be mutable)
    # Is this a good idea?
    # typeof(x)(_map(l -> _getfield(x, l), field_indices(typeof(x))))

    # or punt on that and just return a tuple:
    _getindex(Dims(1), x, field_indices(typeof(x)))
end
_getfield(x, ::Label{L}) where {L} = getfield(x, L)

@inline function _setindex!(::Dims{1}, x, val, ::Label{L}) where L
    setfield!(x, val, L)
end

# Table

# associatives
@propagate_inbounds function _getindex(::Dims{2}, dict::Associative{K}, i1, i2::K) where {K}
    _getindex(Dims(1), dict[i2], i1)
end
@propagate_inbounds function _getindex(::Dims{2}, dict::Associative{K}, i1, i2::AbstractVector{K}) where {K}
    @boundscheck if length(i2) == 0
        throw(IndexError("Expected at least one index"))
    end
    @inbounds first_ind = i2[1]
    first_data = _getindex(Dims(1), dict[first_ind], i1)
    V = typeof(first_data)
    out = similar(dict, Pair{K, V})
    out[first_ind] = first_data
    for j = 2:length(i2)
        @inbounds ind = i2[j]
        out[ind] = _getindex(Dims(1), dict[ind], i1)
    end
    return out
end
@propagate_inbounds function _getindex(::Dims{2}, dict::Associative{K}, i1, ::Colon) where {K}
    map(kv -> kv.first => _getindex(Dims(1), kv.second, i1), dict)
end

@propagate_inbounds function _setindex!(::Dims{2}, dict::Associative{K}, val, i1, i2::K) where {K}
    _setindex!(Dims(1), dict[i2], val, i1)
end

# vectors
@propagate_inbounds function _getindex(::Dims{2}, v::AbstractVector, i1, i2::Integer)
    _getindex(Dims(1), v[i2], i1)
end
@propagate_inbounds function _getindex(::Dims{2}, v::AbstractVector{T}, i1, i2::AbstractVector{I}) where {T, I<:Integer}
    # Unlike AbstractArray, must preserve index values (not just element values)
    # TODO: some sort-based dictionary would be better for this kind of data,
    #       and might be very efficient for ranges
    # TODO: try avoid intermediate allocation
    vals = map(ind -> _getindex(Dims(1), v[ind], i1), i2)
    Dict(map((k,v) -> k=>v, i2, vals))
end
@propagate_inbounds function _getindex(::Dims{2}, v::AbstractVector, i1, ::Colon)
    map(x -> _getindex(Dims(1), x, i1), v)
end

@propagate_inbounds function _setindex!(::Dims{2}, v::AbstractVector, val, i1, i2::Integer)
    _setindex!(Dims(1), v[i2], val, i1)
end

# matrices
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::Integer, i2::Integer)
    m[i1, i2]
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::AbstractVector{<:Integer}, i2::Integer)
    vals = map(ind -> m[ind, i2], i1)
    Dict(map((k,v) -> k=>v, i1, vals))
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::Colon, i2::Integer)
    m[i1, i2]
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::Integer, i2::AbstractVector{<:Integer})
    vals = map(ind -> m[i1, ind], i2)
    Dict(map((k,v) -> k=>v, i2, vals))
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::AbstractVector{<:Integer}, i2::AbstractVector{<:Integer})
    # Dict of Dicts - needs optimizing
    Dict([j2 => Dict([j1=>m[j1,j2] for j1 ∈ i1]) for j2 ∈ i2])
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, ::Colon, i2::AbstractVector{<:Integer})
    # Dict of vectors
    vals = map(ind -> m[:, ind], i2)
    Dict(map((k,v) -> k=>v, i2, vals))
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::Integer, i2::Colon)
    m[i1, i2]
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::AbstractVector{<:Integer}, ::Colon)
    # Vector of Dicts - needs optimizing
    [Dict([j1=>m[j1,j2] for j1 ∈ i1]) for j2 ∈ indices(m, 2)]
end
@propagate_inbounds function _getindex(::Dims{2}, m::AbstractMatrix, i1::Colon, i2::Colon)
    m[i1, i2]
end

@propagate_inbounds function _setindex!(::Dims{2}, m::AbstractMatrix, val, i1::Integer, i2::Integer)
    m[i1, i2] = val
end

# tuples of pairs
@propagate_inbounds function _getindex(::Dims{2}, t::Tuple{Vararg{Pair}}, other_inds, i)
    _getindex(Dims(1), get_tuple_value(t, i), other_inds)
end
@propagate_inbounds function _getindex(::Dims{2}, t::Tuple{Vararg{Pair}}, other_inds, i::Label)
    _getindex(Dims(1), get_tuple_value(t, i), other_inds)
end
@propagate_inbounds function _getindex(::Dims{2}, t::Tuple{Vararg{Pair}}, other_inds, this_inds::Tuple)
    _map(this_inds) do i
        kv = get_tuple_pair(t, i)
        Pair(kv.first, _getindex(Dims(1), kv.second, other_inds))
    end
end
@propagate_inbounds function _getindex(::Dims{2}, t::Tuple{Vararg{Pair}}, other_inds, this_inds::Tuple{Vararg{Label}})
    _map(this_inds) do i
        kv = get_tuple_pair(t, i)
        Pair(kv.first, _getindex(Dims(1), kv.second, other_inds))
    end
end
@propagate_inbounds function _getindex(::Dims{2}, t::Tuple{Vararg{Pair}}, other_inds, ::Colon)
    _map(kv -> Pair(kv.first, _getindex(Dims(1), kv.second, other_inds)), t)
end

@propagate_inbounds function _setindex!(::Dims{2}, t::Tuple{Vararg{Pair}}, val, other_inds, i)
    _setindex!(Dims(1), get_tuple_value(t, i), val, other_inds)
end

# everything else - structs
@inline function _getindex(::Dims{2}, x, i1, ::Label{L}) where L
    _getindex(Dims(1), getfield(x, L), i1)
end
@inline function _getindex(::Dims{2}, x, i1, inds::Tuple{Vararg{Label}})
    _map(l -> l => _getindex(Dims(1), _getfield(x, l), i1), inds)
end
@inline function _getindex(::Dims{2}, x, i1, ::Colon)
    # punt on trying to return the structs back
    _getindex(Dims(2), x, i1, field_indices(typeof(x)))
end

@inline function _setindex!(::Dims{2}, x, val, i1, ::Label{L}) where L
    _setindex!(Dims(1), getfield(x, L), val, i1)
end

# Tables are collections of rows:
@propagate_inbounds getindex(t::AbstractTable, i) = t[i, :]
@propagate_inbounds setindex!(t::AbstractTable, v, i) = (t[i, :] = v)

# """
#     Mutability(tabular)
#     Mutability(typeof(tabular))

# Returns `CanMutate()` if `setindex!` works with `tabular`, and `CannotMutate()` otherwise.
# In the latter case, `Tabulars.setindex` may be expected as a "functional" alternative to
# `setindex!`.
# """
# abstract type Mutability; end
# struct CanMutate; end
# struct CannotMutate; end

# # Maybe a trait for static indices??
# abstract type IndexType; end
# struct StaticIndices; end
# struct DynamicIndices; end

# """
#     TabularIndex

# An abstract type that represents a collection of indices of a tabular object.
# Each subtype of `Tabular` may define a corresponding `TabularIndex`, which will
# be used internally to construct new tabular objects during indexing and other
# operations.
# """
# abstract type TabularIndex; end


text:
# # ACOPF with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl) using [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports non-linear AC optimal power flow through a deep integration
# with [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl). This example shows a
# single multi-period optimization of economic dispatch with a full representation of
# AC optimal power flow.

# ## Dependencies
# We can use the same RTS data and some of the initialization as in
# [OperationsProblem example](../../notebook/3_PowerSimulations_examples/1_operations_problems.ipynb)
# by sourcing it as a dependency.
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(
    pkgpath,
    "test",
    "3_PowerSimulations_examples",
    "01_operations_problems.jl",
));

# Since we'll be doing non-linear optimization, we need a solver that supports non-linear
# problems. Ipopt is quite good.
using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

# In the [OperationsProblem example](../../notebook/3_PowerSimulations_examples/1_operations_problems.ipynb)
# we defined a unit-commitment problem with a copper plate representation of the network.
# Here, we want do define an economic dispatch (linear generation decisions) with an ACOPF
# network representation.
# So, starting with the network, we can select from _almost_ any of the endpoints on this
# tree:
#nb TypeTree(PSI.PM.AbstractPowerModel,  init_expand = 10, scopesep="\n")

# For now, let's just choose a standard ACOPF formulation.
ed_template = template_economic_dispatch(network = ACPPowerModel)

# Currently  energy budget data isn't stored in the RTS-GMLC dataset.
ed_template.devices[:Hydro] = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)

# Now we can build a 4-hour economic dispatch / ACOPF problem with the RTS data.
problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 4,
    optimizer = solver,
    balance_slack_variables = true,
)

# And solve it ...
solve!(problem)

ext: jl


text:
module RNASplicing

using BioSequences

export splice

include("stochastic.jl")
# include("armando.jl")
include("library.jl")

end # module


text: Pkg.test("FactoredEikonalFastMarching")


text:
@testset "h-for" begin
    F.def_LOCAL_VARS!()
    s = """
        @def list = [1,2,3]
        """ |> fd2html_td
    hs = raw"""
        ABC
        {{for x in list}}
          {{fill x}}
        {{end}}
        """
    tokens = F.find_tokens(hs, F.HTML_TOKENS, F.HTML_1C_TOKENS)
    hblocks, tokens = F.find_all_ocblocks(tokens, F.HTML_OCB)
    qblocks = F.qualify_html_hblocks(hblocks)
    @test qblocks[1] isa F.HFor
    @test qblocks[1].vname == "x"
    @test qblocks[1].iname == "list"
    @test qblocks[2] isa F.HFun
    @test qblocks[3] isa F.HEnd

    content, head, i = F.process_html_for(hs, qblocks, 1)
    @test isapproxstr(content, "1 2 3")
end

@testset "h-for2" begin
    F.def_LOCAL_VARS!()
    s = """
        @def list = ["path/to/badge1.png", "path/to/badge2.png"]
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for x in list}}
            {{fill x}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        path/to/badge1.png
        path/to/badge2.png
        """)
end


text:
#
# This file is a part of Code module. It contains macroses for
# organism's code.
#
# @author DeadbraiN
#
import Helper
import Creature
#
# Creates new unique variable name and returns it's symbol
# @param {Creature.Organism} org Owner of new variable
# @return {Symbol}
#
macro newVar(org)
  :(Symbol("var_", $org.symbolId += 1))
end
#
# Returns variable with specified index. e.g.: var_2
# @param index
# @return {Symbol}
#
#
macro var(index)
  :(Symbol("var_", $index))
end
#
# Creates new unique custom function name.
# @param {Creature.Organism} org Organism we are working with
# @return {String}
#
macro newFunc(org)
  :(Symbol("func_", $org.symbolId += 1))
end
#
# Returns random variable symbol
# @return {Symbol}
#
macro randVar()
  :(Symbol("var_", rand(1:Creature.VAR_AMOUNT)))
end
#
# Returns random variable symbol converted to Int
# @return {Symbol}
#
macro randIntVar()
  :(:(round(Int, $(Symbol("var_", rand(1:Creature.VAR_AMOUNT))))))
end
#
# Returns blocks array
# @param {Creature.Organism} org
# @param {Helper.CodePos} pos
# @return {Array{Organism.Creature.Block, 1}}
#
macro getBlocks(org, pos)
  :($org.funcs[$pos.fnIdx].blocks)
end
#
# Returns variables of specified block of specified
# function.
# @param {Creature.Organism} org
# @param {Helper.CodePos} pos
# @return {Array{Symbol, 1}}
#
macro getVars(org, pos)
  :($org.funcs[$pos.fnIdx].blocks[$pos.blockIdx].vars)
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
export init_global_grid

import MPI

"""
    init_global_grid(nx, ny, nz)
    me, dims, nprocs, coords, comm_cart = init_global_grid(nx, ny, nz; <keyword arguments>)

Initialize a Cartesian grid of MPI processes (and also MPI itself by default) defining implicitely a global grid.

# Arguments
- {`nx`|`ny`|`nz`}`::Integer`: the number of elements of the local grid in dimension {x|y|z}.
- {`dimx`|`dimy`|`dimz`}`::Integer=0`: the desired number of processes in dimension {x|y|z}. By default, (value `0`) the process topology is created as compact as possible with the given constraints. This is handled by the MPI implementation which is installed on your system. For more information, refer to the specifications of `MPI_Dims_create` in the corresponding documentation.
- {`periodx`|`periody`|`periodz`}`::Integer=0`: whether the grid is periodic (`1`) or not (`0`) in dimension {x|y|z}.
- `quiet::Bool=false`: whether to suppress printing information like the size of the global grid (`true`) or not (`false`).
!!! note "Advanced keyword arguments"
    - {`overlapx`|`overlapy`|`overlapz`}`::Integer=2`: the number of elements adjacent local grids overlap in dimension {x|y|z}. By default (value `2`), an array `A` of size (`nx`, `ny`, `nz`) on process 1 (`A_1`) overlaps the corresponding array `A` on process 2 (`A_2`) by `2` indices if the two processes are adjacent. E.g., if `overlapx=2` and process 2 is the right neighbor of process 1 in dimension x, then `A_1[end-1:end,:,:]` overlaps `A_2[1:2,:,:]`. That means, after every call `update_halo!(A)`, we have `all(A_1[end-1:end,:,:] .== A_2[1:2,:,:])` (`A_1[end,:,:]` is the halo of process 1 and `A_2[1,:,:]` is the halo of process 2). The analog applies for the dimensions y and z.
    - `disp::Integer=1`:  the displacement argument to `MPI.Cart_shift` in order to determine the neighbors.
    - `reorder::Integer=1`: the reorder argument to `MPI.Cart_create` in order to create the Cartesian process topology.
    - `comm::MPI.Comm=MPI.COMM_WORLD`: the input communicator argument to `MPI.Cart_create` in order to create the Cartesian process topology.
    - `init_MPI::Bool=true`: whether to initialize MPI (`true`) or not (`false`).
    For more information, refer to the documentation of MPI.jl / MPI.

# Return values
- `me`: the MPI rank of the process.
- `dims`: the number of processes in each dimension.
- `nprocs`: the number of processes.
- `coords`: the Cartesian coordinates of the process.
- `comm_cart`: the MPI communicator of the created Cartesian process topology.

# Typical use cases
    init_global_grid(nx, ny, nz)                  # Basic call (no optional in and output arguments).
    me, = init_global_grid(nx, ny, nz)            # Capture 'me' (note the ','!).
    me, dims = init_global_grid(nx, ny, nz)       # Capture 'me' and 'dims'.
    init_global_grid(nx, ny, nz; dimx=2, dimy=2)  # Fix the number of processes in the dimensions x and y of the Cartesian grid of MPI processes to 2 (the number of processes can vary only in the dimension z).
    init_global_grid(nx, ny, nz; periodz=1)       # Make the boundaries in dimension z periodic.

See also: [`finalize_global_grid`](@ref)
"""
function init_global_grid(nx::Integer, ny::Integer, nz::Integer; dimx::Integer=0, dimy::Integer=0, dimz::Integer=0, periodx::Integer=0, periody::Integer=0, periodz::Integer=0, overlapx::Integer=2, overlapy::Integer=2, overlapz::Integer=2, disp::Integer=1, reorder::Integer=1, comm::MPI.Comm=MPI.COMM_WORLD, init_MPI::Bool=true, quiet::Bool=false)
    nxyz          = [nx, ny, nz];
    dims          = [dimx, dimy, dimz];
    periods       = [periodx, periody, periodz];
    overlaps      = [overlapx, overlapy, overlapz];
    cudaaware_MPI = [false, false, false]
    if haskey(ENV, "IGG_CUDAAWARE_MPI") cudaaware_MPI .= (parse(Int64, ENV["IGG_CUDAAWARE_MPI"]) > 0); end
    if none(cudaaware_MPI)
        if haskey(ENV, "IGG_CUDAAWARE_MPI_DIMX") cudaaware_MPI[1] = (parse(Int64, ENV["IGG_CUDAAWARE_MPI_DIMX"]) > 0); end
        if haskey(ENV, "IGG_CUDAAWARE_MPI_DIMY") cudaaware_MPI[2] = (parse(Int64, ENV["IGG_CUDAAWARE_MPI_DIMY"]) > 0); end
        if haskey(ENV, "IGG_CUDAAWARE_MPI_DIMZ") cudaaware_MPI[3] = (parse(Int64, ENV["IGG_CUDAAWARE_MPI_DIMZ"]) > 0); end
    end
    if (nx==1) error("Invalid arguments: nx can never be 1.") end
    if (ny==1 && nz>1) error("Invalid arguments: ny cannot be 1 if nz is greater than 1.") end
    if (any((nxyz .== 1) .& (dims .>1 ))) error("Incoherent arguments: if nx, ny, or nz is 1, then the corresponding dimx, dimy or dimz must not be set (or set 0 or 1)."); end
    if (any((nxyz .< 2 .* overlaps .- 1) .& (periods .> 0))) error("Incoherent arguments: if nx, ny, or nz is smaller than 2*overlapx-1, 2*overlapy-1 or 2*overlapz-1, respectively, then the corresponding periodx, periody or periodz must not be set (or set 0)."); end
    dims[(nxyz.==1).&(dims.==0)] .= 1;   # Setting any of nxyz to 1, means that the corresponding dimension must also be 1 in the global grid. Thus, the corresponding dims entry must be 1.
    if (init_MPI)  # NOTE: init MPI only, once the input arguments have been checked.
        if (MPI.Initialized()) error("MPI is already initialized. Set the argument 'init_MPI=false'."); end
        MPI.Init();
    else
        if (!MPI.Initialized()) error("MPI has not been initialized beforehand. Remove the argument 'init_MPI=false'."); end  # Ensure that MPI is always initialized after init_global_grid().
    end
    nprocs    = MPI.Comm_size(comm);
    MPI.Dims_create!(nprocs, dims);
    comm_cart = MPI.Cart_create(comm, dims, periods, reorder);
    me        = MPI.Comm_rank(comm_cart);
    coords    = MPI.Cart_coords(comm_cart);
    neighbors = fill(MPI.MPI_PROC_NULL, NNEIGHBORS_PER_DIM, NDIMS_MPI);
    for i = 1:NDIMS_MPI
        neighbors[:,i] .= MPI.Cart_shift(comm_cart, i-1, disp);
    end
    nxyz_g = dims.*(nxyz.-overlaps) .+ overlaps.*(periods.==0); # E.g. for dimension x with ol=2 and periodx=0: dimx*(nx-2)+2
    set_global_grid(GlobalGrid(nxyz_g, nxyz, dims, overlaps, nprocs, me, coords, neighbors, periods, disp, reorder, comm_cart, cudaaware_MPI, quiet));
    if (!quiet && me==0) println("Global grid: $(nxyz_g[1])x$(nxyz_g[2])x$(nxyz_g[3]) (nprocs: $nprocs, dims: $(dims[1])x$(dims[2])x$(dims[3]))"); end
    init_timing_functions();
    return me, dims, nprocs, coords, comm_cart; # The typical use case requires only these variables; the remaining can be obtained calling get_global_grid() if needed.
end

# Make sure that timing functions which must be fast at the first user call are already compiled now.
function init_timing_functions()
    tic();
    toc();
end


text:
using SubmodularMaximization
using PyPlot

close()

num_robots = 4
steps = 10
horizon = SubmodularMaximization.default_horizon

configs = MultiRobotTargetTrackingConfigs(num_robots,
                                          horizon=horizon)

data = target_tracking_experiment(steps=steps, num_robots=num_robots,
                                  configs=configs)

target_states = cat(map(x->x.target_states, data)..., dims=2)
robot_states = cat(map(x->x.robot_states, data)..., dims=2)

plot_state_space(configs.grid)
xlim([0, configs.grid.width+1])
ylim([0, configs.grid.height+1])

for ii = 1:size(robot_states, 1)
  let robot_states = robot_states[ii,:]
    if !continuous(configs.grid, robot_states)
      @show robot_states
      error("Robot states are not continuous")
    end

    plot_trajectory(robot_states, color=:blue)
  end
end

for ii = 1:size(target_states, 1)
  let target_states = target_states[ii,:]
    if !continuous(configs.grid, target_states)
      @show target_states
      error("Target states are not continuous")
    end

    plot_trajectory(target_states)
  end
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
import .AST: ASTMDP, ASTMetrics

"""
Stacked figure with metrics over episodes:
    - Running miss distance mean
    - Minimum miss distance
    - Cumulative number of failure events
"""
episodic_figures(mdp::ASTMDP; kwargs...) = episodic_figures(mdp.metrics; kwargs...)
function episodic_figures(metrics::ASTMetrics; gui::Bool=true, fillstd::Bool=false)
    miss_distances = metrics.miss_distance
    max_iters = length(miss_distances)

    PyPlot.pygui(gui) # Plot with GUI window (if true)
    fig = figure(figsize=(7,7))

    handles = []

    # Font size changes
    plt.rc("axes", titlesize=15, labelsize=13)
    plt.rc("legend", fontsize=12)

    ## Plot 1: Running miss distance mean
    ax = fig.add_subplot(3,1,1)
    title("Running Miss Distance Mean")

    ax.axhline(y=0, color="black", linestyle="--", linewidth=1.0)

    rolling_mean = []
    d_sum = 0
    for i in 1:max_iters
        d_sum += miss_distances[i]
        push!(rolling_mean, d_sum/i)
    end
    # [mean(miss_distances[1:i]) for i in 1:max_iters]
    ax.plot(rolling_mean, color="darkcyan", zorder=2)
    if fillstd
        miss_std_below = [mean(miss_distances[1:i])-std(miss_distances[1:i]) for i in 1:max_iters]
        miss_std_above = [mean(miss_distances[1:i])+std(miss_distances[1:i]) for i in 1:max_iters]
        ax.fill_between(1:max_iters, miss_std_below, miss_std_above, color="darkcyan", alpha=0.1)
    end

    ylabel("Miss Distance")
    ax.tick_params(labelbottom=false)

    ## Plot 2: Minimum miss distance
    ax = fig.add_subplot(3,1,2)
    title("Minimum Miss Distance")
    pl0 = ax.axhline(y=0, color="black", linestyle="--", linewidth=1.0)
    rolling_min = []
    current_min = Inf
    for i in 1:max_iters
        if miss_distances[i] < current_min
            current_min = miss_distances[i]
        end
        push!(rolling_min, current_min)
    end
    pl1 = ax.plot(rolling_min, color="darkcyan", label="AST")
    ylabel("Miss Distance")
    handles = [pl0, pl1[1]]

    ax.tick_params(labelbottom=false)

    ## Plot 3: Cumulative failures
    ax = fig.add_subplot(3,1,3)
    E = metrics.event
    max_iters = length(E)

    title("Cumulative Number of Failure Events")
    ax.plot(cumsum(E[1:max_iters]), color="darkcyan")
    xlabel("Episode")
    ylabel("Number of Events")

    yscale("log")

    fig.legend(handles, ["Event Horizon", "AST"],
               columnspacing=0.8, loc="lower center", bbox_to_anchor=(0.52, 0), fancybox=true, shadow=false, ncol=5)
    plt.tight_layout()
    fig.subplots_adjust(bottom=0.13) # <-- Change the 0.02 to work for your plot.
end

"""
Stacked figure with distributions:
    - Miss distance distribution
    - Log-likelihood distribution
"""
distribution_figures(mdp::ASTMDP; kwargs...) = distribution_figures(mdp.metrics; kwargs...)
function distribution_figures(metrics; gui=true)
    PyPlot.pygui(gui) # Plot with GUI window (if true)

    fig = figure(figsize=(6,5))

    # Font size changes
    plt.rc("axes", titlesize=12, labelsize=11)
    plt.rc("legend", fontsize=11)

    ## Plot 1: Miss distanace distribution
    subplot(2,1,1)

    Seaborn.kdeplot(-metrics.miss_distance, bw=0.1, color="darkcyan", shade=true, cut=200) # cut=2000
    plt.axvline(x=0, color="black", linestyle="--", linewidth=1.0)
    xl = plt.xlim()
    plt.xlim([xl[1], xl[2]])

    legend(["Event Horizon", "AST"], ncol=2, columnspacing=0.5, numpoints=2)
    title("Miss Distance Distribution")
    xlabel(L"-d")
    ylabel("Density");
    yl=ylim()
    ylim([yl[1], yl[2]+0.001]) # y-buffer for legend fit and placement.

    ## Plot 2: Log-likelihood distribution
    subplot(2,1,2)
    Seaborn.kdeplot(metrics.logprob[findall(metrics.event)], bw=0.1, color="darkcyan", shade=true, cut=100)
    legend(["AST"], loc="upper left")
    title("Log-Likelihood Distribution: Failure Events")
    xlabel(L"\log p")
    ylabel("Density");
    fig.tight_layout()
end


text:
# # Exact Optimization with Rational Arithmetic

# This example can be found in section 4.3 [in the paper](https://arxiv.org/pdf/2104.06675.pdf).
# The package allows for exact optimization with rational arithmetic. For this, it suffices to set up the LMO
# to be rational and choose an appropriate step-size rule as detailed below. For the LMOs included in the
# package, this simply means initializing the radius with a rational-compatible element type, e.g., `1`, rather
# than a floating-point number, e.g., `1.0`. Given that numerators and denominators can become quite large in
# rational arithmetic, it is strongly advised to base the used rationals on extended-precision integer types such
# as `BigInt`, i.e., we use `Rational{BigInt}`.

# The second requirement ensuring that the computation runs in rational arithmetic is
# a rational-compatible step-size rule. The most basic step-size rule compatible with rational optimization is
# the agnostic step-size rule with ``\gamma_t = 2/(2 + t)``. With this step-size rule, the gradient does not even need to
# be rational as long as the atom computed by the LMO is of a rational type. Assuming these requirements are
# met, all iterates and the computed solution will then be rational.

using FrankWolfe
using LinearAlgebra

n = 100
k = n

x = fill(big(1)//100, n)

f(x) = dot(x, x)
function grad!(storage, x)
    @. storage = 2 * x
end

# pick feasible region
# radius needs to be integer or rational
lmo = FrankWolfe.ProbabilitySimplexOracle{Rational{BigInt}}(1)

# compute some initial vertex
x0 = FrankWolfe.compute_extreme_point(lmo, zeros(n));

x, v, primal, dual_gap, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Agnostic(),
    print_iter=k / 10,
    verbose=true,
    emphasis=FrankWolfe.blas,
);

println("\nOutput type of solution: ", eltype(x))

# Another possible step-size rule is `rationalshortstep` which computes the step size by minimizing the
# smoothness inequality as ``\gamma_t=\frac{\langle \nabla f(x_t),x_t-v_t\rangle}{2L||x_t-v_t||^2}``. However, as this step size depends on an upper bound on the
# Lipschitz constant ``L`` as well as the inner product with the gradient ``\nabla f(x_t)``, both have to be of a rational type.

@time x, v, primal, dual_gap, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.RationalShortstep(),
    L=2,
    print_iter=k / 10,
    verbose=true,
    emphasis=FrankWolfe.blas,
);

# Note: at the last step, we exactly close the gap, finding the solution 1//n * ones(n)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct MarketBidCost <: OperationalCost
        variable::Union{Nothing, IS.TimeSeriesKey}
        no_load::Float64
        start_up::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}
        shut_down::Float64
        ancillary_services::Vector{Service}
    end

Data Structure Operational Cost to reflect market bids of energy and ancilliary services.
Compatible with most US Market bidding mechanisms

# Arguments
- `variable::Union{Nothing, IS.TimeSeriesKey}`: Variable Cost TimeSeriesKey
- `no_load::Float64`: no load cost
- `start_up::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}`: start-up cost at different stages of the thermal cycle. Warm is also refered as intermediate in some markets
- `shut_down::Float64`: shut-down cost, validation range: `(0, nothing)`, action if invalid: `warn`
- `ancillary_services::Vector{Service}`: Bids for the ancillary services
"""
mutable struct MarketBidCost <: OperationalCost
    "Variable Cost TimeSeriesKey"
    variable::Union{Nothing, IS.TimeSeriesKey}
    "no load cost"
    no_load::Float64
    "start-up cost at different stages of the thermal cycle. Warm is also refered as intermediate in some markets"
    start_up::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}
    "shut-down cost"
    shut_down::Float64
    "Bids for the ancillary services"
    ancillary_services::Vector{Service}
end

function MarketBidCost(; variable, no_load, start_up, shut_down, ancillary_services=Vector{Service}(), )
    MarketBidCost(variable, no_load, start_up, shut_down, ancillary_services, )
end

# Constructor for demo purposes; non-functional.
function MarketBidCost(::Nothing)
    MarketBidCost(;
        variable=nothing,
        no_load=0.0,
        start_up=(hot = START_COST, warm = START_COST, cold = START_COST),
        shut_down=0.0,
        ancillary_services=Vector{Service}(),
    )
end

"""Get [`MarketBidCost`](@ref) `variable`."""
get_variable(value::MarketBidCost) = value.variable
"""Get [`MarketBidCost`](@ref) `no_load`."""
get_no_load(value::MarketBidCost) = value.no_load
"""Get [`MarketBidCost`](@ref) `start_up`."""
get_start_up(value::MarketBidCost) = value.start_up
"""Get [`MarketBidCost`](@ref) `shut_down`."""
get_shut_down(value::MarketBidCost) = value.shut_down
"""Get [`MarketBidCost`](@ref) `ancillary_services`."""
get_ancillary_services(value::MarketBidCost) = value.ancillary_services

"""Set [`MarketBidCost`](@ref) `variable`."""
set_variable!(value::MarketBidCost, val) = value.variable = val
"""Set [`MarketBidCost`](@ref) `no_load`."""
set_no_load!(value::MarketBidCost, val) = value.no_load = val
"""Set [`MarketBidCost`](@ref) `start_up`."""
set_start_up!(value::MarketBidCost, val) = value.start_up = val
"""Set [`MarketBidCost`](@ref) `shut_down`."""
set_shut_down!(value::MarketBidCost, val) = value.shut_down = val
"""Set [`MarketBidCost`](@ref) `ancillary_services`."""
set_ancillary_services!(value::MarketBidCost, val) = value.ancillary_services = val


text:
using FSimROS
using Test

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


text: @deprecate turnprogress(progress::Bool) setprogress!(progress)


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV

#Set trainingSize
#trainingSizeInput = parse(Int64, ARGS[1])
trainingSize = 240

#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/ML/Lasso_Test"
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/ML"

@everywhere include("ParallelModelGeneration.jl")
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Set parameters
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

modelMatrix = zeros(industriesTotal, possibilities)
testModel = [1 0 0 1 1]
for i=1:industriesTotal
    modelMatrix[i, :] = testModel
end

#Initialise and load data
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)

riskAversions = linspace(0, 4, 16)
XArrays, YArrays = generateXandYs(industries, modelMatrix)
gamma = 0.1 #regularization term in LASSO
nRows = size(XArrays[1])[1]
w1N = repeat([0.1], outer = 11) #1/N weights
return1NMatrix = zeros(nRows-trainingSize)
returnPPDMatrix = zeros(nRows-trainingSize)

weightsPPD = zeros(nRows-trainingSize, 11)
forecastRows = zeros(nRows-trainingSize, 11)

expectedReturnMatrix = zeros(nRows-trainingSize, 11)
forecastErrors = zeros(nRows-trainingSize, 11)

rfRates = loadRiskFreeRate("NoDur", path)
rfRates = rfRates[:,1]

#Generate perfect results with risk free rate
returnPerfectMatrix = zeros(nRows-trainingSize)
weightsPerfect = zeros(nRows-trainingSize, 11)
fileName = "/zhome/9f/d/88706/SpecialeCode/Results/MV/PointPrediction/"
for g = 1:10
    fileName = "PerfectRFR"
    gammaRisk = riskAversions[g]
    for t=1:840
        println("time $t / 840, gammaRisk $g / 10 ")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)
        valY = zeros(11)
        for i = 1:10
            valY[i] = validationY[i][1]
        end
        valY[11] = rfRates[t+trainingSize]
        rfRatesVec = rfRates[t:(t+trainingSize-1)]
        trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
        weightsPerfect[t, :], returnPerfectMatrix[t] = findPerfectRFRResults(trainX, valY, valY, gammaRisk)
    end
    fileName = fileName*"_train"*string(trainingSize)*"_"*string(gammaRisk)
    writedlm(fileName*"Weights.csv", weightsPerfect,",")
    writedlm(fileName*"Returns.csv", returnPerfectMatrix,",")
end

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/IndexData/LassoTest/"

for g = 1:16
    fileName = "/zhome/9f/d/88706/SpecialeCode/Results/MV/PointPrediction/"
    gammaRisk = riskAversions[g] #riskAversion in MV optimization
    total = nRows-trainingSize-1
    for t = 1:840
        println("time $t / $total, gammaRisk $g / 10 ")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)
        expectedReturns = zeros(industriesTotal+1)
        #Load best dataset based on R^2 value
        for i = 1:industriesTotal
            summary = CSV.read(path*industries[i]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)

            maxR2, fileNameIndex = findmax(summary[:,3]) # Find file name for  max R^2

            bestFileName = summary[fileNameIndex,1]
            fileIndex = summary[fileNameIndex,7] # Find gamma

            estimate = CSV.read(path*industries[i]*"/"*string(trainingSize)*"-1/"*"240_"*bestFileName*"_predicted.CSV",nullable=false)
            expectedReturns[i] = estimate[:,3+fileIndex][t]
        end

        expectedReturns[11] = rfRates[t+trainingSize]

        expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100
        expectedReturnMatrix[t, 11]   = exp(expectedReturns[11])-1

        rfRatesVec = rfRates[t:(t+trainingSize-1)]
        trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
        Sigma =  cov(trainX)

        #A=U^(T)U where U is upper triangular with real positive diagonal entries
        F = lufact(Sigma)
        U = F[:U]  #Cholesky factorization of Sigma

        #getting the actual Y values for each industry
        valY = zeros(11)
        for i = 1:10
            valY[i] = validationY[i][1]
        end
        valY[11] = rfRates[t+trainingSize]
        return1N, returnPPD, wPPD, forecastRow = performMVOptimizationRISK(expectedReturns, U, gammaRisk, valY, valY)
        weightsPPD[t, 1:11]    = wPPD
        forecastErrors[t, 1:11]  = forecastRow-expectedReturnMatrix[t,:]
        return1NMatrix[t]      = return1N
        returnPPDMatrix[t]     = returnPPD
    end
    fileName = fileName*"_train"*string(trainingSize)*"_"*string(gammaRisk)
    writedlm(fileName*"ppdRFRWeights.csv", weightsPPD,",")
    writedlm(fileName*"ppdRFRReturns.csv", returnPPDMatrix,",")
    writedlm(fileName*"ppdRFRErrors.csv", forecastErrors,",")
    writedlm(fileName*"ppdRFRForecasts.csv", expectedReturnMatrix,",")
end

println("Finished")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/LassoMVRate.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/LassoMVRate.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## low-level cre2 interface to RE2 ##

module CRE2
using Libdl

include("cre2_h.jl")

const RE2LIB = "libcre2"

function __init__()
    push!(DL_LOAD_PATH, "/usr/local/lib")
end

function new_options(options::UInt32, mem::Integer=DEFAULT_MAX_MEM)
    copt = ccall((:cre2_opt_new, RE2LIB), Ptr{Cvoid}, ())
    set_options(copt, options, mem)
    copt
end

function free_options(copt::Ptr{Cvoid})
    if copt != Ptr{Cvoid}(0)
        ccall((:cre2_opt_delete, RE2LIB), Cvoid, (Ptr{Cvoid},), copt)
    end
end

function set_options(copt::Ptr{Cvoid}, options::UInt32, mem)
    vencoding = options & LATIN1 != 0 ? Latin1 : UTF8
    vlongest = options & LONGEST_MATCH == 0 ? 0 : 1
    vlogerrors = options & LOG_ERRORS == 0 ? 0 : 1
    vliteral = options & LITERAL == 0 ? 0 : 1
    vnevernl = options & NEVER_NL == 0 ? 0 : 1
    vdotnl = options & DOTALL == 0 ? 0 : 1
    vnevercap = options & NEVER_CAPTURE == 0 ? 0 : 1
    vcase = options & CASELESS == 0 ? 1 : 0
    vposix = options & POSIX_SYNTAX == 0 ? 0 : 1
    vperlclass = options & PERL_CLASSES != 0 ? 1 : 0
    vwordbound = options & WORD_BOUNDARY != 0 ? 1 : 0
    voneline = options & MULTILINE != 0 ? 0 : 1

    ccall((:cre2_opt_set_encoding, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vencoding)
    ccall((:cre2_opt_set_longest_match, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vlongest)
    ccall((:cre2_opt_set_log_errors, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vlogerrors)
    ccall((:cre2_opt_set_literal, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vliteral)
    ccall((:cre2_opt_set_never_nl, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vnevernl)
    ccall((:cre2_opt_set_dot_nl, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vdotnl)
    ccall((:cre2_opt_set_never_capture, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vnevercap)
    ccall((:cre2_opt_set_case_sensitive, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vcase)
    ccall((:cre2_opt_set_posix_syntax, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vposix)
    ccall((:cre2_opt_perl_classes, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vperlclass)
    ccall((:cre2_opt_word_boundary, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, vwordbound)
    ccall((:cre2_opt_one_line, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, voneline)

    ccall((:cre2_opt_set_max_mem, RE2LIB), Cvoid, (Ptr{Nothing}, Cint), copt, mem)
    copt
end

function get_options(copt::Ptr{Cvoid})
    vencoding = ccall((:cre2_opt_encoding, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vlongest = ccall((:cre2_opt_longest_match, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vlogerrors = ccall((:cre2_opt_log_errors, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vliteral = ccall((:cre2_opt_literal, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vnevernl = ccall((:cre2_opt_never_nl, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vdotnl = ccall((:cre2_opt_dot_nl, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vnevercap = ccall((:cre2_opt_never_capture, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vcase = ccall((:cre2_opt_case_sensitive, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vposix = ccall((:cre2_opt_posix_syntax, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vperlclass = ccall((:cre2_opt_perl_classes, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    vwordbound = ccall((:cre2_opt_word_boundary, RE2LIB), Cuint, (Ptr{Nothing},), copt)
    voneline = ccall((:cre2_opt_one_line, RE2LIB), Cuint, (Ptr{Nothing},), copt)

    op = UInt32(0)
    vencoding == Latin1 && (op |= LATIN1)
    vlongest != 0 && (op |= LONGEST_MATCH)
    vlogerrors != 0 && (op |= LOG_ERRORS)
    vliteral != 0 && (op |= LITERAL)
    vnevernl != 0 && (op |= NEVER_NL)
    vdotnl != 0 && (op |= DOTALL)
    vnevercap != 0 && (op |= NEVER_CAPTURE)
    vcase == 0 && (op |= CASELESS)
    vposix != 0 && (op |= POSIX_SYNTAX)
    vperlclass != 0 && (op |= PERL_CLASSES)
    vwordbound != 0 && (op |= WORD_BOUNDARY)
    voneline == 0 && (op |= MULTILINE)
    op
end

function get_maxmem(copt::Ptr{Cvoid})
    ccall((:cre2_opt_max_mem, RE2LIB), Cuint, (Ptr{Cvoid},), copt)
end

function get_program_size(regex::Ptr{Cvoid})
    ccall((:cre2_program_size, RE2LIB), Cuint, (Ptr{Cvoid},), regex)
end

# supported options for different use cases

function error_code(regex::Ptr{Cvoid})
    regex == C_NULL && error("NULL regex object")
    res = ccall((:cre2_error_code, RE2LIB), Cint, (Ptr{Cvoid},), regex)
    Int(res)
end

function error_string(regex::Ptr{Cvoid})
    regex == C_NULL && error("NULL regex object")
    res = ccall((:cre2_error_string, RE2LIB), Cstring, (Ptr{Cvoid},), regex)
    unsafe_string(res)
end

function get_ovec(match_data)
    ptr = ccall((:pcre2_get_ovector_pointer_8, RE2LIB), Ptr{Csize_t},
                (Ptr{Cvoid},), match_data)
    n = ccall((:pcre2_get_ovector_count_8, RE2LIB), UInt32,
              (Ptr{Cvoid},), match_data)
    unsafe_wrap(Array, ptr, 2n, own = false)
end

function compile(pattern::AbstractString, options::Ptr{Cvoid})
    errno = Ref{Cint}(0)
    erroff = Ref{Csize_t}(0)
    re_ptr = ccall((:cre2_new, RE2LIB), Ptr{Cvoid},
                   (Ptr{UInt8}, Csize_t, Ptr{Cvoid}),
                   pattern, sizeof(pattern), options)
    if re_ptr == C_NULL
        error("RE2 compilation error: memory allocation failed")
    end
    errno = error_code(re_ptr)
    if errno != NO_ERROR
        es = error_string(re_ptr)
        error("RE2 compilation error $errno: $es")
    end
    re_ptr
end

free_re(re) = ccall((:cre2_delete, RE2LIB), Cvoid, (Ptr{Cvoid},), re)

function exec(re, text, offset, options, match_data, nmatch_data::Integer=-1)
    text_length::Cint = sizeof(text)
    text = pointer(text)
    start_pos::Cint = offset
    end_pos::Cint = min(text_length, typemax(Cint))
    nmatch_data::Cint = nmatch_data < 0 ? length(match_data) ÷ 2 : Cint(nmatch_data)
    match_data = pointer(match_data)

    rc = ccall((:cre2_match, RE2LIB), Cint,
            (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Cint, Cuint, Ptr{Cvoid}, Cint),
            re, text, text_length, start_pos, end_pos, options, match_data, nmatch_data)

    rc != 0
end

function num_capturing_groups(re)
    n = ccall((:cre2_num_capturing_groups, RE2LIB), Cint, (Ptr{Cvoid},), re)
    n < 0 && error("no capturing groups - invalid regex")
    n
end

function create_match_data(re)
    n = num_capturing_groups(re)
    fill(C_NULL, n * 2 + 2)
end

function group_number_from_name(re, name)
  ccall((:cre2_find_named_capturing_groups, RE2LIB), Cint, (Ptr{Cvoid}, Cstring), re, name)
end

function substring_length_bynumber(match_data, number)
    n = length(match_data) ÷ 2 - 1
    0 <= number <= n || error("no data for substring $number")
    len = match_data[number*2+2]
    convert(Int, len)
end

function substring_copy_bynumber(match_data, number, buf, buf_size)
    n = length(match_data) ÷ 2 - 1
    0 <= number <= n || error("no data for substring $number")
    p = Ptr{UInt8}(match_data[number*2+1])
    unsafe_copyto!(buf, p, buf_size)
    convert(Int, buf_size-1)
end

function capture_names(re::Ptr{Cvoid}, pat::AbstractString)
    names = Dict{Int,SubString}()
    i = start(pat)
    n = 0
    while !done(pat, i)
        c, i = next(pat, i)
        if c == '<'
            n = i
        elseif c == '>' && n != 0
            name = SubString(pat, n, prevind(pat, i, 2))
            ix = group_number_from_name(re, name)
            if ix >= 0
                names[ix] = name
            end
            n = 0
        end
    end
    names
end

end # module


text:
using DataFrames

# when dataframes get written, by default they output
# comma delimited with a header, and strings are quoted.
# these things need to be rectified.
"""
    write_bedfile

    Writes a dataframe into a format that is acceptable as a bed file, or can at least be used by bedtools
      i)   Sorts the dataframe - needed for things like tabix indexing
      ii)  Writes the dataframe to disk - bedfiles don't have headers, and are tab delimited
      iii) removes the quotes around strings.
"""
function write_bedfile(output_path,df::DataFrame;sort=true,seq_id_colname=:seq_id,start_pos_colname=:start_pos)
    if sort
        sort!(df,cols=[seq_id_colname,start_pos_colname])
    end
    writetable( output_path, df, separator='\t', header=false)
    cmd = `sed -i -s 's/\"//g' $output_path` #"
    readstring(cmd)
end

# cross over two dataframes using basic bedtools intersect
function bedtools_intersect(A::DataFrame, B::DataFrame; negate=false,sort=false)

    tmpdir = tempdir()
    tmpa=joinpath( tmpdir, "A")
    tmpb=joinpath( tmpdir, "B")
    tmpc=joinpath( tmpdir, "C")
    write_bedfile(tmpa,A,sort=sort)
    write_bedfile(tmpb,B,sort=sort)
    bedtools_cmd="bedtools"
    bedtools_args = ["intersect","-wo","-a",tmpa, "-b",tmpb]
    if negate
     # splices in between 2nd and 1st element with out removing any elements
     splice!(bedtools_args,2:1,["-v"])
    end
    run_cmd(bedtools_cmd, bedtools_args, tmpc)
    df=readtable(tmpc,separator='\t',header=false)
end

function bedtools_intersect(pathA::AbstractString, B::DataFrame;negate=false,sort=false)

    tmpdir = tempdir()
    tmpb=joinpath( tmpdir, "B")
    tmpc=joinpath( tmpdir, "C")
    write_bedfile(tmpb,B,sort=sort)

    bedtools_cmd="bedtools"
    bedtools_args = ["intersect","-wo","-a",pathA, "-b",tmpb]
    if negate
       # splices in between 2nd and 1st element with out removing any elements
       splice!(bedtools_args,2:1,["-v"])
    end
    run_cmd( bedtools_cmd, bedtools_args,tmpc)
    #cmd=`bedtools intersect $minusv -wo -a $pathA -b $tmpb` |> tmpc
    #cmd=`$bedtools_cmd $bedtools_args` |> tmpc
    #readall(cmd)
    df=readtable(tmpc,separator='\t',header=false)
end

function bedtools_intersect(pathA::AbstractString, pathB::AbstractString;negate=false)
    tmpdir = tempdir()
    tmpc=joinpath( tmpdir, "C")
    bedtools_cmd="bedtools"
    bedtools_args = ["intersect","-wo","-a",pathA, "-b",pathB]
    if negate
       # splices in between 2nd and 1st element with out removing any elements
       splice!(bedtools_args,2:1,["-v"])
    end
    #cmd=`bedtools intersect $minusv -wo -a $pathA -b $tmpb` |> tmpc
    #cmd=`$bedtools_cmd $bedtools_args` |> tmpc
    #readall(cmd)
    run_cmd( bedtools_cmd, bedtools_args,tmpc)
    df=readtable(tmpc,separator='\t',header=false)
end

# converts files which have
# chr1 start strand ... to
# chr1 start start strand
function convert_methpipebed_to_bedfile(pathA)
     tmppath = tempname()
     cmd = `zcat $pathA` |> `perl -F"\t" -lane 'print join("\t",$F[0],$F[1],$F[1],$F[2],$F[3],$F[4],$F[5])'` |> tmppath
     readall(cmd)
     return tmppath
end

# -- NOT TESTED --
# -- functions that avoid creating temporary tables
function query_features_to_query_str(df::DataFrame)
   n,r = size(df)
   rows = ASCIIString[]
   for i in 1:n
        row=Any[]
        println("row $i")
        for j in 1:r
            push!(row,  df[i,j] )
        end
        push!(rows, join(row,"\t")  )
    end
   return join(rows,"\n")
end

function query_features_from_bed_file(bedfile::String, query_feature_df::DataFrame)
    df = DataFrame()
    sort!(query_feature_df)
    query_str = query_features_to_query_str( query_feature_df)
    println("Query String: $query_str :")
    bedtools_args = ["-sorted", "-wo","-a", "-", "-b", bedfile]
    bedtools_cmd = `echo -e $query_str` |> `bedtools intersect $bedtools_args`
    println(bedtools_cmd)
    # open iterator
    lines = eachline( bedtools_cmd  )
    firstline = readdlm( IOBuffer( first(lines) ) )

    # get the first line - this will define how many fields are in the dataframe
    println("first line: $firstline")
    for (num,field) in enumerate(firstline)
       fieldname=symbol("x$num") #note just using df[1],df[2]..and so on would make the header x1,x2..automatically poss customise this in future
       df[fieldname]=field
    end

    for line in lines
       fields = readdlm( IOBuffer(line))
       println("FIELDS: $fields")
       push!(df,fields)
    end
    return df
end


text:
using Documenter, GlobalSearchRegression

makedocs(
    # options
    format = Documenter.HTML(assets = ["assets/logo.png"]),
    modules = [GlobalSearchRegression],
    doctest = false,
    clean = false,
    sitename = "GlobalSearchRegression.jl",
    pages = Any[
        "Introduction" => "index.md",
        "API" => "api.md"
    ]
)

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    repo = "github.com/ParallelGSReg/GlobalSearchRegression.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)


text:
using DifferenceEquations
using Distributions
using Optim
using Random
using LinearAlgebra
using Test

# function ar2_transition(u, p, t) # f
#     A = [p[1] p[2]; 1 0]
#     return A * u
# end

# function ar2_noise(u, p, t) # g
#     return [p[3], 0.0]
# end

# function ar2_observation(u, p, t) # h
#     return [1 0] * u
# end

# Random.seed!(1)

# p = [0.5, -0.25]
# tspan = (1, 100)

# phi1 = 0.5
# phi2 = -0.25
# sigma = 1
# true_theta = [phi1, phi2, sigma]

# Y = zeros(tspan[2])
# n = randn(tspan[2])
# Y += sigma .* n

# for t in 3:length(Y)
#     Y[t] += phi1 * Y[t-1] + phi2 * Y[t-2]
# end

# u0 = [Y[2]; Y[1]]

# @testset "AR2" begin
#     @testset "Provided observations" begin
#         prob = StateSpaceProblem(
#             ar2_transition, # f
#             ar2_noise, # g
#             ar2_observation, # h
#             n,
#             [0.0, 0.0], # u0
#             tspan, # timespan
#             nothing, # observation noise
#             Y # observables
#         )

#         trans_x(x) = [x[1], x[2], exp(x[3])]
#         target(x) = -DifferenceEquations._solve(prob, ConditionalGaussian(), trans_x(x)).likelihood
#         res = optimize(target, [0.0, 0.0, 1.0], Optim.Options(iterations=1000))
#         display(res)
#         println(-res.minimum)
#         println(trans_x(res.minimizer))

#         # Calculate whether we inferred phi1 and phi2
#         @test isapprox(res.minimizer[1], phi1, atol=0.01)
#         @test isapprox(res.minimizer[2], phi2, atol=0.01)
#         @test isapprox(exp(res.minimizer[3]), sigma, atol=0.01)
#     end

#     @testset "No observables" begin
#         prob = StateSpaceProblem(
#             ar2_transition, # f
#             ar2_noise, # g
#             ar2_observation, # h
#             n,
#             [0.0, 0.0], # u0
#             tspan, # timespan
#             nothing,#n, # observation noise
#             nothing#Y # observables
#         )

#         sol1 = DifferenceEquations._solve(prob, ConditionalGaussian(), true_theta)
#         sol2 = DifferenceEquations._solve(prob, ConditionalGaussian(), [0.0, 0.0, 1.0])
#         @test sol1.likelihood > sol2.likelihood
#     end
# end

@testset "Linear model" begin

    A = [0.8 0.0; 0.1 0.7]
    B = Diagonal([0.1, 0.5])
    C = [0.5 0.5] # one observable
    R = [0.01]

    # Simulate data
    T = 10
    u0 =[0.0, 0.1]
    tspan = (1, T)

    prob1 = LinearStateSpaceProblem(A, B, C, u0, tspan, R=R)
    sol1 = solve(prob1, LinearGaussian())

    prob2 = LinearStateSpaceProblem(A, B, C, u0, tspan, R=R, observables=sol1.z)
    sol2 = solve(prob2, LinearGaussian())

    prob3 = LinearStateSpaceProblem(A, B, C, u0, tspan, R=R, observables=sol2.z, noise=DefinedNoise(sol2.n))
    sol3 = solve(prob3, LinearGaussian())

    @test sol2.n == sol3.n
end


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "MHD Multicomponent" begin

  @trixi_testset "elixir_mhdmulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [4.32946856e-02, 4.32790096e-02, 2.58822718e-02, 1.62024804e-01, 1.77548775e-02, 1.77564324e-02, 2.69452591e-02, 3.45826315e-03, 1.20451485e-02, 2.40902969e-02],
      linf = [3.28937103e-01, 3.31201677e-01, 1.86307657e-01, 1.03382043e+00, 8.77373400e-02, 8.85611301e-02, 1.46918707e-01, 4.84055914e-02, 8.80607765e-02, 1.76121553e-01])
  end

  @trixi_testset "elixir_mhdmulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
      l2   = [4.25336086e-02, 4.25277388e-02, 2.38503184e-02, 1.15505865e-01, 1.63985643e-02, 1.63986403e-02, 2.58201977e-02, 3.18799432e-04, 1.08086571e-02, 2.16173142e-02],
      linf = [2.34445200e-01, 2.34641311e-01, 1.18747959e-01, 5.33819926e-01, 6.20443752e-02, 6.20664068e-02, 9.73946982e-02, 4.75640645e-03, 6.21579619e-02, 1.24315924e-01])
  end

  @trixi_testset "elixir_mhdmulti_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_convergence.jl"),
      l2   = [2.59813133e-04, 2.59813133e-04, 4.80423374e-04, 3.28657614e-04, 3.67440987e-04, 3.67440987e-04, 4.99348315e-04, 3.88655564e-04, 5.37357512e-05, 1.07471502e-04, 2.14943005e-04],
      linf = [1.54109844e-03, 1.54109844e-03, 2.39709057e-03, 1.01030991e-03, 1.88632092e-03, 1.88632092e-03, 2.59329609e-03, 1.08934177e-03, 2.13125874e-04, 4.26251747e-04, 8.52503495e-04])
  end

  @trixi_testset "elixir_mhdmulti_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_rotor.jl"),
      l2   = [7.23675501e-01, 7.46495099e-01, 0.00000000e+00, 8.17974304e-01, 1.30738332e-01, 9.21666349e-02, 0.00000000e+00, 1.61636223e-02, 1.97849351e-01, 9.89246753e-02],
      linf = [1.07823615e+01, 1.11367261e+01, 0.00000000e+00, 1.22971727e+01, 2.28094366e+00, 2.93333365e+00, 0.00000000e+00, 4.07011671e-01, 3.44365149e+00, 1.72182574e+00],
      tspan = (0.0, 0.01))
end

end

end # module


text:
using Revise
using Distributed
addprocs(1)

@everywhere using Pkg
@everywhere Pkg.activate("/home/evwaco/.julia/dev/FLICameras.jl/")
@everywhere import FLICameras as FC

interface = FC.FLIDOMAIN_USB
deviceType = FC.FLIDEVICE_CAMERA

camList = FC.CameraList(interface, deviceType)
camera = FC.Camera(camList, 1)

dev = FC.create(FC.SharedCamera)
shcam = FC.attach(FC.SharedCamera, dev.shmid)

FC.register(shcam,camera)
dims = (800,800)    # image ROI
remcam = FC.RemoteCamera{UInt8}(shcam, dims)

#--- listening
# 1. broadcasting shmid of cmds, state, img, imgBuftime, remote camera monitor
img_shmid = FC.get_shmid(remcam.img)
imgTime_shmid = FC.get_shmid(remcam.imgTime)
cmds_shmid = FC.get_shmid(remcam.cmds)
shmids = [img_shmid,imgTime_shmid,cmds_shmid]
FC.broadcast_shmids(shmids)

## 2. initialize the server
RemoteCameraEngine = FC.listening(shcam, remcam)

## This rest of the code is supposed to be done by a remote client
begin
    remcam.cmds[1] = FC._to_Cint(SC.CMD_INIT)
    notify(remcam.no_cmds)
end

## 3. configure camera
# update ImageConfigContext in shared camera
new_conf = FC.ImageConfigContext()
#  milli exposure time
new_conf.exposuretime = 50.0
# ROI
new_conf.width = 800
new_conf.height = 800
new_conf.offsetX = (1600-new_conf.width )/2
new_conf.offsetY = (1200-new_conf.height)/2

## configure
FC.set_img_config(shcam,new_conf)

begin
    remcam.cmds[1] = SC._to_Cint(SC.CMD_CONFIG)
    notify(remcam.no_cmds)
end

## 4. start acquisition
begin
    remcam.cmds[1] = SC._to_Cint(SC.CMD_WORK)
    notify(remcam.no_cmds)
end

# 5. stop acquisition
begin
    remcam.cmds[1] = SC._to_Cint(SC.CMD_STOP)
    notify(remcam.no_cmds)
end

#6. update and restart acquisition
begin
    remcam.cmds[1] = SC._to_Cint(SC.CMD_UPDATE)
    notify(remcam.no_cmds)
end


text:
"""
    showpower(p::Point, power::P;
        maxpower = 10.0kW,
        fs = EM,
        symbsize = 1.5fs,
        backgroundcolor = colorant"white",
        labelpower::Bool = true) where {P <: Power}

p                   bottom left corner of indicator.
power               to be shown.
maxpower = 10.0kW   the symbol will be coloured to indicate power / maxpower
fs = EM             font size
symbsize = 1.5EM    lightning bolt size.
backgroundcolor     is used to pick a color for an outline.
labelpower = true   text or no text
"""
function showpower(p::Point, power::Power;
                   maxpower = 15.0kW,
                   fs = EM,
                   symbsize = 1.5fs,
                   backgroundcolor = colorant"white",
                   labelpower::Bool = true)

    fractionpower = power / ( maxpower|> typeof(power))
    rounded = round(unit(power), power, digits = 1)
    symbsize = 1.5fs
    # Prepare color selection
    luminback = lumin(backgroundcolor)
    luminfront = get_current_lumin()
    deltalumin = luminback - luminfront
    avglumin = luminfront + 0.5deltalumin
    curcol = get_current_RGB()
    avgcol = color_with_lumin(curcol, luminfront + 0.5 * deltalumin)
    contrastcol = color_with_lumin(curcol, luminfront + 2 * deltalumin)

    gsave()
    fontsize(symbsize)
    fontface("DejaVu Sans")
    str = "⚡"
    wi, he = textextents(str)[3:4]
    # Outline symbol
    sethue(avgcol)
    fontface("DejaVu Sans Bold")
    fontsize(1.0989010989010988 * symbsize)
    text(str, p - fs.*(0.03, -0.04))
    # Shadow text
    fontface("DejaVu Sans")
    fontsize(fs)
    labelpower && text(string(rounded), p + (1.5wi, 0) - fs.*(0.03, -0.04))

    # Inside symbol
    sethue(color_with_lumin(color_from_palette("red"), luminback))
    fontsize(symbsize)
    text(str, p)
    # Front text
    fontsize(fs)
    labelpower && text(string(rounded), p + (1.5wi, 0))

    # Make a mask of the symbol
    clipreset()
    newpath()
    translate(p)
    fontsize(symbsize)
    textpath(str)
    fillpreserve()
    clip()
    translate(-p)

    # fill the symbol from bottom according to relative power
    sethue(curcol)
    rect(p, wi , -fractionpower * he, :fill)
    sethue(contrastcol)
    setline(PT)
    rect(p, wi , -fractionpower * he, :stroke)

    # Clean up
    clipreset()
    grestore()
    nothing
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using StaticArrays
import SymPy
using SymPy: symbols, simplify
using ACE.SphericalHarmonics: RSHBasis

# These are the expressions we use to test the real SH implementation
#  Y_l^m    =      1/√2 (Y_{lm} - i Y_{l,-m})
#  Y_l^{-m} = (-1)^m/√2 (Y_{lm} + i Y_{l,-m})
# cYt[i_p] = (1/sqrt(2)) * (rY[i_p] - im * rY[i_m])
# cYt[i_m] = (-1)^m * (1/sqrt(2)) * (rY[i_p] + im * rY[i_m])

"""
core of the C->R conversion, uses SymPy to express a product of
complex spherical harmonics as a product of real spherical harmonics, and
then computes the new coefficients.
"""
function convert_c2r_1b(ll, mm, c; verbose=false)
   # length of the tuple / length of the product
   n = length(ll)
   verbose && println("mm = $mm")

   # create symbols representing the real spherical harmonics
   # by analogy with exp(ikx) = cos(kx) + i sin(kx)
   #   So Cj ≡ Y_{lj, |mj|} and Sj ≡ Y_{lj, -|mj|}
   #   and Y_{lj}^{±|mj|} can be expressed in terms of Cj, Sj.
   C = [ symbols("C$j", real=true) for j in 1:length(ll) ]
   S = [ symbols("S$j", real=true) for j in 1:length(ll) ]

   # now express the complex spherical harmonics in terms of the real ones
   # (note we don't need to refer to the cY list since it stores the same
   #  symbols as the C, S lists)
   cY = Vector{Any}(undef, n)
   for j in 1:n
      if mm[j] == 0
         cY[j] = C[j]
      elseif mm[j] > 0
         cY[j] = (1/SymPy.sqrt(2)) * (C[j] - im * S[j])
      else
         cY[j] = ((-1)^mm[j]/SymPy.sqrt(2)) * (C[j] + im * S[j])
      end
    end

   # finally create the symbols for the prefactor
   a = symbols("a", real=true)
   b = symbols("b", real=true)
   coeff = a + im * b

   # let Sympy multiply and simplify the expression
   expr = (real(coeff * prod(cY))).expand()
   verbose && println(expr)

   # next, we need to extract the prefactors
   # to get these we loop over all possible {C,S}1*{C,S}2*... combinations
   # below, CS is to conveniently access the relevant symbols
   # while signs is used to decide which symbol correspond to a +|m| or -|m|
   # basis function
   CS = [S'; C']
   signs = [- ones(Int, n)'; ones(Int, n)']
   real_basis = Any[]

   verbose && println("mm = $mm")
   for ii in CartesianIndices(ntuple(_ -> 2, n))
      ivec = SVector(ii.I...)
      term = prod(CS[ivec[α], α] for α = 1:n)
      pref = (expr.coeff(term)).subs(a, real(c)).subs(b, imag(c))
      if pref == 0; continue; end
      mm_i = [ signs[ivec[α], α] * abs(mm[α])  for α = 1:n ]
      verbose && println(term, " -> (", mm_i, ", ", pref, ")")
      push!(real_basis, (mm = mm_i, c = SymPy.N(pref)))
   end
   return [b for b in real_basis]
end

function convertc2r(cship::SHIP{T, NZ}) where {T, NZ}
   RSH = RSHBasis(cship.SH.maxL, T)
   alists =  ntuple(iz -> deepcopy(cship.alists[iz]), NZ)
   aalists = ntuple(iz -> deepcopy(cship.aalists[iz]), NZ)
   rcoeffs = ntuple(iz -> _convert_c2r_inner(cship.coeffs[iz],
                                             alists[iz],
                                             aalists[iz]), NZ )
   return RSHIP(cship.J, RSH, deepcopy(cship.zlist),
                alists, aalists, rcoeffs)
end

"""
This is the main conversion loop, which is called from the glue code
`convertc2r` and which itself calls the core function
`convert_c2r_1b` which converts a single basis function using SymPy.
"""
function _convert_c2r_inner(ccoeffs::AbstractVector{Complex{T}},
                            alist, aalist) where {T}
   rcoeffs = zeros(T, length(ccoeffs))
   missed = Any[]
   for iAA = 1:length(aalist)
      N = aalist.len[iAA]    # number of As to be multiplied
      # get the (z, k, l, m) infor the the terms in the product
      As = [ alist[aalist.i2Aidx[iAA, α]]  for α = 1:N ]
      izz = Int16[ As[α].z for α = 1:N ]
      kk  = IntS[ As[α].k for α = 1:N ]
      ll  = IntS[ As[α].l for α = 1:N ]
      mm  = IntS[ As[α].m for α = 1:N ]
      # now we convert the ∏ Y_{lα}^{mα} to several ∏ Y_{lα, mα},
      # i.e. we get a LIST of several real basis functions
      rbasis = convert_c2r_1b(ll, mm, ccoeffs[iAA])
      # so now we need to loop through the rbasis and add the new
      # coefficients to rcoeffs as specified by  the basis functions
      for b in rbasis
         # b -> b.mm, b.c -> (ll, b.mm) define a product of real SHs
         # and b.c is the new coefficient; but we first need to identify
         # the correct index in `rcoeffs` which is the same as the
         # corresponding index in `aalist`
         b_mm = IntS.(b.mm)
         # if the izklm key doesn't exist it means we have to add a
         # new AA-basis function to the aalist
         if !haskey(aalist.zklm2i, (izz, kk, ll, b_mm))
            push!(aalist, (izz, kk, ll, b_mm), alist)
            push!(rcoeffs, 0.0)
         end
         i_rAA = aalist.zklm2i[(izz, kk, ll, b_mm)]
         rcoeffs[i_rAA] += b.c
      end
   end
   if length(rcoeffs) > length(ccoeffs)
      @info("Added $(length(rcoeffs) - length(ccoeffs)) basis function(s).")
   end
   return rcoeffs
end


text:
using JuMP
using Dates

const OPTIMAL = JuMP.MathOptInterface.OPTIMAL
const INFEASIBLE = JuMP.MathOptInterface.INFEASIBLE
const UNBOUNDED = JuMP.MathOptInterface.DUAL_INFEASIBLE

function resolvePlne(model, verbose=1, nameOfPL="")

	optimize!(model)

	if verbose == 3
		println(solution_summary(model, verbose=true))
		return model
	end

	status = termination_status(model)

	# if status == JuMP.MathOptInterface.OPTIMAL
	# 	println("Valeur optimale = ", objective_value(model))
	# end

	if verbose >= 2

		last = ""
		for var in all_variables(model)

			which = split(JuMP.name(var), "[")[1]
			if which != last
				last = which
				println("\nVariables $(last)")
			end
			println(var, ": ", value(var))

		end

	end

	if verbose >= 1
		println("\nObjetive value: ", objective_value(model))
	end

	if ! isempty(nameOfPL)
		write_to_file(model, nameOfPL*"_model_"*string(Dates.now())*".mps")
	end

	return model

end

function getModelVariables(model)

end


text:
#=

Author: Shunsuke Hori

=#
using LS
using CES
using BGP

"""
Time Varying Example
"""
PP_seq_time = LS.Planners_Allocation_Sequential(CES.M_time_example) #solve sequential problem

sHist_h = [1, 2, 3, 4, 6, 6, 6]
sHist_l = [1, 2, 3, 5, 6, 6, 6]

sim_seq_h = LS.simulate(PP_seq_time, 1.0, 1, 7, sHist_h)
sim_seq_l = LS.simulate(PP_seq_time, 1.0, 1, 7, sHist_l)

using PyPlot

plt[:figure](figsize=[14, 10])
plt[:subplot](3, 2, 1)
plt[:title]("Consumption")
plt[:plot](sim_seq_l[1], "-ok")
plt[:plot](sim_seq_h[1], "-or")
plt[:subplot](3, 2, 2)
plt[:title]("Labor Supply")
plt[:plot](sim_seq_l[2], "-ok")
plt[:plot](sim_seq_h[2], "-or")
plt[:subplot](3, 2, 3)
plt[:title]("Government Debt")
plt[:plot](sim_seq_l[3], "-ok")
plt[:plot](sim_seq_h[3], "-or")
plt[:subplot](3, 2, 4)
plt[:title]("Tax Rate")
plt[:plot](sim_seq_l[4], "-ok")
plt[:plot](sim_seq_h[4], "-or")
plt[:subplot](3, 2, 5)
plt[:title]("Government Spending")
plt[:plot](CES.M_time_example.G[sHist_l], "-ok")
plt[:plot](CES.M_time_example.G[sHist_h], "-or")
plt[:subplot](3, 2, 6)
plt[:title]("Output")
plt[:plot](CES.M_time_example.Theta[sHist_l].*sim_seq_l[2], "-ok")
plt[:plot](CES.M_time_example.Theta[sHist_h].*sim_seq_h[2], "-or")

plt[:tight_layout]()
plt[:savefig]("TaxSequence_time_varying.png")

plt[:figure](figsize=[8, 5])
plt[:title]("Gross Interest Rate")
plt[:plot](sim_seq_l[end], "-ok")
plt[:plot](sim_seq_h[end], "-or")
plt[:tight_layout]()
plt[:savefig]("InterestRate_time_varying.png")

"""
Time 0 example
"""
PP_seq_time0 = LS.Planners_Allocation_Sequential(CES.M2) #solve sequential problem

B_vec = linspace(-1.5, 1.0, 100)
taxpolicy = hcat([LS.simulate(PP_seq_time0, B_, 1, 2)[4] for B_ in B_vec]...)'
interest_rate = hcat([LS.simulate(PP_seq_time0, B_, 1, 3)[end] for B_ in B_vec]...)'

plt[:figure](figsize=(14, 6))
plt[:subplot](211)
plt[:plot](B_vec,taxpolicy[:, 1], linewidth=2.)
plt[:plot](B_vec,taxpolicy[:, 2], linewidth=2.)

plt[:title]("Tax Rate")
plt[:legend]((latexstring("Time ", "t=0"),
        latexstring("Time ", L"t \geq 1")), loc=2, shadow=true)
plt[:subplot](212)
plt[:title]("Gross Interest Rate")
plt[:plot](B_vec,interest_rate[:, 1], linewidth=2.)
plt[:plot](B_vec,interest_rate[:, 2], linewidth=2.)
plt[:xlabel]("Initial Government Debt")
plt[:tight_layout]()

plt[:savefig]("Time0_taxpolicy.png")

#compute the debt entered with at time 1
B1_vec = hcat([LS.simulate(PP_seq_time0, B_, 1, 2)[3][2] for B_ in B_vec]...)'
#now compute the optimal policy if the government could reset
tau1_reset = hcat([LS.simulate(PP_seq_time0, B1, 1, 1)[4] for B1 in B1_vec]...)'

plt[:figure](figsize=[10, 6])
plt[:plot](B_vec, taxpolicy[:, 2], linewidth=2.)
plt[:plot](B_vec, tau1_reset, linewidth=2.)
plt[:xlabel]("Initial Government Debt")
plt[:title]("Tax Rate")
plt[:legend]((L"\tau_1", L"\tau_1^R"), loc=2, shadow=true)
plt[:tight_layout]()

plt[:savefig]("Time0_inconsistent.png")

"""
BGP Example
"""
#initialize mugrid for value function iteration
muvec = linspace(-0.6, 0.0, 200)

PP_seq = LS.Planners_Allocation_Sequential(BGP.M1) #solve sequential problem

PP_bel = LS.Planners_Allocation_Bellman(BGP.M1,muvec) #solve recursive problem

T = 20
#sHist = utilities.simulate_markov(M1.Pi,0,T)
sHist = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1]

#simulate
sim_seq = LS.simulate(PP_seq, 0.5, 1, T, sHist)
sim_bel = LS.simulate(PP_bel, 0.5, 1, T, sHist)

#plot policies
plt[:figure](figsize=[14, 10])
plt[:subplot](3, 2, 1)
plt[:title]("Consumption")
plt[:plot](sim_seq[1], "-ok")
plt[:plot](sim_bel[1], "-xk")
plt[:legend](("Sequential", "Recursive"), loc="best")
plt[:subplot](3, 2, 2)
plt[:title]("Labor Supply")
plt[:plot](sim_seq[2], "-ok")
plt[:plot](sim_bel[2], "-xk")
plt[:subplot](3, 2, 3)
plt[:title]("Government Debt")
plt[:plot](sim_seq[3], "-ok")
plt[:plot](sim_bel[3], "-xk")
plt[:subplot](3, 2, 4)
plt[:title]("Tax Rate")
plt[:plot](sim_seq[4], "-ok")
plt[:plot](sim_bel[4], "-xk")
plt[:subplot](3, 2, 5)
plt[:title]("Government Spending")
plt[:plot](BGP.M1.G[sHist], "-ok")
plt[:plot](BGP.M1.G[sHist], "-xk")
plt[:plot](BGP.M1.G[sHist], "-^k")
plt[:subplot](3, 2, 6)
plt[:title]("Output")
plt[:plot](BGP.M1.Theta[sHist].*sim_seq[2], "-ok")
plt[:plot](BGP.M1.Theta[sHist].*sim_bel[2], "-xk")

plt[:tight_layout]()
plt[:savefig]("TaxSequence_LS.png")

plt[:figure](figsize=[8, 5])
plt[:title]("Gross Interest Rate")
plt[:plot](sim_seq[end], "-ok")
plt[:plot](sim_bel[end], "-xk")
plt[:legend](("Sequential", "Recursive"), loc="best")
plt[:tight_layout]()


text:
"""
Package: Forecast

    ar(x::DataFrame, or, constant = true;
                     alpha = 1.0, dΦ0 = nothing, dΦ = nothing)

    ar(x::AbstractArray, or, constant = true;
                         alpha = false, dΦ0 = nothing, dΦ = nothing, varnames = nothing)

Fit a multivariate autoregressive series model.

The fitting is done via Ordinary Least Squared and implements the following model:

```math
Xt = \\Phi_0 + \\sum_{i=1}^p \\Phi_i \\cdot X_{t-i} + \\mathcal{N}(\\vec{0},\\Sigma)
```

# Arguments
- `x`: Multivariate series each column containing a dimension and ordered by time ascending rows.
- `or`: Number of parameters Φ to be estimated.
- `constant`: If `true` `ar` estimates Φ0 otherwise it is assume to be zero.
- `alpha`: fixes to zero all coefficients which significance is below its value. It defaults to one.
- `dΦ0`: Tuple containing two Φ0 objects, the first one will act as an original reference to the second one and different values will be fixed in the fitting process to the values in the second Φ0.
- `dΦ`:  Equivalent to dΦ0 but for Φ.
- `varnames`: Names of the dimensions (by default xi where `i` is an integer)

# Returns
An AR object containing the model coefficients, the error sigma matrix, residuals and a collection of information criteria

# Examples
```julia-repl
julia> ar(rand(100,2),2)
AR([...])
"""
function ar(df::DataFrame,
            order::Integer = 1,
            constant::Bool = true;
            alpha::Real = 1.0,
            dΦ0::Tuple = d_ar_dΦ0(size(df,2)-1,constant,Float64),
            dΦ::Tuple = d_ar_dΦ(size(df,2)-1,order,Float64))

    #TODO promote dΦ0 and dΦ types based on df and control for Integer values on x
    ttype = [Date,DateTime,Day,Month,Week]
    ar_df = eltype(df[:,1]) in ttype ? Array{Float64}(df[:,2:end]) : Array{Float64}(df)
    names_df = eltype(df[:,1]) in ttype ? names(df[:,2:end]) : names(df)
    xar = ar(ar_df, order, constant; alpha, dΦ0, dΦ, varnames = names_df)
    xar.x = df

    return xar

end

# Default ar values functions
d_ar_dΦ0(m::Integer,constant::Bool,T::Type) = (repeat([T(1)],m),  repeat([T(constant ? 1 : 0)],m))
d_ar_dΦ(m::Integer,order::Integer,T::Type) =
    (r11=reshape(repeat([T(1)],m*m*order),m,m,order); (r11,r11))
d_ar_varnames(m::Integer)::Vector{String} = ["x"*string(i) for i in 1:m]

function ar(x::AbstractArray{T},
            order::Integer = 1,
            constant::Bool = true;
            alpha::Real = 1.0,
            dΦ::Tuple = d_ar_dΦ(size(x,2),order,T),
            dΦ0::Tuple = d_ar_dΦ0(size(x,2),constant,T),
            varnames::AbstractVector{String} = d_ar_varnames(size(x,2))) where T<:Real

    @assert 0.0 < alpha <= 1.0

    m = size(x,2)
    np = order

    dΦ  = ( expand(dΦ[1],(m,m,np)), expand(dΦ[2],(m,m,np)) )
    dΦ0 = ( expand(dΦ0[1],(m,)),    expand(dΦ0[2],(m,)) )

    xar = ar_ols(x, order, constant; dΦ0, dΦ, varnames)

    alpha == 1.0 && return xar

    m, np = xar.ndims, xar.order

    Φ = expand(xar.Φ,(m,m,np))
    Φ0 = expand(xar.Φ0,(m,))
    Φpv = expand(xar.Φpv,(m,m,np))
    Φ0pv = expand(xar.Φ0pv,(m,))

    dΦs = (Φpv .<= alpha)
    dΦ0s = (Φ0pv .<= alpha)

    all(dΦs) && all(dΦ0s) && return xar

    if isnothing(dΦ)
        dΦ = (Φ, Φ .* dΦs)
        dΦ0 = (Φ0, Φ0 .* dΦ0s)
    else
        dΦ = (dΦ[1], dΦ[2] .* dΦs)
        dΦ0 = (dΦ0[1], dΦ0[2] .* dΦ0s)
    end

    return ar_ols(x, order, constant; dΦ0, dΦ, varnames)

end

function ar_ols(x::AbstractArray{T},
                or::Integer,
                constant::Bool;
                dΦ::Tuple,
                dΦ0::Tuple,
                varnames::AbstractVector{String}) where T<:Real

    @assert 1 <= ndims(x) <= 2
    @assert 1 <= or < size(x,1)-1

    n = size(x,1)
    m = size(x,2)
    np = m*m*or + m

    x = x[end:-1:1,:]
    nx = x
    for i in 1:or
        nx = hcat(x,nx[vcat(2:end,1),:])
    end
    nx = nx[1:end-or,:]

    Y = nx[:,1:m]
    X = nx[:,m+1:end]
    X = hcat(repeat([1.0],size(X,1)),X)

    # Fixing parameters
    W = Array{T,2}(undef,(or*m+1,m))
    for i in 1:m
        Xi,Yi = fixΦ(X,Y,i,dΦ0,dΦ)
        dW = Xi\Yi
        W[:,i] = fixW(dW,i,dΦ0,dΦ)
    end

    # Fitted values
    fitted = X*W

    # Prediction error
    e = Y - fitted

    Φ0 = W[1,:]
    Φ = W[2:end,:]
    Φ = Array(reshape(Φ',m,m,or))

    # Maximum Likelihood noise covariance
    k = or*m*m
    @assert n-k > 0 "Insufficient data for the model"
    Σ2 = 1/(n-k)*e'*e

    # ML parameters std. error.
    Φse = sqrt.(abs.(diag(kron(Σ2, (X'*X)^-1))))
    Φse = fixΦse(Φse,dΦ0,dΦ)

    rΦse = reshape(Φse,:,m)'
    Φ0se = rΦse[:,1]
    Φse  = reshape(rΦse[:,2:end],m,m,:)

    #fix p
    np = fixnp(dΦ0,dΦ)

    # TODO check weher to use 'or' or 'np' and if 'or' then a fixor function is needed
    # Information Criteria
    lΣ2   = log(norm(Σ2))
    ic = Dict([("AIC",  lΣ2 + 2*np*m^2/n),
               ("AICC", lΣ2 + 2*(np*m^2+1)/(n-(np*m^2+2))),
               ("BIC",  n*log(norm(Σ2)+m)+(m^2*np+m*(m+1)/2)*log(n)),
               ("H&Q",  lΣ2 + 2*log(log(n))*np*m^2/n)])

    # Statistics
    SStot = sum((Y .- mean(Y,dims=1)).^2,dims=1)
    SSres = sum(e .^ 2,dims=1)
    R2 = vec(1 .- (SSres ./ SStot))

    pvf(mu,se) = se == 0 ? 1.0 : cdf(Normal(abs(mu),se),0) #1 to make log(1) = 0
    Φpv = pvf.(Φ,Φse)
    Φ0pv = pvf.(Φ0,Φ0se)
    slpv = vec(reshape(sum(log.(vcat(reshape(Φ0pv,:,m),reshape(Φpv,:,m))),dims=1),:,1))

    stats = Dict([(" Variable", varnames),
                  ("R2",    R2),
                  ("R2adj", 1 .- (1 .- R2) * (n-1)/(n-(np-1)/m-1)),
                  ("Fisher's p-test", 1 .- cdf(Chisq(2*np/m),-2*slpv))])

    Σ = sqrt.(diag(Σ2))

    call = "ar(X, order="*string(or)*
        ", constant="*string(constant)*")"

    AR(varnames,
       or, m,
       compact.([Φ,Φ0,Σ2,Σ,x[end:-1:1,:],fitted,e])...,
       ic,stats,
       compact.([Φse,Φ0se,Φpv,Φ0pv])...,
       call)
end

"""
Package: Forecast

    fixΦ(X,Y,dΦ0,dΦ)

For a given X and Y OLS matrices returns the X and Y resulting from fixing parameters given dΦ0 and dΦ
"""
function fixΦ(X::AbstractMatrix{T},
              Y::AbstractMatrix{T},
              i::Integer,
              dΦ0::Tuple{AbstractArray{T},AbstractArray{T}},
              dΦ::Tuple{AbstractArray{T},AbstractArray{T}}) where T<:Real

    Φ0, fΦ0 = dΦ0
    Φ, fΦ = dΦ
    (Φ0 == fΦ0) & (Φ == fΦ) && return((X,Y[:,i]))

    rΦ0 = reshape(Φ0,:,1)
    rΦ  = reshape(Φ,size(Φ,1),:)
    rΦ  = hcat(rΦ0,rΦ)

    rfΦ0 = reshape(fΦ0,:,1)
    rfΦ = reshape(fΦ,size(fΦ,1),:)
    rfΦ = hcat(rfΦ0,rfΦ)

    dc = findall(rΦ[i,:] .!== rfΦ[i,:])

    @assert length(rΦ) > length(dc) "No degrees of freedom"

    # Create new Yi
    Yi = Y[:,i]
    for c in dc
        Yi = Yi .- rfΦ[i,c]*X[:,c]
    end

    # Create new X
    Xi = drop(X,c=dc)

    (Xi,Yi)

end

"""
Package: Forecast

    fixW(W,dΦ0,dΦ)

For a given Weight matrix returns a version with fixed values based on dΦ0 and dΦ
"""
function fixW(dWi::AbstractVector{T},
              i::Integer,
              dΦ0::Tuple{AbstractArray{T},AbstractArray{T}},
              dΦ::Tuple{AbstractArray{T},AbstractArray{T}}) where T<:Real

    Φ0, fΦ0 = dΦ0
    Φ, fΦ = dΦ
    (Φ0 == fΦ0) & (Φ == fΦ) && return(dWi)

    rΦ0 = reshape(Φ0,:,1)
    rΦ  = reshape(Φ,size(Φ,1),:)
    rΦ  = hcat(rΦ0,rΦ)

    rfΦ0 = reshape(fΦ0,:,1)
    rfΦ = reshape(fΦ,size(fΦ,1),:)
    rfΦ = hcat(rfΦ0,rfΦ)

    dc = findall(rΦ[i,:] .!== rfΦ[i,:])

    for c in dc
        insert!(dWi,c,rfΦ[i,c])
    end

    dWi

end

"""
Package: Forecast

    fixnp(dΦ0,dΦ)

return the number of free parameters
"""
function fixnp(dΦ0::Tuple{AbstractArray{T},AbstractArray{T}},
               dΦ::Tuple{AbstractArray{T},AbstractArray{T}}) where T<:Real

    Φ0, fΦ0 = dΦ0
    Φ, fΦ = dΦ
    (Φ0 == fΦ0) & (Φ == fΦ) && return(length(Φ) + length(Φ0))

    rΦ0 = reshape(Φ0,:,1)
    rΦ  = reshape(Φ,size(Φ,1),:)
    rΦ  = hcat(rΦ0,rΦ)

    rfΦ0 = reshape(fΦ0,:,1)
    rfΦ = reshape(fΦ,size(fΦ,1),:)
    rfΦ = hcat(rfΦ0,rfΦ)

    dc = findall(rΦ .!== rfΦ)

    length(Φ) + length(Φ0) - length(dc)
end

"""
Package: Forecast

    fixΦse(M,dΦ0,dΦ)

For a given `se` matrix returns an version with zeroes based on dΦ0 and dΦ
"""
function fixΦse(Φse::AbstractVector{T},
                dΦ0::Tuple{AbstractArray{T},AbstractArray{T}},
                dΦ::Tuple{AbstractArray{T},AbstractArray{T}}) where T<:Real

    Φ0, fΦ0 = dΦ0
    Φ, fΦ = dΦ
    (Φ0 == fΦ0) & (Φ == fΦ) && return(Φse)

    _,m,np = size(Φ)

    rΦ  = reshape(hcat(Φ0 ,reshape(Φ ,m,m*np))',m*(m*np+1),1)
    rfΦ = reshape(hcat(fΦ0,reshape(fΦ,m,m*np))',m*(m*np+1),1)

    dc = findall(rΦ .!== rfΦ)

    fΦse = Vector{T}(undef, length(Φse))
    fΦse[:] = Φse
    for ci in dc
        c = Tuple(ci)
        fΦse[c[1]] = 0.0
    end

    fΦse

end


text:
#!/usr/bin/env julia6
using LESFilter, ReadGlobal, InplaceRealFFT

function filter(inputfile::String,N::Int,fil::String)
    nx,ny,nz,xs,ys,zs = getdimsize()
    dtype, padded = checkinput(inputfile,nx,ny,nz)
    field = PaddedArray{dtype}(inputfile,(nx,ny,nz),padded)

    # Ugly fix in case dtype is Float32
    if eltype(field) == Complex{Float32}
        field2 = field
        field = PaddedArray{Float64}(size(real(field2)))
        real(field) .= real(field2)
    end

    boxdim = N*2*zs*π/nz
    isfile("fftw_wisdom") && FFTW.import_wisdom("fftw_wisdom")
    info("Filtering File $inputfile")
    lesfilter!(field, fil, boxdim, (xs,ys,zs))
    info("Done")

    if fil == "G"
        pathoutput = "./Filtered_Fields/Gaussian/N$N"
    elseif fil =="C"
        pathoutput = "./Filtered_Fields/CutOff/N$N"
    elseif fil =="B"
        pathoutput = "./Filtered_Fields/Box/N$N"
    end

    isdir(pathoutput) || mkpath(pathoutput)
    outfilename = "$(pathoutput)/$(fil)$(split(split(inputfile,"/")[end],".")[1])_N$(N)"
    info("Saving filtered field in $outfilename")
    write(outfilename,field.r)
    FFTW.export_wisdom("fftw_wisdom")
end

function filter(inputfile::String, N1::Int, N2::Int, fil::String)
    nx,ny,nz,xs,ys,zs = getdimsize()
    dtype, padded = checkinput(inputfile,nx,ny,nz)
    field = PaddedArray{dtype}(inputfile,(nx,ny,nz),padded)

    # Ugly fix in case dtype is Float32
    if eltype(field) == Complex{Float32}
        field2 = field
        field = PaddedArray{Float64}(size(real(field2)))
        real(field) .= real(field2)
    end

    boxdimxy = N1*2*zs*π/nz
    boxdimz = N2*2*zs*π/nz
    isfile("fftw_wisdom") && FFTW.import_wisdom("fftw_wisdom")
    info("Filtering File $inputfile")
    lesfilter!(field, fil, (boxdimxy, boxdimz), (xs,ys,zs))
    info("Done")

    if fil == "G"
        pathoutput = "./Filtered_Fields/Anisotropic/Gaussian/N$(N2)N$(N1)"
    elseif fil =="C"
        pathoutput = "./Filtered_Fields/Anisotropic/CutOff/N$(N2)N$(N1)"
    end

    isdir(pathoutput) || mkpath(pathoutput)
    outfilename = "$(pathoutput)/$(fil)$(split(split(inputfile,"/")[end],".")[1])_N$(N2)N$(N1)"
    info("Saving filtered field in $outfilename")
    write(outfilename,field.r)
    FFTW.export_wisdom("fftw_wisdom")
end

if length(ARGS) == 3
    inputfile = ARGS[1]
    N = parse(Int,ARGS[2])
    fil = ARGS[3]

    filter(inputfile,N,fil)
else
    inputfile = ARGS[1]
    N1 = parse(Int,ARGS[2])
    N2 = parse(Int,ARGS[3])
    fil = ARGS[4]

    filter(inputfile,N1,N2,fil)
end


text:
struct SequentialBackend <: AbstractBackend end

const sequential = SequentialBackend()

function get_part_ids(b::SequentialBackend,nparts::Integer)
  parts = [ part for part in 1:nparts ]
  SequentialData(parts)
end

function get_part_ids(b::SequentialBackend,nparts::Tuple)
  parts = collect(LinearIndices(nparts))
  SequentialData(parts)
end

function prun_debug(driver::Function,b::SequentialBackend,nparts)
  prun(driver,b,nparts)
end

struct SequentialData{T,N} <: AbstractPData{T,N}
  parts::Array{T,N}
end

Base.size(a::SequentialData) = size(a.parts)

i_am_main(a::SequentialData) = true

get_backend(a::SequentialData) = sequential

function Base.iterate(a::SequentialData)
  next = map_parts(iterate,a)
  if eltype(next.parts) == Nothing || any(i->i==Nothing,next.parts)
    return nothing
  end
  item = map_parts(first,next)
  state = map_parts(_second,next)
  item, state
end

_second(a) = a[2]

function Base.iterate(a::SequentialData,state)
  next = map_parts(iterate,a,state)
  if eltype(next.parts) == Nothing || any(i->i==Nothing,next.parts)
    return nothing
  end
  item = map_parts(first,next)
  state = map_parts(_second,next)
  item, state
end

function map_parts(task,args::SequentialData...)
  @assert length(args) > 0
  @assert all(a->length(a.parts)==length(first(args).parts),args)
  parts_in = map(a->a.parts,args)
  parts_out = map(task,parts_in...)
  SequentialData(parts_out)
end

function Base.show(io::IO,k::MIME"text/plain",data::SequentialData)
  for part in 1:num_parts(data)
    if part != 1
      println(io,"")
    end
    println(io,"On part $part of $(num_parts(data)):")
    show(io,k,data.parts[part])
  end
end

get_part(a::SequentialData,part::Integer) = a.parts[part]
get_part(a::SequentialData) = get_main_part(a)

function gather!(rcv::SequentialData,snd::SequentialData)
  @assert num_parts(rcv) == num_parts(snd)
  @assert length(rcv.parts[MAIN]) == num_parts(snd)
  for part in 1:num_parts(snd)
    rcv.parts[MAIN][part] = snd.parts[part]
  end
  rcv
end

function gather!(rcv::SequentialData{<:Table},snd::SequentialData)
  @assert num_parts(rcv) == num_parts(snd)
  @assert length(rcv.parts[MAIN]) == num_parts(snd)
  for part in 1:num_parts(snd)
    offset = rcv.parts[MAIN].ptrs[part]-1
    for i in 1:length(snd.parts[part])
      rcv.parts[MAIN].data[i+offset] = snd.parts[part][i]
    end
  end
  rcv
end

function gather_all!(rcv::SequentialData,snd::SequentialData)
  @assert num_parts(rcv) == num_parts(snd)
  for part_rcv in 1:num_parts(rcv)
    @assert length(rcv.parts[part_rcv]) == num_parts(snd)
    for part_snd in 1:num_parts(snd)
      rcv.parts[part_rcv][part_snd] = snd.parts[part_snd]
    end
  end
  rcv
end

function gather_all!(rcv::SequentialData{<:Table},snd::SequentialData)
  @assert num_parts(rcv) == num_parts(snd)
  for part_rcv in 1:num_parts(rcv)
    @assert length(rcv.parts[part_rcv]) == num_parts(snd)
    for part_snd in 1:num_parts(snd)
      offset = rcv.parts[part_rcv].ptrs[part_snd]-1
      for i in 1:length(snd.parts[part_snd])
        rcv.parts[part_rcv].data[i+offset] = snd.parts[part_snd][i]
      end
    end
  end
  rcv
end

function scatter(snd::SequentialData)
  @assert length(snd.parts[MAIN]) == num_parts(snd)
  parts = similar(snd.parts,eltype(snd.parts[MAIN]),size(snd.parts))
  copyto!(parts,snd.parts[MAIN])
  SequentialData(parts)
end

function async_exchange!(
  data_rcv::SequentialData,
  data_snd::SequentialData,
  parts_rcv::SequentialData,
  parts_snd::SequentialData,
  t_in::AbstractPData)

  @check num_parts(parts_rcv) == num_parts(data_rcv)
  @check num_parts(parts_rcv) == num_parts(data_snd)
  @check num_parts(parts_rcv) == num_parts(t_in)

  map_parts(schedule,t_in)
  map_parts(wait,t_in)

  @boundscheck _check_rcv_and_snd_match(parts_rcv,parts_snd)
  for part_rcv in 1:num_parts(parts_rcv)
    for (i, part_snd) in enumerate(parts_rcv.parts[part_rcv])
      j = first(findall(k->k==part_rcv,parts_snd.parts[part_snd]))
      data_rcv.parts[part_rcv][i] = data_snd.parts[part_snd][j]
    end
  end

  t_out = map_parts(data_rcv) do data_rcv
    @task nothing
  end
  t_out
end

function _check_rcv_and_snd_match(parts_rcv::SequentialData,parts_snd::SequentialData)
  @check num_parts(parts_rcv) == num_parts(parts_snd)
  for part in 1:num_parts(parts_rcv)
    for i in parts_rcv.parts[part]
      @check length(findall(k->k==part,parts_snd.parts[i])) == 1
    end
    for i in parts_snd.parts[part]
      @check length(findall(k->k==part,parts_rcv.parts[i])) == 1
    end
  end
  true
end

function async_exchange!(
  data_rcv::SequentialData{<:Table},
  data_snd::SequentialData{<:Table},
  parts_rcv::SequentialData,
  parts_snd::SequentialData,
  t_in::SequentialData)

  @check num_parts(parts_rcv) == num_parts(data_rcv)
  @check num_parts(parts_rcv) == num_parts(data_snd)
  @check num_parts(parts_rcv) == num_parts(t_in)

  map_parts(schedule,t_in)
  map_parts(wait,t_in)

  @boundscheck _check_rcv_and_snd_match(parts_rcv,parts_snd)
  for part_rcv in 1:num_parts(parts_rcv)
    for (i, part_snd) in enumerate(parts_rcv.parts[part_rcv])
      j = first(findall(k->k==part_rcv,parts_snd.parts[part_snd]))
      ptrs_rcv = data_rcv.parts[part_rcv].ptrs
      ptrs_snd = data_snd.parts[part_snd].ptrs
      @check ptrs_rcv[i+1]-ptrs_rcv[i] == ptrs_snd[j+1]-ptrs_snd[j]
      for p in 1:(ptrs_rcv[i+1]-ptrs_rcv[i])
        p_rcv = p+ptrs_rcv[i]-1
        p_snd = p+ptrs_snd[j]-1
        data_rcv.parts[part_rcv].data[p_rcv] = data_snd.parts[part_snd].data[p_snd]
      end
    end
  end

  t_out = map_parts(data_rcv) do data_rcv
    @task nothing
  end
  t_out
end


text:
AXIS = ["x","y"]
NEV = 5
N = 400
ΔN = 10

ϕ = linspace(0,0.0005,100)
I = (1,1) #(sqrt(2),1,0).*4
G = (1,2) #(1,1,0)
L = 20
α = 1/6
σ = -0.2 - 2*α^2
#@everywhere
#N = 600 ### What is a discretization of the box

@everywhere include("exciton-hamiltonian.jl")
@everywhere function energies(α,I,G,p0,L,N,NEV,σoff)

    σ = σoff - maximum(I)^2/2 - maximum(I)*maximum(p0)

    println("σ = $σ")
    println("p0 = $p0")
    run(`hostname`)
    BLAS.set_num_threads(1)

    function EXTRAPOLATOR(x,y)
        x1,x2 = x[1],x[2]
        y1,y2 = y[1],y[2]

        A = (x1^3*y2 - x2^3*y1)/(x1^3-x2^3)
        B = (y1-y2)/(x1^3-x2^3)

        return A
    end

    Hm = Hamiltonian(α,I,G,p0,L,N[1])
    @time E1 = eigs(Hm,nev=NEV,tol=1e-6,maxiter=10000,which=:LM,sigma=σ)[1]
    E1 = sort(real(E1))

    Hm = Hamiltonian(α,I,G,p0,L,N[2])
    @time E2 = eigs(Hm,nev=NEV,tol=1e-6,maxiter=10000,which=:LM,sigma=σ)[1]
    E2 = sort(real(E2))

    Eres = EXTRAPOLATOR(1./N,(E1,E2))

    return Dict(:EVALS=>Eres,:E1=>E1,:E2=>E2,:N1=>N[1],:N2=>N[2]) #real(solution[1]),spin
end

knotsx = Iterators.product((ϕ,0:0,0:0)...)
knotsy = Iterators.product((0:0,ϕ,0:0)...)
knotsz = Iterators.product((0:0,0:0,ϕ)...)
knotsxy = zip(((ϕ, ϕ, zeros(ϕ)) ./ sqrt(2))...)
knotsyz = zip(((zeros(ϕ), ϕ, ϕ) ./ sqrt(2))...)
knotszx = zip(((ϕ, zeros(ϕ), ϕ) ./ sqrt(2))...)

knots = []
for i in AXIS
    vi = eval(parse("knots$i"))
    push!(knots,collect(vi))
end

knots = vcat([i[:] for i in knots]...)

results = pmap(x->energies(α,I[1:2],G[1:2],x[1:2],L,(N-ΔN,N+ΔN),NEV,σ),knots)

parlist = []
Lright=0
for i in AXIS
    push!(parlist,"knots$i")
    push!(parlist,eval(parse("knots$i")))
    push!(parlist,"results$i")

    Li = length(eval(parse("knots$i")))
    Lright += Li
    Lleft = 1+Lright-Li
    push!(parlist,results[Lleft:Lright])
end

using JLD
save("outdir/columbenergies-$N $ΔN $(I[1:2]) $α.jld","knots",knots,"results",results,parlist...)

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc