text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

""" Datatypes and post-processing for raw flow graphs.
"""
module RawFlowGraphs
export RawNode, RawPort, RawNodeAnnotationKind,
  FunctionAnnotation, ConstructAnnotation, SlotAnnotation,
  read_raw_graph, rem_literals!, rem_unused_ports

using AutoHashEquals, Parameters
using Nullables

using Catlab.WiringDiagrams

@enum(RawNodeAnnotationKind,
  FunctionAnnotation = 0,
  ConstructAnnotation = 1,
  SlotAnnotation = 2)

function Base.convert(::Type{RawNodeAnnotationKind}, s::String)
  if (s == "function") FunctionAnnotation
  elseif (s == "construct") ConstructAnnotation
  elseif (s == "slot") SlotAnnotation
  else error("Unknown annotation kind \"$s\"") end
end

@with_kw struct RawNode
  language::Dict{String,Any} = Dict{String,Any}()
  annotation::Nullable{String} = Nullable{String}()
  annotation_index::Nullable{Int} = Nullable()
  annotation_kind::RawNodeAnnotationKind = FunctionAnnotation
end

function Base.:(==)(n1::RawNode, n2::RawNode)
  n1.language == n2.language &&
  isequal(n1.annotation, n2.annotation) &&
  isequal(n1.annotation_index, n2.annotation_index) &&
  n1.annotation_kind == n2.annotation_kind
end

@with_kw struct RawPort
  language::Dict{String,Any} = Dict{String,Any}()
  annotation::Nullable{String} = Nullable{String}()
  annotation_index::Nullable{Int} = Nullable()
  value::Nullable = Nullable()
end

function Base.:(==)(p1::RawPort, p2::RawPort)
  p1.language == p2.language &&
  isequal(p1.annotation, p2.annotation) &&
  isequal(p1.annotation_index, p2.annotation_index) &&
  isequal(p1.value, p2.value)
end

# Graph post-processing
#######################

""" Remove literals from raw flow graph.

Removes all nodes that are literal value constructors. (Currently, such nodes
occur in raw flow graphs for R, but not Python.)
"""
function rem_literals!(d::WiringDiagram)
  literals = filter(box_ids(d)) do v
    kind = get(box(d,v).value.language, "kind", "function")
    kind == "literal"
  end
  rem_boxes!(d, literals)
  d
end

""" Remove input and output ports with no connecting wires.

This simplification is practically necessary to visualize raw flow graphs
because scientific computing functions often have dozens of keyword arguments
(which manifest as input ports).
"""
function rem_unused_ports(diagram::WiringDiagram)
  result = WiringDiagram(input_ports(diagram), output_ports(diagram))
  for v in box_ids(diagram)
    # Note: To ensure that port numbers on wires remain valid, we only remove
    # unused ports beyond the last used port.
    b = box(diagram, v)
    last_used_input = maximum([0; [wire.target.port for wire in in_wires(diagram, v)]])
    last_used_output = maximum([0; [wire.source.port for wire in out_wires(diagram, v)]])
    unused_inputs = input_ports(b)[1:last_used_input]
    unused_outputs = output_ports(b)[1:last_used_output]
    @assert add_box!(result, Box(b.value, unused_inputs, unused_outputs)) == v
  end
  add_wires!(result, wires(diagram))
  result
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 10349518-03ca-11eb-09b2-69c80c4662ac
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	using JuMP, Tulip, GLPK, LinearAlgebra
	using Distributions, Plots, StatsPlots, LaTeXStrings, Measures
end

# ╔═╡ f7f3a256-03c6-11eb-2c1e-83cc62bf55e6
md"""
# Linear programming
We will be using [JuMP](https://jump.dev/JuMP.jl/stable/) as a general framework, combined with [Tulip](https://github.com/ds4dm/Tulip.jl) or [GLPK](https://github.com/jump-dev/GLPK.jl) as a solver.
"""

