text:
#
# @author DeadbraiN
#
module Helper
  export CodePos
  export Point
  export RetObj

  export toBytes
  export fRand
  export fastRand
  export fastZRand
  export getProbIndex
  export isopen
  export done
  export info
  export warn
  export error
  export emptyFn
  #
  # Shorthand for Float16 type
  #
  const f16 = Float16::DataType
  #
  # These constants were gotten from julia/base/stream.jl file. They
  # are used for optimization of isopen() function.
  #
  const StatusUninit = 0 # handle is allocated, but not initialized
  const StatusInit   = 1 # handle is valid, but not connected/active
  const StatusClosed = 6 # handle is closed
  const StatusEOF    = 7 # handle is a TTY that has seen an EOF event
  #
  # Describes organism's code position. This position is used
  # for mutations (add, change, del,...)
  #
  type CodePos
    #
    # Index of function in Organism.Creature.funcs array
    #
    fnIdx::Int
    #
    # Index of block inside Creature.Func function
    #
    blockIdx::Int
    #
    # Index of code line inside Creature.Block.lines
    #
    lineIdx::Int
  end
  #
  # One point in the world. Is described by two coordinates.
  #
  type Point
    x::Int
    y::Int
  end
  #
  # Universal structure for returning a value from event handlers.
  # See "beforeclone", "getenergy" and other events for details.
  #
  type RetObj
    #
    # Return value
    #
    ret::Any
    #
    # Position in a world
    #
    pos::Helper.Point
    #
    # ctor
    #
    RetObj(r = nothing, p = nothing) = (x = new(r); p === nothing ? x : (x.pos = p;x))
  end
  #
  # Converts any type to bytes array
  # @param x Value we have to convert
  # @return {Vector{UInt8}}
  #
  function toBytes(x::Any)
     local sz::Int = sizeof(x)
     local ba::Vector{UInt8} = Vector{UInt8}(sz)
     local ptr::Ptr{UInt8} = convert(Ptr{UInt8}, pointer_from_objref(x))
     unsafe_copy!(pointer(ba), ptr, sz)

     ba
  end
  #
  # Prints warning message with white color
  # @param msg Text message to show
  # OPT: console output is heavy operation
  function warn(msg::AbstractString)
    print_with_color(:yellow, "WARNING: ", msg, "\n")
  end
  #
  # Print error message with red color and don't stop an application
  # @param msg Text message to show
  # OPT: console output is heavy operation
  function error(msg::AbstractString)
    print_with_color(:red, "ERROR: ", msg, "\n")
  end
  #
  # Prints info message with green color
  # @param msg Text message to show
  # @param newLine Add new line character at the end of message
  #
  function info(msg::AbstractString, newLine::Bool = true)
    print_with_color(:green, string("INFO: ", msg, newLine ? "\n" : ""))
  end
  #
  # Prints "done" after info message (see info() function)
  # @param msg Text message to show
  #
  function done(msg::AbstractString = "done")
    print_with_color(:green, string(msg, "\n"))
  end
  #
  # This method is an analog of isopen() from julia/base/stream.jl.
  # We created this, because it original version is slow, because
  # of exception inside it. In case if some constants or it's
  # implementation will be changed, we have to update it.
  # @param sock Socket we are checking
  # @return {Bool}
  #
  function isopen{SocketType}(sock::SocketType)
      if sock.status == StatusUninit || sock.status == StatusInit return false end
      sock.status != StatusClosed && sock.status != StatusEOF
  end
  #
  # Returns lines array from AST by specified indexes. e.g.:
  # [2,1] means ex.args[2].args[1].args
  # @param exp Expression we have to start with
  # @return {Array{Int, 1}}
  #
  function getLines(exp::Expr, indexes::Array{Int, 1})
    getArg(exp, indexes).args
  end
  #
  # Returns head of specified AST node. An array of indexes
  # is used like a node position. For example:
  # [2,1] means ex.args[2].args[1].head
  # @param exp Expression we have to start with
  # @return {Array{Int, 1}}
  #
  function getHead(exp::Expr, indexes::Array{Int, 1})
    getArg(exp, indexes).head
  end
  #
  # Returns argument of specified AST node. An array of indexes
  # is used like a node position. For example:
  # [2,1] means ex.args[2].args[1]
  # @param exp Expression we have to start with
  # @return {Array{Int, 1}}
  #
  function getArg(exp::Expr, indexes::Array{Int, 1})
    local args::Array{Any, 1} = exp.args
    local len::Int = length(indexes)

    for i = 1:(len-1)
      args = args[indexes[i]].args
    end

    args[indexes[len]]
  end
  #
  # Generates random Float16 number in all available range -XX.XX...XX.XX
  # @return {Float16}
  #
  function fRand() f16((rand() * 2.0 - 1.0) * 6.55e4) end
  #
  # Fast version of rand(n::Range) function. Generates random Int number in
  # range 1:n
  # @param n Right number value in a range
  #
  function fastRand(n::Int) trunc(Int, rand() * n) + 1 end
  #
  # Fast version of rand(n::Range) function. Generates random UInt number in
  # range 1:n
  # @param n Right number value in a range
  #
  function fastRand(n::UInt) trunc(UInt, rand() * n) + UInt(1) end
  #
  # Fast version of rand(0::Range) function. Generates random Int number in
  # range 0:n
  # @param n Right number value in a range
  #
  function fastZRand(n::Int) trunc(Int, rand() * (n + 1)) end
  #
  # It calculates probability index from variable amount of components.
  # Let's imagine we have three actions: one, two and three. We want
  # these actions to be called randomly, but with different probabilities.
  # For example it may be [3,2]. It means that one should be called
  # in half cases, two in 1/3 cases and three in 1/6 cases. Probabilities
  # should be greated then -1.
  # @param {Array{Int, 1}} prob Probabilities array. e.g.: [3,2] or [1,3]
  # @return {Int} 0 Means that index is invalid
  #
  function getProbIndex(prob::Array{Int, 1})
    local len::Int = length(prob)
    if len < 1 return 0 end
    local s::Int = sum(prob)
    if s < 1 return 0 end
    local num::Int = fastRand(s)
    local i::Int = 1
    #
    # This is small optimization algorithm. if random number in
    # a left part of all numbers sum, the we have to go to it from
    # left to right, if not - then from right to left. Otherwise,
    # going every time from left to right will be a little bit
    # slower then this aproach.
    #
    if num < div(s, 2)
      s = 0
      for i = 1:len if num <= (s += prob[i]) break end end
    else
      for i = len:-1:1 if num > (s -= prob[i]) break end end
    end

    i
  end
  #
  # The same like getProbIndex(prob::Array{Int, 1}), but for UInt type
  # @param {Array{UInt, 1}} prob Probabilities array. e.g.: [3,2] or [1,3]
  # @return {Int} 0 Means that index is invalid
  #
  function getProbIndex(prob::Array{UInt, 1})
    local len::Int = length(prob)
    if len < 1 return 0 end
    local s::UInt = sum(prob)
    if s < 1 return 0 end
    local num::UInt = UInt(fastRand(s))
    local i::Int = 1
    #
    # This is small optimization algorithm. if random number in
    # a left part of all numbers sum, the we have to go to it from
    # left to right, if not - then from right to left. Otherwise,
    # going every time from left to right will be a little bit
    # slower then this aproach.
    #
    if num < div(s, UInt(2))
      s = UInt(0)
      for i = 1:len if num <= (s += prob[i]) break end end
    else
      for i = len:-1:1 if num > (s -= prob[i]) break end end
    end

    i
  end
  #
  # Returns default empty map with supported Julia types
  # @return {Dict{DataType, Array{Symbol, 1}}}
  #
  function getTypesMap()
    Dict{DataType, Array{Symbol, 1}}(f16 => Symbol[])
  end
  #
  # Checks if position is empty. x == y == 0 - this is empty
  # @param pos Position to check
  #
  function empty(pos::Point)
    pos.x === pos.y === 0
  end
  #
  # Saves custom data into the file. If file exists, it will
  # be overrided. It's only rewrites existing file and not
  # append it. It also doesn't work with native types, in sense
  # that you can't save native values into the file without#
  # meta information. For example, you can't store ascii string
  # in a file without special prefic before it.
  # @param data Data to save
  # @param file File name
  # @param asText true if data need to be saved as a text
  # @return {Bool} saving result
  #
  function save(data::Any, file::String = "backup.data", asText::Bool = false)
    local io  = nothing
    local ret = true

    try
      io = open(file, "w")
      if asText write(io, data)
      else serialize(io, data)
      end
    catch(e)
      warn("Helper.save(): $e")
      showerror(STDOUT, e, catch_backtrace())
      ret = false
    finally
      if io !== nothing close(io) end
    end

    ret
  end
  #
  # Loads custom data from the file
  # @param file File name
  # @return {Any|nothing} loading result or nothing
  #
  function load(file::String = "backup.data")
    local io  = nothing
    local ret = nothing

    try
      io  = open(file)
      ret = deserialize(io)
    catch(e)
      warn("Helper.load(): $e")
      showerror(STDOUT, e, catch_backtrace())
      ret = nothing
    finally
      if io !== nothing close(io) end
    end

    ret
  end
  #
  # Just stub function, which is used as a placeholder for saving a backup
  #
  function emptyFn() end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
module TrackingTimers

export TrackingTimer

using Printf
using Distributed # for myid()
using Tables, PrettyTables

#####
##### Implementation
#####

const TrackingTimerElType = @NamedTuple{name::String,time::Float64,gctime::Float64,
                                        n_allocs::Int,bytes::Int,thread_id::Int,pid::Int}

"""
    TrackingTimer

Stores the results of `@timed` calls in a `RemoteChannel` to provide a
distributed and thread-friendly way to collect timing results. Construct
a `TrackingTimer` by `t = TrackingTimer()`. Populate it by
[`TrackingTimers.@timeit`](@ref) or call `t` on a function (optionally providing a name)
to obtain an `InstrumentedFunction`, which automatically populates `t` with
timing results whenever it is called.

TrackingTimer's support the Tables.jl row table interface. Call `Tables.rows(t)`
to obtain a `Vector{$(TrackingTimerElType)}` of the timing results obtained so far.
Note that this uses a lock (via [`synchronize!`](@ref)) so it should be thread-safe
albeit may cause contention if called from multiple threads simultaneously.

## Example

```julia
julia> using Distributed

julia> addprocs(2);

julia> @everywhere using TrackingTimers

julia> t = TrackingTimer()
TrackingTimer: 1.09 s since creation (0% measured).
No entries.

julia> @everywhere function f(i)
           v = 1:(1000*myid())
           return sum( v .^ (1/π))
       end

julia> f_inst = t(f) # instrument `f` with TrackingTimer `t`
(::TrackingTimers.InstrumentedFunction{typeof(f)}) (generic function with 1 method)

julia> pmap(f_inst, 1:10)
10-element Vector{Float64}:
 17056.850202253918
 29106.968991882364
 29106.968991882364
 17056.850202253918
 17056.850202253918
 29106.968991882364
 17056.850202253918
 29106.968991882364
 17056.850202253918
 29106.968991882364

julia> t
TrackingTimer: 2.54 s since creation (0% measured).
 name   time   gcfraction  n_allocs    allocs    thread ID  proc ID
────────────────────────────────────────────────────────────────
 f     0.00 s      0%         2  23.516 KiB          1        3
 f     0.00 s      0%         1  15.750 KiB          1        2
 f     0.00 s      0%         2  23.516 KiB          1        3
 f     0.00 s      0%         2  23.516 KiB          1        3
 f     0.00 s      0%         2  23.516 KiB          1        3
 f     0.00 s      0%         2  23.516 KiB          1        3
 f     0.00 s      0%         1  15.750 KiB          1        2
 f     0.00 s      0%         1  15.750 KiB          1        2
 f     0.00 s      0%         1  15.750 KiB          1        2
 f     0.00 s      0%         1  15.750 KiB          1        2
```
"""
struct TrackingTimer
    start_nanosecond::Int64
    chan::RemoteChannel{Channel{TrackingTimerElType}}
    results::Vector{TrackingTimerElType}
    results_lock::ReentrantLock
    function TrackingTimer()
        return new(time_ns(),
                   RemoteChannel(() -> Channel{TrackingTimerElType}(typemax(Int))),
                   TrackingTimerElType[], ReentrantLock())
    end
end

Base.put!(t::TrackingTimer, val) = put!(t.chan, val)

"""
    synchronize!(t::TrackingTimer)

Populates `t.results` with any timing results collected so far.
This uses a lock so it is safe to call from multiple threads, but
may cause contention. Called automatically by `Tables.rows(::TrackingTimer)`.
"""
function synchronize!(t::TrackingTimer)
    lock(t.results_lock)
    try
        while isready(t.chan)
            push!(t.results, take!(t.chan))
        end
    finally
        unlock(t.results_lock)
    end
    return nothing
end

"""
    @timeit(t::TrackingTimer, name, expr)

Evaluates `expr` under `@timed`, storing the results in the [`TrackingTimer`](@ref) `t`.
"""
macro timeit(t, name, expr)
    quote
        local results = @timed $(esc(expr))
        put!($(esc(t)),
             (; name=$(esc(name)), results.time, results.gctime,
              n_allocs=Base.gc_alloc_count(results.gcstats), results.bytes,
              thread_id=Threads.threadid(), pid=myid()))
        results.value
    end
