text:
function test_episode(e::Episode, controller::Function)
    for i in 1:10
        inputs = get_inputs(e)
        @test length(inputs) == 10
        @test all(inputs .>= 0.0)
        @test all(inputs .<= 1.0)
        output = controller(inputs)
        println(RoboGrid.ACTIONS[output])
        @time step!(e, output)
        @test e.robot.y >= 1
        @test e.robot.y <= size(e.grid.cells, 1)
        @test e.robot.x >= 1
        @test e.robot.x <= size(e.grid.cells, 2)
    end
end

@testset "Inputs" begin
    g = Grid("maps/empty.yaml")
    e = Episode(g)
    obj = rand()
    g.cells[e.robot.y, e.robot.x].obj = obj
    inputs = get_inputs(e)
    orig_inputs = copy(inputs)
    @test inputs[1] == 0.0
    @test inputs[3] == RoboGrid.type_to_float("wall")
    @test inputs[6] == obj
    RoboGrid.turn_left!(e.robot)
    inputs = get_inputs(e)
    @test inputs[1] == 0.0
    @test inputs[3] == RoboGrid.type_to_float("wall")
    @test inputs[6] == obj
    RoboGrid.turn_left!(e.robot)
    inputs = get_inputs(e)
    @test inputs[1] == 0.0
    @test inputs[3] == RoboGrid.type_to_float("empty")
    @test inputs[6] == obj
    RoboGrid.turn_left!(e.robot)
    inputs = get_inputs(e)
    @test inputs[1] == 0.0
    @test inputs[3] == RoboGrid.type_to_float("empty")
    @test inputs[6] == obj
    RoboGrid.turn_left!(e.robot)
    inputs = get_inputs(e)
    @test all(inputs .== orig_inputs)
end

@testset "Random controller" begin
    e = Episode(Grid())
    c(x::Array{Float64}) = rand(1:length(RoboGrid.ACTIONS))
    test_episode(e, c)
end

@testset "Run and terminate" begin
    e = Episode(Grid(); meta=Dict{String, Any}("max_step"=>10))
    c(x::Array{Float64}) = rand(1:length(RoboGrid.ACTIONS))
    total_reward = run!(e, c)
    @test total_reward == 0.0
    @test e.meta["step"] <= e.meta["max_step"]
    @test e.meta["time"] <= e.meta["max_time"]
    @test e.meta["memory"] <= e.meta["max_memory"]
end


text:
"""
k: an index from subspace
indices: array of indices from comspace
"""
# generic U, single state
@inline function subspace_mul_kernel_generic!(S::AbstractVector{T}, y, indices, U, k::Int, offset::Int) where T
    @inbounds for i in axes(U, 1)
        y[i] = zero(T)
        for j in axes(U, 2)
            idx_j = k + indices[j] + offset
            y[i] += U[i, j] * S[idx_j]
        end
    end

    @inbounds for i in axes(U, 1)
        idx_i = k + indices[i] + offset
        S[idx_i] = y[i]
    end
    return
end

# generic U, multiple states
@inline function subspace_mul_kernel_generic!(S::AbstractMatrix{T}, y, indices, U, k::Int, b::Int, offset::Int) where T
    @inbounds for i in axes(U, 1)
        y[i] = zero(T)
        for j in axes(U, 2)
            idx_j = k + indices[j] + offset
            y[i] += U[i, j] * S[b, idx_j]
        end
    end

    @inbounds for i in axes(U, 1)
        idx_i = k + indices[i] + offset
        S[b, idx_i] = y[i]
    end
    return
end

# Diagonal U, single state
@inline function subspace_mul_kernel_generic!(S::AbstractVector{T}, y, indices, U::Diagonal, k::Int, offset::Int) where T
    @inbounds for i in axes(U, 1)
        idx_i = k + indices[i] + offset
        S[idx_i] = U[i, i] * S[idx_i]
    end
    return
end

# Diagonal U, multiple states
@inline function subspace_mul_kernel_generic!(S::AbstractMatrix{T}, y, indices, U::Diagonal, k::Int, b::Int, offset::Int) where T
    @inbounds for i in axes(U, 1)
        idx_i = k + indices[i] + offset
        S[b, idx_i] = U[i, i] * S[b, idx_i]
    end
    return
end

@inline function subspace_mul_kernel_generic_4x4!(S::AbstractVector{T}, y, indices, U, k::Int, offset::Int) where T
    @inbounds begin
        Base.Cartesian.@nexprs 4 i -> begin
            y[i] = zero(T)
            Base.Cartesian.@nexprs 4 j -> begin
                idx_j = k + indices[j] + offset
                y[i] += U[i, j] * S[idx_j]
            end
        end

        Base.Cartesian.@nexprs 4 i -> begin
            idx_i = k + indices[i] + offset
            S[idx_i] = y[i]
        end
    end
    return
end

@inline function subspace_mul_kernel_generic_4x4!(S::AbstractMatrix{T}, y, indices, U, k::Int, b::Int, offset::Int) where T
    @inbounds begin
        Base.Cartesian.@nexprs 4 i -> begin
            y[i] = zero(T)
            Base.Cartesian.@nexprs 4 j -> begin
                idx_j = k + indices[j] + offset
                y[i] += U[i, j] * S[b, idx_j]
            end
        end

        Base.Cartesian.@nexprs 4 i -> begin
            idx_i = k + indices[i] + offset
            S[b, idx_i] = y[i]
        end
    end
    return
end

@inline function subspace_mul_kernel_generic_4x4!(S::AbstractVector{T}, y, indices, U::Diagonal, k::Int, offset::Int) where T
    @inbounds begin
        Base.Cartesian.@nexprs 4 i -> begin
            idx_i = k + indices[i] + offset
            S[idx_i] = U[i, i] * S[idx_i]
        end
    end
    return
end

@inline function subspace_mul_kernel_generic_4x4!(S::AbstractMatrix{T}, y, indices, U::Diagonal, k::Int, b::Int, offset::Int) where T
    @inbounds begin
        Base.Cartesian.@nexprs 4 i -> begin
            idx_i = k + indices[i] + offset
            S[b, idx_i] = U[i, i] * S[b, idx_i]
        end
    end
    return
end

ext: jl


text:
using Pkg; Pkg.add("Documenter")
using Documenter, SeisMain

makedocs(
	modules = [SeisMain],
	sitename = "SeisMain.jl",
	format = Documenter.HTML(),
	pages = [
	       "Home" => "index.md",
	       "Manual" => Any[
	       		"Guide" => "man/guide.md"],
		"Library" => Any[
			  "Public" => "lib/public.md",
				],
		],
)

deploydocs(
    repo = "github.com/SeismicJulia/SeisMain.jl.git",
)


text:
__precompile__(true)
module Mehs

import Base: *, <, ==, !=, <=, !, +, -, ^, /, &, |, xor
using Compat: AbstractRange

export ismeh, meh, mehs, Meh, MehException, levels

"""
    Meh

A type with no fields whose singleton instance [`meh`](@ref) is used
to represent meh values.
"""
struct Meh end

"""
    meh

The singleton instance of type [`Meh`](@ref) representing a meh value.
"""
const meh = Meh()

Base.show(io::IO, x::Meh) = print(io, "meh")

"""
    MehException(msg)

Exception thrown when a [`meh`](@ref) value is encountered in a situation
where it is not supported. The error message, in the `msg` field
may provide more specific details.
"""
struct MehException <: Exception
    msg::AbstractString
end

Base.showerror(io::IO, ex::MehException) =
    print(io, "MehException: ", ex.msg)

T(::Type{Union{T1, Meh}}) where {T1} = T1
T(::Type{Meh}) = Union{}
T(::Type{T1}) where {T1} = T1
T(::Type{Any}) = Any

ismeh(::Any) = false
ismeh(::Meh) = true

# vector constructors
mehs(dims...) = fill(meh, dims)
mehs(::Type{T}, dims...) where {T >: Meh} = fill!(Array{T}(dims), meh)
mehs(::Type{T}, dims...) where {T} = fill!(Array{Union{T, Meh}}(dims), meh)

Base.promote_rule(::Type{T}, ::Type{Meh}) where {T} = Union{T, Meh}
Base.promote_rule(::Type{T}, ::Type{Union{S,Meh}}) where {T,S} = Union{promote_type(T, S), Meh}
Base.promote_rule(::Type{T}, ::Type{Any}) where {T} = Any
Base.promote_rule(::Type{Any}, ::Type{Meh}) = Any
Base.promote_rule(::Type{Meh}, ::Type{Any}) = Any
Base.promote_rule(::Type{Meh}, ::Type{Meh}) = Meh

Base.convert(::Type{Union{T, Meh}}, x) where {T} = convert(T, x)

# Comparison operators
==(::Meh, ::Meh) = meh
==(::Meh, b) = meh
==(a, ::Meh) = meh
# != must be defined explicitly since fallback expects a Bool
!=(::Meh, ::Meh) = meh
!=(::Meh, b) = meh
!=(a, ::Meh) = meh
Base.isequal(::Meh, ::Meh) = true
Base.isequal(::Meh, b) = false
Base.isequal(a, ::Meh) = false
<(::Meh, ::Meh) = meh
<(::Meh, b) = meh
<(a, ::Meh) = meh
Base.isless(::Meh, ::Meh) = false
Base.isless(::Meh, b) = false
Base.isless(a, ::Meh) = true
if VERSION < v"0.7.0-DEV.300"
    <=(::Meh, ::Meh) = meh
    <=(::Meh, b) = meh
    <=(a, ::Meh) = meh
end

# Unary operators/functions
for f in (:(!), :(+), :(-), :(Base.identity), :(Base.zero),
          :(Base.abs), :(Base.abs2), :(Base.sign),
          :(Base.acos), :(Base.acosh), :(Base.asin), :(Base.asinh), :(Base.atan), :(Base.atanh),
          :(Base.sin), :(Base.sinh), :(Base.cos), :(Base.cosh), :(Base.tan), :(Base.tanh),
          :(Base.exp), :(Base.exp2), :(Base.expm1), :(Base.log), :(Base.log10), :(Base.log1p),
          :(Base.log2), :(Base.exponent), :(Base.sqrt), :(Base.gamma), :(Base.lgamma),
          :(Base.iseven), :(Base.ispow2), :(Base.isfinite), :(Base.isinf), :(Base.isodd),
          :(Base.isinteger), :(Base.isreal), :(Base.isimag), :(Base.isnan), :(Base.isempty),
          :(Base.iszero), :(Base.transpose), :(Base.ctranspose), :(Base.float))
    @eval $(f)(d::Meh) = meh
end

Base.zero(::Type{Union{T, Meh}}) where {T <: Number} = zero(T)
Base.zero(::Type{Union{T, Meh}}) where {T <: Base.Dates.Period} = zero(T)

# Binary operators/functions
for f in (:(+), :(-), :(*), :(/), :(^),
          :(Base.div), :(Base.mod), :(Base.fld), :(Base.rem), :(Base.min), :(Base.max))
    @eval begin
        # Scalar with meh
        ($f)(::Meh, ::Meh) = meh
        ($f)(d::Meh, x::Number) = meh
        ($f)(d::Number, x::Meh) = meh
    end
end

# Rounding and related functions
for f in (:(Base.ceil), :(Base.floor), :(Base.round), :(Base.trunc))
    @eval begin
        ($f)(::Meh, digits::Integer=0, base::Integer=0) = meh
        ($f)(::Type{>:Meh}, ::Meh) = meh
        ($f)(::Type{T}, ::Meh) where {T} =
            throw(MehException("cannot convert a meh value to type $T"))
    end
end

# to avoid ambiguity warnings
(^)(::Meh, ::Integer) = meh

# Bit operators
(&)(::Meh, ::Meh) = meh
(&)(a::Meh, b::Bool) = ifelse(b, meh, false)
(&)(b::Bool, a::Meh) = ifelse(b, meh, false)
(&)(::Meh, ::Integer) = meh
(&)(::Integer, ::Meh) = meh
(|)(::Meh, ::Meh) = meh
(|)(a::Meh, b::Bool) = ifelse(b, true, meh)
(|)(b::Bool, a::Meh) = ifelse(b, true, meh)
(|)(::Meh, ::Integer) = meh
(|)(::Integer, ::Meh) = meh
xor(::Meh, ::Meh) = meh
xor(a::Meh, b::Bool) = meh
xor(b::Bool, a::Meh) = meh
xor(::Meh, ::Integer) = meh
xor(::Integer, ::Meh) = meh

# String functions
*(d::Meh, x::AbstractString) = meh
*(d::AbstractString, x::Meh) = meh