# ╔═╡ 7910d53e-03ca-11eb-1536-7fc4c236e10a
md"""
## Application - Employee planning
We manage a crew of call center employees and want to optimise our shifts in order to reduce the total payroll cost. Employees have to work for five consecutive days and are then given two days off. The current policy is simple: each day gets the same amount of employees (currently we have 5 persons per shift, which leads to 25 persons on any given day).

We have some historical data that gives us the minimum amount of calls we can expect: Mon: 22, Tue: 17, Wed:13, Thu:14, Fri: 15, Sat: 18, Sun: 24

Employees are payed € 96 per day of work. This lead to the current payroll cost of 25x7x96 = € 16.800. You need to optimize employee planning to reduce the payroll cost.

| Schedule | Days worked | Attibuted Pers | Mon | Tue | Wed | Thu | Fri | Sat | Sun |
|----------|-------------|----------------|-----|-----|-----|-----|-----|-----|-----|
| A | Mon-Fri | 5 | W | W | W | W | W | O | O |
| B | Tue-Sat | 5 | O | W | W | W | W | W | O |
| C | Wed-Sun | 5 | O | O | W | W | W | W | W |
| D | Thu-Mon | 5 | W | O | O | W | W | W | W |
| E | Fri-Tue | 5 | W | W | O | O | W | W | W |
| F | Sat-Wed | 5 | W | W | W | O | O | W | W |
| G | Sun-Thu | 5 | W | W | W | W | O | O | W |
| Totals: | - | 35 | 5 | 5 | 5 | 5 | 5 | 5 | 5 |
| Required: | - | - | 22 | 17 | 13 | 14 | 15 | 18 | 24 |

### Mathematical formulation
We need do formaly define our decision variables, constraints and objective function.
* decision variables: the amount of persons attributed to each schedule ( ``Y = [y_1,y_2,\dots,y_7]^{\intercal} ``)
* objective function: the payroll cost

  Suppose the matrix ``A`` is the matrix indicating the workload for each schedule (in practice ``W=1`` and ``O=0``):
```math
A = \begin{bmatrix}
W & W & W & W & W & O & O \\
O & W & W & W & W & W & O \\
O & O & W & W & W & W & W \\
W & O & O & W & W & W & W \\
W & W & O & O & W & W & W \\
W & W & W & O & O & W & W \\
W & W & W & W & O & O & W 	\\
\end{bmatrix}
```

  Now $$Y*A'$$ gives us a vector indicating the amount of employees working on a given day. Suppose we also use the vector $$c$$ to indicate the salary for a given day (in this case $$c = [96,96,96,\dots,96]$$).

We are now able to write our objective function:
```math
\min Z = c^\intercal A^\intercal Y
```

* constraints (1): each day we need at least enough employees to cover all incoming calls. Suppose we use the vector $$b$$ to indicate the amount of incoming calls for a given day. We are able to write the constraints in a compact way:

```math
\text{subject to } A^\intercal Y  \ge b
```

* constraints (2): we also want to avoid a negative amount of attributed employees on any given day, since this would lead to a negative payroll cost:
```math
\text{and }Y \ge 0
```

* $\forall Y : Y \in \mathbb{N}$
### Implementation
"""

# ╔═╡ 572d14c4-03cb-11eb-2b68-234b3d7e9e8e
begin
	# basic data
	A = ones(Bool,7,7) - diagm(-1=>ones(Bool,6), -2=> ones(Bool,5), 5=>ones(Bool,2), 6=>ones(Bool,1))
	Y = [5,5,5,5,5,5,5]
	B = [22,17,13,14,15,18,24]
	C = [96,96,96,96,96,96,96];
	A
end

# ╔═╡ 4ba4687e-087f-11eb-1cb0-2581818cbd93
# A' * Y .> B
C' * A' * Y

# ╔═╡ 79351ee4-087f-11eb-147d-d389de200857
let
	model = Model(GLPK.Optimizer)
	@variable(model, Y[1:7] >= 0, Int)
	@constraint(model, A' * Y .>= B)
	@objective(model, Min, C' * A' * Y)
	optimize!(model)
	termination_status(model)
	objective_value(model)
	#A' * value.(Y) .> B
	sum(value.(Y))
end

# ╔═╡ 3bafd152-6724-4847-a79d-f2340e2ab5e4

# ╔═╡ 4bc4c2aa-04b4-11eb-2b6e-452d4ecc258a
md"""
### Adding uncertainty
Up to now, we have had constant numbers the minimum number of employees needed per day. In reality these quantities are uncertain. The actual number of calls will fluctuate each day. For simplicity's sake will we use a [lognormal distribution](https://en.wikipedia.org/wiki/Log-normal_distribution#Occurrence_and_applications) for the amount of calls (using their initial value as mean and a standard deviation of two). Working this way, we avoid having negative calls.
"""

# ╔═╡ 6e48e306-04b4-11eb-2561-0151a5e0a908
B

# ╔═╡ 7a54f9b4-04b4-11eb-3a7c-8d90eb026392
begin
	# generating the distributions
	B_u = Distributions.LogNormal.(B,2) # array with distributions

	# quick sample to illustrate amount of calls being randomized
	log.(rand.(B_u))
end

# ╔═╡ 9d3ceafe-0880-11eb-36f0-0f8530d6b285
begin
	cost = Float64[]
	for _ in 1:10000
		let cout=cost
			model = Model(GLPK.Optimizer)
			@variable(model, Y[1:7] >= 0, Int)
			@constraint(model, A' * Y .>= log.(rand.(B_u)))
			@objective(model, Min, C' * A' * Y)
			optimize!(model)
			push!(cout, objective_value(model))
		end
	end
	cost
end

# ╔═╡ 14a32a20-0881-11eb-1a20-715b99417266
StatsPlots.histogram(cost, yscale=:log10)

# ╔═╡ 259bbbc6-04b5-11eb-1ad4-c567e45ba4b6
md"""
### Small variant: adding a commission
Suppose each worker receives extra pay for the amount of calls that have been treated. We can easily include this in our model
"""

# ╔═╡ 83130656-0881-11eb-2c87-a3ae5f81f04b
begin
	cost_c = Float64[]
	commission = 20
	for _ in 1:1000
		let cout=cost_c
			model = Model(GLPK.Optimizer)
			appels = log.(rand.(B_u))
			@variable(model, Y[1:7] >= 0, Int)
			@constraint(model, A' * Y .>= appels)
			@objective(model, Min, C' * A' * Y + sum(appels) * commission)
			optimize!(model)
			push!(cout, objective_value(model))
		end
	end
	StatsPlots.histogram(cost_c)
end

# ╔═╡ e1e41ea6-04b5-11eb-174e-1d43f601a07c
md"""
#### Playing it safe
The above has given us some information on what the distributions of the payroll cost may be, however in reality, you would want to make sure that the clients calling to center are taken care off. To realise this, one might say that for any given day, you want to make sure that 90% of all calls can be treated by the specific capacity.
"""

# ╔═╡ 62b516be-0882-11eb-094b-d183f8d00ab8
log.(quantile.(B_u, 0.90))

# ╔═╡ 6274be66-0882-11eb-11be-05163cda6633
 let
	model = Model(GLPK.Optimizer)
	@variable(model, Y[1:7] >= 0, Int)
	@constraint(model, A' * Y .>= log.(quantile.(B_u, 0.99)))
	@objective(model, Min, C' * A' * Y)
	optimize!(model)
	termination_status(model)
	objective_value(model)
end

# ╔═╡ 1f177098-04b6-11eb-2508-6bd8d7e1e996
md"""
### Additional questions
* The example we have treated so far has very traditional working patterns for the employees. How woud you deal with modern working patterns (e.g. 4/5 or parttime working)?
* We took a look at the stochastic nature of the amount of calls, however, the personnel might not show up for various reasons. How would you describe the possible influence? Hint: make a discrete event model of this setting, using the optimal design and controlling for employees showing up or not.
"""

# ╔═╡ Cell order:
# ╟─f7f3a256-03c6-11eb-2c1e-83cc62bf55e6
# ╠═10349518-03ca-11eb-09b2-69c80c4662ac
# ╟─7910d53e-03ca-11eb-1536-7fc4c236e10a
# ╠═572d14c4-03cb-11eb-2b68-234b3d7e9e8e
# ╠═4ba4687e-087f-11eb-1cb0-2581818cbd93
# ╠═79351ee4-087f-11eb-147d-d389de200857
# ╠═3bafd152-6724-4847-a79d-f2340e2ab5e4
# ╟─4bc4c2aa-04b4-11eb-2b6e-452d4ecc258a
# ╠═6e48e306-04b4-11eb-2561-0151a5e0a908
# ╠═7a54f9b4-04b4-11eb-3a7c-8d90eb026392
# ╠═9d3ceafe-0880-11eb-36f0-0f8530d6b285
# ╠═14a32a20-0881-11eb-1a20-715b99417266
# ╟─259bbbc6-04b5-11eb-1ad4-c567e45ba4b6
# ╠═83130656-0881-11eb-2c87-a3ae5f81f04b
# ╟─e1e41ea6-04b5-11eb-174e-1d43f601a07c
# ╠═62b516be-0882-11eb-094b-d183f8d00ab8
# ╠═6274be66-0882-11eb-11be-05163cda6633
# ╟─1f177098-04b6-11eb-2508-6bd8d7e1e996

ext: jl

max_stars_repo_path: Exercises/PS05 - Linear programming I.jl


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@doc raw"""
    luminosity(c::RGB{T}) :: Float64

Return the best average luminosity of a color ``c = (R_i, G_i, B_i)``, through
the Shirley & Morley proposal:

```math
l_i = \frac{\max\{R_i, G_i, B_i \} + \min\{R_i, G_i, B_i \}}{2}
```
"""
luminosity(c::RGB{T}) where {T} = (max(c.r, c.g, c.b) + min(c.r, c.g, c.b))/2.

"""
    lum_max(img::HDRimage) :: Float64

Return the maximum luminosity of the given `img` according to the
`luminosity` function of a color.

See also: [`HDRimage`](@ref), [`luminosity`](@ref)
"""
function lum_max(img::HDRimage)
    lum_max=0.0
    for pix in img.rgb_m
        (lum_max>luminosity(pix)) || (lum_max=luminosity(pix))
    end
    lum_max
end

"""
    avg_lum(img::HDRimage, δ::Number=1e-10) :: Float64

Return the average luminosity of the `img` according to the
`luminosity` function of a color.
The `δ` parameter is used to prevent  numerical problems for
under-illuminated pixels.

See also: [`HDRimage`](@ref), [`luminosity`](@ref)
"""
function avg_lum(img::HDRimage, δ::Number=1e-10)
    cumsum=0.0
    for pix in img.rgb_m
        cumsum += log10(δ + luminosity(pix))
    end
    10^(cumsum/(img.width*img.height))
end

@doc raw"""
    normalize_image!(
            img::HDRimage,
            a::Float64 = 0.18,
            lum::Union{Number, Nothing} = nothing,
            δ::Number = 1e-10
            )

Normalize the `img` colors through the following formula:

```math
\forall \; \mathrm{colors} \, c =(R_i, G_i, B_i) \;
\mathrm{of \, the \, image}   \\\\

X_i \rightarrow \frac{a}{\langle l \rangle} X_i
    \quad , \quad
\forall X_i = R_i, G_i, B_i
```
where ``\langle l \rangle`` is the average luminosity returned
by the `avg_lum` function.

See also: [`HDRimage`](@ref), [`avg_lum`](@ref)
"""
function normalize_image!(
            img::HDRimage,
            a::Float64=0.18,
            lum::Union{Number, Nothing}=nothing,
            δ::Number=1e-10
            )

    (isnothing(lum) || lum≈0.0 ) && (lum = avg_lum(img, δ))
    img.rgb_m .= img.rgb_m .* a .* (1.0/lum)
    nothing
end

@doc raw"""
    clamp(x::Number) :: Float64

Clamping function:

```math
x \rightarrow x/(x+1)
```
"""
clamp(x::Number) = x/(x+1)

"""
    clamp_image!(img::HDRimage)

Adjust the color levels of the brightest pixels in the `img`
through the `clamp` function.

See also: [`HDRimage`](@ref), [`clamp`](@ref)
"""
function clamp_image!(img::HDRimage)
    h=img.height
    w=img.width
    for y in h-1:-1:0, x in 0:w-1
        col = get_pixel(img, x, y)
        T = typeof(col).parameters[1]
        new_col = RGB{T}( clamp(col.r), clamp(col.g), clamp(col.b) )
        set_pixel(img, x,y, new_col)
    end
    nothing
end

"""
    γ_correction!(img::HDRimage, γ::Float64=1.0, k::Float64=1.0)

Corrects the image using the `γ` factor, assuming a potential dependence
between the input and output signals of a monitor/screen.

As third optional argument, you can pass the maximum value 'k' of the range
you want the RGB colors may have. The default value is 'k=1.0', so the range
RGB colors can span is '[0.0, 1.0]'

See also: [`HDRimage`](@ref)
"""
function γ_correction!(img::HDRimage, γ::Float64=1.0, k::Float64=1.0)
    h=img.height
    w=img.width
    for y in h-1:-1:0, x in 0:w-1
        cur_color = get_pixel(img, x, y)
        T = typeof(cur_color).parameters[1]
        new_col = RGB{T}( floor(255 * cur_color.r^(1/γ)),
                          floor(255 * cur_color.g^(1/γ)),
                          floor(255 * cur_color.b^(1/γ))
        )
        set_pixel(img, x, y, k/255.0*new_col)
    end
    nothing
end

##########################################################################################92

function tone_mapping(x::(Pair{T1,T2} where {T1,T2})...)
	tone_mapping( parse_tonemapping_settings(  Dict( pair for pair in [x...]) )... )
end

function tone_mapping(
            infile::String,
            outfile::String,
            a::Float64 = 0.18,
            γ::Float64 = 1.0,
            lum::Union{Number, Nothing} = nothing,
            ONLY_FOR_TESTS::Bool=false
            )
    (ONLY_FOR_TESTS==false) || (return nothing)
    tone_mapping(["$(infile)", "$(outfile)", "$a", "$γ", "$lum"])
end

function tone_mapping(args_lum::Vector{String})
    correct_usage =
        "\ncorrect usage of tone mapping function for vector of string arguments:\n"*
        "julia>  tonemapping([\"infile\",\"outfile\" ])\n"*
        "julia>  tonemapping([\"infile\",\"outfile\", \"a\"])\n"*
        "julia>  tonemapping([\"infile\",\"outfile\",  \"a\",  \"γ\" ])\n\n"*
        "julia>  tonemapping([\"infile\",\"outfile\",  \"a\",  \"γ\", \"lum\" ])\n\n"*
        "default values are a=0.18, γ=1.0, lum=nothing\n\n"
    if isempty(args_lum) || length(args_lum)==1 || length(args_lum)>5
        throw(ArgumentError(correct_usage))
		return nothing
    end

	parameters = nothing
	try
		parameters =  Parameters(parse_command_line(args_lum[begin:end-1])...)
	catch e
		println("Error: ", e)
        println(correct_usage)
		return nothing
	end

	img = open(parameters.infile, "r") do inpf; read(inpf, HDRimage); end

	println("\nfile $(parameters.infile) has been read from disk.\n")

	normalize_image!(img, parameters.a, string2positive(args_lum[end]))
	clamp_image!(img)
	Raytracing.γ_correction!(img, parameters.γ)

	#println(img, 3)

	matrix = get_matrix(img)
	Images.save(parameters.outfile, matrix)

	println("\nFile $(parameters.outfile) has been written into the disk.\n")
end

"""
$(SIGNATURES)

Tone-map the given input pfm file `infile` with luminosity
normalisation `a`, gamma factor `γ` and average luminosity `lum`.
Return a file with the input outfile name `outfile` and
of the specified LDR format, if possible.

In order to do the tone-mapping, this function relies on the
following three function, called in the presented order:

- `normalize_image!` : normalize the image colors
- `clamp_image!` : adjust the color levels of the brightest pixels
- `γ_correction!` : corrects the image using the `γ` factor
- `avg_lum` : if not specified an average luminosity `lum` (or if it is =0),
  it is calculated through this function

See also: [`normalize_image!`](@ref), [`clamp_image!`](@ref),
[`γ_correction!`](@ref), [`avg_lum`](@ref)
"""
tone_mapping


text:
# Query.jl is an embarrassment when it comes to group by Performance

# using StatsBase, RCall, Query, DataFrames, DataFramesMeta, Distributions, PooledArrays, IndexedTables
#
# function query_benches(df::DataFrame)
#
#     # timings
#     ti = Dict()
#
#     ti[:sum1] = @elapsed @from i in df begin
#                              @group i by i.id1 into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum2] = @elapsed @from i in df begin
#                              @group i by i.id1 into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum3] = @elapsed @from i in df begin
#                              @group i by (i.id1,i.id2) into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum4] = @elapsed @from i in df begin
#                              @group i by (i.id1,i.id2) into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum_mean1] = @elapsed @from i in df begin
#                              @group i by i.id3 into g
#                              @select {s=sum(g..v1),m=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:sum_mean2] = @elapsed @from i in df begin
#                              @group i by i.id3 into g
#                              @select {s=sum(g..v1),m=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:mean7_9_by_id4_1] = @elapsed @from i in df begin
#                              @group i by i.id4 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:mean7_9_by_id4_2] = @elapsed @from i in df begin
#                              @group i by i.id4 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:sum7_9_by_id6_1] = @elapsed @from i in df begin
#                              @group i by i.id6 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:sum7_9_by_id6_2] = @elapsed @from i in df begin
#                              @group i by i.id6 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     return ti
# end
#
# function query_benches(df::IndexedTable)
#
#     # timings
#     ti = Dict()
#
#     ti[:sum1] = @elapsed @from i in df begin
#                              @group i by i.id1 into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum2] = @elapsed @from i in df begin
#                              @group i by i.id1 into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum3] = @elapsed @from i in df begin
#                              @group i by (i.id1,i.id2) into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum4] = @elapsed @from i in df begin
#                              @group i by (i.id1,i.id2) into g
#                              @select {r=sum(g..v1)}
#                              @collect DataFrame
#                          end
#     ti[:sum_mean1] = @elapsed @from i in df begin
#                              @group i by i.id3 into g
#                              @select {s=sum(g..v1),m=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:sum_mean2] = @elapsed @from i in df begin
#                              @group i by i.id3 into g
#                              @select {s=sum(g..v1),m=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:mean7_9_by_id4_1] = @elapsed @from i in df begin
#                              @group i by i.id4 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:mean7_9_by_id4_2] = @elapsed @from i in df begin
#                              @group i by i.id4 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:sum7_9_by_id6_1] = @elapsed @from i in df begin
#                              @group i by i.id6 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     ti[:sum7_9_by_id6_2] = @elapsed @from i in df begin
#                              @group i by i.id6 into g
#                              @select {m7=mean(g..v1),m8=mean(g..v2),m9=mean(g..v3)}
#                              @collect DataFrame
#                          end
#     return ti
# end
#
#
# function run_query_benches(N=1_000_000;K=100)
#     # get small data for JIT warmup
#     d_ = createSynDataFrame(10,3)
#     # dt_ = createDataTable(10,3)
#
#
#     # warm up julia benchmarks
#     query_benches(d_);
#     # benches(dt_);
#     benches(di_);
#     DfMeta_benches(d_);
#
#     # get real data
#     d = createDataFrame(N,K)
#     # measure DataFrames
#     query_df = benches(d)
#     meta = DfMeta_benches(d)
#     d = 0
#     gc()
#     # measure DataTables
#     # dt = createDataTable(N,K)
#     # query_dt = benches(dt)
#     # dt = 0
#     # gc()
#     # measure IndexedTables
#     di = createIndexedTable(N,K)
#     query_di = benches(di)
#     di = 0
#     gc()
#
#
#     # get R time
#     R = R_bench(N,K)
#
#     # get
#     k = collect(keys(query_df))
#     #out = DataFrame(bench = k,Query_DF = [query_df[kk] for kk in k],Query_DT= [query_dt[kk] for kk in k],Query_idxT= [query_di[kk] for kk in k],DataFramesMeta=[meta[kk] for kk in k], Rdatatable=[R[kk] for kk in k])
#     out = DataFrame(bench = k,Query_DF = [query_df[kk] for kk in k] ,Query_idxT= [query_di[kk] for kk in k],DataFramesMeta=[meta[kk] for kk in k], Rdatatable=[R[kk] for kk in k])
#     sort!(out,cols=:bench)
#     rel = deepcopy(out)
#     #rel = @transform(rel,Query_DF = :Query_DF./:Rdatatable, Query_DT = :Query_DT./:Rdatatable, Query_idxT = :Query_idxT./:Rdatatable, DataFramesMeta = :DataFramesMeta ./ :Rdatatable, Rdatatable = 1.0)
#     rel = @transform(rel,Query_DF = :Query_DF./:Rdatatable, Query_idxT = :Query_idxT./:Rdatatable, DataFramesMeta = :DataFramesMeta ./ :Rdatatable, Rdatatable = 1.0)
#     return (out,rel)
# end
#
# function run_all()
#     d=Dict()
#     for n in [1_000_000, 2_000_000]
#         d[n] = run_benches(n)
#         gc()
#         info("results with $n rows:")
#         info("-- absolute time in seconds:")
#         println(d[n][1])
#         info("-- time relative to R data.table:")
#         println(d[n][2])
#         println()
#         println()
#     end
#     return d
# end


text:
# Ported to Julia by Andy Ferris, 2016, and re-released under an MIT license.
#/**
# * Copyright (c) Charles Karney (2008-2015) <charles@karney.com> and licensed
# * under the MIT/X11 License.  For more information, see
# * http://geographiclib.sourceforge.net/
# **********************************************************************/

"""
    (x, y, gamma, k) = polarst_fwd(northpole::Bool, k0::Float64, tm::TransverseMercator, lat, lon)

Perform polar-stereographic projection of `lat` and `lon` with respect to north
or south pole `northpole` and horizontal scaling `k0` (`= 0.994` for UPS).
`γ` and `k` are the local convergence and scaling factors, respectively.
"""
function polarst_fwd(northp::Bool, k0::Float64, tm::TransverseMercator, lat, lon) # k0 is scale factor...
    lat = LatFix(lat) * (northp ? 1 : -1)

    tau = tand(lat)
    secphi = hypot(1.0, tau)
    taup = taupf(tau, tm.e2) # TODO revert to C++ es here?
    rho = hypot(1.0, taup) + abs(taup)
    rho = (taup >= 0 ? (lat != 90 ? 1/rho : 0.0) : rho)
    rho *= 2 * k0 * tm.a / tm.c
    k = (lat != 90 ? (rho / tm.a) * secphi * sqrt(tm.e2m + tm.e2 / (secphi*secphi)) : k0)
    x = sind(lon)
    y = cosd(lon)
    x *= rho
    y *= (northp ? -rho : rho)
    gamma = AngNormalize(northp ? lon : -lon)

    return (x, y, gamma, k)
end

"""
    (lat, lon, gamma, k) = polarst_inv(northp::Bool, k0::Float64, tm::TransverseMercator, x, y)

Invert polar-stereographic projection of `x` and `y` with respect to north
or south pole `northpole` and horizontal scaling `k0` (`= 0.994` for UPS).
`γ` and `k` are the local convergence and scaling factors, respectively.
"""
function polarst_inv(northp::Bool, k0::Float64, tm::TransverseMercator, x, y) # k0 is scale factor...
    rho = hypot(x, y)
    t = (rho > 0 ? rho / (2 * k0 * tm.a / tm.c) : eps(Float64)*eps(Float64))
    taup = (1 / t - t) / 2
    tau = tauf(taup, tm.e2) # TODO revert to C++ es here?
    secphi = hypot(1.0, tau)
    k = (rho > 0 ? (rho / tm.a) * secphi * sqrt(tm.e2m + tm.e2 / (secphi*secphi)) : k0)
    lat = (northp ? 1 : -1) * atand(tau)
    lon = atan2(x, northp ? -y : y ) * 180/pi
    gamma = AngNormalize(northp ? lon : -lon)

    return (lat, lon, gamma, k)
end


text:
const HAS_ATOMIC = VERSION ≥ v"1.7"

macro if_has_atomic(ex)
    if HAS_ATOMIC
        esc(ex)
    end
end

mutable struct StackNode{T}
    object::T
    next::Union{Nothing,StackNode{T}}
end

function _CentralizedCache end


text:
typeof_model_elements(model::StateSpaceModel) = eltype(model.system.y)
num_states(model::StateSpaceModel) = num_states(system(model))
num_series(model::StateSpaceModel) = num_series(system(model))
system(model::StateSpaceModel) = model.system
isunivariate(model::StateSpaceModel) = isa(model.system.y, Vector)
model_name(model::StateSpaceModel) = "$(typeof(model))"
length_observations(model::StateSpaceModel) = length(model.system.y)
observations(model::StateSpaceModel) = model.system.y

function lagmat(y::Vector{Fl}, k::Int) where Fl
    X = Matrix{Fl}(undef, length(y) - k, k)
    for i in 1:k
        X[:, i] = lag(y, i)[k + 1:end]
    end
    return X
end

function variance_of_valid_observations(y::Vector{Fl}) where Fl
    return var(filter(!isnan, y))
end

function mean_of_valid_observations(y::Vector{Fl}) where Fl
    return mean(filter(!isnan, y))
end

function assert_zero_missing_values(y::Vector)
    for el in y
        if isnan(el)
            return error("This model does not support missing values.")
        end
    end
    return nothing
end


text:
using NESSie.Format
using LightXML: parse_string, root, name, child_elements, attribute, content

@testset "writevtk (surface model)" begin
    for T in testtypes
        # empty system
        model = Model{T, Triangle{T}}()
        (fname, fh) = mktemp()
        try
            xroot = root(parse_string(readback(fh -> writevtk(fh, model))))
            @test name(xroot) == "VTKFile"
            @test attribute(xroot, "type") == "PolyData"
            xgrid = collect(child_elements(xroot))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "PolyData"
            xgrid = collect(child_elements(xgrid[1]))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "Piece"
            @test attribute(xgrid[1], "NumberOfPoints") == "0"
            @test attribute(xgrid[1], "NumberOfVerts") == "0"
            @test attribute(xgrid[1], "NumberOfLines") == "0"
            @test attribute(xgrid[1], "NumberOfStrips") == "0"
            @test attribute(xgrid[1], "NumberOfPolys") == "0"
            xpiece = collect(child_elements(xgrid[1]))
            @test length(xpiece) == 2
            @test name(xpiece[1]) == "Points"
            xdata = collect(child_elements(xpiece[1]))
            @test length(xdata) == 1
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "$T"
            @test attribute(xdata[1], "NumberOfComponents") == "3"
            @test attribute(xdata[1], "format") == "ascii"
            @test strip(content(xdata[1])) == ""
            @test name(xpiece[2]) == "Polys"
            xdata = collect(child_elements(xpiece[2]))
            @test length(xdata) == 2
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "Int32"
            @test attribute(xdata[1], "Name") == "connectivity"
            @test attribute(xdata[1], "format") == "ascii"
            @test strip(content(xdata[1])) == ""
            @test name(xdata[2]) == "DataArray"
            @test attribute(xdata[2], "type") == "Int32"
            @test attribute(xdata[2], "Name") == "offsets"
            @test attribute(xdata[2], "format") == "ascii"
            @test strip(content(xdata[2])) == ""
        finally
            close(fh)
            rm(fname)
        end

        # small system
        nodes = Vector{T}[[0, 0, 0], [0, 0, 3], [0, 3, 0], [1, -3, 3]]
        model = Model(
            nodes,
            Triangle{T}[
                Triangle(nodes[1], nodes[2], nodes[3]),
                Triangle(nodes[1], nodes[4], nodes[2])
            ]
        )
        (fname, fh) = mktemp()
        try
            xroot = root(parse_string(readback(fh -> writevtk(fh, model))))
            @test name(xroot) == "VTKFile"
            @test attribute(xroot, "type") == "PolyData"
            xgrid = collect(child_elements(xroot))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "PolyData"
            xgrid = collect(child_elements(xgrid[1]))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "Piece"
            @test attribute(xgrid[1], "NumberOfPoints") == "4"
            @test attribute(xgrid[1], "NumberOfVerts") == "0"
            @test attribute(xgrid[1], "NumberOfLines") == "0"
            @test attribute(xgrid[1], "NumberOfStrips") == "0"
            @test attribute(xgrid[1], "NumberOfPolys") == "2"
            xpiece = collect(child_elements(xgrid[1]))
            @test length(xpiece) == 2
            @test name(xpiece[1]) == "Points"
            xdata = collect(child_elements(xpiece[1]))
            @test length(xdata) == 1
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "$T"
            @test attribute(xdata[1], "NumberOfComponents") == "3"
            @test attribute(xdata[1], "format") == "ascii"
            @test [parse(T, e) for e in split(strip(content(xdata[1])))] == T[0, 0, 0, 0, 0, 3, 0, 3, 0, 1, -3, 3]
            @test name(xpiece[2]) == "Polys"
            xdata = collect(child_elements(xpiece[2]))
            @test length(xdata) == 2
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "Int32"
            @test attribute(xdata[1], "Name") == "connectivity"
            @test attribute(xdata[1], "format") == "ascii"
            @test [parse(Int, e) for e in split(strip(content(xdata[1])))] == [0, 1, 2, 0, 3, 1]
            @test name(xdata[2]) == "DataArray"
            @test attribute(xdata[2], "type") == "Int32"
            @test attribute(xdata[2], "Name") == "offsets"
            @test attribute(xdata[2], "format") == "ascii"
            @test [parse(Int, e) for e in split(strip(content(xdata[2])))] == [3, 6]
        finally
            close(fh)
            rm(fname)
        end
    end
end

@testset "writevtk (volume model)" begin
    for T in testtypes
        # empty system
        model = Model{T, Tetrahedron{T}}()
        (fname, fh) = mktemp()
        try
            xroot = root(parse_string(readback(fh -> writevtk(fh, model))))
            @test name(xroot) == "VTKFile"
            @test attribute(xroot, "type") == "UnstructuredGrid"
            xgrid = collect(child_elements(xroot))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "UnstructuredGrid"
            xgrid = collect(child_elements(xgrid[1]))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "Piece"
            @test attribute(xgrid[1], "NumberOfPoints") == "0"
            @test attribute(xgrid[1], "NumberOfCells") == "0"
            xpiece = collect(child_elements(xgrid[1]))
            @test length(xpiece) == 2
            @test name(xpiece[1]) == "Points"
            xdata = collect(child_elements(xpiece[1]))
            @test length(xdata) == 1
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "$T"
            @test attribute(xdata[1], "NumberOfComponents") == "3"
            @test attribute(xdata[1], "format") == "ascii"
            @test strip(content(xdata[1])) == ""
            @test name(xpiece[2]) == "Cells"
            xdata = collect(child_elements(xpiece[2]))
            @test length(xdata) == 3
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "Int32"
            @test attribute(xdata[1], "Name") == "connectivity"
            @test attribute(xdata[1], "format") == "ascii"
            @test strip(content(xdata[1])) == ""
            @test name(xdata[2]) == "DataArray"
            @test attribute(xdata[2], "type") == "Int32"
            @test attribute(xdata[2], "Name") == "offsets"
            @test attribute(xdata[2], "format") == "ascii"
            @test strip(content(xdata[2])) == ""
            @test name(xdata[3]) == "DataArray"
            @test attribute(xdata[3], "type") == "UInt8"
            @test attribute(xdata[3], "Name") == "types"
            @test attribute(xdata[3], "format") == "ascii"
            @test strip(content(xdata[3])) == ""
        finally
            close(fh)
            rm(fname)
        end

        # small system
        nodes = Vector{T}[[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [-2, 0, 0]]
        model = Model(
            nodes,
            Tetrahedron{T}[
                Tetrahedron(nodes[1], nodes[2], nodes[3], nodes[4]),
                Tetrahedron(nodes[1], nodes[2], nodes[3], nodes[5])
            ]
        )
        (fname, fh) = mktemp()
        try
            xroot = root(parse_string(readback(fh -> writevtk(fh, model))))
            @test name(xroot) == "VTKFile"
            @test attribute(xroot, "type") == "UnstructuredGrid"
            xgrid = collect(child_elements(xroot))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "UnstructuredGrid"
            xgrid = collect(child_elements(xgrid[1]))
            @test length(xgrid) == 1
            @test name(xgrid[1]) == "Piece"
            @test attribute(xgrid[1], "NumberOfPoints") == "5"
            @test attribute(xgrid[1], "NumberOfCells") == "2"
            xpiece = collect(child_elements(xgrid[1]))
            @test length(xpiece) == 2
            @test name(xpiece[1]) == "Points"
            xdata = collect(child_elements(xpiece[1]))
            @test length(xdata) == 1
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "$T"
            @test attribute(xdata[1], "NumberOfComponents") == "3"
            @test attribute(xdata[1], "format") == "ascii"
            @test [parse(T, e) for e in split(strip(content(xdata[1])))] == T[0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, -2, 0, 0]
            @test name(xpiece[2]) == "Cells"
            xdata = collect(child_elements(xpiece[2]))
            @test length(xdata) == 3
            @test name(xdata[1]) == "DataArray"
            @test attribute(xdata[1], "type") == "Int32"
            @test attribute(xdata[1], "Name") == "connectivity"
            @test attribute(xdata[1], "format") == "ascii"
            @test [parse(Int, e) for e in split(strip(content(xdata[1])))] == [0, 1, 2, 3, 0, 1, 2, 4]
            @test name(xdata[2]) == "DataArray"
            @test attribute(xdata[2], "type") == "Int32"
            @test attribute(xdata[2], "Name") == "offsets"
            @test attribute(xdata[2], "format") == "ascii"
            @test [parse(Int, e) for e in split(strip(content(xdata[2])))] == [4, 8]
            @test name(xdata[3]) == "DataArray"
            @test attribute(xdata[3], "type") == "UInt8"
            @test attribute(xdata[3], "Name") == "types"
            @test attribute(xdata[3], "format") == "ascii"
            @test [parse(Int, e) for e in split(strip(content(xdata[3])))] == [10, 10]
        finally
            close(fh)
            rm(fname)
        end
    end
end