end

#####
##### Tables interface
#####

Tables.istable(::TrackingTimer) = true
function Tables.rows(t::TrackingTimer)
    synchronize!(t)
    return t.results
end
Tables.rowaccess(::Type{TrackingTimer}) = true
Tables.columnaccess(::Type{TrackingTimer}) = false

#####
##### Instrumented functions
#####

# implemented as a callable type to get better error messages
# (i.e. you see `F` explictly, which might be `typeof(f)` telling
# you that `f` is involved).
struct InstrumentedFunction{F} <: Function
    func::F
    t::TrackingTimer
    name::String
end

InstrumentedFunction(f, t) = InstrumentedFunction(f, t, string(repr(f)))

function (inst::InstrumentedFunction)(args...; kwargs...)
    @timeit inst.t inst.name inst.func(args...; kwargs...)
end

"""
    (t::TrackingTimer)(f, name=string(repr(f))) -> InstrumentedFunction

Instruments `f` by the [`TrackingTimer`](@ref) `t` returning an `InstrumentedFunction`.
This function can be used just like `f`, but whenever it is called it stores timing
results in `t`.
"""
(t::TrackingTimer)(f, name=string(repr(f))) = InstrumentedFunction(f, t, name)

#####
##### Display
#####

nt_keys(::Type{NamedTuple{K,V}}) where {K,V} = K

function formatter(tbl)
    return function (v, i, j)
        col = nt_keys(TrackingTimerElType)[j]
        col == :time && return @sprintf("%.2f s", v)
        col == :gctime && return string(@sprintf("%.0f", (v / tbl[i].time) * 100), "%")
        col == :bytes && return Base.format_bytes(v)
        return v
    end
end

one_line_show(io::IO, ::TrackingTimer) = print(io, "TrackingTimer(…)")

Base.show(io::IO, t::TrackingTimer) = one_line_show(io, t)

function Base.show(io::IO, ::MIME"text/plain", t::TrackingTimer)
    get(io, :compact, false) === true && return one_line_show(io, t)
    total_time_in_seconds = (time_ns() - t.start_nanosecond) * 1e-9
    tbl = Tables.rows(t)
    tot_measured_time = isempty(tbl) ? 0.0 : sum(x.time for x in tbl)
    println(io, "TrackingTimer: ", @sprintf("%.2f", total_time_in_seconds),
            " s since creation (",
            @sprintf("%.0f", tot_measured_time / total_time_in_seconds * 100),
            "% measured).")
    if isempty(tbl)
        print(io, "No entries.")
    else
        tbl_sorted = sort(tbl; by=x -> x.time, rev=true)
        pretty_table(io, tbl_sorted,
                     ["name", "time", "gctime", "n_allocs", "allocs", "thread ID",
                      "proc ID"]; newline_at_end=false, formatters=formatter(tbl_sorted),
                     hlines=[1], vlines=[],
                     alignment=[:l, (:r for _ in 1:(length(tbl_sorted[1]) - 1))...],
                     header_alignment=:c)
    end
    return nothing
end

end # module

ext: jl

max_stars_repo_name: ericphanson/TrackingTimers

max_issues_repo_name: ericphanson/TrackingTimers

max_forks_repo_name: ericphanson/TrackingTimers


text:
type LinePlot <: Graph
	dataGroups::Vector{DataGroup}
end

LinePlot(dg::DataGroup) = LinePlot([dg])

function plotData(p::PlotFrame, g::LinePlot, showLegend::Bool)
	global dgCount

	for x in 1:length(g.dataGroups)
		dg = g.dataGroups[x]
		if dg.lineColor == "blank"
			dg.lineColor = goodColors[mod(x, length(goodColors))]
		end

		if dg.markerColor == "blank"
			dg.markerColor = goodColors[mod(x, length(goodColors))]
		end

		if dg.markerType == "blank"
			dg.markerType = goodSymbols[mod(x, length(goodSymbols))]
		end

		plotCurve(p, dg)
		if dg.plotPoints
			plotPoints(p, dg)
		end

		if showLegend
			dgCount += 1
			if dg.legend.label == "Data Group "
				dg.legend.label = dg.legend.label*string(dgCount)
			end
		end
	end
end

function lineplot(dgs::Vector{DataGroup})
	plot(PlotFrame(LinePlot(dgs)))
end

function lineplot(dg::DataGroup)
	plot(PlotFrame(LinePlot(dg)))
end

function lineplot(dfs::Vector{DataFrame})
	dgs = DataGroup[]
	for df in dfs
		push!(dgs, DataGroup(df))
	end
	lineplot(dgs)
end

function lineplot(df::DataFrame)
	lineplot(DataGroup(df))
end