# Iterators
"""
    Mehs.replace(itr, replacement)

Return an iterator wrapping iterable `itr` which replaces [`meh`](@ref) values with
`replacement`. When applicable, the size of `itr` is preserved.
If the type of `replacement` differs from the element type of `itr`,
it will be converted.

See also: [`Mehs.skip`](@ref), [`Mehs.fail`](@ref)

# Examples
```jldoctest
julia> collect(Mehs.replace([1, meh, 2], 0))
3-element Array{Int64,1}:
 1
 0
 2

julia> collect(Mehs.replace([1 meh; 2 meh], 0))
2×2 Array{Int64,2}:
 1  0
 2  0

```
"""
replace(itr, replacement) = EachReplaceMeh(itr, convert(eltype(itr), replacement))
struct EachReplaceMeh{T, U}
    x::T
    replacement::U
end
Base.iteratorsize(::Type{<:EachReplaceMeh{T}}) where {T} =
    Base.iteratorsize(T)
Base.iteratoreltype(::Type{<:EachReplaceMeh{T}}) where {T} =
    Base.iteratoreltype(T)
Base.length(itr::EachReplaceMeh) = length(itr.x)
Base.size(itr::EachReplaceMeh) = size(itr.x)
Base.start(itr::EachReplaceMeh) = start(itr.x)
Base.done(itr::EachReplaceMeh, state) = done(itr.x, state)
Base.eltype(itr::EachReplaceMeh) = Mehs.T(eltype(itr.x))
@inline function Base.next(itr::EachReplaceMeh, state)
    v, s = next(itr.x, state)
    (v isa Meh ? itr.replacement : v, s)
end

"""
    Mehs.skip(itr)

Return an iterator wrapping iterable `itr` which skips [`meh`](@ref) values.

Use [`collect`](@ref) to obtain an `Array` containing the non-`meh` values in
`itr`. Note that even if `itr` is a multidimensional array, the result will always
be a `Vector` since it is not possible to remove mehs while preserving dimensions
of the input.

See also: [`Mehs.replace`](@ref), [`Mehs.fail`](@ref)

# Examples
```jldoctest
julia> collect(Mehs.skip([1, meh, 2]))
2-element Array{Int64,1}:
 1
 2

julia> collect(Mehs.skip([1 meh; 2 meh]))
2-element Array{Int64,1}:
 1
 2

```
"""
skip(itr) = EachSkipMeh(itr)
struct EachSkipMeh{T}
    x::T
end
Base.iteratorsize(::Type{<:EachSkipMeh}) =
    Base.SizeUnknown()
Base.iteratoreltype(::Type{EachSkipMeh{T}}) where {T} =
    Base.iteratoreltype(T)
Base.eltype(itr::EachSkipMeh) = Mehs.T(eltype(itr.x))
# Fallback implementation for general iterables: we cannot access a value twice,
# so after finding the next non-meh element in start() or next(), we have to
# pass it in the iterator state, which introduces a type instability since the value
# is meh if the input does not contain any non-meh element.
# As of Julia 0.6 and early 0.7, this instability kills performance.
@inline function Base.start(itr::EachSkipMeh)
    s = start(itr.x)
    v = meh
    @inbounds while !done(itr.x, s) && v isa Meh
        v, s = next(itr.x, s)
    end
    (v, s)
end
@inline Base.done(itr::EachSkipMeh, state) = ismeh(state[1]) && done(itr.x, state[2])
@inline function Base.next(itr::EachSkipMeh, state)
    v1, s = state
    v2 = meh
    @inbounds while !done(itr.x, s) && v2 isa Meh
        v2, s = next(itr.x, s)
    end
    (v1, (v2, s))
end
# Optimized implementation for AbstractArray, relying on the ability to access x[i] twice:
# once in done() to find the next non-meh entry, and once in next() to return it.
# This works around the type instability problem of the generic fallback.
@inline function _next_nonmeh_ind(x::AbstractArray, s)
    idx = eachindex(x)
    @inbounds while !done(idx, s)
        i, new_s = next(idx, s)
        x[i] isa Meh || break
        s = new_s
    end
    s
end
@inline Base.start(itr::EachSkipMeh{<:AbstractArray}) =
    _next_nonmeh_ind(itr.x, start(eachindex(itr.x)))
@inline Base.done(itr::EachSkipMeh{<:AbstractArray}, state) =
    done(eachindex(itr.x), state)
@inline function Base.next(itr::EachSkipMeh{<:AbstractArray}, state)
    i, state = next(eachindex(itr.x), state)
    @inbounds v = itr.x[i]::eltype(itr)
    (v, _next_nonmeh_ind(itr.x, state))
end

"""
    Mehs.fail(itr)

Return an iterator wrapping iterable `itr` which will throw a [`MehException`](@ref)
if a [`meh`](@ref) value is found.

Use [`collect`](@ref) to obtain an `Array` containing the resulting values.
If `itr` is an array, the resulting array will have the same dimensions.

See also: [`Mehs.skip`](@ref), [`Mehs.replace`](@ref)

# Examples
```jldoctest
julia> collect(Mehs.fail([1 2; 3 4]))
2×2 Array{Int64,2}:
 1  2
 3  4

julia> collect(Mehs.fail([1, meh, 2]))
ERROR: MehException: meh value encountered by Mehs.fail
[...]
```
"""
fail(itr) = EachFailMeh(itr)
struct EachFailMeh{T}
    x::T
end
Base.iteratorsize(::Type{EachFailMeh{T}}) where {T} =
    Base.iteratorsize(T)
Base.iteratoreltype(::Type{EachFailMeh{T}}) where {T} =
    Base.iteratoreltype(T)
Base.length(itr::EachFailMeh) = length(itr.x)
Base.size(itr::EachFailMeh) = size(itr.x)
Base.start(itr::EachFailMeh) = start(itr.x)
Base.done(itr::EachFailMeh, state) = done(itr.x, state)
Base.eltype(itr::EachFailMeh) = Mehs.T(eltype(itr.x))
@inline function Base.next(itr::EachFailMeh, state)
    v, s = next(itr.x, state)
    # NOTE: v isa Meh currently gives incorrect code, cf. JuliaLang/julia#24177
    ismeh(v) && throw(MehException("meh value encountered by Mehs.fail"))
    (v::eltype(itr), s)
end

"""
    coalesce(x, y...)

Return the first non-`meh` value in the arguments, or `meh` if all arguments are `meh`.

In its broadcasted form, this function can be used to replace all meh values
in an array with a given value (see examples).

# Examples

```jldoctest
julia> coalesce(meh, 1)
1

julia> coalesce(1, meh)
1

julia> coalesce(meh, meh)
meh

julia> coalesce.([meh, 1, meh], 0)
3-element Array{$Int,1}:
 0
 1
 0

julia> coalesce.([meh, 1, meh], [0, 10, 5])
3-element Array{$Int,1}:
 0
 1
 5

```
"""
coalesce(x) = x
coalesce(x, y...) = ifelse(x !== meh, x, coalesce(y...))

"""
    levels(x)

Return a vector of unique values which occur or could occur in collection `x`,
omitting `meh` even if present. Values are returned in the preferred order
for the collection, with the result of [`sort`](@ref) as a default.

Contrary to [`unique`](@ref), this function may return values which do not
actually occur in the data, and does not preserve their order of appearance in `x`.
"""
function levels(x)
    T = Mehs.T(eltype(x))
    levs = convert(AbstractArray{T}, filter!(!ismeh, unique(x)))
    if method_exists(isless, Tuple{T, T})
        try; sort!(levs); end
    end
    levs
end

# AbstractArray{>:Meh} functions

function ==(A::AbstractArray{>:Meh}, B::AbstractArray)
    if indices(A) != indices(B)
        return false
    end
    if isa(A,AbstractRange) != isa(B,AbstractRange)
        return false
    end
    anymeh = false
    @inbounds for (a, b) in zip(A, B)
        eq = (a == b)
        if eq === false
            return false
        else
            anymeh |= ismeh(eq)
        end
    end
    return anymeh ? meh : true
end

==(A::AbstractArray, B::AbstractArray{>:Meh}) = (B == A)
==(A::AbstractArray{>:Meh}, B::AbstractArray{>:Meh}) =
    invoke(==, Tuple{AbstractArray{>:Meh}, AbstractArray}, A, B)

!=(x::AbstractArray{>:Meh}, y::AbstractArray) = !(x == y)
!=(x::AbstractArray, y::AbstractArray{>:Meh}) = !(x == y)
!=(x::AbstractArray{>:Meh}, y::AbstractArray{>:Meh}) = !(x == y)

function Base.any(f, A::AbstractArray{>:Meh})
    anymeh = false
    @inbounds for x in A
        v = f(x)
        if v === true
            return true
        else
            anymeh |= ismeh(v)
        end
    end
    return anymeh ? meh : false
end

function Base.all(f, A::AbstractArray{>:Meh})
    anymeh = false
    @inbounds for x in A
        v = f(x)
        if v === false
            return false
        else
            anymeh |= ismeh(v)
        end
    end
    return anymeh ? meh : true
end

function Base.float(A::AbstractArray{Union{T, Meh}}) where {T}
    U = typeof(float(zero(T)))
    convert(AbstractArray{Union{U, Meh}}, A)
end
Base.float(A::AbstractArray{Meh}) = A

end # module


text:
__precompile__()

"""
$(DocStringExtensions.README)
"""
module ParameterizedFunctions

using DocStringExtensions
  using DataStructures, DiffEqBase, Latexify
  using Reexport
  @reexport using ModelingToolkit
  using ModelingToolkit: Sym, FnType, tosymbol

  import LinearAlgebra

  import Base: getindex

  include("ode_def_opts.jl")
  include("ode_findrep.jl")
  include("utils.jl")
  include("dict_build.jl")
  include("macros.jl")
  include("latexify.jl")

  export @ode_def,ode_def_opts,@ode_def_bare, @ode_def_all
end # module

##### Extra

# Jacobian Factorization

#=
local fsymjac_L
local fsymjac_U
local Jex_L
local Jex_U
try
  # Factorize the Jacobian
  fsymjac = lufact(symjac)
  fsymjac_L = fsymjac[:L]
  fsymjac_U = fsymjac[:U]
  Jex_L,Jex_U = build_fjac_func(fsymjac_L,fsymjac_U,indvar_dict,param_dict,inline_dict)
catch
  fsymjac_L = Matrix{SymEngine.Basic}(0,0)
  fsymjac_U = Matrix{SymEngine.Basic}(0,0)
end

"""
Builds the LU-factorized Jacobian functions
"""
function build_fjac_func(symjac_L,symjac_U,indvar_dict,param_dict,inline_dict)
  # Lower Triangle
  Jex = :()
  for i in 1:size(symjac_L,1)
    for j in 1:i
      ex = parse(string(symjac_L[i,j]))
      if typeof(ex) <: Expr
        ode_findreplace(ex,ex,indvar_dict,param_dict,inline_dict)
      else
        ex = ode_symbol_findreplace(ex,indvar_dict,param_dict,inline_dict)
      end
      push!(Jex.args,:(J[$i,$j] = $ex))
    end
  end
  Jex.head = :block
  push!(Jex.args,nothing)
  Jex_L = :(jac = (t,u,J)->$Jex)

  # Upper Triangle
  Jex = :()
  for j in 1:size(symjac_U,2)
    for i in 1:j
      ex = parse(string(symjac_U[i,j]))
      if typeof(ex) <: Expr
        ode_findreplace(ex,ex,indvar_dict,param_dict,inline_dict)
      else
        ex = ode_symbol_findreplace(ex,indvar_dict,param_dict,inline_dict)
      end
      push!(Jex.args,:(J[$i,$j] = $ex))
    end
  end
  Jex.head = :block
  push!(Jex.args,nothing)
  Jex_U = :(jac = (t,u,J)->$Jex)

  Jex_L,Jex_U
end

=#


text:
#Code written by Skipper af Rosenborg and Esben Bager
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")

#Set path
#Esben's path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexData"

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

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

#Define datasets to load
modelMatrix = zeros(industriesTotal, possibilities)
noDurModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

#Initialize arrays
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)

#Fill arrays with data
XArrays, YArrays = generateXandYs(industries, modelMatrix)

standY = YArrays[1]
nRows = size(standY)[1]

#Initialization of parameters
gamma = 0 #risk aversion
returnPerfectMatrix = zeros(nRows-trainingSize)
weightsPerfect = zeros(nRows-trainingSize, 10)

#Establishing perfect results in order to avoid doing same mean-variance calculation over and over
for t=1:(nRows-trainingSize-2)
    trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)
    valY = zeros(10)
    for i = 1:10
        valY[i] = validationY[i][1]
    end
    rfRatesVec = rfRates[t:(t+trainingSize-1)]
    trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
    weightsPerfect[t, :], returnPerfectMatrix[t] = findPerfectResults(trainX, valY, valY, gamma)
end

#Output results
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"
writedlm(path*"weightsPerfect.csv", weightsPerfect, ",")
writedlm(path*"returnPerfectMatrix.csv", returnPerfectMatrix, ",")

max_stars_repo_name: skipperafrosenborg/Thesis

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

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
# code to printEdges, Nodes and addHybridizationUpdate to avoid steps
# in tests_5taxon.jl
# Claudia Septiembre 2014

#printEdges(net)
#printNodes(net)
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(net);
#printEdges(net)
#printNodes(net)


text:
#
# Propagators for CR3BP orbits.
#

"""
Dynamics for ideal Circular Restricted Three-body numerical integration.
Compatable with `DifferentialEquations` solvers!
"""
function CR3BPTic!(∂u, u, p, t=0)
    ∂u.r =  u.v
    accel!(∂u.v, u.r, u.v, p.μ)
    return nothing
end

"""
Returns dynamics tic for combined Halo iterative solver state vector.

__Arguments:__
- `∂u`: Derivative of state `u`.
- `u`: State vector: `[x, y, z, ẋ, ẏ, ż, Φ₁, Φ₂, Φ₃, Φ₄, Φ₅, Φ₆]`.
- `p`: Parameters (contains non-dimensional mass parameter `μ`, positions `x₁`, `x₂`, and configuration).
- `t`: Time in seconds.

__Outputs:__
- None (sets derivative `∂u` in-place).

__References:__
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/).
"""
function CR3BPSTMTic!(∂u, u, p, t)
    # Cartesian state
    ∂u.r .=  u.v
    accel!(∂u.v, u.r, u.v, p.μ)

    # State transition matrix
    ∂Φ     = state_transition_dynamics(p.μ, u.r) * SMatrix{6,6}(transpose(hcat(u.Φ₁, u.Φ₂, u.Φ₃, u.Φ₄, u.Φ₅, u.Φ₆)))
    ∂u.Φ₁ .= ∂Φ[1,:]
    ∂u.Φ₂ .= ∂Φ[2,:]
    ∂u.Φ₃ .= ∂Φ[3,:]
    ∂u.Φ₄ .= ∂Φ[4,:]
    ∂u.Φ₅ .= ∂Φ[5,:]
    ∂u.Φ₆ .= ∂Φ[6,:]

    return nothing
end

"""
Given a `CartesianOrbit`, return a `DifferentialEquations.ODEProblem` instance
that describes the Restricted Two-body Problem.
"""
function SciMLBase.ODEProblem(orbit::NormalizedSynodicCR3BPOrbit, Δt::Real)

    return ODEProblem(
        AstrodynamicalModels.CR3BPVectorField,
        Vector(orbit.state.rv),
        (epoch(orbit.state), Δt),
        [normalized_mass_parameter(orbit.system)]
    )

end

"""
Given a `CartesianOrbit`, return a `DifferentialEquations.ODEProblem` instance
that describes the Restricted Two-body Problem.
"""
function SciMLBase.ODEProblem(orbit::NormalizedSynodicSTMCR3BPOrbit, Δt::Real)

    return ODEProblem(
        AstrodynamicalModels.CR3BPWithSTMVectorField,
        vcat(orbit.state.cart.rv, vec(orbit.state.stm)),
        (epoch(orbit.state), Δt),
        [normalized_mass_parameter(orbit.system)]
    )

end

"""
Given a `CartesianOrbit`, return a `DifferentialEquations.ODEProblem` instance
that describes the Restricted Two-body Problem.
"""
SciMLBase.ODEProblem(orbit::NormalizedSynodicCR3BPOrbit, Δt::Unitful.Time) = ODEProblem(orbit, ustrip(timeunit(orbit.state), Δt))

"""
Given a `CartesianOrbit`, return a `DifferentialEquations.ODEProblem` instance
that describes the Restricted Two-body Problem.
"""
SciMLBase.ODEProblem(orbit::NormalizedSynodicSTMCR3BPOrbit, Δt::Unitful.Time) = ODEProblem(orbit, ustrip(timeunit(orbit.state), Δt))

"""
Uses `DifferentialEquations` solvers to propagate normalized,
Synodic CR3BP states Δt into the future. All keyword arguments are passed
directly to `DifferentialEquations` solvers.
"""
function propagate(orbit::NormalizedSynodicCR3BPOrbit, Δt::Real; kwargs...)

    # Set default kwargs (modified from [3])
    defaults = (;  reltol=1e-14, abstol=1e-14)
    options = merge(defaults, kwargs)

    # Initial conditions
    problem = ODEProblem(orbit, Δt)

    # Numerically integrate!
    sols = solve(problem; options...)

    # Return Trajectory structure
    return [
        let
            cart = CartesianState(sols.u[i][1:3], sols.u[i][4:6], sols.t[i], Synodic;
                                  lengthunit = NormalizedLengthUnit(), timeunit = NormalizedTimeUnit())
            CircularRestrictedThreeBodyOrbit(cart, orbit.system)
        end
        for i ∈ 1:length(sols.t)
    ]
end

"""
Uses `DifferentialEquations` solvers to propagate normalized,
Synodic CR3BP states Δt into the future. All keyword arguments are passed
directly to `DifferentialEquations` solvers.
"""
function propagate(orbit::NormalizedSynodicSTMCR3BPOrbit, Δt::Real; kwargs...)

    # Set default kwargs (modified from [3])
    defaults = (;  reltol=1e-14, abstol=1e-14)
    options = merge(defaults, kwargs)

    # Initial conditions
    problem = ODEProblem(orbit, Δt)

    # Numerically integrate!
    sols = solve(problem; options...)

    # Return Trajectory structure
    return [
        let
            cart = CartesianState(sols.u[i][1:3], sols.u[i][4:6], sols.t[i], Synodic;
                                  lengthunit = NormalizedLengthUnit(), timeunit = NormalizedTimeUnit())
            state = SynodicCartesianSTMState(cart, SMatrix{6,6}(sols.u[i][7:end]...))
            CircularRestrictedThreeBodyOrbit(state, orbit.system)
        end
        for i ∈ 1:length(sols.t)
    ]
end

"""
Returns the partial derivative matrix of potential `U`.

__Arguments:__
- `μ`: Non-dimensional mass parameter for the CR3BP system.
- `r`: Non-dimensional position vector for the spacecraft.

__Outputs:__
- Partial derivative matrix of potential `U`.

__References:__
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/)
"""
potential_energy_hessian = include("PotentialEnergyHessian.jl")

"""
Returns the derivative mapping of CR3BP state transition matrix, `F`.

__Arguments:__
- `μ`: Non-dimensional mass parameter for the CR3BP system.
- `r`: Non-dimensional position vector for the spacecraft.

__Outputs:__
- Linear mapping from Φ to Φ̇, `F`.

__References:__
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/)
"""
function state_transition_dynamics(μ, r)

    return SMatrix{6,6}(vcat(
        hcat(zeros((3,3)), I(3)),
        hcat(potential_energy_hessian(r..., μ), [0 2 0; -2 0 0; 0 0 0])
    ))

end


text:
using Documenter
using ThreadedScans

makedocs(
    sitename = "ThreadedScans",
    format = Documenter.HTML(),
    modules = [ThreadedScans]
)

deploydocs(
    repo = "github.com/JuliaFolds/ThreadedScans.jl",
    push_preview = true,
)


text:
using Documenter, InteractiveFiction

makedocs(;
    modules=[InteractiveFiction],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/PPKFS/InteractiveFiction.jl/blob/{commit}{path}#L{line}",
    sitename="InteractiveFiction.jl",
    authors="Avery",
    assets=String[],
)


text:
include("../utils.jl")

@testset "Group wrapped in ValidationManifold" begin
    G = SpecialOrthogonal(3)
    M = Rotations(3)
    AG = ValidationManifold(G)
    @test is_group_manifold(AG)

    eg = Matrix{Float64}(I, 3, 3)
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0]]
    p, q = [exp(M, eg, hat(M, eg, ωi)) for ωi in ω]
    X = hat(M, eg, [-1.0, 2.0, 0.5])

    e = Identity(G)
    p2, q2 = ValidationMPoint(p), ValidationMPoint(q)
    @test q2 === Manifolds.array_point(q2) # test that double wraps are avoided.
    X2 = ValidationTVector(X)

    @test Identity(G) isa Identity

    eg = allocate(p)
    identity_element!(G, eg)
    eg2 = allocate(p2)
    identity_element!(AG, eg2)
    @test eg2.value == eg
    eg = allocate(p2)

    @test inv(AG, p2) isa ValidationMPoint
    @test isapprox(G, inv(AG, p2).value, inv(G, p))
    @test inv(AG, e) isa ValidationMPoint{<:Identity}
    @test inv(AG, Identity(AG)) == inv(AG, e)

    pinvq = allocate(p2)
    inv!(AG, pinvq, p2)
    @test isapprox(G, pinvq.value, inv(G, p))
    eg = allocate(p2)
    inv!(AG, eg, e)
    @test isapprox(G, eg.value, Manifolds.array_value(e))

    @test compose(AG, p2, q2) isa ValidationMPoint
    @test isapprox(G, compose(AG, p2, q2).value, compose(G, p, q))
    @test compose(AG, p2, e) === p2
    @test compose(AG, e, p2) === p2

    pq = allocate(p2)
    compose!(AG, pq, p2, q2)
    @test isapprox(G, pq.value, compose(G, p, q))
    pq = allocate(p2)
    compose!(AG, pq, e, q2)
    @test isapprox(G, pq.value, compose(G, e, q))
    pq = allocate(p2)
    compose!(AG, pq, p2, e)
    @test isapprox(G, pq.value, compose(G, p, e))

    @test exp_lie(AG, X2) isa ValidationMPoint
    @test isapprox(G, exp_lie(AG, X2).value, exp_lie(G, X))
    expX = allocate(p2)
    exp_lie!(AG, expX, X2)
    @test isapprox(G, expX.value, exp_lie(G, X))

    @test log_lie(AG, p2) isa ValidationTVector
    @test isapprox(G, e, log_lie(AG, p2).value, log_lie(G, p))
    logp = allocate(X2)
    log_lie!(AG, logp, p2)
    @test isapprox(G, e, logp.value, log_lie(G, p))

    @test lie_bracket(AG, X, X2) isa ValidationTVector
    Xlb = allocate(X2)
    lie_bracket!(AG, Xlb, X, X2)
    @test isapprox(G, e, Xlb.value, lie_bracket(G, X, X2.value))

    @test adjoint_action(AG, p, X) isa ValidationTVector
    Xaa = allocate(X2)
    adjoint_action!(AG, Xaa, p, X2)
    @test isapprox(G, e, Xaa.value, adjoint_action(G, p, X2.value))

    for conv in (LeftAction(), RightAction())
        @test translate(AG, p2, q2, conv) isa ValidationMPoint
        @test isapprox(G, translate(AG, p2, q2, conv).value, translate(G, p, q, conv))

        pq = allocate(p2)
        translate!(AG, pq, p2, q2, conv)
        @test isapprox(G, pq.value, translate(G, p, q, conv))

        @test inverse_translate(AG, p2, q2, conv) isa ValidationMPoint
        @test isapprox(
            G,
            inverse_translate(AG, p2, q2, conv).value,
            inverse_translate(G, p, q, conv),
        )

        pinvq = allocate(p2)
        inverse_translate!(AG, pinvq, p2, q2, conv)
        @test isapprox(G, pinvq.value, inverse_translate(G, p, q, conv))
    end

    for conv in (LeftAction(), RightAction())
        @test translate_diff(AG, q2, p2, X2, conv; atol=1e-10) isa ValidationTVector
        @test isapprox(
            G,
            translate_diff(AG, q2, p2, X2, conv; atol=1e-10).value,
            translate_diff(G, q, p, X, conv),
        )

        Y = allocate(X2)
        translate_diff!(AG, Y, q2, p2, X2, conv; atol=1e-10)
        @test isapprox(Y.value, translate_diff(G, q, p, X, conv))

        @test inverse_translate_diff(AG, q2, p2, X2, conv; atol=1e-10) isa ValidationTVector
        @test isapprox(
            G,
            inverse_translate_diff(AG, q2, p2, X2, conv; atol=1e-10).value,
            inverse_translate_diff(G, q, p, X, conv),
        )

        Y = allocate(X2)
        inverse_translate_diff!(AG, Y, q2, p2, X2, conv; atol=1e-10)
        @test isapprox(Y.value, inverse_translate_diff(G, q, p, X, conv))
    end
end


text:
### Quadrado: área e perímetro
# Faça um programa em que receba (entrada de dados)
# o valor correspondente ao lado de um quadrado,
# calcule e imprima (saída de dados) seu perímetro e sua área.

println("informe o valor do lado do quadrado")
square = parse(UInt8, readline())

println("Perímetro = ", (square*4) )
println("Área = ", (square^2) )

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
#=
BFSParentsTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#
using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BFSParents.jl")
matrix =  GBMatrix(sparse([[0, 0, 0, 1, 0, 0, 0] [1, 0, 0, 0, 0, 0, 0] [0, 0, 0, 1, 0, 1, 1] [1, 0, 0, 0, 0, 0, 1] [0, 1, 0, 0, 0, 0, 1] [0, 0, 1, 0, 1, 0, 0] [0, 1, 0, 0, 0, 0, 0]]))
bfs_parents(matrix, 1, 7)
@test bfs_parents(matrix, 1, 7) == GBVector([0, 1, 4, 1, 2, 3, 2])


text:
module DataStorage

## Imports
import Base: size #must import to add a definition to size

##Exports
export DataContainer, PairedDataContainer
export size
export get_data, get_inputs, get_outputs

## Objects
"""
    DataContainer{FT <: Real}

Container to store data samples as columns in an array.
"""
struct DataContainer{FT <: Real}
    #stored data, each piece of data is a column [data dimension × number samples]
    stored_data::Array{FT, 2}
    #constructor with 2D arrays
    function DataContainer(stored_data::Array{FT, 2}; data_are_columns = true) where {FT <: Real}

        if data_are_columns
            new{FT}(stored_data)
        else
            new{FT}(permutedims(stored_data, (2, 1)))
        end
    end
end

"""
    PairedDataContainer{FT <: Real}

stores input - output pairs as data containers, there must be an equal number of inputs and outputs
"""
struct PairedDataContainer{FT <: Real}
    # container for inputs and ouputs, each Container holds an array
    # size [data/parameter dimension × number samples]
    inputs::DataContainer{FT}
    outputs::DataContainer{FT}

    #constructor with 2D Arrays
    function PairedDataContainer(
        inputs::Array{FT, 2},
        outputs::Array{FT, 2};
        data_are_columns = true,
    ) where {FT <: Real}

        sample_dim = data_are_columns ? 2 : 1
        if !(size(inputs, sample_dim) == size(outputs, sample_dim))
            throw(DimensionMismatch("There must be the same number of samples of both inputs and outputs"))
        end

        stored_inputs = DataContainer(inputs; data_are_columns = data_are_columns)
        stored_outputs = DataContainer(outputs; data_are_columns = data_are_columns)
        new{FT}(stored_inputs, stored_outputs)

    end
    #constructor with DataContainers
    function PairedDataContainer(inputs::DataContainer, outputs::DataContainer)

        if !(size(inputs, 2) == size(outputs, 2))
            throw(DimensionMismatch("There must be the same number of samples of both inputs and outputs."))
        else
            FT = eltype(get_data(inputs))
            new{FT}(inputs, outputs)
        end
    end

end

## functions
"""
    size(dc::DataContainer, idx::IT) where {IT <: Integer}

Returns the size of the stored data. If `idx` provided, it returns the size along dimension `idx`.
"""
size(dc::DataContainer) = size(dc.stored_data)

size(dc::DataContainer, idx::IT) where {IT <: Integer} = size(dc.stored_data, idx)

size(pdc::PairedDataContainer) = size(pdc.inputs), size(pdc.outputs)

"""
    size(pdc::PairedDataContainer, idx::IT) where {IT <: Integer}

Returns the sizes of the inputs and ouputs along dimension `idx` (if provided)
"""
size(pdc::PairedDataContainer, idx::IT) where {IT <: Integer} = size(pdc.inputs, idx), size(pdc.outputs, idx)

get_data(dc::DataContainer) = deepcopy(dc.stored_data)

get_data(pdc::PairedDataContainer) = get_inputs(pdc), get_outputs(pdc)

get_inputs(pdc::PairedDataContainer) = get_data(pdc.inputs)

get_outputs(pdc::PairedDataContainer) = get_data(pdc.outputs)

end # module


text:
using Printf
using Random

using GLFW
using ModernGL
using CSyntax
using CSyntax.CStatic

using CImGui
using CImGui.LibCImGui
using CImGui.GLFWBackend
using CImGui.OpenGLBackend

using ImPlot
import CImGui.LibCImGui: ImGuiCond_Always

using Formatting
import DataStructures.CircularBuffer

#using .World

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    tmp_glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    tmp_glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    GLFW.WindowHint(GLFW.FOCUSED,GL_TRUE)
end
const glsl_version=tmp_glsl_version

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

struct GlyphRange
	first::UInt16
	last::UInt16
	zero::UInt16
	function GlyphRange(first::UInt16,last::UInt16)
		new(first,last,0x0)
	end
end
greek=GlyphRange(0x0300, 0x03ff)

font_config=ImFontConfig_ImFontConfig()  #not used
#let offset=0
#	for i in eachindex(Base.datatype_fieldtypes(ImFontConfig))
#		fieldnames(ImFontConfig)[i] == :MergeMode ? break : nothing
#		offset+=sizeof(Base.datatype_fieldtypes(ImFontConfig)[i])
#	end
#	unsafe_store!(Ptr{Bool}(font_config+offset),true,1)
#end

mutable struct Display
    isOpen::Bool
    window::GLFW.Window
    ctx::Ptr{ImGuiContext}
    ctxp::Ptr{ImPlot.LibCImPlot.ImPlotContext}
    ranges::Ptr{ImVector_ImWchar}
    builder::Ptr{ImFontGlyphRangesBuilder}
    config::Ptr{ImFontConfig}
    img_width::Int
    img_height::Int
    image_id::Int
    world_image::Array{GLubyte,3}
    image_cache::Array{Int,2}
    x_values::CircularBuffer{Int}
    y_org::CircularBuffer{Int}
    y_res::CircularBuffer{Int}
    filename_buffer::Vector{UInt8}
    runTime::Float64
    showLoad::Bool
    showReset::Bool
end

mutable struct Action
    stop::Bool
    exit::Bool
    save::Bool
    load::Bool
    reset::Bool
    function Action()
        new(false,false,false,false,false)
    end
    function Action(stop::Bool)
        new(stop,false,false,false,false)
    end
end

function display_initialize()
    # create window
    window = GLFW.CreateWindow(1100, 720, "PlayGround")

    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()
    ctxp = ImPlot.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsClassic()

    #fonts_dir = joinpath(pathof(CImGui), "..","..","fonts")
    fonts_dir = raw"c:\Windows\Fonts";
    fonts = CImGui.GetIO().Fonts

    builder=ImFontGlyphRangesBuilder_ImFontGlyphRangesBuilder()
    ImFontGlyphRangesBuilder_AddRanges(builder, ImFontAtlas_GetGlyphRangesDefault(fonts))
    greek_ptr=Core.Ptr{ImWchar}(Base.pointer_from_objref(Ref(greek)))
    ImFontGlyphRangesBuilder_AddRanges(builder, greek_ptr)

    ranges=ImVector_ImWchar_create()
    ImVector_ImWchar_Init(ranges)
    ImFontGlyphRangesBuilder_BuildRanges(builder, ranges)

    r = unsafe_wrap(Vector{ImVector_ImWchar}, ranges, 1)
    CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "arial.ttf"), 16, C_NULL, r[1].Data )

    #CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "arial.ttf"), 16, fonts_config, ImFontAtlas_GetGlyphRangesDefault(fonts) )
    #greek=GlyphRange(0x0300, 0x03ff)
    #greek_ptr=Core.Ptr{UInt16}(Base.pointer_from_objref(Ref(greek)))
    #CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "arial.ttf"), 16, fonts_config, greek_ptr )

    CImGui.Build(fonts)
    default_font = CImGui.AddFontDefault(fonts)
    #@assert default_font != C_NULL

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)

    img_width=400
    img_height=400
    (world_image,image_cache) = new_display_image(img_width, img_height)
    image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)

    y_org=CircularBuffer{Int}(1000)
    y_res=CircularBuffer{Int}(1000)
    x_values=CircularBuffer{Int}(1000)
    fill!(y_org,0)
    #fill!(y_res,0)
    for i in 1:1000
        push!(x_values,i)
    end
    default_filename="ws.serialize"
    filename_buffer=resize!(UInt8.(collect(default_filename)),100)
    filename_buffer[length(default_filename)+1]=UInt8('\0')
    return Display(true,window,ctx,ctxp,ranges,builder,font_config,img_width,img_height,image_id,world_image,image_cache,
        x_values,y_org,y_res,
        filename_buffer,
        0,
        true,true
    )
end

function new_display_image(img_width, img_height)
    display_image=zeros(GLubyte, 4, img_width, img_height)
    display_image[4,:,:].=255
    image_cache=zeros(Int,img_width, img_height)
    return (display_image,image_cache)
end

function is_visible(px,py,img_width,img_height)
    return px>0 && px<=img_width && py>0 && py<=img_height
end

#function display_update(display::Display,ws::World.State,action::Action)
function display_update(display::Display,action::Action,ws)
    try
        clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
        if display.isOpen
            if !GLFW.WindowShouldClose(display.window)
                GLFW.PollEvents()

                # start the Dear ImGui frame
                ImGui_ImplOpenGL3_NewFrame()
                ImGui_ImplGlfw_NewFrame()
                CImGui.NewFrame()
                CImGui.SetNextWindowSize(CImGui.ImVec2(350.0,270.0),CImGui.ImGuiCond_Once);

                #CImGui.Begin("Special strings")
                #CImGui.Text("ÄÜÖ π ϕ θ ♐ ☽ ")
                #CImGui.End()

                CImGui.Begin("Status")

                CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))
                CImGui.Text(@sprintf("Time running (s): %.3f",display.runTime))
                CImGui.Text(@sprintf("Current time: %.20e",World.get(ws,World.currentTime)))
                CImGui.Text(@sprintf("Current steps: %i of %i",World.get_step(ws),World.get_maxsteps(ws)))
                CImGui.Text(@sprintf("Steps per second: %.5f", Float64(World.get_step(ws))/display.runTime ))
                CImGui.BeginChild("Org",CImGui.ImVec2(0, CImGui.GetTextLineHeightWithSpacing() * 2) )
                CImGui.Text(@sprintf("Ancestor count: %i",World.get(ws,org_ancestor_count)))
                CImGui.Text(@sprintf("Max generations: %i",World.get(ws,max_generations)))
                CImGui.EndChild()
                if CImGui.IsItemHovered()
                    CImGui.BeginTooltip()
                    if World.get(ws,max_generations) > 0
                        CImGui.Text(@sprintf("trigger_split_expression: %s",ws.cur_max_ancestor[1].trigger_split_expression))
                        CImGui.Text(@sprintf("trigger_split_expression parameters: %s",ws.cur_max_ancestor[1].trigger_split_expression_parameter))
                        CImGui.Text(@sprintf("housekeeping_expression: %s",ws.cur_max_ancestor[1].housekeeping_expression))
                        CImGui.Text(@sprintf("housekeeping_expression parameters: %s",ws.cur_max_ancestor[1].housekeeping_expression_parameter))
                        CImGui.Text(@sprintf("energy_expression: %s",ws.cur_max_ancestor[1].energy_expression))
                        CImGui.Text(@sprintf("energy_expression parameters: %s",ws.cur_max_ancestor[1].energy_expression_parameter))
                        CImGui.Text(@sprintf("ϕ expression: %s",ws.cur_max_ancestor[1].ϕ_expression))
                        CImGui.Text(@sprintf("ϕ expression parameters: %s",ws.cur_max_ancestor[1].ϕ_expression_parameter))
                        CImGui.Text(@sprintf("θ expression: %s",ws.cur_max_ancestor[1].θ_expression))
                        CImGui.Text(@sprintf("θ expression parameters: %s",ws.cur_max_ancestor[1].θ_expression_parameter))
                    end
                    CImGui.EndTooltip()
                end

                CImGui.Text(@sprintf("Free energy: %i",World.get(ws,World.freeEnergy)))
                CImGui.Text(@sprintf("# Resources: %i",length(ws.resources)))
                CImGui.Text(@sprintf("# organisms: %i",length(ws.organisms)))
                CImGui.Text(@sprintf("Max life time: %.20e",World.get(ws,World.maxLifeTime)))
                CImGui.End()

                CImGui.Begin("Action")
                if !action.stop
                    CImGui.Button("Stop") && (action.stop = true)
                else
                    CImGui.Button("Run") && (action.stop = false)
                end
                CImGui.SameLine()
                CImGui.Button("Exit") && (action.exit = true)
                if display.showReset
                    CImGui.SameLine()
                    CImGui.Button("Reset") && (action.reset = true)
                end
                CImGui.Button("Save") && (action.save = true)
                if display.showLoad
                    CImGui.SameLine()
                    CImGui.Button("Load") && (action.load = true)
                end

                CImGui.SameLine()
                CImGui.InputText("",display.filename_buffer,length(display.filename_buffer)-1)
                CImGui.End()

                CImGui.Begin("Parameter")
                delta=Cfloat[get(ws,randomizeDelta)]
                #CImGui.DragFloat("Randomize delta:",pointer(delta),0.01)
                CImGui.SliderFloat("Randomize delta:",pointer(delta),0.0,1.0)
                set!(ws,randomizeDelta,Float64(delta[1]))
                CImGui.End()

                CImGui.Begin("Print")
                CImGui.BeginChild("Child")

                size=CImGui.GetWindowSize()
                cur_img_width=floor(Int,size.x)
                cur_img_height=floor(Int,size.y)
                if cur_img_width != display.img_width || cur_img_height != display.img_height
                    display.img_width=cur_img_width
                    display.img_height=cur_img_height
                    (display.world_image,display.image_cache) = new_display_image(display.img_width, display.img_height)
                    ImGui_ImplOpenGL3_DestroyImageTexture(display.image_id)
                    display.image_id = ImGui_ImplOpenGL3_CreateImageTexture(display.img_width, display.img_height)
                    #empty!(ws.helper.res_index_positions2resource)
                    draw_all(display,ws)
                    #wipe_resources_incremental(display,ws)
                    #wipe_org_incremental(display,ws)
                    #draw_resource_incremental(display,ws)
                    #draw_org_incremental(display,ws)
                end
                ImGui_ImplOpenGL3_UpdateImageTexture(display.image_id, display.world_image, display.img_width, display.img_height)
                CImGui.Image(Ptr{Cvoid}(display.image_id), (display.img_width, display.img_height))

                io = CImGui.GetIO()
                pos = CImGui.GetCursorScreenPos()
                rel_mouse_x=io.MousePos.x-pos.x+1
                rel_mouse_y=display.img_height-(pos.y-io.MousePos.y-5)
                if CImGui.IsItemHovered()
                    CImGui.BeginTooltip()
                    CImGui.Text(@sprintf("Image: (%.2f, %.2f)", display.img_width, display.img_height))
                    CImGui.Text(@sprintf("Region: (%.2f, %.2f)", pos.x, pos.y))
                    CImGui.Text(@sprintf("Mouse: (%.2f, %.2f)",rel_mouse_x,rel_mouse_y))
                    CImGui.Text(@sprintf("pos: %i %i",floor(Int,rel_mouse_x),floor(Int,rel_mouse_y)))
                    CImGui.Text(@sprintf("Resources:"))
                    c=World.Coordinates2dIndex(floor(Int,rel_mouse_x),floor(Int,rel_mouse_y))
                    area_size=5
                    for cx in -area_size:area_size
                        for cy in -area_size:area_size
                            ctmp=World.Coordinates2dIndex(c.x+cx,c.y+cy)
                            #if haskey(ws.helper.res_index_positions2resource,ctmp)
                            #    for res in ws.helper.res_index_positions2resource[ctmp]
                            #        cr=World.polar_to_2d_index(res.position,display.img_width, display.img_height)
                            #        CImGui.Text(@sprintf("Resource: %i %i",cr.x,cr.y))
                            #        for restype in res.resources
                            #            CImGui.Text(@sprintf("  ResourceType: %s %i %i",String(Symbol(restype.name)),restype.energy,restype.free_energy))
                            #        end
                            #    end
                            #end
                        end
                    end
                    CImGui.EndTooltip()
                end

                CImGui.EndChild()
                CImGui.End()

                CImGui.Begin("Plot")
                size=CImGui.GetWindowSize()
                cur_img_width=floor(Int,size.x)
                cur_img_height=floor(Int,size.y)
                if length(display.y_res) == 0
                    fill!(display.y_res,length(ws.resources))
                end
                if !action.stop
                    push!(display.y_org,length(ws.organisms))
                    push!(display.y_res,length(ws.resources))
                    push!(display.x_values,display.x_values[1000]+1)
                end
                if length(display.y_org) > 0
                    x_min_index=1
                    x_max_index=1000
                    y_org_max = maximum(display.y_org)
                    y_res_min = Float64(minimum(display.y_res))
                    y_res_min -= 0.01*y_res_min
                    y_res_max = Float64(maximum(display.y_res))
                    y_res_max += 0.01*y_res_max

                    nlabels=5
                    fmt = "%10.0f"
                    fmtrfunc = generate_formatter( fmt )

                    yticks=collect(range(0,y_org_max,length=nlabels))
                    #ylabels=join.(setindex!.(collect.(fmtrfunc.(yticks)),'.',1))
                    #ylabels=replace.(fmtrfunc.(yticks)," " => "_")
                    ylabels=fmtrfunc.(yticks)
                    ylabels[1]=replace(ylabels[1]," " => "_")

                    ImPlot.SetNextPlotTicksY(yticks,nlabels;labels=ylabels,show_default=false,y_axis=0)
                    ImPlot.SetNextPlotLimits(display.x_values[1], display.x_values[1000], 0.0, y_org_max, ImGuiCond_Always)
                    if ImPlot.BeginPlot("##line1", C_NULL, "org", CImGui.ImVec2(-1,cur_img_height/2-20);
                        x_flags = ImPlotAxisFlags_NoTickLabels
                        )
                        ImPlot.PlotLine(display.x_values,display.y_org)
                        ImPlot.EndPlot()
                    end

                    yticks=collect(range(y_res_min,y_res_max,length=nlabels))
                    #ylabels=join.(setindex!.(collect.(fmtrfunc.(yticks)),'.',1))
                    #ylabels=replace.(fmtrfunc.(yticks)," " => "_")
                    ylabels=fmtrfunc.(yticks)
                    ylabels[1]=replace(ylabels[1]," " => "_")

                    ImPlot.SetNextPlotTicksY(yticks,nlabels;labels=ylabels,show_default=false,y_axis=0)
                    ImPlot.SetNextPlotLimits(display.x_values[1], display.x_values[1000], y_res_min, y_res_max, ImGuiCond_Always)
                    if ImPlot.BeginPlot("##line2", "steps", "res", CImGui.ImVec2(-1,cur_img_height/2-20) )
                        ImPlot.PlotLine(display.x_values,display.y_res)
                        ImPlot.EndPlot()
                    end
                end

                #y1 = rand(1000)
                #y2 = rand(1000).+2.0
                #ImPlot.SetNextPlotLimits(0.0,1000,0.0,4.0,ImGuiCond_Always)
                #ImPlot.SetNextPlotLimitsY(0.0, 1.0, ImGuiCond_Always,1)
                #ImPlot.SetNextPlotLimitsY(0.0, 4.0, ImGuiCond_Always,2)
                #if ImPlot.BeginPlot("##line", "x", "y", CImGui.ImVec2(-1,-1); flags = ImPlot.ImPlotFlags_YAxis2 )
                #    ImPlot.PlotLine(y1)
                #    ImPlot.SetPlotYAxis(2)
                #    ImPlot.PlotLine(y2)
                #    ImPlot.EndPlot()
                #end

                CImGui.End()

                # rendering
                CImGui.Render()
                GLFW.MakeContextCurrent(display.window)
                display_w, display_h = GLFW.GetFramebufferSize(display.window)
                glViewport(0, 0, display_w, display_h)
                glClearColor(clear_color...)
                glClear(GL_COLOR_BUFFER_BIT)
                ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

                GLFW.MakeContextCurrent(display.window)
                GLFW.SwapBuffers(display.window)
            else
                display=display_close(display)
            end
        end
    catch e
        @error "Error in renderloop!" exception=e
        Base.show_backtrace(stderr, catch_backtrace())
        display_close(display)
    end
    return display
end

function draw_resource(display::Display,res::World.Resource)
    c=World.polar_to_2d_index(res.position,display.img_width,display.img_height)
    color=[0,0,0,255]
    if res.type.name == res_type_a
        color[1]=255
    elseif res.type.name == res_type_b
        color[2]=255
    end
    if res.free_energy>0
        xindex=c.x
        yindex=c.y
        if xindex>=0 && xindex<=display.img_width && yindex>=0 && yindex<=display.img_height
            for rgb in 1:4
                display.world_image[rgb,xindex,yindex]=color[rgb]
            end
            display.image_cache[xindex,yindex]+=1
        end
    end
end

function wipe_resource(display::Display,res::World.Resource)
    c=World.polar_to_2d_index(res.position,display.img_width,display.img_height)
    xindex=c.x
    yindex=c.y
    color=[0,0,0,255]
    if xindex>=0 && xindex<=display.img_width && yindex>=0 && yindex<=display.img_height
        if display.image_cache[xindex,yindex] == 1
            for rgb in 1:4
                display.world_image[rgb,xindex,yindex]=color[rgb]
            end
        end
        display.image_cache[xindex,yindex]-=1
        if display.image_cache[xindex,yindex]<0
            println("res: ",display.image_cache[xindex,yindex])
        end
    end
end

#function wipe_resources_incremental(display::Display,ws::World.State)
#    for res in ws.helper.resource2wipe
#        wipe_resources(display,res)
#    end
#    empty!(ws.helper.resource2wipe)
#end

#function wipe_org_incremental(display::Display,ws::World.State)
#    for org in ws.organisms
#        if org.needs_wipe
#            wipe_org(display,org)
#        end
#    end
#end

function wipe_org(display::Display,org::World.Organism)
    cDel=World.polar_to_2d_index(org.current_position,display.img_width,display.img_height)
    xindexDel=cDel.x
    yindexDel=cDel.y
    colorDel=[0,0,0,255]
    for rx in -3:3
        for ry in [-3,3]
            pxDel=xindexDel+rx
            pyDel=yindexDel+ry
            if is_visible(pxDel,pyDel,display.img_width,display.img_height)
                if display.image_cache[pxDel,pyDel] == 1
                    for rgb in 1:4
                        display.world_image[rgb,pxDel,pyDel]=colorDel[rgb]
                    end
                end
                display.image_cache[pxDel,pyDel]-=1
                if display.image_cache[pxDel,pyDel]<0
                    println("org: ",display.image_cache[pxDel,pyDel])
                end
            end
        end
    end
    for rx in [-3,3]
        for ry in -2:2
            pxDel=xindexDel+rx
            pyDel=yindexDel+ry
            if is_visible(pxDel,pyDel,display.img_width,display.img_height)
                if display.image_cache[pxDel,pyDel] == 1
                    for rgb in 1:4
                        display.world_image[rgb,pxDel,pyDel]=colorDel[rgb]
                    end
                end
                display.image_cache[pxDel,pyDel]-=1
                if display.image_cache[pxDel,pyDel]<0
                    println("org: ",display.image_cache[pxDel,pyDel])
                end
            end
        end
    end
end

function draw_org(display::Display,org::World.Organism)
    c=World.polar_to_2d_index(org.current_position,display.img_width,display.img_height)
    xindex=c.x
    yindex=c.y
    color=[255,255,255,255]
    if org.is_sibling
        color=[0,255,0,255]
    end
    for rx in -3:3
        for ry in [-3,3]
            px=xindex+rx
            py=yindex+ry
            if is_visible(px,py,display.img_width,display.img_height)
                for rgb in 1:4
                    display.world_image[rgb,px,py]=color[rgb]
                end
                display.image_cache[px,py]+=1
            end
        end
    end
    for rx in [-3,3]
        for ry in -2:2
            px=xindex+rx
            py=yindex+ry
            if is_visible(px,py,display.img_width,display.img_height)
                for rgb in 1:4
                    display.world_image[rgb,px,py]=color[rgb]
                end
                display.image_cache[px,py]+=1
            end
        end
    end
end

function wipe_all(display::Display,ws)
    colorDel=[0,0,0,255]
    for xindex in 1:1:display.img_width
        for yindex in 1:1:display.img_height
            display.image_cache[xindex,yindex] = 0
            for rgb in 1:4
                display.world_image[rgb,xindex,yindex]=colorDel[rgb]
            end
        end
    end
end

function draw_all(display::Display,ws)
    wipe_all(display,ws)
    for res in ws.resources
        #draw_resource(display,ws,res)
        draw_resource(display,res)
    end
    #empty!(ws.helper.resource2draw)
    #empty!(ws.helper.resource2wipe)
    for org in ws.organisms
        #c=World.polar_to_2d_index(org.current_position,display.img_width,display.img_height)
        #if haskey(ws.helper.org_index_positions2organism,c)
        #    ws.helper.org_index_positions2organism[c]=vcat(ws.helper.org_index_positions2organism[c],org)
        #else
        #    ws.helper.org_index_positions2organism[c]=[org]
        #end
        #if org.energy>0
        draw_org(display,org)
        #end
        #org.last_displayed_position=org.current_position
        #org.needs_draw=false
        #org.needs_wipe=false
    end
end

#function draw_resource_incremental(display::Display,ws::World.State)
#    for res in ws.helper.resource2draw
#        draw_resource(display,ws,res)
#    end
#    empty!(ws.helper.resource2draw)
#end

#function draw_org_incremental(display::Display,ws::World.State)
#    for org in ws.organisms
#        if org.needs_draw
#            draw_org(display,org)
#        end
#    end
#end

function display_close(display::Display)
    if display.isOpen
        display.isOpen = false
        ImGui_ImplOpenGL3_Shutdown()
        ImGui_ImplGlfw_Shutdown()
        ImVector_ImWchar_destroy(display.ranges)
        ImFontGlyphRangesBuilder_destroy(display.builder)
        ImFontConfig_destroy(display.config)
        CImGui.DestroyContext(display.ctx)
        ImPlot.DestroyContext(display.ctxp)
        GLFW.DestroyWindow(display.window)
    end
    return display
end


text:
module ApproximateComputations

include("FittingFunctionApproximation.jl")
include("ASTReplacementApproximation.jl")
include("LoopPerforation.jl")
include("Memoisation.jl")

### Fitting functions exports
# Exporting the approximation wrapper type
export Approximation, Get

# Exporting function generation results types
export GeneratedFunctionType, ErrorResultsContainer

# Exporting getters for function result type
export GetAbsoluteError, GetMeanDifference, GetMedianDifference, GetMedianBenchmarkTime, GetMeanBenchmarkTime

# Exporting function selection functions
export GetFastestAcceptable, FilterFunctionList, GetFunctionName

# Export function visualisation plotting functions
export PlotApproximationFunctionResults, PlotApproximationFunctionDiff, PlotApproximationFunctionDiffHist, PlotMedianError, PlotMedianRuntime, PlotRuntimeErrorPair

# Export main function for generating replacement functions
export GenerateAllApproximationFunctions

# Exporting Loop Perforation function and helpers
export LoopPerforation
export OnlyEvenIterations, OnlyOddIterations, ClipFrontAndBack, OnlyFirstHalf, OnlySecondHalf

# Exporting Memoisation
export TrendingMemoisation, ApproximateHashingMemoise

### AST Replacement Exports
export Operator, Variable, TreeMember, ResetGlobalID, GetGlobalID
export EmulateTree, FullUnwrap, UnwrapTree, WrapTree, ReplaceSubTree, GetAllTrees, UpdateEnvironmentForFunction
export printtree, SetSymbolValue, ClearSymbolDict
export GetAllLeaves, GetAllSymbols, GetAllSymbolsList
export GetOverrides, GetConstructionFunction, @BuildOverrideFromArray
export ReplaceConstantsWithVariables, ReplaceTypeOfSpecifiedVariable, ReplaceAllVariablesOfType
export InArray, GetOperatorIDs, GetOperators
export GetResultForID, SetResultForID
export GetErrorInTree
export GetSubTree
export ToString
export PlotASTError
export TreeToFunction, TreeToFunctionLeaf
export PrintTreeComparisonError, TreeComparison
end


text:
@coeff function K_n(reactor::Reactor)
  K_G / K_beta * sqrt(P_F/K_F)
end


text:
"""
prototyping tools to help with boxy & friends
"""

using Images, FileIO           # load elevation from PNG
using Interpolations
using Gadfly, Colors
using Optim
using Test

include("/home/pvt/workspace/libraries/Caesar.jl/examples/dev/scalar/Sequences.jl")

# some Gadfly stuff
Gadfly.set_default_plot_size(35cm,25cm)
st=style( major_label_font="CMU Serif",
          minor_label_font="CMU Serif",
          major_label_font_size=20pt,
          minor_label_font_size=20pt)

img = load("/home/pvt/workspace/libraries/Caesar.jl/examples/dev/scalar/dem.png")
img = Float64.(img)
terrain = 1e3*Float64.(@view img[300,:])

## 0. parameters
n = 150
sigma_y = 0.1

## 1. test w/ globally-referenced sequences

s1 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1:n]+sigma_y*randn(n))
s2 = Sequences.MeasurementSequence(collect(LinRange(5,20,n)), terrain[1+Int(n/3):1+Int(n/3)+n-1]+sigma_y*randn(n))

Sequences.ssd(s1, s2) # should be low

p = plot(st)
push!(p, layer(x=s1.x, y=s1.y,Geom.line,color=[colorant"orange"]))
push!(p, layer(x=s2.x, y=s2.y,Geom.line,color=[colorant"hotpink"]))

f(x) = Sequences.ssd(s1, Sequences.displace(s2, x))
r = optimize(f, -10.0, 10.0)

@test abs(r.minimizer)<0.1 # valid for low noise

# show results
p = plot(st)
push!(p, layer(x=s1.x, y=s1.y,Geom.line,color=[colorant"orange"]))
s3 = Sequences.displace(s2, r.minimizer)
push!(p, layer(x=s3.x, y=s3.y,Geom.line,color=[colorant"hotpink"]))

## 2. test w/ locally-referenced sequences

s1 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1:n]+sigma_y*randn(n))
s2 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1+Int(n/3):1+Int(n/3)+n-1]+sigma_y*randn(n))

Sequences.ssd(s1, s2) # should be low

p = plot(st)
push!(p, layer(x=s1.x, y=s1.y,Geom.line,color=[colorant"orange"]))
push!(p, layer(x=s2.x, y=s2.y,Geom.line,color=[colorant"hotpink"]))

f(x) = Sequences.ssd(s1, Sequences.displace(s2, x))
r = optimize(f, -10.0, 10.0)

@test abs(r.minimizer - 5.0)<0.1

# show results
p = plot(st)
push!(p, layer(x=s1.x, y=s1.y,Geom.line,color=[colorant"orange"]))
s3 = Sequences.displace(s2, r.minimizer)
push!(p, layer(x=s3.x, y=s3.y,Geom.line,color=[colorant"hotpink"]))

## 3. slightly more practical use - map ssd to probability

s1 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1:n]+sigma_y*randn(n))
s2 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1+Int(n/3):1+Int(n/3)+n-1]+sigma_y*randn(n))

p = plot(st)
push!(p, layer(x=s1.x, y=s1.y, Geom.line,color=[colorant"orange"]))
push!(p, layer(x=s2.x, y=s2.y, Geom.line,color=[colorant"hotpink"]))

f(x) = Sequences.ssd(s1, Sequences.displace(s2, x))
qr = collect(LinRange(-10.0,10.0,10001))
v = f.(qr)

k = abs(maximum(v))
pv = v./k
pv = exp.(-pv)
pv ./= sum(pv)
plot(x=qr,y=pv,Geom.line,Scale.y_log)
qr[argmax(pv)]


text:
module FiniteMarkovDecisionProcesses

include("_0_mdp.jl")
include("_1_mdp_utils.jl")
include("_2_simulation.jl")
include("_3_dynamic_programming.jl")
include("_4_monte_karlo.jl")
include("_5_temporal_difference.jl")

end # module


text:
using Plots

plotsdir(args...) = joinpath(@__DIR__, "output", args...)
isdir(plotsdir()) && rm(plotsdir(), recursive=true)
mkpath(plotsdir())

function mnerror(prec_frac, prec_abs; kwargs...)
    term = 1:length(prec_frac)
    m1 = @. !iszero(prec_abs)
    prec_abs = prec_abs[m1]
    prec_frac = filter(!iszero, prec_frac)

    plot(term[m1], prec_abs; xlabel="term", ylabel="abs. err", label="", yscale=:log10, rightmargin=70Plots.px, kwargs...)
    plot!(twinx(), prec_frac; yscale=:log10, ylabel="frac. err", label="", ygrid=false, kwargs...)
end


text:
function do_optimal_routing(input_data)

    spatial_data, spatial_data_matrices, cost_data, equipment_data = prepare_spatial_data(input_data)
    c_tot, optimal_path, ac_dc, ac_cab, dc_cab = optimize_route(spatial_data, spatial_data_matrices, cost_data, equipment_data, input_data)
    route_impedance, route_length = determine_impedance_parameters(input_data, spatial_data, spatial_data_matrices, optimal_path, equipment_data)
    return spatial_data, spatial_data_matrices, cost_data, equipment_data, c_tot, optimal_path, ac_dc, ac_cab, dc_cab, route_impedance, route_length
end


text:
using Documenter, LPIP

makedocs(
    modules = [LPIP],
    doctest  = false,
    clean = true,
    # assets = ["assets/favicon.ico"],
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "LPIP.jl",
    authors = "Guilherme Bodin, and contributors",
    pages = [
        "Home" => "index.md",
        "manual.md",
        "examples.md",
        "reference.md"
    ]
)

deploydocs(
    repo = "github.com/guilhermebodin/LPIP.jl.git",
)


text:
using Statistics
using TimeZones

testfile = joinpath(dir, "Example_Drive01.d7d")

f = DWDataReader.File(testfile)
@test string(f) == "DWDataReader.File(\"$(testfile)\"): 20 channels\n"
@test DWDataReader.getname(f.channels[1]) == "GPSvel"
@test f.channels[1].name == "GPSvel"
@test f.nchannels == 20

g = DWDataReader.File(testfile)
@test f.readerid != g.readerid

@test size(DWDataReader.scaled(f.channels[1]))[1] == 9580
@test size(DWDataReader.scaled(f.channels[1]))[2] == 2

@test size(DWDataReader.scaled(g.channels[1]))[1] == 9580
@test size(DWDataReader.scaled(g.channels[1]))[2] == 2

@test size(DWDataReader.scaled(f[:ENG_RPM]))[1] == 4791

c = DWDataReader.scaled(f[:ENG_RPM])
@test abs(mean(c[5.0 .<= c[:, 1] .<= 5.5, :][:, 2]) - 3098.5) < 1

@test f.info.start_store_time == ZonedDateTime(DateTime("2003-10-09T21:27:46.812"), tz"UTC")
@test f.info.sample_rate == 100.0
@test f.info.duration == 95.8
@test string(f.info) == "2003-10-09T21:27:46.812+00:00 100.0 Hz 95.8 s\n"


text:
##### Weighted var & std

## var

"""
    varm(x, wv::WeightVec, m, [dim])

Return the variance of a real-valued array `x` with a known mean `m`, optionally
over a dimension `dim`. The weighting vector `wv` specifies frequency weights
(also called case weights) for the result.

This function differs from its counterpart in Base in that Bessel's correction
is not used. That is, here the denominator for the variance is `sum(wv)`,
whereas it's `length(x)-1` in `Base.varm`. The impact is that this is not a
weighted estimate of the population variance based on the sample; it's the weighted
variance of the sample.
"""
Base.varm(v::RealArray, wv::WeightVec, m::Real) = _moment2(v, wv, m)

"""
    var(x, wv::WeightVec, [dim]; mean=nothing)

Return the variance of a real-valued array `x`, optionally over a dimension `dim`.
The weighting vector `wv` specifies frequency weights (also called case weights)
for the estimate.

This function differs from its counterpart in Base in that Bessel's correction
is not used. That is, here the denominator for the variance is `sum(wv)`,
whereas it's `length(x)-1` in `Base.var`. The impact is that this is not a
weighted estimate of the population variance based on the sample; it's the weighted
variance of the sample.
"""
function Base.var(v::RealArray, wv::WeightVec; mean=nothing)
    mean == 0 ? Base.varm(v, wv, 0) :
    mean == nothing ? varm(v, wv, Base.mean(v, wv)) :
    varm(v, wv, mean)
end

## var along dim

Base.varm!(R::AbstractArray, A::RealArray, wv::WeightVec, M::RealArray, dim::Int) =
    scale!(_wsum_centralize!(R, @functorize(abs2), A, values(wv), M, dim, true), inv(sum(wv)))

function var!(R::AbstractArray, A::RealArray, wv::WeightVec, dim::Int; mean=nothing)
    if mean == 0
        Base.varm!(R, A, wv,
            Base.reducedim_initarray(A, dim, 0, eltype(R)), dim)
    elseif mean == nothing
        Base.varm!(R, A, wv, Base.mean(A, wv, dim), dim)
    else
        # check size of mean
        for i = 1:ndims(A)
            dA = size(A,i)
            dM = size(mean,i)
            if i == dim
                dM == 1 || throw(DimensionMismatch("Incorrect size of mean."))
            else
                dM == dA || throw(DimensionMismatch("Incorrect size of mean."))
            end
        end
        Base.varm!(R, A, wv, mean, dim)
    end
end

Base.varm(A::RealArray, wv::WeightVec, M::RealArray, dim::Int) =
    @static if VERSION < v"0.6.0-dev.1121"
        Base.varm!(similar(A, Float64, Base.reduced_dims(size(A), dim)), A, wv, M, dim)
    else
        Base.varm!(similar(A, Float64, Base.reduced_indices(indices(A), dim)), A, wv, M, dim)
    end

Base.var(A::RealArray, wv::WeightVec, dim::Int; mean=nothing) =
    @static if VERSION < v"0.6.0-dev.1121"
        var!(similar(A, Float64, Base.reduced_dims(size(A), dim)), A, wv, dim; mean=mean)
    else
        var!(similar(A, Float64, Base.reduced_indices(indices(A), dim)), A, wv, dim; mean=mean)
    end

## std
"""
    stdm(v, wv::WeightVec, m, [dim])

Return the standard deviation of a real-valued array `v` with a known mean `m`,
optionally over a dimension `dim`. The weighting vector `wv` specifies frequency
weights (also called case weights) for the estimate.
"""
Base.stdm(v::RealArray, wv::WeightVec, m::Real) = sqrt(varm(v, wv, m))

"""
    std(v, wv::WeightVec, [dim]; mean=nothing)

Return the standard deviation of a real-valued array `v`, optionally over a
dimension `dim`. The weighting vector `wv` specifies frequency weights (also
called case weights) for the estimate.
"""
Base.std(v::RealArray, wv::WeightVec; mean=nothing) = sqrt.(var(v, wv; mean=mean))

Base.stdm(v::RealArray, m::RealArray, dim::Int) = Base.sqrt!(varm(v, m, dim))
Base.stdm(v::RealArray, wv::WeightVec, m::RealArray, dim::Int) = sqrt.(varm(v, wv, m, dim))
Base.std(v::RealArray, wv::WeightVec, dim::Int; mean=nothing) = sqrt.(var(v, wv, dim; mean=mean))

##### Fused statistics
"""
    mean_and_var(x, [wv::WeightVec], [dim]) -> (mean, var)

Return the mean and variance of a real-valued array `x`, optionally over a dimension
`dim`, as a tuple. A weighting vector `wv` can be specified to weight the estimates.
The weights are assumed to be frequency weights, also called case weights.
"""
mean_and_var(A::RealArray) = (m = mean(A); (m, varm(A, m)))

"""
    mean_and_std(x, [wv::WeightVec], [dim]) -> (mean, std)

Return the mean and standard deviation of a real-valued array `x`, optionally
over a dimension `dim`, as a tuple. A weighting vector `wv` can be specified
to weight the estimates. The weights are assumed to be frequency weights, also
called case weights.
"""
mean_and_std(A::RealArray) = (m = mean(A); (m, stdm(A, m)))

mean_and_var(A::RealArray, wv::WeightVec) = (m = mean(A, wv); (m, varm(A, wv, m)))
mean_and_std(A::RealArray, wv::WeightVec) = (m = mean(A, wv); (m, stdm(A, wv, m)))

mean_and_var(A::RealArray, dim::Int) = (m = mean(A, dim); (m, varm(A, m, dim)))
mean_and_std(A::RealArray, dim::Int) = (m = mean(A, dim); (m, stdm(A, m, dim)))

mean_and_var(A::RealArray, wv::WeightVec, dim::Int) = (m = mean(A, wv, dim); (m, varm(A, wv, m, dim)))
mean_and_std(A::RealArray, wv::WeightVec, dim::Int) = (m = mean(A, wv, dim); (m, stdm(A, wv, m, dim)))

##### General central moment

function _moment2(v::RealArray, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += z * z
    end
    s / n
end

function _moment2(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += (z * z) * w[i]
    end
    s / sum(wv)
end

function _moment3(v::RealArray, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += z * z * z
    end
    s / n
end

function _moment3(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += (z * z * z) * w[i]
    end
    s / sum(wv)
end

function _moment4(v::RealArray, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += abs2(z * z)
    end
    s / n
end

function _moment4(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += abs2(z * z) * w[i]
    end
    s / sum(wv)
end

function _momentk(v::RealArray, k::Int, m::Real)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += (z ^ k)
    end
    s / n
end

function _momentk(v::RealArray, k::Int, wv::WeightVec, m::Real)
    n = length(v)
    s = 0.0
    w = values(wv)
    for i = 1:n
        @inbounds z = v[i] - m
        @inbounds s += (z ^ k) * w[i]
    end
    s / sum(wv)
end

"""
    moment(v, k, [wv::WeightVec], m=mean(v))

Return the `k`th order central moment of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""
function moment(v::RealArray, k::Int, m::Real)
    k == 2 ? _moment2(v, m) :
    k == 3 ? _moment3(v, m) :
    k == 4 ? _moment4(v, m) :
    _momentk(v, k, m)
end

function moment(v::RealArray, k::Int, wv::WeightVec, m::Real)
    k == 2 ? _moment2(v, wv, m) :
    k == 3 ? _moment3(v, wv, m) :
    k == 4 ? _moment4(v, wv, m) :
    _momentk(v, k, wv, m)
end

moment(v::RealArray, k::Int) = moment(v, k, mean(v))
moment(v::RealArray, k::Int, wv::WeightVec) = moment(v, k, wv, mean(v, wv))

##### Skewness and Kurtosis

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    skewness(v, [wv::WeightVec], m=mean(v))

Compute the standardized skewness of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""
function skewness(v::RealArray, m::Real)
    n = length(v)
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    for i = 1:n
        @inbounds z = v[i] - m
        z2 = z * z

        cm2 += z2
        cm3 += z2 * z
    end
    cm3 /= n
    cm2 /= n
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

function skewness(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    w = values(wv)

    @inbounds for i = 1:n
        x_i = v[i]
        w_i = w[i]
        z = x_i - m
        z2w = z * z * w_i
        cm2 += z2w
        cm3 += z2w * z
    end
    sw = sum(wv)
    cm3 /= sw
    cm2 /= sw
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

skewness(v::RealArray) = skewness(v, mean(v))
skewness(v::RealArray, wv::WeightVec) = skewness(v, wv, mean(v, wv))

# (excessive) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    kurtosis(v, [wv::WeightVec], m=mean(v))

Compute the excess kurtosis of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""
function kurtosis(v::RealArray, m::Real)
    n = length(v)
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    for i = 1:n
        @inbounds z = v[i] - m
        z2 = z * z
        cm2 += z2
        cm4 += z2 * z2
    end
    cm4 /= n
    cm2 /= n
    return (cm4 / (cm2 * cm2)) - 3.0
end

function kurtosis(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    w = values(wv)

    @inbounds for i = 1 : n
        x_i = v[i]
        w_i = w[i]
        z = x_i - m
        z2 = z * z
        z2w = z2 * w_i
        cm2 += z2w
        cm4 += z2w * z2
    end
    sw = sum(wv)
    cm4 /= sw
    cm2 /= sw
    return (cm4 / (cm2 * cm2)) - 3.0
end

kurtosis(v::RealArray) = kurtosis(v, mean(v))
kurtosis(v::RealArray, wv::WeightVec) = kurtosis(v, wv, mean(v, wv))


text:
H = Hilbert(10)
@test H[10,10] == 1//(10+10-1)
@test inv(H) == InverseHilbert(10)


text:
module mION

JUPL_PATH = ARGS[1] #get path from command line

push!(LOAD_PATH, "$JUPL_PATH/init")
push!(LOAD_PATH, "$JUPL_PATH/src")
using Glob

include("$JUPL_PATH/src/io.jl")
output_files = glob("output/*.csv")
for file in output_files
    rm(file)
end
jpg_files = glob("output/*.jpg")
for file in jpg_files
    rm(file)
end

using Constants
using Get_Sim_Params

include("$JUPL_PATH/init/process_params.jl")
const dr, mi_g, mi, qi, nquiet, cs, CFL, dtm,
        lm, nsmooth, vsmooth, tsmooth = process_params()

include("$JUPL_PATH/src/types.jl")
hydro = hydro_vars(nz, nspec, neqi)

include("$JUPL_PATH/init/init_sim.jl")
const r = init_spacegrid(dr)
const eps_visc = init_art_viscosity(geom)
if !restart
    hydro = init_variables(hydro)
end
hydro = do_smoothing(hydro)
hydro = init_fluxes(hydro)
init_predictor_corrector(hydro)

erf_table = get_erf_integral(JUPL_PATH)
write_sim_parameters()
open_files()
write_data(hydro)
close_files()

include("$JUPL_PATH/src/utils.jl")
include("$JUPL_PATH/src/boundary.jl")
include("$JUPL_PATH/src/numerics.jl")
include("$JUPL_PATH/src/collisions.jl")

nq = 0
j = 0
while tm <= maxTime
    j += 1
    tm += dtm
    is_print = is_print_time(tm, dt_print)

    nq, time1, time2 = is_quiet_time!(time1, time2, tm, tm_quiet, geom, nz, nq, dt_print, dtm)

    apply_BC!(hydro)

    #----- predictor -----------
    #---------------------------
    predictor!(hydro, nq, j)
    update_variables!(hydro.U1D_p, hydro)
    hydro, k_DT, ke, Qextra = calculate_collisions!(hydro, erf_table)
    source_terms!(hydro, k_DT, ke, Qextra, nq)
      # #----start debug
      # if j==1
      #     open_files()
      #     write_data(hydro)
      #     close_files()
      #     write_all_data(hydro)
      #     hydro, k_DT, ke, Qextra = calculate_collisions!(hydro, erf_table)
      #     systemerror(0)
      # end
      # #----end debug

    #----- corrector -----------
    #---------------------------
    corrector!(hydro, nq, j)

    #----- final step ----------
    hydro.U1D[2:nq-1,1:neqi+1] = 0.5 * (hydro.U1D[2:nq-1,1:neqi+1] + hydro.U1D_c[2:nq-1,1:neqi+1]) +
        eps_visc[1:nq-2,1:neqi+1] .*
        (hydro.U1D[3:nq,1:neqi+1] - 2 * hydro.U1D[2:nq-1,1:neqi+1] + hydro.U1D[1:nq-2,1:neqi+1])

    update_variables!(hydro.U1D, hydro)
    hydro, k_DT, ke, Qextra = calculate_collisions!(hydro, erf_table)
    source_terms!(hydro, k_DT, ke, Qextra, nq)

    if is_print
        open_files()
        write_data(hydro)
        close_files()
        is_print = false
    end

end

println("---Simulation ended---")

close_files()

end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
# defining some constants for use in the tests
const tol = 1e-6

const ρinit = [1.0 0.0]' * [1.0 0.0 + 0im]
const ρfin = [0.0 1.0]' * [0.0 1.0 + 0im]

const Sρinit = SArray{Tuple{2,2},ComplexF64}(ρinit)
const Sρfin = SArray{Tuple{2,2},ComplexF64}(ρfin)

const Sx = [0.0 1.0; 1.0 0.0+0.0im] / 2
const Sy = [0.0 -1.0im; 1.0im 0.0+0.0im] / 2
const Sz = [1.0 0.0; 0.0 -1.0+0.0im] / 2

const SSx = SArray{Tuple{2,2},ComplexF64}(Sx)
const SSy = SArray{Tuple{2,2},ComplexF64}(Sy)
const SSz = SArray{Tuple{2,2},ComplexF64}(Sz)

const Uinit = Array{ComplexF64,2}((I(2)))
const Ufin = [0 1; 1 -0.0+0.0im]

const ρinitS = SMatrix{2,2,ComplexF64}([1.0 0.0]' * [1.0 0.0 + 0im])
const ρfinS = SMatrix{2,2,ComplexF64}([0.0 1.0]' * [0.0 1.0 + 0im])

const UinitS = SMatrix{2,2,ComplexF64}(Array{ComplexF64,2}((I(2))))
const UfinS = SMatrix{2,2,ComplexF64}([0 1; 1 -0.0+0.0im])

A_gens = k -> (k - 2.5) / 2.5 * Sz * 5
B_gens = k -> [Sx, Sy]

A_gens_static = k -> (k - 2.5) / 2.5 * SSz * 5
B_gens_static = k -> [SSx, SSy]

function odd_switch(x)
    if Bool(x % 2)
        return ρfin
    else
        return ρinit
    end
end

X_init_gens = k -> ρinit
X_target_gens = k -> odd_switch(k)

X_init_gens_static = k -> ρinitS
X_target_gens_static = k -> SMatrix{2,2,ComplexF64}(odd_switch(k))

function odd_switch_unitary(x)
    if Bool(x % 2)
        return Ufin
    else
        return Uinit
    end
end

U_init_gens = k -> Uinit
U_target_gens = k -> odd_switch_unitary(k)

U_init_gens_static = k -> UinitS
U_target_gens_static = k -> SMatrix{2,2,ComplexF64}(odd_switch_unitary(k))


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Test Queue" begin

    # Initial a Queue with data '0'
    T_e = Int64
    T = DS.Queue{T_e}
    data= T_e[0]
    queue = T(0)
    @test queue.data == data

    # Type test
    @test queue isa T

    # Push test: Push data into queue 20 times by 1:20
    passed_check_data = []
    for i = 1:20
        push!(queue, i)
        push!(data, i)
        push!(passed_check_data, queue.data == data)
    end
    @test all(passed_check_data)

    # Pop test: Pop data from queue 5 times
    passed_check_data = []
    for i = 1:5
        pop!(queue)
        popfirst!(data)
        push!(passed_check_data, queue.data == data)
    end
    @test all(passed_check_data)

    # Show Test
    @test repr(queue) == "Queue(["*
        "5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]"*
    ")"

end


text:
# "Fundamental" literals
crawl(ex::String) = :(JSTerminal($ex))
crawl(b::Bool) = crawl(Symbol(b))

# Symbols
# We sub-dispatch on symbols to allow for special handling of certain keywords.
crawl(sym::Symbol) = crawl_symbol(Val(sym))
function crawl_symbol(sym::Val{S}) where {S}
    return :(JSTerminal($(QuoteNode(S))))
end
deparse(sym::JSTerminal)::JSString = sym.s

# Convert the nothing symbol into null.
crawl_symbol(::Val{:nothing}) = :(JSTerminal(nothing))


text:
"""
Compute preimages of monotonic sequences
"""

using IntervalArithmetic
using .Contractors

"""
Type used to represent a "branch" of a dynamic. The branch is represented by a monotonic map `f` with domain `X=(a,b)` with a≤b (where typically a,b are intervals).
`Y=(f(a),f(b))` and `increasing` may be provided (for instance if we know that `Y=(0,1)`), otherwise they are computed automatically.
"""
struct Branch{T,S}
    f::T
    X::Tuple{S, S}
    Y::Tuple{S, S}
    increasing::Bool
end
Branch(f, X, Y=(f(Interval(X[1])), f(Interval(X[2]))), increasing=unique_increasing(Y[1], Y[2])) = Branch{typeof(f), typeof(interval(X[1]))}(f, X, Y, increasing)

"""
Return Branches for a given PwMap, in an iterable.

TODO: in future, maybe it is a better idea to replace the type PwMap directly with an array of branches, since that's all we need
"""
function branches(D::PwMap)
    return [Branch(D.Ts[k], (D.endpoints[k], D.endpoints[k+1]), (D.y_endpoints[k,1], D.y_endpoints[k,2]), D.increasing[k]) for k in 1:length(D.Ts)]
end

"""
Smallest possible i such that a is in the semi-open interval [y[i], y[i+1]).

This should work properly even if `a, y` are intervals; in this case it returns the *smallest* possible value of i over all possible "assignments" of a, y inside those intervals.
Assumes y is sorted, i.e., map(y, x->Interval(x).lo) and map(y, x->Interval(x).hi) are sorted.
"""
function first_overlapping(y, a)
    if iszero(a) # avoids -0 crap
        a = zero(a)
    end
    searchsortedlast(y, Interval(a).lo, by=x->Interval(x).hi)
end

"""
Largest possible j such that a-ε is in the semi-open interval [y[j], y[j+1]).

This should work properly even if `a, y` are intervals; in this case it returns the *largest* possible value of i over all possible "assignments" of a, y inside those intervals.
Assumes y is sorted, i.e., map(y, x->Interval(x).lo) and map(y, x->Interval(x).hi) are sorted.
"""
function last_overlapping(y, a)
    if iszero(a) # avoids -0 crap
        a = zero(a)
    end
    searchsortedfirst(y, Interval(a).hi, by=x->Interval(x).lo) - 1
end

"""
Construct preimages of an increasing array y under a monotonic branch defined on X = (a, b), propagating additional labels `ylabel`

The sequence y subdivides the y-axis into semi-open intervals [y[l], y[l+1]); each of them is identified by the label `ylabel[l]`. We construct an increasing sequence
x that splits X (in the x-axis) into semi-open intervals, each of them with f([x[k], x[k+1]) ⊂ [y[l], y[l+1]) for a certain l.
We set xlabel[k] = ylabel[l], and return the pair (x, xlabel).

In the simplest case where D is full-branch, the points in x are preimages of the points in y, but in the general case they can also include D.endpoints:
in general, there may be a certain number of points in y that have no preimage at the beginning and the end of the sequence, because
they fall out of the range R = [f(a), f(b)]. In the worst case, no point has a preimage, because y[i] < R < y[i+1] for some
i (or vice versa with orientations), and in this case we just return the 1-element vectors x = [branch.X[1]] and xlabel = [i].

x[begin] always coincides with branch.X[1], while branch.X[2] is "the point after x[end]", and is not stored explicitly in x, for easier composing.
In this way x and xlabel have the same length.

This function fills the array by using a bisection strategy to save computations: if y ∈ [a,b], then f⁻¹(y) ∈ [f⁻¹(a),f⁻¹(b)] (paying attention to orientation).
So we can fill v by filling in first entries `v[k+1]` with higher dyadic valuation of k.

For a dynamic with multiple branches, preimages(y, D) is simply the concatenation of x, xlabel for b in all branches. These values still form an increasing sequence that
splits X into intervals, each of which is mapped into a different semi-open interval [y[k], y[k+1]).
"""
function preimages(y, br::Branch, ylabel = 1:length(y), ϵ = 0.0)

    if br.increasing
        i = first_overlapping(y, br.Y[1])  # smallest possible i such that a is in the semi-open interval [y[i], y[i+1]).
        j = last_overlapping(y, br.Y[2]) # largest possible j such that b-ε is in the semi-open interval [y[j], y[j+1]).
        n = j - i + 1
        x = fill((-∞..∞)::typeof(Interval(br.X[1])), n)
        xlabel = collect(ylabel[i:j]) # we collect to avoid potential type instability, since this may be an UnitRange while in the other branch we could have a StepRange
        x[1] = br.X[1]
        if n == 1
            return (x, xlabel)
        end
        # the bisection strategy: fill the array in "strides" of length `stride`, then halve the stride and repeat
        # for instance, if the array is 1..13 (with x[1] filled in already), we first take stride=8 and fill in x[9],
        # then stride=4 and fill in x[5], x[13], (note that the distance is `2stride`, since x[9], and in general all the even multiples of `stride`, is already filled in)
        # then stride=2 and fill in x[3], x[7], x[11],
        # then stride=1 and fill in x[2], x[4], x[6], x[8], x[10], x[12]
        # at each step we have bracketed the preimage in a "search range" given by already-computed preimages x[k-stride] and x[k+stride].
        stride = prevpow(2, n-1)
        while stride >= 1
            # fill in v[i] using x[i-stride].lo and x[i+stride].hi as range for the preimage search
            for k = 1+stride:2*stride:n
                search_range = Interval(x[k-stride].lo, (k+stride <= n ? x[k+stride] : Interval(br.X[2])).hi)
                x[k] = preimage(y[i-1+k], br.f, search_range, ϵ)
            end
            stride = stride ÷ 2
        end
    else # branch decreasing
        i = last_overlapping(y, br.Y[1]) # largest possible j such that b-ε is in the semi-open interval [y[j], y[j+1]).
        j = first_overlapping(y, br.Y[2]) # smallest possible i such that a is in the semi-open interval [y[i], y[i+1]).
        n = i - j + 1
        x = fill((-∞..∞)::typeof(Interval(br.X[1])), n)
        xlabel = collect(ylabel[i:-1:j])
        x[1] = br.X[1]
        if n == 1
            return (x, xlabel)
        end
        stride = prevpow(2, n-1)
        while stride >= 1
            # fill in v[i] using x[i-stride].lo and x[i+stride].hi as range for the preimage search
            for k = 1+stride:2*stride:n
                search_range = Interval(x[k-stride].lo, (k+stride <= n ? x[k+stride] : Interval(br.X[2])).hi)
                x[k] = preimage(y[i+2-k], br.f, search_range, ϵ)
            end
            stride = stride ÷ 2
        end
    end
    return (x, xlabel)
end

function preimages(y, D::Dynamic, ylabel = 1:length(y), ϵ = 0.0)
    results = collect(preimages(y, b, ylabel, ϵ) for b in branches(D))
    x = vcat((result[1] for result in results)...)
    xlabel = vcat((result[2] for result in results)...)
    return x, xlabel
end

"""
Compute preimages of D *and* the derivatives f'(x) in each point.

Returns: x, xlabel, x′

Assumes that the dynamic is full-branch, because otherwise things may compose the wrong way.
This is not restrictive because we'll need it only for the Hat assembler (at the moment)

We combine them in a single function because there are avenues to optimize by recycling some computations (not completely exploited for now)
"""
function preimages_and_derivatives(y, br::Branch, ylabel = 1:length(y), ϵ = 0.0)
    x, xlabel = preimages(y, br, ylabel, ϵ)
    f′ = Contractors.derivative(br.f)
    x′ = f′.(x)
    return x, xlabel, x′
end
function preimages_and_derivatives(y, D::Dynamic, ylabel = 1:length(y), ϵ = 0.0)
    @assert is_full_branch(D)
    results = collect(preimages_and_derivatives(y, b, ylabel, ϵ) for b in branches(D))
    x = vcat((result[1] for result in results)...)
    xlabel = vcat((result[2] for result in results)...)
    x′ = vcat((result[3] for result in results)...)
    return x, xlabel, x′
end

"""
Composed map D1 ∘ D2 ∘ D3. We store with [D1, D2, D3] in this order.

We overwrite ∘ in base, so one can simply write D1 ∘ D2 or ∘(D1, D2, D3) to construct them.
"""
struct ComposedDynamic <: Dynamic
    dyns::Tuple{Vararg{Dynamic}}
end
Base.:∘(d::Dynamic...) = ComposedDynamic(d)

"""
Utility function to return the domain of a dynamic
"""
domain(D::PwMap) = (D.endpoints[begin], D.endpoints[end])
domain(D::ComposedDynamic) = domain(D.dyns[end])

function preimages(z, Ds::ComposedDynamic, zlabel = 1:length(z), ϵ = 0.0)
    for d in Ds.dyns
        z, zlabel = preimages(z, d, zlabel, ϵ)
    end
    return z, zlabel
end
function preimages_and_derivatives(z, Ds::ComposedDynamic, zlabel = 1:length(z), ϵ = 0.0)
    derivatives = fill(1, length(z))
    for d in Ds.dyns
        z, zindex, z′ = preimages_and_derivatives(z, d, 1:length(z), ϵ)
        zlabel = zlabel[zindex]
        derivatives = derivatives[zindex] .* z′
    end
    return z, zlabel, derivatives
end

"""
Replacement of DualComposedWithDynamic.
"""
abstract type Dual end

struct UlamDual <: Dual
    x::Vector{Interval} #TODO: a more generic type may be needed in future
    xlabel::Vector{Int}
    lastpoint::Interval
end
Dual(B::Ulam, D, ϵ) = UlamDual(preimages(B.p, D, 1:length(B.p)-1, ϵ)..., domain(D)[end])

Base.length(dual::UlamDual) = length(dual.x)
Base.eltype(dual::UlamDual) = Tuple{eltype(dual.xlabel), Tuple{eltype(dual.x), eltype(dual.x)}}
function iterate(dual::UlamDual, state = 1)
    n = length(dual.x)
    if state < n
        return (dual.xlabel[state], (dual.x[state], dual.x[state+1])), state+1
    elseif state == n
        return (dual.xlabel[n], (dual.x[n], dual.lastpoint)), state+1
    else
        return nothing
    end
end

struct HatDual <: Dual
    x::Vector{Interval} #TODO: a more generic type may be needed in future
    xlabel::Vector{Int}
    x′::Vector{Interval}
end

Dual(B::Hat, D, ϵ) = HatDual(preimages_and_derivatives(B.p, D, 1:length(B.p)-1, ϵ)...)
Base.length(dual::HatDual) = length(dual.x)
Base.eltype(dual::HatDual) = Tuple{eltype(dual.xlabel), Tuple{eltype(dual.x), eltype(dual.x′)}}
function iterate(dual::HatDual, state=1)
    if state <= length(dual.x)
        return ((dual.xlabel[state], (dual.x[state], abs(dual.x′[state]))), state+1)
    else
        return nothing
    end
end

# Variants of assemble and DiscretizedOperator; the code is repeated here for easier comparison with the older algorithm
function assemble2(B, D, ϵ=0.0; T = Float64)
	I = Int64[]
	J = Int64[]
	nzvals = Interval{T}[]
	n = length(B)

	# TODO: reasonable size hint?

	for (i, dual_element) in Dual(B, D, ϵ)
		if !is_dual_element_empty(B, dual_element)
			for (j, x) in ProjectDualElement(B, dual_element)
				push!(I, i)
				push!(J, mod(j,1:n))
				push!(nzvals, x)
			end
		end
	end

	return sparse(I, J, nzvals, n, n)
end

function DiscretizedOperator2(B, D, ϵ=0.0; T = Float64)
	L = assemble2(B, D, ϵ; T)
	if is_integral_preserving(B)
		return IntegralPreservingDiscretizedOperator(L)
	else
		f = integral_covector(B)
		e = one_vector(B)
		w = f - f*L #will use interval arithmetic when L is an interval matrix
		return NonIntegralPreservingDiscretizedOperator(L, e, w)
	end
end