function lineplot(ys::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(y, 1)
		push!(dfs, DataFrame(x = collect(1:length(y)), y = y[r,:]))
	end
	lineplot(dfs)
end

function lineplot(xs::Array{Float64, 2}, ys::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(y, 1)
		push!(dfs, DataFrame(x = x[r,:], y = y[r,:]))
	end
	lineplot(dfs)
end

function lineplot(xs::Array{Float64, 2}, ys::Array{Float64, 2}, yerrs::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(ys, 1)
		push!(dfs, DataFrame(x = vec(xs[r,:]), y = vec(ys[r,:]), yerr = vec(yerrs[r,:])))
	end
	lineplot(dfs)
end


text:
module TestCLI
using Test
using FileIO
using DataStructures

include("utils.jl")

# include source of bin file here to avoid world age issues
include("../bin/qxrun.jl")

@testset "Test prepare rqc input cli script" begin
    ghz_example_dir = joinpath(dirname(@__DIR__), "examples", "ghz")
    dsl_input = joinpath(ghz_example_dir, "ghz_5.qx")
    # create empty temporary directory
    mktempdir() do path
        output_fname = joinpath(path, "out.jld2")
        args = ["-d", dsl_input,
                "-o", output_fname]
        main(args)
        @test isfile(output_fname)

        output = load(output_fname, "results")
        @test all([x[2] ≈ ghz_results[x[1]] for x in zip(output...)])
    end

    mktempdir() do path
        output_fname = joinpath(path, "out.jld2")
        args = ["-d", dsl_input,
                "-o", output_fname,
                "--number-amplitudes", "1"]
        main(args)
        @test isfile(output_fname)

        output = load(output_fname, "results")
        @test length(output[1]) == 1
        @test output[2][1] ≈ ghz_results["11001"]
    end

    mktempdir() do path
        output_fname = joinpath(path, "out.jld2")
        args = ["-d", dsl_input,
                "-o", output_fname,
                "--number-amplitudes", "2",
                "--number-slices", "1"]
        main(args)
        @test isfile(output_fname)

        output = load(output_fname, "results")
        @test length(output) == 2
        @test all([output[2][x] ≈ ghz_results[output[1][x]] for x in 1:2])
    end

end
end


text:
using Pkg
Pkg.activate(".")
using PlotlyJS
using Dash, DashHtmlComponents, DashCoreComponents

# Data for solving (dx/dt,dy/dt) = F(x,y), where F=field!
# (x0, y0) is the initial state and
# (tmin, tmax) the time interval for the solution
mutable struct ODEData
    field!::Function
    # Field function formulas
    dxdt::String
    dydt::String
    # Initial conditions
    x0::Float64
    y0::Float64
    # Time interval
    tmin::Float64
    tmax::Float64
end

# Evaluates function in (x,y) from formula
function fxy(x, y, formula::String)
    ex = quote
        x = $x
        y = $y
        x + y
    end
    ex.args[6] = Meta.parse(formula)
    eval(ex)
end

# Parsing ODE data from field formulas
function parseODEData(dxdt, dydt, x0, y0, t0, t1)
    function field!(du, u, p, t)
        du[1] = fxy(u[1], u[2], dxdt)
        du[2] = fxy(u[1], u[2], dydt)
    end
    ODEData(field!, dxdt,dydt,x0, y0, t0, t1)
end

# Solves ODE explicitely using RK4
function solveODE(data::ODEData)
    f = data.field!
    n = 250 # steps count
    zero = [0.0;0.0]
    h = (data.tmax - data.tmin)/n # time step
    k1 = [0.0;0.0]
    k2 = [0.0;0.0]
    k3 = [0.0;0.0]
    k4 = [0.0;0.0]
    # Setting Initial conditions
    ps = [[0.0;0.0] for i = 1:(n+1)]
    ps[1] = [data.x0;data.y0]
    for i = 2:(n+1)
        p = ps[i-1]
        f(k1,p,zero,0.0)
        f(k2,p + (h/2.0)*k1,zero,0.0)
        f(k3,p + (h/2.0)*k2,zero,0.0)
        f(k4,p + h*k3,zero,0.0)
        ps[i] = p + (h/6.0)*(k1 + 2*k2 +2*k3 + k4)
    end
    ps
end

# Gets the Plotly Scatter plot from the solution to the ODE
function trace(sol)
    x = [u[1] for u in sol]
    y = [u[2] for u in sol]
    scatter(;x=x, y=y, mode="lines", name="Solution")
end

examples = Dict(
    "Spring"=>parseODEData("y","-x",3.0,0.0,0.0,2*3.1417),
    "DampedSpring" => parseODEData("y", "-x-y/3", 0.0, 4.0, 0.0, 30.0),
    "Pendulum"=>parseODEData("y", "-sin(x)", 3.0, 0.0, 0.0, 16),
    "Predator-Prey"=>parseODEData("x-x*y","x*y-y",2.0,3.0,0.0,8.0)
    )

default_example="Spring"
layout = Layout(;
    title = "",
    xaxis_title="x",
    yaxis_title="y",
    xaxis_range=[-5, 5],
    yaxis_range=[-5, 5], width=550,height=500)
#
# Creation of Dash App starts here
#
#app = dash(external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"])

mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js?config=TeX-MML-AM_CHTML"
app = dash(external_scripts=[mathjax])
#app = dash()
app.title = "DashODE"
# Utility fuction to retrieve inputs
function odeinput(id, label; value="")
    html_div(
    [dcc_input(id=id, type="text", value=value),html_span("  "*label)])
end

# Utility function to do n html breaks
function html_break(n)
    html_div([html_br() for i = 1:n])
end

odeoptions = [(label = k, value = k) for k in keys(examples)]

function ode_graph(sol)
    dcc_graph(id="odedash", figure=plot([trace(sol)], layout))
end

function ode_all_inputs()
    e = examples[default_example]
    [html_break(3),
    odeinput("dxdt",raw"\(dx/dt\)",value=e.dxdt),
    odeinput("dydt",raw"\(dy/dt\)",value=e.dydt),
    odeinput("x0",raw"\(x_0\)",value=string(e.x0)),
    odeinput("y0",raw"\(y_0\)",value=string(e.y0)),
    odeinput("time",raw"\(t_0\)",value=string(e.tmax)),
    html_break(1),
    html_button("Solve", id="solve-button", n_clicks=0),
    html_break(1),
    html_h6("Load ODE"),
    html_div(style=Dict("max-width" => "200px"),
    dcc_dropdown(id="examples", options=odeoptions, value=default_example))]
end

title = "Solving Differential Equations Interactively"

description = raw"
This [Dash Julia](https://dash-julia.plotly.com/) application is designed to interactively solve
[differential equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation) in two variables  $x=x(t)$ and $y=y(t)$ during a time interval $T$, and with initial conditions
$$x(0)=x_0, y(0)=y_0.$$
In particular, we will consider the following examples:

  1. ***Mass on a spring***. In this case $x$ stands for the displacement of the mass $m$ from its equilibrium position, and for $y$ the velocity.
  The equations modelling this oscillating system are
  $$ \frac{dx}{dt} = y, m \frac{dy}{dt} = -k x,$$
  where $k$ stands for the stiffness of the spring.
  2. ***Mass on a spring with damping***. If we add to the above equations a damping force that depends linearly on velocity, we obtain
  $$ \frac{dx}{dt} = y, m \frac{dy}{dt} = -k x - \alpha y.$$
  3. The [pendulum equation](https://en.wikipedia.org/wiki/Pendulum_(mechanics)) describes the motion of a mass suspended from a fixed support and under the influence of gravity. In this case
  $x$ is the angle with respect to vertical direction and $y$ is the angular velocity. The nonlinear differential equation associated to this system is
  $$ \frac{dx}{dt} = y, \frac{dy}{dt} = -(g/l)\sin(x),$$
  where $g$ is the acceleration of gravity and $l$ is the distance from the point of attachment to the location of the mass.
  4. The [predator-prey equations](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations) model the interaction of two species (whose populations are measued by $x$ and $y$) in which one is a predator and the other is a prey.
  The corresponding equations are
  $$ \frac{dx}{dt} = \alpha x - \beta x y, \frac{dy}{dt} = \delta x y - \gamma y$$

"
app_instructions = raw"To explore the above equations interactively use the drop down menu below to select one of them, and then press the solve button. You can
modify these equations manually; try changing their parameters to visualize the effect on their solution. You can also create your own differential equations from scratch by entering the necessary data in the corresponding input fields"
# Setting the app's layout: ODE graph to the left and controls to the right.

mathjax = html_script(type="text/javascript", async=true, id="MathJax-script", src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js")

app.layout = html_div(
    [html_center([html_h1("DashODE"),
    html_h2(title)]),
    html_div(style=Dict("max-width" => "1000px"),dcc_markdown(description)),
    html_h2("Interactive application"),
    html_div(style=Dict("max-width" => "1000px"),dcc_markdown(app_instructions)),
    html_div([
        html_div([ode_graph(solveODE(examples[default_example])),
        ], className="six columns"),
        html_div(ode_all_inputs(), className="six columns"),
    ], className="row"), html_break(length(examples)) # Making space for dropdown
])

# Update ODE data from dropdown
callback!(app,Output("dxdt","value"), Output("dydt","value"),
              Output("x0","value"), Output("y0","value"), Output("time", "value"),
              Input("examples", "value")) do value
    e = examples[value]
    return (e.dxdt, e.dydt, string(e.x0), string(e.y0), string(e.tmax))
end

# Update figure fro ODE data
callback!(
    app,
    Output("odedash", "figure"),
    Input("solve-button", "n_clicks"),
    State("x0", "value"),
    State("y0", "value"),
    State("dxdt", "value"),
    State("dydt", "value"),
    State("time", "value"),
) do clicks, x0, y0, dxdt, dydt, time
    x0f = parse(Float64,x0)
    y0f = parse(Float64,y0)
    odedata = parseODEData(dxdt, dydt, x0f, y0f , 0.0, parse(Float64,time))
    sol = solveODE(odedata)
    p0G = scatter(;x=[x0f], y=[y0f], mode="markers",name="(x0,y0)")
    eq = raw"$\frac{dx}{dt} ="*dxdt*raw", "*raw"\frac{dy}{dt} ="*dydt*raw"$"
    eqG = scatter(;x=0, y = [-4.4], mode="text", text=eq, name="Equation", textfont_size=15)
    return plot([trace(sol),p0G,eqG], layout)
end
run_server(app, "0.0.0.0", debug=true)


text:
"A `Literal` is given by a symbol and a depth > 0, defaulting to 1."
immutable Literal
    symbol::Symbol
    depth::Int
end
Literal(tup::Tuple{Symbol,Int}) = Literal(tup[1], tup[2])
Literal(sym::Symbol) = Literal(sym, 1)

function Base.string(literal::Literal)
    if literal.depth == 1
        return string(literal.symbol)
    else
        return string(literal.symbol) * string(literal.depth)
    end
end

Base.isless(a::Literal, b::Literal) =
    (a.symbol<b.symbol) || ((a.symbol==b.symbol) && (a.depth<b.depth))

"A `PathKey` is a `Vector` of `Literal`s."
immutable PathKey
    literals::Vector{Literal}
    PathKey() = new(Literal[])
    PathKey(pathkey::PathKey) = pathkey
    PathKey(tuple_args::Tuple) = PathKey(tuple_args...)
    PathKey(literals::Vector{Literal}) = new(literals)
    function PathKey(args...)
        literals = Literal[]
        idarg = 1
        nargs = length(args)
        while idarg <= nargs
            if (idarg+1)<=nargs &&
                    isa(args[idarg], Symbol) && isa(args[idarg+1], Int)
                push!(literals, Literal(args[idarg], args[idarg+1]))
                idarg += 2
            else
                push!(literals, Literal(args[idarg]))
                idarg += 1
            end
        end
        new(literals)
    end
end

Base.:(==)(x::PathKey, y::PathKey) = (x.literals == y.literals)

Base.convert(::Type{PathKey}, sym::Symbol) = PathKey(sym)
Base.convert(::Type{PathKey}, tup::Tuple) = PathKey(tup...)

function Base.isless(a::PathKey, b::PathKey)
    for (x, y) in zip(a.literals, b.literals)
        if (x.symbol > y.symbol) || (x.depth > y.depth)
            return false
        end
    end
    if length(a.literals) > length(b.literals)
        return false
    end
    return a.literals != b.literals
end

function Base.string(pathkey::PathKey)
    return join([string(literal) for literal in pathkey.literals], "_")
end

Base.Symbol(pathkey::PathKey) = Symbol(string(pathkey))

prepend(prefix, pathkey::PathKey) =
    PathKey(vcat(Literal(prefix), pathkey.literals))

"""A `Path` is a sorted dictionary whose keys are `PathKey`s and whose
values are integer indices."""
immutable Path
    sdict::DataStructures.SortedDict{PathKey,Int}
    function Path(pairs...)
        sdict = DataStructures.SortedDict(Pair{PathKey,Int}[])
        for pair in pairs
            push!(sdict, PathKey(pair.first) => pair.second)
        end
        new(sdict)
    end
end

Base.:(==)(x::Path, y::Path) = reduce(&, map(==, x.sdict, y.sdict))

function Base.isless(a::Path, b::Path)
    for (x, y) in zip(a.sdict, b.sdict)
        if (x.first < y.first)
            return true
        elseif (x.first > y.first)
            return false
        else
            if (x.second < y.second)
                return true
            elseif (x.second > y.second)
                return false
            end
        end
    end
    return length(a.sdict) < length(b.sdict)
end

function scattering_order(path::Path, suffix::PathKey)
    pathkeys = collect(keys(path.sdict))
    order = 0
    for idkey in eachindex(pathkeys)
        if pathkeys[idkey].literals[2:end] == suffix.literals
            order = max(order, pathkeys[idkey].literals[1].depth)
        end
    end
    return order
end

"""A `PathRange` is a sorted dictionary whose keys are `PathKey`'s
and whose values are integer ranges."""
immutable PathRange
    sdict::DataStructures.SortedDict{PathKey,StepRange{Int,Int}}
    function PathRange(pairs...)
        sdict = DataStructures.SortedDict(Pair{PathKey,StepRange{Int,Int}}[])
        for pair in pairs
            rhs = pair.second
            if isa(rhs, Int)
                push!(sdict, PathKey(pair.first) => rhs:1:rhs)
            elseif isa(pair.second, UnitRange{Int})
                push!(sdict, PathKey(pair.first) => rhs.start:1:rhs.stop)
            elseif isa(pair.second, StepRange{Int,Int})
                push!(sdict, PathKey(pair.first) => rhs)
            end
        end
        new(sdict)
    end
end

Base.:(==)(x::PathRange, y::PathRange) = reduce(&, map(==, x.sdict, y.sdict))


text:
# Fractals
# TODO: super-fat Cantor sets, asymmetric Cantor sets...

export cantor, thincantor, thinnercantor, thinnestcantor, smithvolterracantor

# Standard Cantor set removes the middle third at every level

for set in (:cantor,:thincantor,:thinnercantor,:thinnestcantor)
    @eval begin
        $set(d::Domain,n::Int) = $set(d,n,3)

        function $set(d::Circle{T,V},n::Int,α::Number) where {T,V}
            c,r = d.center,d.radius
            α = convert(promote_type(real(T),V,typeof(α)),α)
            if n == 0
                return d
            else
                C = $set(Segment{T}(),n,α)
                return UnionDomain(map(d->Arc(c,r,(d.a+1/2α)π,(d.b+1/2α)π),
                                       components(C)[1:(ncomponents(C)÷2)])) ∪
                        UnionDomain(map(d->Arc(c,r,(d.a-1/2α)π,(d.b-1/2α)π),
                                        components(C)[(ncomponents(C)÷2)+1:ncomponents(C)]))
            end
        end

    end
end

# α is width, n is number of levels

function cantor(d::Segment{T},n::Int,α::Number) where T
    a,b = d.a,d.b
    if n == 0
        return d
    else
        C = Segment{T}(zero(T),one(T))
        for k=n:-1:1
            C = C/α ∪ (α-1+C)/α
        end
        return a+(b-a)*C
    end
end

# Thin Cantor set removes the middle n/(n+2)th at the nth level

function thincantor(d::Segment{T},n::Int,α::Number) where T
    a,b = d.a,d.b
    if n == 0
        return d
    else
        C = Segment{T}(zero(T),one(T))
        for k=n:-1:1
            C = C/(α+k-1) ∪ (α+k-2+C)/(α+k-1)
        end
        return a+(b-a)*C
    end
end

function thinnercantor(d::Segment{T},n::Int,α::Number) where T
    a,b = d.a,d.b
    if n == 0
        return d
    else
        C = Segment{T}(zero(T),one(T))
        for k=n:-1:1
            C = C/α^k ∪ (α^k-1+C)/α^k
        end
        return a+(b-a)*C
    end
end

function thinnestcantor(d::Segment{T},n::Int,α::Number) where T
    a,b = d.a,d.b
    if n == 0
        return d
    else
        C = Segment{T}(zero(T),one(T))
        for k=n:-1:1
            C = C/α^(2^(k-1)) ∪ (α^(2^(k-1))-1+C)/α^(2^(k-1))
        end
        return a+(b-a)*C
    end
end

smithvolterracantor(d::Domain{T},n::Int) where {T} = smithvolterracantor(d,n,one(T)/4,one(T)/4)

function smithvolterracantor(d::Segment{T},n::Int,α::Number,β::Number) where {T}
    a,b = d.a,d.b
    if n == 0
        return d
    else
        C = Segment{T}(zero(T),one(T))
        ln = (one(T)-β*(one(T)-(2α)^n)/(one(T)-2α))/2^n
        lnm1 = (2ln+β*α^(n-1))
        for k=n:-1:1
            C = (C*ln ∪ (lnm1+(C-1)*ln))/lnm1
            ln,lnm1 = lnm1,(2lnm1+β*α^(k-2))
        end
        return a+(b-a)*C
    end
end


text:
function device_linear_rateofchange(ps_m::CanonicalModel,
                                    rate_data::Tuple{Vector{String}, Vector{UpDown}},
                                    initial_conditions::Vector{InitialCondition},
                                    time_range::UnitRange{Int64},
                                    cons_name::Symbol,
                                    var_name::Symbol)

    up_name = Symbol(cons_name,:_up)
    down_name = Symbol(cons_name,:_down)

    set_name = rate_data[1]

    ps_m.constraints[up_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[down_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

    for (ix, name) in enumerate(rate_data[1])
        ps_m.constraints[up_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, 1] - initial_conditions[ix].value <= rate_data[2][ix].up)
        ps_m.constraints[down_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, initial_conditions[ix].value - ps_m.variables[var_name][name, 1] <= rate_data[2][ix].down)
    end

    for t in time_range[2:end], (ix, name) in enumerate(rate_data[1])
        ps_m.constraints[up_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, t-1] - ps_m.variables[var_name][name, t] <= rate_data[2][ix].up)
        ps_m.constraints[down_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, t] - ps_m.variables[var_name][name, t-1] <= rate_data[2][ix].down)
    end

    return

end

function device_mixedinteger_rateofchange(ps_m::CanonicalModel,
                                            rate_data::Tuple{Vector{String}, Vector{UpDown}, Vector{MinMax}},
                                            initial_conditions::Vector{InitialCondition},
                                            time_range::UnitRange{Int64},
                                            cons_name::Symbol,
                                            var_names::Tuple{Symbol,Symbol,Symbol})

    up_name = Symbol(cons_name,:_up)
    down_name = Symbol(cons_name,:_down)

    set_name = rate_data[1]

    ps_m.constraints[up_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[down_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

    for (ix, name) in enumerate(rate_data[1])
        ps_m.constraints[up_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[1]][name, 1] - initial_conditions[ix].value <= rate_data[2][ix].up + rate_data[3][ix].max*ps_m.variables[var_names[2]][name, 1])
        ps_m.constraints[down_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, initial_conditions[ix].value - ps_m.variables[var_names[1]][name, 1] <= rate_data[2][ix].down + rate_data[3][ix].min*ps_m.variables[var_names[3]][name, 1])
    end

    for t in time_range[2:end], (ix, name) in enumerate(rate_data[1])
        ps_m.constraints[up_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[1]][name, t-1] - ps_m.variables[var_names[1]][name, t] <= rate_data[2][ix].up + rate_data[3][ix].max*ps_m.variables[var_names[2]][name, t])
        ps_m.constraints[down_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[1]][name, t] - ps_m.variables[var_names[1]][name, t-1] <= rate_data[2][ix].down + rate_data[3][ix].min*ps_m.variables[var_names[3]][name, t])
    end

    return

end

#Old implementation of initial_conditions
#=
function device_linear_rateofchange(ps_m::CanonicalModel,
                                    rate_data::Vector{UpDown},
                                    initial_conditions::Vector{Float64},
                                    time_range::UnitRange{Int64},
                                    cons_name::Symbol,
                                    var_name::Symbol)

    up_name = Symbol(cons_name,:_up)
    down_name = Symbol(cons_name,:_down)

    set_name = [name for r in rate_data]

    ps_m.constraints[up_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[down_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

    for (ix,r) in enumerate(rate_data)
        ps_m.constraints[up_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, 1] - initial_conditions[ix] <= rate_data[2][ix].up)
        ps_m.constraints[down_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, initial_conditions[ix] - ps_m.variables[var_name][name, 1] <= rate_data[2][ix].down)
    end

    for t in time_range[2:end], r in rate_data
        ps_m.constraints[up_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, t-1] - ps_m.variables[var_name][name, t] <= rate_data[2][ix].up)
        ps_m.constraints[down_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, t] - ps_m.variables[var_name][name, t-1] <= rate_data[2][ix].down)
    end

    return

end

function device_mixedinteger_rateofchange(ps_m::CanonicalModel,
                                            rate_data::Array{Tuple{String,NamedTuple{(:up, :down),Tuple{Float64,Float64}},NamedTuple{(:min, :max),Tuple{Float64,Float64}}},1},
                                            initial_conditions::Vector{Float64},
                                            time_range::UnitRange{Int64},
                                            cons_name::Symbol,
                                            var_names::Tuple{Symbol,Symbol,Symbol})

    up_name = Symbol(cons_name,:_up)
    down_name = Symbol(cons_name,:_down)

    set_name = [name for r in rate_data]

    ps_m.constraints[up_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[down_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

    for (ix,r) in enumerate(rate_data)
        ps_m.constraints[up_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[1]][name, 1] - initial_conditions[ix] <= rate_data[2][ix].up + rate_data[3][ix].max*ps_m.variables[var_names[2]][name, 1])
        ps_m.constraints[down_name][name, 1] = JuMP.@constraint(ps_m.JuMPmodel, initial_conditions[ix] - ps_m.variables[var_names[1]][name, 1] <= rate_data[2][ix].down + rate_data[3][ix].min*ps_m.variables[var_names[3]][name, 1])
    end

    for t in time_range[2:end], r in rate_data
        ps_m.constraints[up_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[1]][name, t-1] - ps_m.variables[var_names[1]][name, t] <= rate_data[2][ix].up + rate_data[3][ix].max*ps_m.variables[var_names[2]][name, t])
        ps_m.constraints[down_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[1]][name, t] - ps_m.variables[var_names[1]][name, t-1] <= rate_data[2][ix].down + rate_data[3][ix].min*ps_m.variables[var_names[3]][name, t])
    end

    return

end
=#


text:
module TestSimplePass

const SCRIPT = joinpath(@__DIR__, "../../samples/simplepass/runtests.jl")

test() = include(SCRIPT)

end  # module


text:
module TestSetExtensions

export ExtendedTestSet, @includetests

using Compat.Test
import Compat.Test: record, finish
using Compat.Test: DefaultTestSet, AbstractTestSet
using Compat.Test: get_testset_depth, scrub_backtrace
using Compat.Test: Result, Pass, Fail, Error

using Base: @deprecate
@deprecate DottedTestSet ExtendedTestSet

using DeepDiffs

"""
Includes the given test files, given as a list without their ".jl" extensions.
If none are given it will scan the directory of the calling file and include all
the julia files.
"""
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = @__FILE__
        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(f->endswith(f, ".jl") && f!= basename(rootfile), tests)
        else
            tests = map(f->string(f, ".jl"), tests)
        end
        println();
        for test in tests
            print(splitext(test)[1], ": ")
            include(test)
            println()
        end
    end
end

struct ExtendedTestSet{T<:AbstractTestSet} <: AbstractTestSet
    wrapped::T

    ExtendedTestSet{T}(desc) where {T} = new(T(desc))
end

struct FailDiff <: Result
    result::Fail
end

function ExtendedTestSet(desc; wrap=DefaultTestSet)
    ExtendedTestSet{wrap}(desc)
end

function record(ts::ExtendedTestSet, res::Fail)
    if myid() == 1
        println("\n=====================================================")
        print_with_color(:white, ts.wrapped.description, ": ")
        if res.test_type == :test && isa(res.data,Expr) && res.data.head == :comparison
            dd = deepdiff(res.data.args[1], res.data.args[3])
            if !isa(dd, DeepDiffs.SimpleDiff)
                # The test was an comparison between things we can diff,
                # so display the diff
                print_with_color(Base.error_color(), "Test Failed\n"; bold = true)
                println("  Expression: ", res.orig_expr)
                print_with_color(Base.info_color(), "\nDiff:\n")
                display(dd)
                println()
            else
                # fallback to the default printing if we don't have a pretty diff
                print(res)
            end
        else
            # fallback to the default printing for non-comparisons
            print(res)
        end
        Base.show_backtrace(STDOUT, scrub_backtrace(backtrace()))
        # show_backtrace doesn't print a trailing newline
        println("\n=====================================================")
    end
    push!(ts.wrapped.results, res)
    res, backtrace()
end

function record(ts::ExtendedTestSet, res::Error)
    println("\n=====================================================")
    record(ts.wrapped, res)
    print("=====================================================\n")
end

function record(ts::ExtendedTestSet, res::Pass)
    print_with_color(:green, ".")
    record(ts.wrapped, res)
    res
end

record(ts::ExtendedTestSet, res) = record(ts.wrapped, res)

function finish(ts::ExtendedTestSet)
    get_testset_depth() == 0 && print("\n\n")
    finish(ts.wrapped)
end

end # module


text:
using VegaLite, DataFrames, Query, VegaDatasets

cars = dataset("cars")

cars |> @select(:Acceleration, :Name) |> collect

function foo(data, origin)
    df = data |> @filter(_.Origin==origin) |> DataFrame

    return df |> @vlplot(:point, :Acceleration, :Miles_per_Gallon)
end

p = foo(cars, "USA")

p |> save("foo.png")


text:
import StaticArrays: SUnitRange

# Computing Integrals Involving the Matrix Exponential (elements of Van Loan '78)
function _integrate_expAt_B(A::StaticMatrix{N,N}, B::StaticMatrix{N,M}, dt::T) where {N,M,T}
    S = promote_type(eltype(A), eltype(B), T)
    X = exp([[A B]; zeros(SMatrix{M,N+M,S})]*dt)
    (eᴬᵗ=X[SUnitRange(1,N),SUnitRange(1,N)], ∫eᴬᵗB=X[SUnitRange(1,N),SUnitRange(N+1,N+M)])
end
function integrate_expAt_B(A::StaticMatrix{N,N}, B::StaticVector{N}, dt) where {N}
    ans = _integrate_expAt_B(A, SMatrix{N,1}(B), dt)
    (eᴬᵗ=ans.eᴬᵗ, ∫eᴬᵗB=SVector(ans.∫eᴬᵗB))
end
function integrate_expAt_B(A::StaticMatrix{N,N}, B::StaticMatrix{N,M}, dt) where {N,M}
    eᴬᵗ, ∫eᴬᵗB₁ = integrate_expAt_B(A, B[:,1], dt)
    (eᴬᵗ=eᴬᵗ, ∫eᴬᵗB=hcat(∫eᴬᵗB₁, ntuple(i -> integrate_expAt_B(A, B[:,i+1], dt).∫eᴬᵗB, Val(M-1))...))
end

function _integrate_expAt_Bt_dtinv(A::StaticMatrix{N,N}, B::StaticMatrix{N,M}, dt::T) where {N,M,T}
    S = promote_type(eltype(A), eltype(B), T)
    X = exp([[A B zeros(typeof(B))]*dt;
             [zeros(SMatrix{M,N+M,S}) SMatrix{M,M,S}(I)];
             zeros(SMatrix{M,N+2*M,S})])
    ∫eᴬᵗB=X[SUnitRange(1,N),SUnitRange(N+1,N+M)]
    (eᴬᵗ=X[SUnitRange(1,N),SUnitRange(1,N)], ∫eᴬᵗB=∫eᴬᵗB, ∫eᴬᵗBtdt⁻¹=∫eᴬᵗB - X[SUnitRange(1,N),SUnitRange(N+M+1,N+2*M)])
end
function integrate_expAt_Bt_dtinv(A::StaticMatrix{N,N}, B::StaticVector{N}, dt) where {N}
    ans = _integrate_expAt_Bt_dtinv(A, SMatrix{N,1}(B), dt)
    (eᴬᵗ=ans.eᴬᵗ, ∫eᴬᵗB=SVector(ans.∫eᴬᵗB), ∫eᴬᵗBtdt⁻¹=SVector(ans.∫eᴬᵗBtdt⁻¹))
end
function integrate_expAt_Bt_dtinv(A::StaticMatrix{N,N}, B::StaticMatrix{N,M}, dt) where {N,M}
    eᴬᵗ, ∫eᴬᵗB₁, ∫eᴬᵗB₁tdt⁻¹ = integrate_expAt_Bt_dtinv(A, B[:,1], dt)
    X₁ = [∫eᴬᵗB₁; ∫eᴬᵗB₁tdt⁻¹]
    X = hcat(X₁, ntuple(i -> (ans = integrate_expAt_Bt_dtinv(A, B[:,i+1], dt); [ans.∫eᴬᵗB; ans.∫eᴬᵗBtdt⁻¹]), Val(M-1))...)
    (eᴬᵗ=eᴬᵗ, ∫eᴬᵗB=X[SUnitRange(1,N),:], ∫eᴬᵗBtdt⁻¹=X[SUnitRange(N+1,2*N),:])
end

function integrate_expAt_B_expATt(A::StaticMatrix{N,N}, B::StaticMatrix{N,N}, dt) where {N}
    X = exp([[-A B]; [zeros(typeof(A)) A']]*dt)
    X[SUnitRange(N+1,2*N),SUnitRange(N+1,2*N)]'*X[SUnitRange(1,N),SUnitRange(N+1,2*N)]
end

# Univariate Optimization (@inline here and in optimal_time necessary to eliminate allocations)
@inline function golden_section(f, a, b; ϵ=1e-3)
    T = typeof((a + b)/2)
    gr = T((1 + sqrt(5))/2)
    fa, fb = f(a), f(b)
    while fa < fb; a /= 2; fa = f(a); end
    m1 = b - (b - a)/gr
    m2 = a + (b - a)/gr
    while abs(m1 - m2) > ϵ
        if f(m1) < f(m2)
            b = m2
        else
            a = m1
        end
        m1 = b - (b - a)/gr
        m2 = a + (b - a)/gr
    end
    (a + b)/2
end

@inline function bisection(f, a, b; ϵ=1e-3)
    # Bisection
    # Broken for Float32 triple integrator with R = [1f-3],
    # q0 = SVector(-6.094263f0, 0.014560595f0, -0.6846263f0)
    # qf = SVector(-6.0940447f0, 0.0f0, 0.0f0)
    f(b) < 0 && return b
    fa = f(a)
    while fa > 0
        a /= 2
        fa = f(a)
        (a == 0 || isnan(fa) || isinf(fa)) && return nothing    # workaround for above
    end
    m = (a + b)/2
    fm = f(m)
    while abs(fm) > ϵ && abs(a - b) > ϵ
        fm > 0 ? b = m : a = m
        m = (a + b)/2
        fm = f(m)
    end
    m
end

@inline function false_position(f, a, b; ϵ=1e-3)
    f(b) < 0 && return b
    fa = f(a)
    while fa > 0
        a /= 2
        fa = f(a)
        (a == 0 || isnan(fa) || isinf(fa)) && return nothing
    end
    fb = f(b)
    x = (a*fb - b*fa)/(fb - fa)
    fx = f(x)
    while abs(fx) > ϵ && abs(a - b) > ϵ
        fx > 0 ? b = x : a = x
        x = (a*fb - b*fa)/(fb - fa)
        fx = f(x)
    end
    x
end

@inline function newton(f, f′, a, b; ϵ=1e-6)
    # Bisection / Newton's method combo
    f(b) < 0 && return b
    fa = f(a)
    while fa > 0
        a /= 2
        fa = f(a)
        (a == 0 || isnan(fa) || isinf(fa)) && return nothing
    end
    x = (a + b)/2
    fx = f(x)
    while abs(fx) > ϵ && abs(a - b) > ϵ
        x = x - fx / f′(x)
        (x < a || x > b) && (x = (a + b)/2)
        fx = f(x)
        fx > 0 ? b = x : a = x
    end
    x
end


text:
# Similar to JuMP/test/solvers.jl

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

grb = false && try_import(:Gurobi)
isgrb(solver) = occursin("GurobiSolver", string(typeof(solver)))
cpx = false && try_import(:CPLEX)
iscpx(solver) = occursin("CplexSolver", string(typeof(solver)))
xpr = false && try_import(:Xpress)
clp = false && try_import(:Clp)
isclp(solver) = occursin("ClpSolver", string(typeof(solver)))
glp = try_import(:GLPKMathProgInterface)
msk = false && try_import(:Mosek)
ismsk(solver) = occursin("MosekSolver", string(typeof(solver)))

lp_solvers = Any[]
grb && push!(lp_solvers, Gurobi.GurobiSolver(OutputFlag=0))
cpx && push!(lp_solvers, CPLEX.CplexSolver(CPX_PARAM_SCRIND=0, CPX_PARAM_REDUCE=0))
xpr && push!(lp_solvers, Xpress.XpressSolver(OUTPUTLOG=0))
clp && push!(lp_solvers, Clp.ClpSolver())
glp && push!(lp_solvers, GLPKMathProgInterface.GLPKSolverLP())
msk && push!(lp_solvers, Mosek.MosekSolver(QUIET=true))


text:
using NeidArchive
using Dates
last_day_retreived_solar_data = Date(2021,3,19)

include("password.jl")  # Needs to set user_nexsci and passwd_nexsci
cookie_nexsci = "./neidadmincookie.txt"
 query_result_file = "./criteria.csv"

NeidArchive.login(userid=user_nexsci, password=passwd_nexsci, cookiepath=cookie_nexsci)

param = Dict{String,String}()
 param["datalevel"] = "solarl1"
 param["piname"] = "Mahadevan"
 param["object"] = "Sun"
 param["datetime"] = NeidArchive.datetime_range_after(last_day_retreived_solar_data)

NeidArchive.query(param, cookiepath=cookie_nexsci, outpath=query_result_file)
num_lines = countlines(query_result_file) - 1
println("# Query resulted in file with ", num_lines, " entries.")

#= If you want to parse the file
using CSV, DataFrames
if (@isdefined CSV) && (@isdefined DataFrames)
 df = CSV.read(query_result_file, DataFrame)
 nfiles = size(df,1)
 @assert nfiles >= 1
end
=#


text:
using Mdp.PricingProblem
using Mdp.PricingAlgorithmFH

"""
Model of test problem:

3 nodes, 2 edges: graph:      A---B---C
capacity:     1   1
sell end:     3   5
"""
P = Problem(
    [0:1, 0:1],             # Capacity of edges
    [3, 5],                 # Selling period end of edges
    5,                      # Number of timesteps (Start at 1)
    [10, 20],               # Actions (prices)
    [(1,), (2,), (1, 2)],   # Products (seqeuences of edge indeces)
    Dict((1,) => 0.5,       # λ: Dictionary of demand intensities for products
        (2,) => 0.5,
        (1, 2) => 0.4),
)

"""
Model of test problem 2:

2 nodes, 1 edge: graph:      A---B
capacity:     1
sell end:     5
"""
P_2 = Problem(
    [0:1,],                 # Capacity of edges
    [5,],                   # Selling period end of edges
    5,                      # Number of timesteps (Start at 1)
    [10, 20],               # Actions (prices)
    [(1,),],                # Products (seqeuences of edge indeces)
    Dict((1,) => 0.5)       # λ: Dictionary of demand intensities for products
)

@testset "Q, π, V" begin
    @testset "Test Problem 1" begin
        # See docs/FH_mdp_testcase.jpg
        @testset "V_next=0" begin
            V_next = zeros(Float64, map(length, P.C)...)
            k = 3#P.N
            Q = get_Q(k, V_next, P)
            # Q array dimensions are [c1, a, c2]
            @test Q == cat([0.0 1.3333333333333333;
                            0.0 1.3333333333333333;
                            0.0 0.0],
                            [0.8000000000000002 3.2;
                            0.8000000000000002 3.2;
                            0.0 0.0], dims=3)
            @test get_V(Q) == [0.0 0.8000000000000002;
                           1.3333333333333333 3.2]
            @test get_π(Q) == [1 1; 1 1]

            k=5 # difference from selling period end for products (1,) and (1,2)
            Q = get_Q(k, V_next, P)
            @test Q == cat([0.0 0.0;
                            0.0 0.0;
                            0.0 0.0],
                            [0.8000000000000002 0.8000000000000002;
                            0.8000000000000002 0.8000000000000002;
                            0.0 0.0], dims=3)
            @test get_V(Q) == [0.0 0.8000000000000002;
                          0.0 0.8000000000000002]
            @test get_π(Q) == [1 1; 1 1]
        end

        @testset "V_next=10" begin
            V_next = ones(Float64, map(length, P.C)...).*10
            k = 3
            #TODO: The action 2 is caused by numerical instability!
            Q = get_Q(k, V_next, P)
            @test Q == cat([10.0 11.333333333333334;
                            10.0 11.333333333333334;
                            10.0 10.0],
                            [10.799999999999999 13.199999999999998;
                            10.8 13.200000000000001;
                            10.0 10.0], dims=3)
            @test get_V(Q) == [0.0 0.8000000000000002;
                          1.3333333333333333 3.200000000000001].+10.
            @test get_π(Q) == [1 2; 1 2]

            k = 5
            Q = get_Q(k, V_next, P)
            @test Q == cat([10.0 10.0;
                            10.0 10.0;
                            10.0 10.0],
                            [10.799999999999999 10.799999999999999;
                            10.8 10.8;
                            10.0 10.0], dims=3)
            @test get_V(Q) == [0.0 0.8000000000000002;
                          0.0 0.8000000000000002].+10.
            @test get_π(Q) == [1 2; 1 2]
        end
    end
    @testset "Test Problem 2" begin
    end
end

@testset "state2id" begin
    @test state2id([0,1])==[1,2]
    @test state2id([0,0])==[1,1]
end

@testset "get_R" begin
    R_s = state_rewards(P.A)
    @test get_R([0,1], [0,1], R_s)==[0.,0.,0.]
    @test get_R([0,0], [0,0], R_s)==[0.,0.,0.]
    @test get_R([0,1], [0,0], R_s)==[P.A...,0.]
    @test get_R([1,0], [0,0], R_s)==[P.A...,0.]
    @test get_R([1,1], [0,1], R_s)==[P.A...,0.]
end

@testset "state_transitions" begin
    """
    These probabilities are used to set up the following testcases.
    """
    k = 3
    @test state_transitions([0,0],k, P.A, P.products, P) == ([[0,0]],
        reshape([1,1,1], 3,1))
    @test state_transitions([0,1],k, P.A, P.products, P) == ([[0,1],[0,0]],
        [0.9199999999999999 0.08000000000000002;
        0.96 0.04000000000000001;
        1.0 0.0])
    @test state_transitions([1,0],k, P.A, P.products, P) == ([[1,0],[0,0]],
        [0.8666666666666667 0.13333333333333333;
        0.9333333333333333 0.06666666666666667;
        1.0 0.0])
    @test state_transitions([1,1],k, P.A, P.products, P) == ([[1,1],[0,1],[1,0],[0,0]],
        [0.6799999999999999 0.13333333333333333 0.08000000000000002 0.10666666666666667;
        0.84 0.06666666666666667 0.04000000000000001 0.05333333333333334;
        1.0 0.0 0.0 0.0])

    k = 4
    @test state_transitions([0,0],k, P.A, P.products, P) == ([[0,0]],
        reshape([1,1,1], 3,1))
    @test state_transitions([0,1],k, P.A, P.products, P) == ([[0,1],[0,0]],
        [0.9199999999999999 0.08000000000000002;
        0.96 0.04000000000000001;
        1.0 0.0])
    @test state_transitions([1,0],k, P.A, P.products, P) == ([[1,0],[0,0]],
        [1. 0.0;
        1. 0.0;
        1.0 0.0])
    @test state_transitions([1,1],k, P.A, P.products, P) == ([[1,1],[0,1],[1,0],[0,0]],
        [0.91999999999999998 0.0 0.08000000000000002 0.0;
        0.95999999999999999 0.0 0.04000000000000001 0.0;
        1.0 0.0 0.0 0.0])
end

@testset "get_next_states" begin
    @test get_next_states([0,0], P.products) == ([()],
        [[0,0]])
    @test get_next_states([0,1], P.products) == ([(), (2,)],
        [[0,1],[0,0]])
    @test get_next_states([1,0], P.products) == ([(), (1,)],
        [[1,0],[0,0]])
    @test get_next_states([1,1], P.products) == ([(), (1,), (2,), (1,2)],
        [[1,1],[0,1],[1,0],[0,0]])
end


text:
function elementpartition2(dmd,t,t2t,nproc,metis)

nve,ne = size(t);

if nproc==1
    i = 1;
    dmd[i].elempart = reshape(1:ne,1,ne);
    dmd[i].elempartpts = reshape([ne],1,1);
    dmd[i].nbsd = reshape([0],1,1);
    #dmd[i].elem2cpu = reshape([0],1,1);
    dmd[i].elemrecv = reshape([0],1,1);
    dmd[i].elemsend = reshape([0],1,1);
    dmd[i].elemrecvpts = reshape([0],1,1);
    dmd[i].elemsendpts = reshape([0],1,1);
    return dmd;
end

elem2cpu,~ = partition(t,ne,nproc,metis);

for i = 1:nproc
    intelem = findall(elem2cpu[:] .== (i-1)); # elements in subdomain i
    elem = neighboringelements(t2t, intelem); # all elements connected to elements in subdomain i
    extelem = sort(setdiff(elem,intelem)); # exterior elements

    elem = neighboringelements(t2t, extelem); # all elements connected to exterior elements
    bndelem = intersect(elem,intelem);  # boundary elements in subdmain i
    outelem = sort(setdiff(elem,vcat(intelem, extelem))); #  elements outside subdmain i

    tmp = [setdiff(intelem,bndelem); bndelem; extelem; outelem];
    dmd[i].elempart = reshape(tmp,1,length(tmp)); # partitions of elements
    dmd[i].elempartpts = [length(intelem)-length(bndelem) length(bndelem) length(extelem) length(outelem)];
    #dmd[i].elem2cpu = elem2cpu[dmd[i].elempart];
    nelem = dmd[i].elempartpts;

    recvelem = [extelem[:]; outelem[:]]; # elements received from neighbors
    tmp = hcat(elem2cpu[recvelem].+1, Vector((sum(nelem[1:2])+1):sum(nelem)));
    dmd[i].elemrecv = hcat(tmp, recvelem);
    ind = sortperm(dmd[i].elemrecv[:,1]);
    dmd[i].elemrecv = dmd[i].elemrecv[ind,:];
    dmd[i].nbsd = (sort(unique(dmd[i].elemrecv[:,1])))'; # neighboring subdomains
end

# store elements sent to neighboring subdomains to assemble the linear system
for k = 1:nproc
    dmd[k].elemsend = zeros(Int,0,1);
end
for i = 1:nproc
    for j = 1:length(dmd[i].nbsd)
        # cpu k sends information to cpu i
        k = dmd[i].nbsd[j];
        ii = findall(dmd[i].elemrecv[:,1] .== k);
        tm = dmd[i].elemrecv[ii,:];
        tm[:,1] .= i;
        tm[:,2] = indexin(tm[:,3],dmd[k].elempart[:]);
        if length(dmd[k].elemsend)==0
            dmd[k].elemsend = tm;
        else
            dmd[k].elemsend = [dmd[k].elemsend; tm];
        end
    end
end

for i = 1:nproc
    dmd[i].elemsendpts = zeros(Int,1, length(dmd[i].nbsd));
    dmd[i].elemrecvpts = zeros(Int,1, length(dmd[i].nbsd));
    for j = 1:length(dmd[i].nbsd)
        dmd[i].elemsendpts[j] = length(findall(dmd[i].elemsend[:,1] .== dmd[i].nbsd[j]));
        dmd[i].elemrecvpts[j] = length(findall(dmd[i].elemrecv[:,1] .== dmd[i].nbsd[j]));
    end
    dmd[i].elemsend = reshape(dmd[i].elemsend[:,2],1,size(dmd[i].elemsend,1));
    dmd[i].elemrecv = reshape(dmd[i].elemrecv[:,2],1,size(dmd[i].elemrecv,1));
end

return dmd; # elem2cpu; #elempart,elempartall

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module WandbMacros

using PyCall

export wandb, @wandbinit, @wandbconfig, @wandbsave, @wandblog, @wandbfinish

const wandb = PyNULL()

function __init__()
    copy!(wandb, pyimport_conda("wandb", "wandb", "conda-forge"))
end

include("wandbinit.jl")
include("wandbconfig.jl")
include("wandbsave.jl")
include("wandblog.jl")
include("wandbfinish.jl")

end


text:
export interpret

function interpret(m::Model{A,B,M}, tilde, ctx0) where {A,B,M}
    theModule = getmodule(m)
    mk_function(theModule, make_body(theModule, m.body, tilde, ctx0))
end

function make_body(M, f, m::AbstractModel)
    make_body(M, body(m))
end

function make_body(M, f, ast::Expr, return_action)
    function go(ex, scope=(bounds = Var[], freevars = Var[], bound_inits = Symbol[]))
        @match ex begin
            :(($x, $l) ~ $rhs) => begin
                varnames = Tuple(locals(l)) # ∪ locals(rhs))
                varvals = Expr(:tuple, varnames...)

                x = unsolve(x)
                l = unsolve(l)
                # q = quote
                #     _vars = NamedTuple{$varnames}($varvals)
                #     # @show _vars
                # end

                # unsolved_lhs = unsolve(lhs)
                # x == unsolved_lhs && delete!(varnames, x)

                sx = static(x)
                # X = to_type(unsolved_lhs)
                M = to_type(unsolve(rhs))

                q = quote
                    __old_x = $l == identity ? nothing : $x
                    ($x, _ctx, _retn) = $tilde($f, $l, $sx, __old_x, $rhs, _cfg, _ctx)
                    _retn isa Tilde.ReturnNow && return _retn.value
                end

                q
            end

            Expr(:scoped, new_scope, ex) => begin
                go(ex, new_scope)
            end

            Expr(head, args...) => Expr(head, map(Base.Fix2(go, scope), args)...)

            x => x
        end
    end

    if return_action isa DropReturn
        ast = drop_return(ast)
    end

    body = go(@q begin
            $(solve_scope(opticize(ast)))
    end) |> unsolve |> MacroTools.flatten

    body
end

function _get_gg_func_body(::RuntimeFn{Args,Kwargs,Body}) where {Args,Kwargs,Body}
    Body
end

function _get_gg_func_body(ex)
    error(ex)
end

struct DropReturn end
struct KeepReturn end

@generated function gg_call(_mc::MC, ::F, _cfg, _ctx, R) where {MC, F}
    _m = type2model(MC)
    M = getmodule(_m)

    _args = argvalstype(MC)
    _obs = obstype(MC)

    body = _m.body |> loadvals(_args, _obs)

    f = MeasureBase.instance(F)
    return_action = MeasureBase.instance(R)
    body = make_body(M, f, body, return_action)

    q = MacroTools.flatten(@q function (_mc, _cfg, _ctx)
            local _retn
            _args = Tilde.argvals(_mc)
            _obs = Tilde.observations(_mc)
            _cfg = merge(_cfg, (args=_args, obs=_obs))
            $body
            _retn
        end)

    from_type(_get_gg_func_body(mk_function(M, q)))
end


text:
export
    Constant

"""
    mutable struct Constant{O} <: Dist{O,Nothing}

`Constant` represents an sfunc that always produces the same value.  It has no
inputs and no parameters.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `logcpdf`
- `bounded_probs`
- `compute_pi`

# Type parameters
- `O`: the output type(s) of the `Constant`
"""
mutable struct Constant{O} <: Dist{O}
    "the constant value to be returned"
    x :: O
end

@impl begin
    struct ConstantSupport end

    function support(sf::Constant{O},
            parranges::NTuple{N,Vector},
            size::Integer,
            curr::Vector{<:O}) where {O,N}
        [sf.x]
    end
end

@impl begin
    struct ConstantSupportQuality end
    function support_quality(::Constant, parranges)
        :CompleteSupport
    end
end

@impl begin
    struct ConstantSample end

    function sample(sf::Constant{O}, i::Tuple{})::O where {O}
        sf.x
    end
end

@impl begin
    struct ConstantLogcpdf end

    function logcpdf(sf::Constant{O}, i::Tuple{}, o::O)::AbstractFloat where {O}
        return o == sf.x ? 0.0 : -Inf
    end
end

@impl begin
    struct ConstantBoundedProbs end

    function bounded_probs(sf::Constant{O},
            range::Vector{<:O},
            parranges::NTuple{N,Vector})::Tuple{Vector{<:AbstractFloat}, Vector{<:AbstractFloat}} where {O,N}

        p = sf.x in range ? [1.0] : [0.0]
        return (p, p)
    end
end

#=
@impl begin
    struct ConstantInitialStats end

    initial_stats(sf::Constant) = nothing
end
=#
@impl begin
    struct ConstantComputePi end

    function compute_pi(sf::Constant{O},
            range::Vector{<:O},
            parranges::NTuple{N,Vector},
            incoming_pis::Tuple)::Dist{<:O} where {O,N}

        ps = [x == sf.x ? 1.0 : 0.0 for x in range]
        Cat(range, ps)
    end
end
#=
@impl begin
    struct ConstantAccumulateStats end

    function accumulate_stats(::Constant, existing_stats, new_stats)
        nothing
    end
end

@impl begin
    struct ConstantExpectedStats end

    function expected_stats(::Constant{O},
            range::Vector{},
            parranges::NTuple{N,Vector},
            pis::NTuple{M,Dist},
            lambda::Score{<:O}) where {O,N,M}

        nothing
    end
end

@impl begin
    struct ConstantMaximizeStats end

    function maximize_stats(::Constant, stats)
        nothing
    end
end
=#


text:
gausspeaks(x, p) = sum([@. A * 1/√(2π * σ^2) * exp(-(x - μ)^2 / (2σ^2)) for (A, μ, σ) in p])

"""
    generate_testdata(grid, peaks, nm::AbstractNoiseModel; seedvalue=1) :: Curve

Return simulated test data (sum of Gauss peaks plus noise).

# Arguments

- `grid`: grid on which simulation shall be calculated
- `peaks`: list of tuples, each holding area, position, and width (standard deviation) of Gauss peaks
- `nm :: AbstractNoiseModel`: model that controls noise in simulation
"""
function generate_testdata(grid, peaks, nm::AbstractNoiseModel; seedvalue=1, baseline=nothing)
    seed!(seedvalue)
    # simulate FTIR spectrum with two bands
    # the true intensity ratio is 1 : 2
    signals = gausspeaks(grid, peaks)
    baseline = baseline === nothing ? zeros(eltype(grid), length(grid)) : baseline
    uc = add_noise(Curve(grid, baseline + signals), nm, 1)
    return get_draw(1, uc)
end

"""
Return test dataset with two peaks with area ratio 1 : 2 and correlated noise
(α = 0.05, λ = 0.5).
"""
testdata_1() = let
    x = collect(0:0.1:100)
    baseline = @. 1.0 + x*1.5e-2 - (x-50)^3*3e-6
    signals = [(1, 15, √0.5), (2, 30, √0.5)]
    generate_testdata(x, signals, MvGaussianNoiseModel(0.1, 0.05, 0.5); baseline=baseline)
end

"""
Return test dataset with two peaks with area ratio 1 : 2 and uncorrelated noise (σ = 0.1).
"""
testdata_2() = let
    x = collect(0:0.1:100)
    baseline = @. 1.0 + 1.5e-4x^2 - 3e-6x^3
    signals = [(1, 15, √0.5), (2, 30, √0.5)]
    generate_testdata(x, signals, GaussianNoiseModel(0.1); baseline=baseline)
end

"""
Return test dataset from testdata_1() with uneven grid spacing.
"""
testdata_3() = let
    x = collect(0:0.1:100)
    x = @. x + 1e-3x^2
    baseline = @. 1.0 + x*1.5e-2 - (x-50)^3*3e-6
    signals = [(1, 15, √0.5), (2, 30, √0.5)]
    generate_testdata(x, signals, MvGaussianNoiseModel(0.1, 0.05, 0.5); baseline=baseline)
end

"""
Return test dataset with four peaks with area ratio 1 : 2 : 3 : 0.5 and correlated noise
(α = 0.05, λ = 0.5).
"""
testdata_4() = let
    x = collect(0:0.1:200)
    baseline = @. 1.0 + x*2.5e-3 - (x-50)^3*8e-8
    signals = [(1, 15, √0.5), (2, 30, √0.5), (3, 60, √0.5), (1, 85, √0.5)]
    generate_testdata(x, signals, MvGaussianNoiseModel(0.1, 0.05, 0.5); baseline=baseline, seedvalue=42)
end


text:
using LinearAlgebra: norm, normalize

set(obj, ::typeof(last), val) = @set obj[lastindex(obj)] = val
set(obj, ::typeof(first), val) = @set obj[firstindex(obj)] = val
delete(obj, ::typeof(last)) = delete(obj, IndexLens((lastindex(obj),)))
delete(obj, ::typeof(first)) = delete(obj, IndexLens((firstindex(obj),)))
insert(obj, ::typeof(last), val) = insert(obj, IndexLens((lastindex(obj) + 1,)), val)
insert(obj, ::typeof(first), val) = insert(obj, IndexLens((firstindex(obj),)), val)

set(obj, ::typeof(identity), val) = val
set(obj, ::typeof(inv), new_inv) = inv(new_inv)

################################################################################
##### eltype
################################################################################
set(obj::Array, ::typeof(eltype), T::Type) = collect(T, obj)
set(obj::Number, ::typeof(eltype), T::Type) = T(obj)
set(obj::Type{<:Number}, ::typeof(eltype), ::Type{T}) where {T} = T
set(obj::Type{<:Array{<:Any,N}}, ::typeof(eltype), ::Type{T}) where {N,T} = Array{T,N}
set(obj::Type{<:Dict}, ::typeof(eltype), ::Type{Pair{K,V}}) where {K,V} = Dict{K,V}
set(obj::Dict, ::typeof(eltype), ::Type{T}) where {T} = set(typeof(obj), eltype, T)(obj)

set(obj::Dict, lens::Union{typeof(keytype),typeof(valtype)}, T::Type) =
    set(typeof(obj), lens, T)(obj)
set(obj::Type{<:Dict{<:Any,V}}, lens::typeof(keytype), ::Type{K}) where {K,V} = Dict{K,V}
set(obj::Type{<:Dict{K}}, lens::typeof(valtype), ::Type{V}) where {K,V} = Dict{K,V}

################################################################################
##### os
################################################################################
set(path, ::typeof(splitext), (stem, ext))     = string(stem, ext)
set(path, ::typeof(splitdir), (dir, last))     = joinpath(dir, last)
set(path, ::typeof(splitdrive), (drive, rest)) = joinpath(drive, rest)
set(path, ::typeof(splitpath), pieces)         = joinpath(pieces...)
set(path, ::typeof(dirname), new_name)         = @set splitdir(path)[1] = new_name
set(path, ::typeof(basename), new_name)        = @set splitdir(path)[2] = new_name
delete(path, ::typeof(basename)) = dirname(path)
delete(path, ::typeof(dirname)) = basename(path)

################################################################################
##### math
################################################################################
set(x::Real, ::typeof(real), y) = y
set(x,       ::typeof(real), y) = y + im*imag(x)
set(x,       ::typeof(imag), y) = real(x) + im*y

set(arr, ::typeof(normalize), val) = norm(arr) * val
set(arr, ::typeof(norm), val)      = val/norm(arr) * arr # should we check val is positive?


text:
## modules.j. Modules used in ROCAnalysis
## (c) 2015 David A. van Leeuwen.
##
## Licensed under the MIT software license, see LICENSE.md

## Loading modules has been moved to the source files that actually use these modules.

## These are just some functions from NumericalFuns that we somehow had to implement ourselves.

logit(x::AbstractFloat) = log(x / (one(x) - x))
logit(x::Real) = logit(float(x))

softplus(x::AbstractFloat) = x ≤ 0 ? log1p(exp(x)) : x + log1p(exp(-x))
softplus(x::Real) = softplus(float(x))

function sigmoid(x::AbstractFloat)
   if x ≤ 0
	   y = exp(x)
	   return y / (y + one(x))
   else
	   return one(x) / (one(x) + exp(-x))
   end
end
sigmoid(x::Real) = sigmoid(float(x))


text:
#=##############################################################################
# DESCRIPTION
    Class of wind turbine geometry.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jul 2018
  * License   : MIT License
=###############################################################################

################################################################################
# WIND FARM
################################################################################

# import JLD
# import CSV
# import Dierckx
# import PyPlot
# using PyCall
#
# const plt = PyPlot

# ------------ FLOW MODULES ----------------------------------------------------
# https://github.com/byuflowlab/GeometricTools.jl
# import GeometricTools
# const gt = GeometricTools

# include("WFVisual_turbine.jl")

using PyCall

"""
`generate_layout(D::Array{T,1}, H::Array{T,1}, N::Array{Int64,1},
                          x::Array{T,1}, y::Array{T,1}, z::Array{T,1},
                          glob_yaw::Array{T,1};
                          # TURBINE GEOMETRY OPTIONS
                          hub::Array{String,1}=String[],
                          tower::Array{String,1}=String[],
                          blade::Array{String,1}=String[],
                          data_path::String=def_data_path,
                          # FILE OPTIONS
                          save_path=nothing, file_name="windfarm",
                          paraview=true
                         ) where{T<:Real}`

  Generates and returns a wind farm layout consisting of a MultiGrid object
containing every turbine at the indicated position `x,y,z`, and the given
geometry `D,H,N` and orientation `glob_yaw`.

  **Arguments**
  * `D::Array{Float64,1}`         : Rotor diameter of every turbine.
  * `H::Array{Float64,1}`         : Tower height of every turbine.
  * `N::Array{Float64,1}`         : Number of blades of every turbine.
  * `x::Array{Float64,1}`         : x-position of the every turbine base.
  * `y::Array{Float64,1}`         : y-position of the every turbine base.
  * `z::Array{Float64,1}`         : z-position of the every turbine base.
  * `glob_yaw::Array{Float64,1}`  : Angle of the plane of rotation relative to
                                    the x-axis of the global coordinate system
                                    IN DEGREES.

"""
function generate_layout(D::Array{T,1}, H::Array{T,1}, N::Array{Int64,1},
                          x::Array{T,1}, y::Array{T,1}, z::Array{T,1},
                          glob_yaw::Array{T,1};
                          # TURBINE GEOMETRY OPTIONS
                          hub::Array{String,1}=String[],
                          tower::Array{String,1}=String[],
                          blade::Array{String,1}=String[],
                          data_path::String=def_data_path,
                          # FILE OPTIONS
                          save_path=nothing, file_name="windfarm",
                          paraview=true
                         ) where{T<:Real}

  nturbines = size(D, 1)        # Number of turbines

  # Default turbine geometry
  if size(hub,1)==0
    hub = ["hub" for i in 1:nturbines]
  end
  if size(tower,1)==0
    tower = ["tower1" for i in 1:nturbines]
  end
  if size(blade,1)==0
    blade = ["NREL5MW" for i in 1:nturbines]
  end

  # Generates layout
  windfarm = gt.MultiGrid(3)

  for i in 1:nturbines
    # Generate wind turbine geometry
    turbine = generate_windturbine(D[i]/2, H[i], blade[i], hub[i], tower[i];
                                    nblades=N[i], data_path=data_path,
                                    save_path=nothing)

    # Places it at the location and orientation
    Oaxis = gt.rotation_matrix(glob_yaw[i], 0, 0)
    gt.lintransform!(turbine, Oaxis, [x[i], y[i], z[i]])

    # Adds it to the farm
    gt.addgrid(windfarm, "turbine$i", turbine)
  end

  if save_path!=nothing
    gt.save(windfarm, file_name; path=save_path)

    if paraview
      strn = ""
      for i in 1:nturbines
        strn *= file_name*"_turbine$(i)_tower.vtk;"
        strn *= file_name*"_turbine$(i)_rotor_hub.vtk;"
        for j in 1:N[i]
          strn *= file_name*"_turbine$(i)_rotor_blade$(j).vtk;"
        end
      end

      run(`paraview --data=$save_path/$strn`)
    end

  end

  return windfarm::gt.MultiGrid
end

"""
`generate_perimetergrid(perimeter::Array{Array{T, 1}, 1},
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min::Real=0, z_max::Real=0,
                                  # SPLINE OPTIONS
                                  verify_spline::Bool=true,
                                  spl_s=0.001, spl_k="automatic",
                                  # FILE OPTIONS
                                  save_path=nothing, file_name="perimeter",
                                  paraview=true
                                )`

  Generates the perimeter grid with `perimeter` the array of points of the
  contour (must be a closed contour), and `NDIVS_` the number of cells in
  each parametric dimension (give it `NDIVSz=0` for a flat surface, otherwise
  it'll generate a volumetric grid between `z_min` and `z_max`).
"""
function generate_perimetergrid(perimeter::Array{Array{T, 1}, 1},
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min::Real=0, z_max::Real=0,
                                  # SPLINE OPTIONS
                                  verify_spline::Bool=true,
                                  spl_s=0.001, spl_k="automatic", atol=1.0,
                                  # FILE OPTIONS
                                  save_path=nothing, file_name="perimeter",
                                  paraview=true
                                ) where{T<:Real}
  # Error cases
  multidiscrtype = Array{Tuple{Float64,Int64,Float64,Bool},1}
  if typeof(NDIVSx)==Int64
    nz = NDIVSz
  elseif typeof(NDIVSz)==multidiscrtype
    nz = 0
    for sec in NDIVSz
      nz += sec[2]
    end
  else
    error("Expected `NDIVSz` to be type $(Int64) or $MultiDiscrType,"*
            " got $(typeof(NDIVSz)).")
  end

  # --------- REPARAMETERIZES THE PERIMETER ---------------------------
  org_x = [p[1] for p in perimeter]
  org_y = [p[2] for p in perimeter]
  # Separate upper and lower sides to make the contour injective in x
  upper, lower = gt.splitcontour(org_x, org_y)

  # # Parameterize both sides independently
  # fun_upper = gt.parameterize(upper[1], upper[2], zeros(upper[1]); inj_var=1,
  #                                                     s=spl_s, kspl=spl_k)
  # fun_lower = gt.parameterize(lower[1], lower[2], zeros(lower[1]); inj_var=1,
  #                                                     s=spl_s, kspl=spl_k)
  # # Discretizes both sides
  # if NDIVSx==multidiscrtype
  #   new_upper = gt.multidiscretize(fun_upper, 0, 1, NDIVSx)
  #   new_lower = gt.multidiscretize(fun_lower, 0, 1, NDIVSx)
  # else
  #   new_upper = gt.discretize(fun_upper, 0, 1, NDIVSx, 1.0)
  #   new_lower = gt.discretize(fun_lower, 0, 1, NDIVSx, 1.0)
  # end

  splt_up = Int(ceil((size(upper[1],1)/2)))
  splt_low = Int(ceil((size(lower[1],1)/2)))

  # Splits the perimeter into four faces
  upper1 = [[x for x in upper[1][1:splt_up]], [y for y in upper[2][1:splt_up]]]
  upper2 = [[x for x in upper[1][splt_up:end]], [y for y in upper[2][splt_up:end]]]
  lower1 = [[x for x in lower[1][1:splt_low]], [y for y in lower[2][1:splt_low]]]
  lower2 = [[x for x in lower[1][splt_low:end]], [y for y in lower[2][splt_low:end]]]

  # Parameterize both sides independently
  fun_upper1 = gt.parameterize(upper1[1], upper1[2], zeros(upper1[1]); inj_var=1,
                                              s=spl_s, kspl=spl_k, atol=atol)
  fun_upper2 = gt.parameterize(upper2[1], upper2[2], zeros(upper2[1]); inj_var=1,
                                              s=spl_s, kspl=spl_k, atol=atol)
  fun_lower1 = gt.parameterize(lower1[1], lower1[2], zeros(lower1[1]); inj_var=1,
                                              s=spl_s, kspl=spl_k, atol=atol)
  fun_lower2 = gt.parameterize(lower2[1], lower2[2], zeros(lower2[1]); inj_var=1,
                                              s=spl_s, kspl=spl_k, atol=atol)
  # Discretizes both sides
  if NDIVSx==multidiscrtype
    new_upper1 = gt.multidiscretize(fun_upper1, 0, 1, NDIVSx)
    new_upper2 = gt.multidiscretize(fun_upper2, 0, 1, NDIVSy)
    new_lower1 = gt.multidiscretize(fun_lower1, 0, 1, NDIVSy)
    new_lower2 = gt.multidiscretize(fun_lower2, 0, 1, NDIVSx)
  else
    new_upper1 = gt.discretize(fun_upper1, 0, 1, NDIVSx, 1.0)
    new_upper2 = gt.discretize(fun_upper2, 0, 1, NDIVSy, 1.0)
    new_lower1 = gt.discretize(fun_lower1, 0, 1, NDIVSy, 1.0)
    new_lower2 = gt.discretize(fun_lower2, 0, 1, NDIVSx, 1.0)
  end

  # ----------------- SPLINE VERIFICATION --------------------------------------

  if verify_spline
  #   new_points = vcat(reverse(new_upper), new_lower)
    # new_x = [p[1] for p in new_points]
    # new_y = [p[2] for p in new_points]
    plt.plot(org_x, org_y, "--^k", label="Original", alpha=0.5)
    # plt.plot(new_x, new_y, ":.b", label="Parameterized")

    plt.plot(upper1[1], upper1[2], "*b", label="upper1", alpha=0.75)
    plt.plot(upper2[1], upper2[2], "^b", label="upper2", alpha=0.75)
    plt.plot(lower1[1], lower1[2], "*r", label="lower1", alpha=0.75)
    plt.plot(lower2[1], lower2[2], "^r", label="lower2", alpha=0.75)

    new_points = vcat(new_upper1, new_upper2)
    new_x = [p[1] for p in new_points]
    new_y = [p[2] for p in new_points]
    plt.plot(new_x, new_y, ":.b", label="Parameterized Upper", alpha=0.75)
    new_points = vcat(new_lower1, new_lower2)
    new_x = [p[1] for p in new_points]
    new_y = [p[2] for p in new_points]
    plt.plot(new_x, new_y, ":.r", label="Parameterized Lower", alpha=0.75)
    plt.xlabel(plt.L"x")
    plt.ylabel(plt.L"y")
    plt.legend(loc="best")
    plt.grid(true, color="0.8", linestyle="--")
  end

  # --------- GRIDS THE INSIDE OF THE PERIMETER ---------------------
  # Parametric grid
  P_min = zeros(3)
  P_max = [1, 1, 1*(nz!=0)]
  param_grid = gt.Grid(P_min, P_max, [NDIVSx, NDIVSy, NDIVSz])

  # function my_space_transform(X, ind)
  #     i = ind[1]                      # Arc length point
  #     w = X[2]                        # Weight
  #     z = z_min + X[3]*(z_max-z_min)  # z-position
  #
  #     Y = new_lower[i] + w*(new_upper[i]-new_lower[i])
  #     Y[3] = z
  #
  #     return Y
  # end

  rev_new_upper1 = reverse(new_upper1) # Left edge
  rev_new_upper2 = reverse(new_upper2) # Upper edge
  rev_new_lower1 = reverse(new_lower1) # Lower edge
  rev_new_lower2 = reverse(new_lower2) # Right edge

  function my_space_transform(X, ind)
      beta = 1e-2*pi/2                 # Warping parameter

      z = z_min + X[3]*(z_max-z_min)   # z-position

      xw = X[2]                        # x weight
      yw = X[1]                        # y weight

      # Linearly interpolate the edges
      x = new_upper1[ind[1]]*xw + new_lower2[ind[1]]*(1-xw) # Left-right edges
      y = rev_new_upper2[ind[2]]*yw + rev_new_lower1[ind[2]]*(1-yw) # Upper-lower edges
      # Y = x + y

      # A clever way of merging the two interpolations
      # First weight
      yw = abs(X[1]-0.5)/0.5
      yw = tan(yw*pi/2 - beta)/tan(pi/2 - beta)
      # Second weight
      xw = abs(X[2]-0.5)/0.5
      xw = tan(xw*pi/2 - beta)/tan(pi/2 - beta)
      # Merged weight
      w = ( (1-xw) + yw ) /2
      # Merged point
      Y = x*(1-w) + y*w

      Y[3] = z

      return Y
  end

  # Applies the space transformation to the parametric grid
  gt.transform3!(param_grid, my_space_transform)

  if save_path!=nothing
    gt.save(param_grid, file_name; path=save_path)

    if paraview
      strn = file_name*".vtk"
      run(`paraview --data=$save_path/$strn`)
    end

  end

  return param_grid::gt.Grid
end

function generate_windfarm(D::Array{T,1}, H::Array{T,1}, N::Array{Int64,1},
                          x::Array{T,1}, y::Array{T,1}, z::Array{T,1},
                          glob_yaw::Array{T,1}, perimeter::Array{T, 2},
                          wake; optargs...
                         ) where{T<:Real}

    _perimeter = M2arr(perimeter)
    return generate_windfarm(D, H, N, x, y, z, glob_yaw, _perimeter, wake;
                                                                    optargs...)
end

function generate_windfarm(D::Array{T,1}, H::Array{T,1}, N::Array{Int64,1},
                          x::Array{T,1}, y::Array{T,1}, z::Array{T,1},
                          glob_yaw::Array{T,1}, perimeter::Array{Array{T, 1}},
                          wake;
                          # TURBINE GEOMETRY OPTIONS
                          hub::Array{String,1}=String[],
                          tower::Array{String,1}=String[],
                          blade::Array{String,1}=String[],
                          data_path::String=def_data_path,
                          # PERIMETER AND FLUID DOMAIN OPTIONS
                          NDIVSx=50, NDIVSy=50, NDIVSz=50,
                          z_min="automatic", z_max="automatic",
                          z_off=0.0,
                          fdom_perimeter::Union{Array{Array{T, 1}}, Void}=nothing,
                          # PERIMETER SPLINE OPTIONS
                          verify_spline::Bool=true,
                          spl_s=0.001, spl_k="automatic",
                          # FILE OPTIONS
                          save_path=nothing, file_name="mywindfarm",
                          paraview=true, num=nothing
                         ) where{T<:Real}

  windfarm = generate_layout(D, H, N, x, y, z, glob_yaw;
                                  hub=hub, tower=tower, blade=blade,
                                  data_path=data_path, save_path=nothing)

  perimeter_grid = generate_perimetergrid(perimeter, NDIVSx, NDIVSy, 0;
                                     z_min=z_off, z_max=z_off,
                                      verify_spline=verify_spline, spl_s=spl_s,
                                      spl_k=spl_k, save_path=nothing)

  _zmin = z_min=="automatic" ? 0 : z_min
  _zmax = z_max=="automatic" ? maximum(H) + 1.25*maximum(D)/2 : z_max
  fdom = generate_perimetergrid(fdom_perimeter != nothing ? fdom_perimeter : perimeter,
                                    NDIVSx, NDIVSy, NDIVSz;
                                    z_min=_zmin+z_off, z_max=_zmax+z_off,
                                    verify_spline=false,
                                    spl_s=spl_s, spl_k=spl_k,
                                    save_path=nothing,
                                  )

  gt.calculate_field(fdom, wake, "wake", "vector", "node")

  if save_path!=nothing
    gt.save(windfarm, file_name; path=save_path, num=num)
    gt.save(perimeter_grid, file_name*"_perimeter"; path=save_path, num=num)
    gt.save(fdom, file_name*"_fdom"; path=save_path, num=num)

    if paraview
      strn = ""
      for i in 1:size(D,1)
        strn *= file_name*"_turbine$(i)_tower.vtk;"
        strn *= file_name*"_turbine$(i)_rotor_hub.vtk;"
        for j in 1:N[i]
          strn *= file_name*"_turbine$(i)_rotor_blade$(j).vtk;"
        end
      end

      strn *= file_name*"_perimeter.vtk;"
      strn *= file_name*"_fdom.vtk;"

      run(`paraview --data=$save_path/$strn`)
    end

  end

  return (windfarm, perimeter, fdom)
end

function M2arr(M::Array{T, 2}) where{T<:Real}
  return [ [p for p in M[i, :]] for i in 1:size(M,1) ]
end

# ------------ END OF WIND FARM ------------------------------------------------

max_issues_repo_name: byuflowlab/wind-farm-utilities

max_forks_repo_name: byuflowlab/wind-farm-utilities


text:
#!/usr/bin/env julia --threads 6
using XAM,Statistics,GenomicFeatures

bName = ARGS[1]

println("Threads: ",Threads.nthreads())

function histQuals(bamName,chrName)
	allquals=Int[]
	bamIn = open(BAM.Reader,bamName,index="$bamName.bai")
	count = 0

	for a in eachoverlap(bamIn,chrName,1:40000000)
		count+=1
    	push!(allquals,BAM.mappingquality(a))
	end

	close(bamIn)

	println(chrName," Count: ",count)
	println("Mean quality: ",mean(allquals))
end

chromList=["chr1","chr2","chr3","chr4","chr5","chr6","chr7","chr8","chr9","chr10","chr11",
			"chr12","chr13","chr14","chr15","chr16","chr17","chr18","chr19","chr20","chr21","chrX"]

Threads.@threads for cName in chromList
	histQuals(bName,cName)
end

# ./averagequal_threads.jl PC-HLA-EG4.bam  16.93s user 0.81s system 342% cpu 5.185 total

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
using Test
using QuantumLattices.Prerequisites.TypeTraits

struct WithTrait{F1,F2,F3}
    f1::F1
    f2::F2
    f3::F3
end
Base.:(==)(fc1::WithTrait,fc2::WithTrait) = ==(efficientoperations,fc1,fc2)
Base.isequal(fc1::WithTrait,fc2::WithTrait) = isequal(efficientoperations,fc1,fc2)
Base.:<(fc1::WithTrait,fc2::WithTrait) = <(efficientoperations,fc1,fc2)
Base.isless(fc1::WithTrait,fc2::WithTrait) = isless(efficientoperations,fc1,fc2)
Base.isapprox(fc1::WithTrait,fc2::WithTrait;atol::Real=10^-5,rtol::Real=10^-5) = isapprox(efficientoperations,Val((:f1,:f2)),fc1,fc2;atol=atol,rtol=rtol)
Base.replace(fc::WithTrait;kwargs...) = replace(efficientoperations,fc;kwargs...)

@testset "efficientoperations" begin
    @test ==(efficientoperations,(),())
    @test isequal(efficientoperations,(),())
    @test ==(efficientoperations,(1,2),(1,2,3))==false
    @test isequal(efficientoperations,(1,2),(1,2,3))==false

    fc1,fc2=WithTrait(1.0,2,3),WithTrait(1,2,3)
    @test fc1==fc2
    @test isequal(fc1,fc2)
    fc1,fc2=WithTrait(1.0,2,3),WithTrait(1,2,4.0)
    @test fc1<fc2
    @test isless(fc1,fc2)

    fc1,fc2=WithTrait(1.0+10^-6,2,3),WithTrait(1,2,3)
    @test fc1≈fc2
    fc1,fc2=WithTrait(1.0,2-10^-6,3),WithTrait(1,2,3)
    @test fc1≈fc2
    fc1,fc2=WithTrait(1.0,2,3-10^-6),WithTrait(1,2,3)
    @test fc1≉ fc2

    @test replace(WithTrait(1,2,3),f1='c')==WithTrait('c',2,3)
end

@testset "forder/corder" begin
    dims=(2,2,2)
    finds=[(1,1,1),(2,1,1),(1,2,1),(2,2,1),(1,1,2),(2,1,2),(1,2,2),(2,2,2)]
    cinds=[(1,1,1),(1,1,2),(1,2,1),(1,2,2),(2,1,1),(2,1,2),(2,2,1),(2,2,2)]
    @test collect(indtosub(dims,i,forder) for i=1:prod(dims))==finds
    @test collect(indtosub(dims,i,corder) for i=1:prod(dims))==cinds
    @test collect(subtoind(dims,inds,forder) for inds in finds)==collect(1:prod(dims))
    @test collect(subtoind(dims,inds,corder) for inds in cinds)==collect(1:prod(dims))
end


text:
"""
    side_indices(c::GridContent)::RowCols{Int}

Indices of the rows / cols for each side
"""
function side_indices(c::GridContent)
    return RowCols(
        c.span.cols.start,
        c.span.cols.stop,
        c.span.rows.start,
        c.span.rows.stop,
    )
end

# These functions tell whether an object in a grid touches the left, top, etc. border
# of the grid. This means that it is relevant for the grid's own protrusion on that side.
ismostin(gc::GridContent, grid, ::Left) = gc.span.cols.start == 1
ismostin(gc::GridContent, grid, ::Right) = gc.span.cols.stop == grid.ncols
ismostin(gc::GridContent, grid, ::Bottom) = gc.span.rows.stop == grid.nrows
ismostin(gc::GridContent, grid, ::Top) = gc.span.rows.start == 1

function protrusion(x::T, side::Side) where T
    protrusions = protrusionsobservable(x)
    @match side begin
        si::Left => protrusions[].left
        si::Right => protrusions[].right
        si::Bottom => protrusions[].bottom
        si::Top => protrusions[].top
        si => error("Can't get a protrusion value for side $(typeof(si)), only
            Left, Right, Bottom, or Top.")
    end
end

function protrusion(gc::GridContent, side::Side)
    prot = @match gc.side begin
        gcside::Inner => protrusion(gc.content, side)
        # gcside::Outer => BBox(l - pl, r + pr, b - pb, t + pt)
        gcside::Union{Left, Right} => @match side begin
            si::typeof(gcside) => determinedirsize(gc.content, Col(), gc.side)
            si => 0.0
        end
        gcside::Union{Top, Bottom} => @match side begin
            si::typeof(gcside) => determinedirsize(gc.content, Row(), gc.side)
            si => 0.0
        end
        gcside::TopLeft => @match side begin
            si::Top => determinedirsize(gc.content, Row(), gc.side)
            si::Left => determinedirsize(gc.content, Col(), gc.side)
            si => 0.0
        end
        gcside::TopRight => @match side begin
            si::Top => determinedirsize(gc.content, Row(), gc.side)
            si::Right => determinedirsize(gc.content, Col(), gc.side)
            si => 0.0
        end
        gcside::BottomLeft => @match side begin
            si::Bottom => determinedirsize(gc.content, Row(), gc.side)
            si::Left => determinedirsize(gc.content, Col(), gc.side)
            si => 0.0
        end
        gcside::BottomRight => @match side begin
            si::Bottom => determinedirsize(gc.content, Row(), gc.side)
            si::Right => determinedirsize(gc.content, Col(), gc.side)
            si => 0.0
        end
        gcside => error("Invalid side $gcside")
    end
    ifnothing(prot, 0.0)
end

getside(m::Mixed, ::Left) = m.sides.left
getside(m::Mixed, ::Right) = m.sides.right
getside(m::Mixed, ::Top) = m.sides.top
getside(m::Mixed, ::Bottom) = m.sides.bottom

function inside_protrusion(gl::GridLayout, side::Side)
    prot = 0.0
    for elem in gl.content
        if ismostin(elem, gl, side)
            # take the max protrusion of all elements that are sticking
            # out at this side
            prot = max(protrusion(elem, side), prot)
        end
    end
    return prot
end

function protrusion(gl::GridLayout, side::Side)
    # when we align with the outside there is by definition no protrusion
    if gl.alignmode[] isa Outside
        return 0.0
    elseif gl.alignmode[] isa Inside
        inside_protrusion(gl, side)
    elseif gl.alignmode[] isa Mixed
        si = getside(gl.alignmode[], side)
        if isnothing(si)
            inside_protrusion(gl, side)
        elseif si isa Protrusion
            si.p
        else
            # Outside alignment
            0.0
        end
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode[]))")
    end
end

function bbox_for_solving_from_side(maxgrid::RowCols, bbox_cell::Rect2f, idx_rect::RowCols, side::Side)
    pl = maxgrid.lefts[idx_rect.lefts]
    pr = maxgrid.rights[idx_rect.rights]
    pt = maxgrid.tops[idx_rect.tops]
    pb = maxgrid.bottoms[idx_rect.bottoms]

    l = left(bbox_cell)
    r = right(bbox_cell)
    b = bottom(bbox_cell)
    t = top(bbox_cell)

    @match side begin
        s::Inner => bbox_cell
        s::Outer => BBox(l - pl, r + pr, b - pb, t + pt)
        s::Left => BBox(l - pl, l, b, t)
        s::Top => BBox(l, r, t, t + pt)
        s::Right => BBox(r, r + pr, b, t)
        s::Bottom => BBox(l, r, b - pb, b)
        s::TopLeft => BBox(l - pl, l, t, t + pt)
        s::TopRight => BBox(r, r + pr, t, t + pt)
        s::BottomRight => BBox(r, r + pr, b - pb, b)
        s::BottomLeft => BBox(l - pl, l, b - pb, b)
        s => error("Invalid side $s")
    end
end

startside(c::Col) = Left()
stopside(c::Col) = Right()
startside(r::Row) = Top()
stopside(r::Row) = Bottom()

getspan(gc::GridContent, dir::Col) = gc.span.cols
getspan(gc::GridContent, dir::Row) = gc.span.rows

"""
Determine the size of a protrusion layout along a dimension. This size is dependent
on the `Side` at which the layout is placed in its parent grid. An `Inside` side
means that the protrusion layout reports its width but not its protrusions. `Left`
means that the layout reports only its full width but not its height, because
an element placed in the left protrusion loses its ability to influence height.
"""
function determinedirsize(content, gdir::GridDir, side::Side)
    reportedsize = reportedsizeobservable(content)
    if gdir isa Row
        @match side begin
            # TODO: is reportedsize the correct thing to return? or plus protrusions depending on the side
            si::Union{Inner, Top, Bottom, TopLeft, TopRight, BottomLeft, BottomRight} =>
                ifnothing(reportedsize[][2], nothing)
            si::Union{Left, Right} => nothing
            si => error("$side not implemented")
        end
    else
        @match side begin
            si::Union{Inner, Left, Right, TopLeft, TopRight, BottomLeft, BottomRight} =>
                ifnothing(reportedsize[][1], nothing)
            si::Union{Top, Bottom} => nothing
            si => error("$side not implemented")
        end
    end
end

function to_ranges(g::GridLayout, rows::Indexables, cols::Indexables)
    if rows isa Int
        rows = rows:rows
    elseif rows isa Colon
        rows = 1:g.nrows
    end
    if cols isa Int
        cols = cols:cols
    elseif cols isa Colon
        cols = 1:g.ncols
    end
    rows, cols
end

function adjust_rows_cols!(g::GridLayout, rows, cols; update = true)
    rows, cols = to_ranges(g, rows, cols)

    if rows.start < 1
        n = 1 - rows.start
        prependrows!(g, n, update = update)
        # adjust rows for the newly prepended ones
        rows = rows .+ n
    end
    if rows.stop > g.nrows
        n = rows.stop - g.nrows
        appendrows!(g, n, update = update)
    end
    if cols.start < 1
        n = 1 - cols.start
        prependcols!(g, n, update = update)
        # adjust cols for the newly prepended ones
        cols = cols .+ n
    end
    if cols.stop > g.ncols
        n = cols.stop - g.ncols
        appendcols!(g, n, update = update)
    end

    rows, cols
end


text:
include(joinpath(pwd(),"src/models/rocket.jl"))
include(joinpath(pwd(),"src/constraints/free_time.jl"))

optimize = true

# Free-time model
model_nom = free_time_model(model_nominal)

function fd(model::RocketNominal, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, u[end])
end

# Horizon
T = 41

# Bounds
ul, uu = control_bounds(model_nom, T, [-5.0; 0.0; 0.01], [5.0; 100.0; 1.0])

# Initial and final states
x1 = [15.0;
	  model_nom.l1 + 10.0;
	  -45.0 * pi / 180.0;
	  -10.0;
	  -10.0;
	  -1.0 * pi / 180.0]

xT = [0.0;
	  model_nom.l1;
	  0.0;
	  0.0;
	  0.0;
	  0.0]

# Radius of landing pad
r_pad = 0.25

# xl <= x <= xl
_xl = -Inf * ones(model_nom.n)
_xl[2] = model_nom.l1
_xu = Inf * ones(model_nom.n)

xl, xu = state_bounds(model_nom, T, _xl, _xu, x1 = x1, xT = xT)

xl[T][1] = -1.0 * r_pad
xu[T][1] = r_pad
xl[T][2] = xT[2] - 0.001
xu[T][2] = xT[2] + 0.001
xl[T][3] = -1.0 * pi / 180.0
xu[T][3] = 1.0 * pi / 180.0
xl[T][4] = -0.001
xu[T][4] = 0.001
xl[T][5] = -0.001
xu[T][5] = 0.001
xl[T][6] = -0.01 * pi / 180.0
xu[T][6] = 0.01 * pi / 180.0

# Objective
obj = quadratic_time_tracking_objective(
	[(t != T ? Diagonal([1.0; 10.0; 1.0; 1.0; 10.0; 1.0])
		: Diagonal([10.0; 100.0; 10.0; 10.0; 100.0; 10.0])) for t = 1:T],
	[Diagonal([1.0e-1*ones(2); 0.0]) for t = 1:T-1],
    [xT for t = 1:T],
	[zeros(model_nom.m) for t = 1:T],
	1.0)

# Time step constraints
con_free_time = free_time_constraints(T)

# Problem
prob_nominal = trajectory_optimization_problem(model_nom,
					obj,
					T,
                    xl = xl,
                    xu = xu,
                    ul = ul,
                    uu = uu,
                    con = con_free_time
                    )

# Trajectory initialization
X0 = linear_interp(x1, xT, T) # linear interpolation on state
U0 = [ones(model_nom.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob_nominal)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
if optimize
    @time Z̄ = solve(prob_nominal, copy(Z0))
    @save joinpath(@__DIR__, "sol.jld2") Z̄
else
    println("Loading solution...")
    @load joinpath(@__DIR__, "sol.jld2") Z̄
end

X, U = unpack(Z̄, prob_nominal.prob)
@show sum([U[t][end] for t = 1:T-1])

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning