text:
# Some useful function for the multiple user occluded crosswalk environment

#### POMDP related stuff #################################################################################
# Dummy type since the action is given by the policy
mutable struct EgoDriver{A} <: DriverModel{A}
    a::A
end

Base.rand(model::EgoDriver) = model.a
AutomotiveDrivingModels.observe!(m::EgoDriver, s::EntityFrame{S,D,I}, roadway::R, egoid::Int64) where {S,D,I,R} = m

"""
    is_crash(scene::Scene)
return true if the ego car is in collision in the given scene, do not check for collisions between
other participants
"""
function is_crash(scene::Scene)
    ego = scene[findfirst(EGO_ID, scene)]
    for veh in scene
        if veh.id != 1
            if collision_checker(ego, veh) #&& !(ego.state.v == 0. && veh.def.class == AgentClass.PEDESTRIAN)
                return true
            end
        end
    end
    return false
end

"""
Remove vehicles that have reached an end scene or that have collided
"""
function clean_scene!(env::OccludedEnv, scene::Scene)
    exit_lanes = get_exit_lanes(env.roadway)
    for veh in scene
        if veh.id == 1
            continue
        end
        s = veh.state.posF.s
        lane = get_lane(env.roadway, veh)
        s_end = get_end(lane)
        if s >= s_end && lane ∈ exit_lanes
            deleteat!(scene, findfirst(veh.id, scene))
        elseif veh.def.class == AgentClass.PEDESTRIAN
            if isapprox(s, 0., atol=0.01) && veh.state.posF.ϕ ≈ float(π)
                deleteat!(scene, findfirst(veh.id, scene))
            elseif s >= s_end && veh.state.posF.ϕ ≈ 0.
                deleteat!(scene, findfirst(veh.id, scene))
            end
        end
    end
    for veh in scene
        for veh_ in scene
            if collision_checker(veh, veh_) &&
               veh.id != 1 && veh_.id != 1 &&
               veh.id != veh_.id &&
               veh_.def.class != AgentClass.PEDESTRIAN && veh.def.class != AgentClass.PEDESTRIAN &&
               !(min(veh_.state.v, veh.state.v) ≈ 0.)
                deleteat!(scene, findfirst(veh_.id, scene))
            end
        end
    end
end

"""
Heuristic to check if a vehicle can be added to a given scene
"""
function can_push(env, scene::Scene, car::Vehicle, delta_headway::Float64 = 2.0)
    car_lane = get_lane(env.roadway, car)
    min_headway = Inf
    headway = min_headway
    for veh in scene
        if is_colliding(car, veh) || car.id == veh.id
            return false
        end
        veh_lane = get_lane(env.roadway, veh)
        if veh_lane == car_lane
            headway = veh.state.posF.s - car.state.posF.s
        end
        if headway <= min_headway
            min_headway = headway
        end
    end
    if min_headway < delta_headway
        return false
    end
    return true
end

"""
return the ego vehicle from a given scene
"""
function get_ego(scene::Scene)
    return scene[findfirst(EGO_ID, scene)]
end

### ROADWAY STUFF

"""
    get_end(lane::Lane)
return the end longitudinal position of a lane
"""
function get_end(lane::Lane)
    s_end = round(lane.curve[end].s)
    return s_end
end

"""
return the list of entrance lanes, i.e. the one that does not have any entrances
"""
function get_start_lanes(roadway::Roadway{T}) where T<:Real
    lanes = Lane{T}[]
    for i=1:length(roadway.segments)
        for lane in roadway.segments[i].lanes
            if isempty(lane.entrances) && !is_crosswalk(lane) && !(lane.tag ∈ [LaneTag(15,1), LaneTag(16, 1), LaneTag(6,1), LaneTag(13, 1), LaneTag(14,1)])
                push!(lanes, lane)
            end
        end
    end
    return lanes
end

"""
    is_crosswalk(lane::Lane)
returns true if a lane is a crosswalk
"""
is_crosswalk(lane::Lane) = lane.width > 3.0  #XXX hack

"""
    get_lanes(roadway::Roadway)
return all the lanes in a roadway that are not crosswalk and not in the ego path
"""
function get_lanes(roadway::Roadway{T}) where T<:Real
    lanes = Lane{T}[]
    for i=1:length(roadway.segments)
        for lane in roadway.segments[i].lanes
            if !is_crosswalk(lane) && !(lane.tag ∈ [LaneTag(15,1), LaneTag(16, 1), LaneTag(6,1), LaneTag(13, 1), LaneTag(14,1)]) #XXX hack
                push!(lanes, lane)
            end
        end
    end
    return lanes
end

"""
    get_all_lanes(roadway::Roadway)
returns a list of all the lanes present in `roadway`
"""
function get_all_lanes(roadway::Roadway{T}) where T<:Real
    lanes = Lane{T}[]
    for i=1:length(roadway.segments)
        for lane in roadway.segments[i].lanes
            push!(lanes, lane)
        end
    end
    return lanes
end

"""
return the list of exit lanes, i.e. the one that does not have any exits
"""
function get_exit_lanes(roadway::Roadway{T}) where T
    lanes = Lane{T}[]
    for i=1:length(roadway.segments)
        for lane in roadway.segments[i].lanes
            if isempty(lane.exits)
                push!(lanes, lane)
            end
        end
    end
    return lanes
end

"""
Given a list of lanes returns the list of exits
"""
function get_exit_lanes(lanes::Vector{Lane{T}}, roadway::Roadway{T}) where T <: Real
    exits = Lane{T}[]
    for lane in lanes
        for exit in lane.exits
            tag = exit.target.tag
            new_exit = roadway[tag]
            if !(new_exit ∈ exits)
                push!(exits, new_exit)
            else
                continue
            end
        end
    end
    return exits
end

"""
    lane_to_segment(lane::Lane, roadway::Roadway)
Converts a lane to a LineSegment.
It only works for straight lanes!
"""
function lane_to_segment(lane::Lane, roadway::Roadway)
    # only works for straight lanes
    lane_a = lane.curve[1].pos #posg(Frenet(lane, 0.), roadway)
    lane_b = lane.curve[end].pos #posg(Frenet(lane, get_end(lane)), roadway)
    return LineSegment(lane_a, lane_b)
end

"""
    get_conflict_lanes(crosswalk::Lane, roadway::Roadway)
Find the lanes in roadway that intersects with a crosswalk
"""
function  get_conflict_lanes(crosswalk::Lane{T}, roadway::Roadway{T}) where T<:Real
    # find lane intersecting with crosswalk
    cw_seg = lane_to_segment(crosswalk, roadway)
    conflict_lanes = Lane{T}[]
    lanes = get_all_lanes(roadway)
    for lane in lanes
        if  is_crosswalk(lane)
            continue
        end
        lane_seg = lane_to_segment(lane, roadway)
        if intersects(lane_seg, cw_seg) && !(lane ∈ conflict_lanes)
            push!(conflict_lanes, lane)
        end
    end
    return conflict_lanes
end

"""
    random_route(rng::AbstractRNG, roadway::Roadway, start_lane::Lane)
Generate a random route starting from start_lane to a random end node
"""
function random_route(rng::AbstractRNG, roadway::Roadway, start_lane::Lane{T}) where T
    lanes = Lane{T}[start_lane]
    cur_lane = start_lane
    while !isempty(cur_lane.exits)
        rand_exit = rand(rng, cur_lane.exits)
        next_lane_tag = rand_exit.target.tag
        next_lane = roadway[next_lane_tag]
        push!(lanes, next_lane)
        cur_lane = next_lane
    end
    return lanes
end

# return +1 if going toward, -1 if going away
function direction_from_center(ped::Vehicle, crosswalk::Lane)
    direction_from_center(ped.state, crosswalk)
end
function direction_from_center(ped::VehicleState, crosswalk::Lane)
    s_ped = ped.posF.s
    Δs = get_end(crosswalk)/2 - s_ped
    return sign(Δs*cos(ped.posF.ϕ))
end

struct LonAccelDirection
    a_lon::Float64
    direction::Int64
end

function Base.hash(a::LonAccelDirection, h::UInt)
    return hash(a.a_lon, hash(a.direction, h))
end
function Base.:(==)(a1::LonAccelDirection, a2::LonAccelDirection)
    return a1.a_lon == a2.a_lon && a1.direction == a2.direction
end

function AutomotiveDrivingModels.propagate(veh::Entity{VehicleState, D, Int}, action::LonAccelDirection,  roadway::Roadway, Δt::Float64) where {D<:Union{VehicleDef, BicycleModel}}
    previousInd = veh.state.posF.roadind
    a_lon = action.a_lon
    v = veh.state.v
    a_lat = 0.

    v = veh.state.v
       ϕ = veh.state.posF.ϕ
      ds = v*cos(ϕ)
       t = veh.state.posF.t
      dt = v*sin(ϕ)

      ΔT = Δt
      ΔT² = ΔT*ΔT
      Δs = ds*ΔT + 0.5*a_lon*ΔT²
      Δs = max(0, Δs) # no backup
      Δt = dt*ΔT + 0.5*a_lat*ΔT²

      ds₂ = ds + a_lon*ΔT
    ds₂ = max(0, ds₂) # no backup
      dt₂ = dt + a_lat*ΔT
      speed₂ = sqrt(dt₂*dt₂ + ds₂*ds₂)
      v₂ = sign(ds₂)*speed₂
      ϕ₂ = atan(dt₂, ds₂) + (v₂ < 0.0)*π # handle negative speeds
      roadind = AutoUrban.move_along_with_direction(veh.state.posF.roadind, roadway, Δs, direction = action.direction)
      footpoint = roadway[roadind]
      posF = Frenet(roadind, roadway, t=t+Δt, ϕ = ϕ₂)
      # posG = convert(VecE2, footpoint.pos) + polar(t + Δt, footpoint.pos.θ + π/2)
      # posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)

    state = VehicleState(posF, roadway, v₂)
    # state = VehicleState(posG, roadway[roadind.tag], roadway, v₂)
    return state
end

function AutomotiveDrivingModels.propagate(veh::VehicleState, action::LonAccelDirection,  roadway::Roadway, Δt::Float64)
    previousInd = veh.posF.roadind
    a_lon = action.a_lon
    v = veh.v
    a_lat = 0.

    v = veh.v
       ϕ = veh.posF.ϕ
      ds = v*cos(ϕ)
       t = veh.posF.t
      dt = v*sin(ϕ)

      ΔT = Δt
      ΔT² = ΔT*ΔT
      Δs = ds*ΔT + 0.5*a_lon*ΔT²
      Δs = max(0, Δs) # no backup
      Δt = dt*ΔT + 0.5*a_lat*ΔT²

      ds₂ = ds + a_lon*ΔT
    ds₂ = max(0, ds₂) # no backup
      dt₂ = dt + a_lat*ΔT
      speed₂ = sqrt(dt₂*dt₂ + ds₂*ds₂)
      v₂ = sign(ds₂)*speed₂
      ϕ₂ = atan(dt₂, ds₂) + (v₂ < 0.0)*π # handle negative speeds
      roadind = AutoUrban.move_along_with_direction(veh.posF.roadind, roadway, Δs, direction = action.direction)
      footpoint = roadway[roadind]
      posF = Frenet(roadind, roadway, t=t+Δt, ϕ = ϕ₂)
      posG = convert(VecE2, footpoint.pos) + polar(t + Δt, footpoint.pos.θ + π/2)
      posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)

    state = VehicleState(posG, posF, v₂)
    # state = VehicleState(posG, roadway[roadind.tag], roadway, v₂)
    return state
end

mutable struct LonAccel
    a_lon::Float64
end

function AutomotiveDrivingModels.propagate(veh::Entity{VehicleState, D, Int}, action::LonAccel,  roadway::Roadway, ΔT::Float64) where {D<:Union{VehicleDef, BicycleModel}}
    a_lat = action.a_lat
   a_lon = action.a_lon

    v = veh.state.v
    ϕ = veh.state.posF.ϕ
   ds = v*cos(ϕ)
    t = veh.state.posF.t
   dt = v*sin(ϕ)

   ΔT² = ΔT*ΔT
   Δs = ds*ΔT + 0.5*a_lon*ΔT²
   Δt = dt*ΔT + 0.5*a_lat*ΔT²

   ds₂ = ds + a_lon*ΔT
   dt₂ = dt + a_lat*ΔT
   speed₂ = sqrt(dt₂*dt₂ + ds₂*ds₂)
   v₂ = sqrt(dt₂*dt₂ + ds₂*ds₂) # v is the magnitude of the velocity vector
   ϕ₂ = atan(dt₂, ds₂)

   roadind = move_along(veh.state.posF.roadind, roadway, Δs)
   footpoint = roadway[roadind]

   posG = convert(VecE2, footpoint.pos) + polar(t + Δt, footpoint.pos.θ + π/2)

   posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)

   return VehicleState(posG, roadway, v₂)
end

"""
Returns a Frenet object of the end of the road
Args: Roadway, LaneTag
Returns: Frenet
"""
function get_end_frenet(roadway::Roadway, tag::LaneTag)
    return Frenet(roadway[tag], get_end(roadway[tag]))
end


text:
"""
    struct TS
      coredata :: DataFrame
    end

`::TS` - A type to hold ordered data with an index.

A TS object is essentially a `DataFrame` with a specific column marked
as an index. The input `DataFrame` is sorted during construction and
is stored under the property `coredata`. The index is stored in the
`Index` column of `coredata`.

Permitted data inputs to the constructors are `DataFrame`, `Vector`,
and 2-dimensional `Array`. If an index is already not present in the
constructor then a sequential integer index is created
automatically.

`TS(coredata::DataFrame)`: Here, the constructor looks for a column
named `Index` in `coredata` as the index column, if this is not found
then the first column of `coredata` is made the index by default. If
`coredata` only has a single column then a new sequential index is
generated.

Since `TS.coredata` is a DataFrame it can be operated upon
independently using methods provided by the DataFrames package
(ex. `transform`, `combine`, etc.).

# Constructors
```julia
TS(coredata::DataFrame, index::Union{String, Symbol, Int})
TS(coredata::DataFrame, index::AbstractVector{T}) where {T<:Union{Int, TimeType}}
TS(coredata::DataFrame)
TS(coredata::DataFrame, index::UnitRange{Int})
TS(coredata::AbstractVector{T}, index::AbstractVector{V}) where {T, V}
TS(coredata::AbstractVector{T}) where {T}
TS(coredata::AbstractArray{T,2}) where {T}
TS(coredata::AbstractArray{T,2}, index::AbstractVector{V}) where {T, V}
```

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;
julia> random(x) = rand(MersenneTwister(123), x);

julia> df = DataFrame(x1 = random(10))
10×1 DataFrame
 Row │ x1
     │ Float64
─────┼───────────
   1 │ 0.768448
   2 │ 0.940515
   3 │ 0.673959
   4 │ 0.395453
   5 │ 0.313244
   6 │ 0.662555
   7 │ 0.586022
   8 │ 0.0521332
   9 │ 0.26864
  10 │ 0.108871

julia> ts = TS(df)   # generates index
(10 x 1) TS with Int64 Index

 Index  x1
 Int64  Float64
──────────────────
     1  0.768448
     2  0.940515
     3  0.673959
     4  0.395453
     5  0.313244
     6  0.662555
     7  0.586022
     8  0.0521332
     9  0.26864
    10  0.108871

# ts.coredata is a DataFrame
julia> combine(ts.coredata, :x1 => Statistics.mean, DataFrames.nrow)
1×2 DataFrame
 Row │ x1_mean  nrow
     │ Float64  Int64
─────┼────────────────
   1 │ 0.49898    418

julia> df = DataFrame(ind = [1, 2, 3], x1 = random(3))
3×2 DataFrame
 Row │ ind    x1
     │ Int64  Float64
─────┼─────────────────
   1 │     1  0.768448
   2 │     2  0.940515
   3 │     3  0.673959

julia> ts = TS(df, 1)        # the first column is index
(3 x 1) TS with Int64 Index

 Index  x1
 Int64  Float64
─────────────────
     1  0.768448
     2  0.940515
     3  0.673959

julia> df = DataFrame(x1 = random(3), x2 = random(3), Index = [1, 2, 3]);
3×3 DataFrame
 Row │ x1        x2        Index
     │ Float64   Float64   Int64
─────┼───────────────────────────
   1 │ 0.768448  0.768448      1
   2 │ 0.940515  0.940515      2
   3 │ 0.673959  0.673959      3

julia> ts = TS(df)   # uses existing `Index` column
(3 x 2) TS with Int64 Index

 Index  x1        x2
 Int64  Float64   Float64
───────────────────────────
     1  0.768448  0.768448
     2  0.940515  0.940515
     3  0.673959  0.673959

julia> dates = collect(Date(2017,1,1):Day(1):Date(2017,1,10));

julia> df = DataFrame(dates = dates, x1 = random(10))
10×2 DataFrame
 Row │ dates       x1
     │ Date        Float64
─────┼───────────────────────
   1 │ 2017-01-01  0.768448
   2 │ 2017-01-02  0.940515
   3 │ 2017-01-03  0.673959
   4 │ 2017-01-04  0.395453
   5 │ 2017-01-05  0.313244
   6 │ 2017-01-06  0.662555
   7 │ 2017-01-07  0.586022
   8 │ 2017-01-08  0.0521332
   9 │ 2017-01-09  0.26864
  10 │ 2017-01-10  0.108871

julia> ts = TS(df, :dates)
(10 x 1) TS with Date Index

 Index       x1
 Date        Float64
───────────────────────
 2017-01-01  0.768448
 2017-01-02  0.940515
 2017-01-03  0.673959
 2017-01-04  0.395453
 2017-01-05  0.313244
 2017-01-06  0.662555
 2017-01-07  0.586022
 2017-01-08  0.0521332
 2017-01-09  0.26864
 2017-01-10  0.108871

julia> ts = TS(DataFrame(x1=random(10)), dates);

julia> ts = TS(random(10))
(10 x 1) TS with Int64 Index

 Index  x1
 Int64  Float64
──────────────────
     1  0.768448
     2  0.940515
     3  0.673959
     4  0.395453
     5  0.313244
     6  0.662555
     7  0.586022
     8  0.0521332
     9  0.26864
    10  0.108871

julia> ts = TS(random(10), dates);

julia> ts = TS([random(10) random(10)], dates) # matrix object
(10 x 2) TS with Date Index

 Index       x1         x2
 Date        Float64    Float64
──────────────────────────────────
 2017-01-01  0.768448   0.768448
 2017-01-02  0.940515   0.940515
 2017-01-03  0.673959   0.673959
 2017-01-04  0.395453   0.395453
 2017-01-05  0.313244   0.313244
 2017-01-06  0.662555   0.662555
 2017-01-07  0.586022   0.586022
 2017-01-08  0.0521332  0.0521332
 2017-01-09  0.26864    0.26864
 2017-01-10  0.108871   0.108871

```
"""
struct TS

    coredata :: DataFrame

    # From DataFrame, index number/name/symbol
    function TS(coredata::DataFrame, index::Union{String, Symbol, Int})
        if (DataFrames.ncol(coredata) == 1)
            TS(coredata, collect(Base.OneTo(DataFrames.nrow(coredata))))
        end

        sorted_cd = sort(coredata, index)
        index_vals = sorted_cd[!, index]

        cd = sorted_cd[:, Not(index)]
        insertcols!(cd, 1, :Index => index_vals, after=false, copycols=true)

        new(cd)
    end

    # From DataFrame, external index
    function TS(coredata::DataFrame, index::AbstractVector{T}) where {T<:Union{Int, TimeType}}
        sorted_index = sort(index)

        cd = copy(coredata)
        insertcols!(cd, 1, :Index => sorted_index, after=false, copycols=true)

        new(cd)
    end

end

####################################
# Constructors
####################################

function TS(coredata::DataFrame)
    if "Index" in names(coredata)
        return TS(coredata, :Index)
    elseif DataFrames.ncol(coredata) == 1
        return TS(coredata, collect(1:DataFrames.nrow(coredata)))
    else
        return TS(coredata, 1)
    end
end

# From DataFrame, index range
function TS(coredata::DataFrame, index::UnitRange{Int})
    index_vals = collect(index)
    cd = copy(coredata)
    insertcols!(cd, 1, :Index => index_vals, after=false, copycols=true)
    TS(cd, :Index)
end

# From AbstractVector
function TS(coredata::AbstractVector{T}, index::AbstractVector{V}) where {T, V}
    df = DataFrame([coredata], :auto)
    TS(df, index)
end

function TS(coredata::AbstractVector{T}) where {T}
    index_vals = collect(Base.OneTo(length(coredata)))
    TS(coredata, index_vals)
end

# From Matrix and meta
# FIXME: use Metadata.jl
function TS(coredata::AbstractArray{T,2}) where {T}
    index_vals = collect(Base.OneTo(size(coredata)[1]))
    df = DataFrame(coredata, :auto, copycols=true)
    TS(df, index_vals)
end

function TS(coredata::AbstractArray{T,2}, index::AbstractVector{V}) where {T, V}
    df = DataFrame(coredata, :auto, copycols=true)
    TS(df, index)
end

ext: jl


text:
module bench

using RustModels
import HiddenMarkovModels.rsm

include("../test/toymodels.jl")

function bench1()
	srand(425)
	models=classicaltoys()
	dirm=models[1]
	theta0=[0.6,0.8]
	coef!(dirm,theta0)
	data=rand(dirm,100,100)

	timings=zeros(4)
	t=1
	for model in models
		tic()
		@time mle(model,data)
		timings[t]=toq()
		t+=1
	end
	timings
end

end #module


text:
"""
This module provides `tick()`, `tock()`, and `tok()` functions.

- `tick()`        ` start a new timer and start counting
- `tock()`        ` stop counting, show total elapsed time in canonical form
- `tok()`         ` stop counting, show and return total elapsed time in seconds
- `peektimer()    ` continue counting, show and return elapsed seconds so far
- `laptimer()     ` continue counting, show elapsed time so far in canonical form
- `alarm(h, m, s) ` set an alarm in h/m/s from now
- `alarm(dt)      ` set an alarm for `dt`

Don't use these for timing code execution! Julia provides much better facilities for measuring performance, ranging from the `@time` and `@elapsed` macros to packages such as [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl).
"""
module TickTock

export tick, tock, tok, peektimer, laptimer, alarm

using Dates

"""
    tick()

Start a timer.

Other functions: `tock()` (stop counting and show canonical), `tok()` (stop
counting and return seconds), `peektimer()` (continue counting, return elapsed
seconds), and `laptimer()` (continue counting, show canonical)
"""
function tick()
    t0 = time_ns()
    task_local_storage(:TIMERS, (t0, get(task_local_storage(), :TIMERS, ())))
    @info " started timer at: $(now())"
end

function printcanonical(tnow)
    canondc = Dates.canonicalize(
        Dates.CompoundPeriod(Dates.Second(floor(tnow)),
        Dates.Millisecond(floor((tnow - floor(tnow)) * 1000))))
    @info "$(lpad(tnow, 20))s: $canondc"
end

function gettimers()
    result = UInt64[]
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("Use `tick()` to start a timer.")
    end
    while timers[2] != ()
        push!(result, timers[1]::UInt64)
        timers = timers[2]
    end
    push!(result, timers[1]::UInt64)
    return result
end

"""
    peektimer()

Return the elapsed seconds counted by the most recent timer, without stopping it.
"""
function peektimer()
    t1 = time_ns()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("Use `tick()` to start a timer.")
    end
    t0 = timers[1]::UInt64
    return (t1 - t0)/1e9
end

"""
    tok()

Return the elapsed seconds counted by the most recent timer, then stop counting.
"""
function tok()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("Use `tick()` to start a timer.")
    end
    t = peektimer()
    task_local_storage(:TIMERS, timers[2])
    return t
end

"""
    tock()

Print the elapsed time, in canonical form, of the most recent timer, then stop counting.
"""
function tock()
    t = tok()
    printcanonical(t)
end

function showtimes(;canonical=true)
    t1 = time_ns()
    timers = gettimers()
    timernumber = length(timers)
    for t0 in timers
        tnow = (t1 - t0)/1e9
        print(rpad(timernumber, 10))
        canonical ? printcanonical(tnow) : @info(tnow)
        timernumber -= 1
    end
end

"""
laptimer()

Print the elapsed time, in canonical form, of the most recent timer, and continue counting.
"""
laptimer() = showtimes(canonical=true)

"""
    alarm(hours, minutes, seconds;
        action = () -> @info("TickTock: time's up"))

Set an alarm, with the option of providing a anonymous function that executes
when alarm fires.

```
using Dates

@async alarm(0, 5, 0, action = ()-> println("TickTock.jl: 5 minutes is up!"))

@async alarm(0, 0, 5, action = () -> run(`say "your alarm is ringing, sir"`), alarmname="tea's up") # uses macOS speech
```
"""
function alarm(hours, minutes, seconds;
        action=() -> @info("TickTock: alarm"),
        alarmname="TickTock alarm")
    tick()
    while true
        sleep(5)
        if peektimer() > hours * 60 * 60 + minutes * 60 + seconds
            @info alarmname
            action()
            tock()
            break
        end
    end
end

"""
    alarm(dt::DateTime;
        action = () -> @info("TickTock: time's up"))

Run an alarm that fires at time `dt`, with the option of providing a function
that executes when alarm fires.

# Examples

```
using Dates

dt = now() + Dates.Minute(1)

@async alarm(dt, action = () -> println("TickTock.jl: Ready!"))
```

```
@async alarm(now() + Dates.Minute(10) + Dates.Second(0), action = () -> run(`say "wake up, 10 minutes is up"`)) # macOS speech command
```

TODO alarms don't appear in timer lists...
"""
function alarm(dt::DateTime;
        action = () -> @info("TickTock: alarm"),
        alarmname = "TickTock alarm")
    p = Dates.Period(dt - now())
    secs = round(p, Dates.Second).value
    m, s = divrem(secs, 60)
    h, m = divrem(m, 60)
    @info "TickTock: \"$(alarmname)\" alarm for $h hours, $m minutes, $s seconds"
    alarm(h, m, s, action=action, alarmname=alarmname)
end

"""
@async alarm(now() + Dates.Second(5), action = () ->  TickTock.alarmnotify("time's up"))

TODO this is macOS only...
"""
function alarmnotify(subtitle="TickTock alarm")
    !Sys.isapple() && exit()
    command = """
    display notification with title "TickTock.jl" subtitle \"$(subtitle)\" sound name "frog"
    """
    chomp(read(`osascript -e $command`, String))
end

end # module


text:
using KitBase, Plots
using KitBase.JLD2
using Base.Threads: @threads
using KitBase.ProgressMeter: @showprogress
cd(@__DIR__)

function up!(ks, ctr, face, dt, p)
    kn_bzm, nm, phi, psi, phipsi = p

    res = zeros(5)
    avg = zeros(5)

    @inbounds @threads for i = 1:ks.ps.nx
        KitBase.step!(
            face[i].fw,
            face[i].ff,
            ctr[i].w,
            ctr[i].prim,
            ctr[i].f,
            face[i+1].fw,
            face[i+1].ff,
            ks.gas.γ,
            kn_bzm,
            nm,
            phi,
            psi,
            phipsi,
            ks.ps.dx[i],
            dt,
            res,
            avg,
            :fsm,
        )

        #=KitBase.step!(
            face[i].fw,
            face[i].ff,
            ctr[i].w,
            ctr[i].prim,
            ctr[i].f,
            face[i+1].fw,
            face[i+1].ff,
            ks.vs.u,
            ks.vs.v,
            ks.vs.w,
            ks.vs.weights,
            ks.gas.γ,
            ks.gas.μᵣ,
            ks.gas.ω,
            ks.gas.Pr,
            ks.ps.dx[i],
            dt,
            res,
            avg,
            :bgk,
        )=#
    end

    return nothing
end

begin
    set = Setup(case = "layer", space = "1d1f3v", maxTime = 0.2, boundary = ["fix", "fix"], cfl = 0.5)
    ps = PSpace1D(-0.5, 0.5, 500, 1)
    vs = VSpace3D(-6.0, 6.0, 28, -6.0, 6.0, 64, -6.0, 6.0, 28)
    gas = Gas(Kn = 5e-3, K = 0.0)
    fw = function(x)
        prim = zeros(5)
        if x <= 0
            prim .= [1.0, 0.0, 1.0, 0.0, 1.0]
        else
            prim .= [1.0, 0.0, -1.0, 0.0, 2.0]
        end

        return prim_conserve(prim, ks.gas.γ)
    end
    ib = IB1F(fw, vs, gas)
    ks = SolverSet(set, ps, vs, gas, ib)
    ctr, face = init_fvm(ks, ks.ps)
end

fsm = fsm_kernel(ks.vs, ks.gas.μᵣ)
τ0 = vhs_collision_time(ctr[1].prim, ks.gas.μᵣ, ks.gas.ω)
tmax = 10τ0#50τ0
t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(tmax ÷ dt)
res = zero(ctr[1].w)

@time @showprogress for iter = 1:nt
    #reconstruct!(ks, ctr)

    @inbounds @threads for i = 1:ks.ps.nx+1
        flux_kfvs!(
            face[i].fw,
            face[i].ff,
            ctr[i-1].f, #.+ 0.5 .* ctr[i-1].sf .* ks.ps.dx[i-1],
            ctr[i].f, #.- 0.5 .* ctr[i].sf .* ks.ps.dx[i],
            ks.vs.u,
            ks.vs.v,
            ks.vs.w,
            ks.vs.weights,
            dt,
            1.0,
            ctr[i-1].sf,
            ctr[i].sf,
        )
    end

    up!(ks, ctr, face, dt, fsm)

    #=global t += dt
    if abs(t - τ0) < dt
        @save "sol_t.jld2" ctr face
    elseif abs(t - 10 * τ0) < dt
        @save "sol_10t.jld2" ctr face
    end=#
end
#@save "sol_50t.jld2" ctr face

# field
sol = zeros(ks.ps.nx, 5)
for i in axes(sol, 1)
    sol[i, :] .= ctr[i].prim
    sol[i, end] = 1 / sol[i, end]
end

plot(ks.ps.x[1:ks.ps.nx], sol[:, 1])

# distribution function
fc = (ctr[end÷2].f + ctr[end÷2+1].f) ./ 2
hc = reduce_distribution(fc, vs.weights[:, 1, :], 2)
plot(ks.vs.v[1, :, 1], hc)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using JuMP, GLPK, Test

function optimal_price_probabilties(
    price, demand, capacity; verbose = true)
    revenue = price .* demand
    model = Model(with_optimizer(GLPK.Optimizer))
    @variable(model, x[1:length(price)], lower_bound=0)
    # Objective: maximize profit
    @objective(model, Max, revenue' * x)
    # Constraint: can carry all
    @constraint(model, sum(x) <= 1)
    #=
    for i in 1:length(prices)
        @constraint(model, x[i] <= demand[i])
    end
    =#
    # Solve problem using MIP solver
    JuMP.optimize!(model)
    if verbose
        println("Objective is: ", JuMP.objective_value(model))
        println("Solution is:")
        for i in 1:length(price)
            print("x[$i] = ", JuMP.value(x[i]))
            println(", revenue = ", price[i] * JuMP.value(x[i]))
        end
    end
    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    JuMP.value.(x)
end

price = [1.99, 2.49, 2.99, 3.49, 3.99, 4.49]
demand = 50 .- 7 * price
capacity = 60

prob = optimal_price_probabilties(price, demand, capacity)
prob |> display
println()

demand2 = [33.74, 31.29, 25.81, 34.68, 31.84, 15.64]
prob2 = optimal_price_probabilties(price, demand2, capacity)
prob2 |> display


text:
using Compat.Test, MathOptInterface
using LinQuadOptInterface

const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const LQOI = LinQuadOptInterface

@testset "LinQuadOptInterface" begin
    solver = LQOI.MockLinQuadOptimizer()
    # TODO(@joaquim): test with solve=true
    config = MOIT.TestConfig(solve=false)

    @testset "Unit Tests" begin
        MOIT.basic_constraint_tests(solver, config)
        MOIT.unittest(solver, config, [
            "solve_affine_interval",
            "solve_qp_edge_cases",
            "solve_qcp_edge_cases",
            "solve_affine_deletion_edge_cases"
        ])
    end

    @testset "Linear tests" begin
        MOIT.contlineartest(solver, config)
    end

    @testset "Quadratic tests" begin
        MOIT.contquadratictest(solver, config)
    end

    @testset "Linear Conic tests" begin
        MOIT.lintest(solver, config)
    end

    @testset "Integer Linear tests" begin
        MOIT.intlineartest(solver, config, ["int2"])
    end

    @testset "ModelLike tests" begin
        @testset "nametest" begin
            MOIT.nametest(LQOI.MockLinQuadOptimizer())
        end
        @testset "validtest" begin
            MOIT.validtest(LQOI.MockLinQuadOptimizer())
        end
        @testset "emptytest" begin
            MOIT.emptytest(LQOI.MockLinQuadOptimizer())
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(LQOI.MockLinQuadOptimizer())
        end
        @testset "copytest" begin
            MOIT.copytest(LQOI.MockLinQuadOptimizer(),
                          LQOI.MockLinQuadOptimizer())
        end
    end
end


text:
# Simple model objects
parameters_simple  = @with_kw (μ = 0.0048,
                                υ = 0.02,
                                θ = 2.1,
                                r = t -> 0.05,
                                ζ = 14.5,
                                ξ = 1.0, π = (t, x) -> 1,
                                x = t -> ζ)

settings_simple = @with_kw (z_ex = unique([range(0., 0.1, length = 20)' range(0.1, 1., length = 20)' range(1., 5., length = 20)']),
                            T = 100.0, iterations = 1000,
                            ts = range(0., 100., length = length(z_ex)),
                            stationary_x0 = default_simple_stationary_x0,
                            ode_solve_algorithm = IDA())

default_simple_stationary_x0(parameters, settings) = [parameters.r(settings.T)/2] # or something like that


text:
module LocalData

using Dates
using CSV, JSON
using DataFrames

basepath = normpath(@__DIR__, "../../")

export alldata_nj, load_nj

########################
###### New Jersey ######
########################

function load_nj(start_date::Date, end_date::Date; data_type::Symbol, los::Int=-1, pct_beds_available::Real=1.0)
	data_raw = _load_nj()
	data = filter_nj(data_raw, start_date, end_date)
	if data_type == :active_patients
		return active_nj(data)
	elseif data_type == :admitted_patients
		return admitted_nj(data)
	elseif data_type == :initial_patients
		return initial_nj(start_date, end_date)
	elseif data_type == :discharged_patients
		_initial = initial_nj(start_date, end_date)
		return discharged_nj(_initial, start_date, end_date, los=los)
	elseif data_type == :beds
		return beds_nj(data, pct_beds_available=pct_beds_available)
	elseif data_type == :hospital_names
		return hospitals_nj(data)
	elseif data_type == :hospital_locations
		return locations_nj(data)
	elseif data_type == :all
		return alldata_nj(start_date, end_date, los=los, pct_beds_available=pct_beds_available)
	else
		error("Invalid data_type for load_nj")
	end
end

function _load_nj()
	fn = joinpath(basepath, "data/hospitalizations/newjersey/nj_hosp_2020_06_19.csv")
	data_raw = DataFrame(CSV.File(fn))
	return data_raw
end

function filter_nj(data_raw, start_date::Date, end_date::Date)
	data = filter(row -> start_date <= row.date <= end_date, data_raw)

	_locations = sort(unique(data.hosp_name))
	_N = length(_locations)
	_T = (end_date - start_date).value + 1

	# fill missing days
	for l in _locations
		rows = data[data.hosp_name .== l,:]
		if size(rows,1) != _T
			missing_dates = setdiff(start_date:Day(1):end_date, rows.date)
			for d in missing_dates
				push!(data, Dict(:hosp_id => rows.hosp_id[1], :hosp_name => l, :date => d), cols=:subset)
			end
		end
	end

	# filter out hospitals that never have data
	bed_cols = [
		:bed_avail_criticalcare, :bed_avail_intensivecare, :bed_avail_surgical, :bed_avail_other,
		:bed_inuse_criticalcare, :bed_inuse_intensivecare, :bed_inuse_surgical, :bed_inuse_other,
	]
	data.total_beds = sum(hcat([replace(data[:,b], missing => 0) for b in bed_cols]...), dims=2)[:]
	_beds = maximum(vcat([d[:,:total_beds]' for d in groupby(data, :hosp_name, sort=true)]...), dims=2)[:]
	filter_mask = _beds .!= 0
	locations   = _locations[filter_mask]
	filter!(row -> row.hosp_name in locations, data)

	sort!(data, [:hosp_name, :date])
	return data
end

function admitted_nj(data)
	data_by_hosp = groupby(data, :hosp_name)
	_admitted_patients = vcat([d[:,:admitted]' for d in data_by_hosp]...)
	_admitted_patients = Int.(replace(_admitted_patients, missing => 0))
	return _admitted_patients
end

function active_nj(data)
	data_by_hosp = groupby(data, :hosp_name)
	_active_confirmed = replace(vcat([d[:,:active_confirmed]' for d in data_by_hosp]...), missing => 0)
	_active_suspected = replace(vcat([d[:,:active_suspected]' for d in data_by_hosp]...), missing => 0)
	_active = _active_confirmed + _active_suspected
	return _active
end

function initial_nj(start_date::Date, end_date::Date)
	data_raw = _load_nj()
	data = filter_nj(data_raw, start_date-Day(1), end_date)

	day0_df = filter(row -> row.date == start_date-Day(1), data)
	sort!(day0_df, :hosp_name)

	_initial_confirmed = replace(day0_df.active_confirmed, missing => 0)
	_initial_suspected = replace(day0_df.active_suspected, missing => 0)
	_initial = _initial_confirmed + _initial_suspected

	return _initial
end

function discharged_nj(initial::Array{<:Real,1}, start_date::Date, end_date::Date; los::Int)
	N = length(initial)
	T = (end_date - start_date).value
	_discharged_amt = max.(1, round.(Int, initial / los))
	_discharged = zeros(Int, N, T)
	for i in 1:_N
		for t in 1:_T
			_discharged[i,t] = min(_discharged_amt[i], initial[i] - sum(_discharged[i,:]))
		end
	end
	return _discharged
end

function beds_nj(data; pct_beds_available::Real=1.0)
	bed_cols = [
		:bed_avail_criticalcare, :bed_avail_intensivecare, :bed_avail_surgical, :bed_avail_other,
		:bed_inuse_criticalcare, :bed_inuse_intensivecare, :bed_inuse_surgical, :bed_inuse_other,
	]
	data.total_beds = sum(hcat([replace(data[:,b], missing => 0) for b in bed_cols]...), dims=2)[:]

	data_by_hosp = groupby(data, :hosp_name)
	_beds = maximum(vcat([d[:,:total_beds]' for d in data_by_hosp]...), dims=2)[:]
	beds = _beds .* pct_beds_available

	return beds
end

function hospitals_nj(data)
	return sort(unique(data.hosp_name))
end

function locations_nj(data)
	fn = joinpath(basepath, "data/hospitalizations/newjersey/nj_hosp_locs.csv")
	open(fn, "w") do f
		locs = JSON.parsefile(f)
	end

	hosp_names = sort(unique(data.hosp_name))
	hosp_locs = [locs[n] for n in hosp_names]

	return hosp_locs
end

function alldata_nj(start_date::Date, end_date::Date; pct_beds_available::Real=1.0, los::Int)
	@assert Date(2020, 4, 5) <= start_date <= end_date <= Date(2020, 6, 19)
	@assert los > 0

	fn = joinpath(basepath, "data/hospitalizations/newjersey/nj_hosp_2020_06_19.csv")
	data_raw = DataFrame(CSV.File(fn))
	data = filter(row -> start_date <= row.date <= end_date, data_raw)

	_locations = sort(unique(data.hosp_name))

	_N = length(_locations)
	_T = (end_date - start_date).value + 1

	for l in _locations
		rows = data[data.hosp_name .== l,:]
		if size(rows,1) != _T
			missing_dates = setdiff(start_date:Day(1):end_date, rows.date)
			for d in missing_dates
				push!(data, Dict(:hosp_id => rows.hosp_id[1], :hosp_name => l, :date => d), cols=:subset)
			end
		end
	end
	sort!(data, [:date, :hosp_name])

	data_by_hosp = groupby(data, :hosp_name)
	_admitted_patients = vcat([d[:,:admitted]' for d in data_by_hosp]...)
	_admitted_patients = Int.(replace(_admitted_patients, missing => 0))

	_active_confirmed = replace(vcat([d[:,:active_confirmed]' for d in data_by_hosp]...), missing => 0)
	_active_suspected = replace(vcat([d[:,:active_suspected]' for d in data_by_hosp]...), missing => 0)
	_active = _active_confirmed + _active_suspected

	day0_df = filter(row -> row.date == start_date-Day(1), data_raw)
	for l in _locations
		if l ∉ day0_df.hosp_name
			push!(day0_df, Dict(:hosp_name => l, :date => start_date-Day(1)), cols=:subset)
		end
	end
	sort!(day0_df, :hosp_name)
	_initial_confirmed = replace(day0_df.active_confirmed, missing => 0)
	_initial_suspected = replace(day0_df.active_suspected, missing => 0)
	_initial = _initial_confirmed + _initial_suspected

	_discharged_amt = max.(1, round.(Int, _initial / los))
	_discharged = zeros(Int, _N, _T)
	for i in 1:_N
		for t in 1:_T
			_discharged[i,t] = min(_discharged_amt[i], _initial[i] - sum(_discharged[i,:]))
		end
	end

	bed_cols = [
		:bed_avail_criticalcare, :bed_avail_intensivecare, :bed_avail_surgical, :bed_avail_other,
		:bed_inuse_criticalcare, :bed_inuse_intensivecare, :bed_inuse_surgical, :bed_inuse_other,
	]
	data.total_beds = sum(hcat([replace(data[:,b], missing => 0) for b in bed_cols]...), dims=2)[:]

	data_by_hosp = groupby(data, :hosp_name)
	_beds = maximum(vcat([d[:,:total_beds]' for d in data_by_hosp]...), dims=2)[:]

	filter_mask         = _beds .!= 0
	locations           = _locations[filter_mask]

	admitted_patients   = Float32.(_admitted_patients[filter_mask,:])
	active_patients     = Float32.(_active[filter_mask,:])
	initial_patients    = Float32.(_initial[filter_mask])
	discharged_patients = Float32.(_discharged[filter_mask,:])
	beds                = Float32.(_beds[filter_mask]) * Float32(pct_beds_available)

	return (
		admitted_patients = admitted_patients,
		active_patients = active_patients,
		initial_patients = initial_patients,
		discharged_patients = discharged_patients,
		beds = beds,
		hospitals = locations,
	)
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
using BlockArrays, Test

@testset "block product" begin
    @testset "Khatri Rao Product (size)" begin
        #A size
        m = 5
        n = 6

        #A block size
        mi = [3, 2]
        ni = [4, 1, 1]
        #B size
        p = 8
        q = 7

        #B block size
        pi = [5, 3]
        qi = [3, 2, 2]

        @assert sum(mi) == m
        @assert sum(ni) == n
        @assert sum(qi) == q
        @assert sum(pi) == p

        A = BlockArray(ones(m, n), mi, ni)
        B = BlockArray(ones(p, q), pi, qi)

        AB = khatri_rao(A, B)

        @test blocksize(AB) == blocksize(A)
        @test blocksize(AB) == blocksize(B)

        #Test: Size of resulting blocks
        for i in blockaxes(AB,1)
            for j in blockaxes(AB,2)
                @test size(AB[i, j]) == (mi[Int(i)]*pi[Int(i)], ni[Int(j)]*qi[Int(j)])
            end
        end
    end

    @testset "Khatri Rao Product (constant blocks)" begin
        #A size
        m = 5
        n = 6

        #A block size
        mi = [3, 2]
        ni = [4, 1, 1]
        #B size
        p = 8
        q = 7

        #B block size
        pi = [5, 3]
        qi = [3, 2, 2]

        @assert sum(mi) == m
        @assert sum(ni) == n
        @assert sum(qi) == q
        @assert sum(pi) == p

        A = BlockArray(ones(m, n), mi, ni)
        B = BlockArray(ones(p, q), pi, qi)

        #Test: Resulting values for a matrix of constant sub-blocks
        for i in blockaxes(A,1)
            for j in blockaxes(A,2)
                A[i, j] .*= Int(i)+Int(j)
            end
        end

        for i in blockaxes(B,1)
            for j in blockaxes(B,2)
                B[i, j] .*= Int(i)+Int(j)+10
            end
        end

        AB = khatri_rao(A, B)

        for i in blockaxes(AB,1)
            for j in blockaxes(AB,2)
                @test AB[i, j] ≈ (Int(i)+Int(j))*(Int(i)+Int(j)+10)*ones(mi[Int(i)]*pi[Int(i)], ni[Int(j)]*qi[Int(j)])
            end
        end

    end

    @testset "Khatri Rao Product (wrong blocksize)" begin

        A = BlockArray(ones(1, 2), [1], [1,1])
        B = BlockArray(ones(2, 2), [1,1], [1,1])

        @test_throws AssertionError khatri_rao(A, B)
    end

    @testset "Khatri Rao Product (Matrix)" begin

        A = ones(1, 2)
        B = ones(2, 2)

        @test khatri_rao(A, B) ≈ kron(A, B)
    end

    @testset "blockkron" begin
        a = [1,2]
        b = [3,4,5]
        k̃ = BlockKron(a,b) # lazy version of blockkron
        k = blockkron(a,b)
        @test k == k̃ == kron(a,b)
        @test blocksize(k) == blocksize(k̃) == size(a)

        @test k[Block(1)] == k̃[Block(1)] == a[1]*b
        @test k[Block(2)] == k̃[Block(2)] == a[2]*b
        c = 6:8
        k̄ = BlockKron(a,b,c)
        @test k̄ == blockkron(a,b,c) == kron(a,b,c)
        @test k̄[Block(1)][Block(1)] == a[1]*b[1]*c
        @test k̄[Block(1)][Block(2)] == a[1]*b[2]*c
        @test k̄[Block(2)][Block(3)] == a[2]*b[3]*c

        A = randn(2,3)
        B = randn(3,4)
        K = blockkron(A,B)
        K̃ = BlockKron(A,B)
        @test K == K̃ == kron(A,B)
        @test blocksize(K) == blocksize(K̃) == size(A)
        @test K[Block(1,1)] == K̃[Block(1),Block(1)] == A[1,1]*B
        @test K[Block(2,3)] == K̃[Block(2),Block(3)] == A[2,3]*B
        C = randn(2,5)
        K̄ = BlockKron(A,B,C)
        @test K̄ == blockkron(A,B,C) == kron(A,B,C)
        @test K̄[Block(1,1)][Block(1,1)] ≈ A[1,1]*B[1,1]*C
        @test K̄[Block(2,3)][Block(3,4)] ≈ A[2,3]*B[3,4]*C

        @test blockkron(a,B) == kron(a,B)
        @test blockkron(A,b) == kron(A,b)
        @test blockkron(A,b,c) == kron(A,b,c)
        @test blockkron(A,b,C) == kron(A,b,C)

        @test_throws MethodError BlockKron()
        @test_throws MethodError BlockKron(a)
    end
end


text:
###This file intends to demonstrate the rendering and basically follows basic.jl
import Cxx
import DifferentialDynamicProgramming
import GLFW
import MuJoCo

using Cxx
using MuJoCo
using MuJoCo.time

include("../mjderivs.jl")

function distance(pnt1, pnt2, norm=2)
    sum = 0
    for i=1:length(pnt1)
        sum += (pnt1[i] - pnt2[i]) ^ norm
    end
    return sum ^ (1 / norm)
end

mj_activate()

reaching = true
N = 180

xmlfile = nothing
if reaching
    xmlfile = joinpath(@__DIR__, "../../reaching.xml")
else
    xmlfile = joinpath(@__DIR__, "../../cartpoleBalance.xml")
end

model = MJModel(xmlfile)
target = "2-1"

sim = MJSimEnv{SimState, AltControls}(model, s -> (0), s -> (false), SimState(), AltControls())
ResetEnv(sim)
for i in 1:50
    mj_step(sim.model, sim.data)
end
x0 = GetSimState(sim)
u0 = randn(sim.controls.nu, N)

x = u = nothing

if ARGS[1] == "--load"
    print("Loading the previous trajectory...")

    #load from file
    x = readcsv(open("x.csv", "r"))
    u = readcsv(open("u.csv", "r"))

    if length(ARGS) > 1
        x = readcsv(open("x$(ARGS[2]).csv", "r"))
        u = readcsv(open("u$(ARGS[2]).csv", "r"))
    end

    N = size(u, 2)
end

window = GLFW.CreateWindow(640, 480, "demo")
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)

renderer = MJRenderer(sim.data)

width = Int32(0)
height = Int32(0)
step = 1
notStepped = true
elapsedTime = 0
ResetEnv(sim)
if (!reaching)
    for i in 1:50
        mj_step(sim.model, sim.data)
    end
end
tic()
while !GLFW.WindowShouldClose(window)

    #simulate
    while (elapsedTime < 1.0 / 60)
        elapsedTime += toq()
        tic()
        if (step < N && notStepped)
            #for i in 1:3
                mj_step1(sim.model, sim.data)
                #u[:, step] = zeros(size(u[:, step]))
                #SetSimState(sim, x[:, step])
                SetSimControls(sim, u[:, step])
                mj_step2(sim.model, sim.data)
                mj_forward(sim.model, sim.data)
            #end
            step += 1
            notStepped = false
            println("Frame:", step)
        elseif (notStepped)
            #mj_step(sim.model, sim.data)
            step += 1
            notStepped = false
        end

        #Render
        width, height = GLFW.GetFramebufferSize(window);
        mjv_updateScene(renderer)
        mjr_render(renderer, width, height)

        # Swap front and back buffers
    	GLFW.SwapBuffers(window)

    	# Poll for and process events
    	GLFW.PollEvents()
    end
    notStepped = true
    elapsedTime = 0
end

GLFW.Terminate()


text:
using Flux
using Glob
using Images
using Base.Iterators: partition
using Flux: onehotbatch,onecold, crossentropy, throttle
using Flux:@treelike
using Statistics:mean

# Using only modern tables
@show("Load image paths")
img_tr = glob("*_t1*.jpg","icdar2019/ICDAR2019_cTDaR/training/TRACKA/ground_truth/")
mask_tr = glob("*_t1*_mask.png","icdar2019/ICDAR2019_cTDaR/training/TRACKA/ground_truth/")

img_ts = glob("*_t1*.jpg","icdar2019/ICDAR2019_cTDaR/test/TRACKA/")
mask_ts = glob("*_t1*_mask.png","icdar2019/ICDAR2019_cTDaR/test/TRACKA/")

function preprocess(imlist)
    imageArray = []
    for image in imlist
        im = load(image)
        im = Gray.(im)  #Convert to grayscale
        im = imresize(im,(256,256))  # Resize to 2^n
        im = convert(Array{Float32},im) # Convert to Float32
        im = reshape(im,256,256,1,:)  # Reshape array
        append!(imageArray,im)
    end
    return Float32.(reshape(imageArray,256,256,1,:))  # Reshape whole array
end

@show("Converting images to arrays")
trainX = preprocess(img_tr)
trainY = preprocess(mask_tr)
testX = preprocess(img_ts)
testY = preprocess(mask_ts)

@show("Create minibatches")
function make_partition_index(X,batch_size)
    idx = partition(1:Int32(length(X)/(256*256*1)),batch_size)
    indices = [(minimum(i),maximum(i)) for i in idx]
    return indices
end

function make_minibatch(X,Y,batch_size)
    typeof(X)
    indices = [i for i in make_partition_index(X,batch_size)]
    minibatch_X = [X[:,:,:,indices[i][1]:indices[i][2]] for i in 1:length(indices)]
    print(typeof(minibatch_X))
    minibatch_Y = [Y[:,:,:,indices[i][1]:indices[i][2]] for i in 1:length(indices)]
    dataset = [(minibatch_X[i],minibatch_Y[i]) for i in 1:length(indices)]
    return dataset
end

train_set = make_minibatch(trainX,trainY,8)
test_set = make_minibatch(testX,testY,8)

@show("Build U-net")

# Neural Network blocks
# Down Block
struct DownBlock
    filter_in
    filter_out
    kernel_size
    padding
    stride
end

function (m::DownBlock)(x)
    c = Conv(m.kernel_size,m.filter_in=>m.filter_out, relu, pad=m.padding, stride=m.stride)(x)
    c = Conv(m.kernel_size,m.filter_out=>m.filter_out, relu, pad=m.padding, stride=m.stride)(c)
    p = MaxPool((2,2))(c)
    return c,p
end
@treelike DownBlock

# Up Block
struct UpBlock
    skip
    filter_in
    filter_out
    kernel_size
    padding
    stride
end

function (m::UpBlock)(x)
    up = ConvTranspose((2,2),m.filter_in=>m.filter_out, relu, stride=(2,2))(x)
    concat = up .+ m.skip
    c = Conv(m.kernel_size,m.filter_out=>m.filter_out,pad=m.padding, stride=m.stride)(concat)
    c = Conv(m.kernel_size,m.filter_out=>m.filter_out,pad=m.padding, stride=m.stride)(c)
    return c
end
@treelike UpBlock

# Bottleneck layer
struct Bottleneck
    filter_in
    filter_out
    kernel_size
    padding
    stride
end

function (m::Bottleneck)(x)
    c = Conv(m.kernel_size,m.filter_in=>m.filter_out, relu, pad=m.padding, stride=m.stride)(x)
    c = Conv(m.kernel_size,m.filter_out=>m.filter_out, relu, pad=m.padding, stride=m.stride)(c)
    return c
end
@treelike Bottleneck

# Unet
struct UNet
    filters
end

function (m::UNet)(x)
    c1,p1 = DownBlock(1,m.filters[1],(3,3),(1,1),(1,1))(x)
    c2,p2 = DownBlock(m.filters[1],m.filters[2],(3,3),(1,1),(1,1))(p1)
    c3,p3 = DownBlock(m.filters[2],m.filters[3],(3,3),(1,1),(1,1))(p2)
    c4,p4 = DownBlock(m.filters[3],m.filters[4],(3,3),(1,1),(1,1))(p3)
    bn = Bottleneck(m.filters[4],m.filters[5],(3,3),(1,1),(1,1))(p4)
    u1 = UpBlock(c4,m.filters[5],m.filters[4],(3,3),(1,1),(1,1))(bn)
    u2 = UpBlock(c3,m.filters[4],m.filters[3],(3,3),(1,1),(1,1))(u1)
    u3 = UpBlock(c2,m.filters[3],m.filters[2],(3,3),(1,1),(1,1))(u2)
    u4 = UpBlock(c1,m.filters[2],m.filters[1],(3,3),(1,1),(1,1))(u3)
    output = Conv((3,3),m.filters[1]=>1,sigmoid)(u4)
    return output
end
@treelike UNet
# Build model
model() = Chain(UNet([64,128,256,512,1024]))

m = model()


text:
macro juliadots(expr::Expr)
	expr = :(print_with_color(:red, " ●");
				print_with_color(:green, "●");
				print_with_color(:blue, "● ");
				print_with_color(:bold, :($($expr))))
	return expr
end

macro boink(expr::Expr)
	expr = :(:($($expr)) + EPSILON)
	return expr
end

type Foo
	x::Int
end

macro bumper(expr::Expr)
	if (expr.head == :.) || (expr.head == :ref)
		expr = :(:($($expr)) += EPSILON)
	elseif expr.head == :(=)
		expr = :(:($($(expr.args[1]))) = EPSILON + :($($(expr.args[2]))))
	end
	return expr
end

macro buf(args...)
	if isa(args[1], Expr)
		expr = args[1]
	else
		b = args[1]
		expr = args[2]
	end

	model = expr.args[1]

	if expr.args[2] == :(:Npsums)
		quoteblock =
		quote
		$(esc(model)).Npsumsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).Npsums[$(esc(b))])
		end

	elseif expr.args[2] == :(:Jpsums)
		quoteblock =
		quote
		$(esc(model)).Jpsumsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).Jpsums[$(esc(b))])
		end

	elseif expr.args[2] == :(:Rpsums)
		quoteblock =
		quote
		$(esc(model)).Rpsumsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).Rpsums[$(esc(b))])
		end

	elseif expr.args[2] == :(:Ypsums)
		quoteblock =
		quote
		$(esc(model)).Ypsumsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).Ypsums[$(esc(b))])
		end

	elseif expr.args[2] == :(:terms)
		quoteblock =
		quote
		$(esc(model)).termsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).terms[$(esc(b))])
		end

	elseif expr.args[2] == :(:counts)
		quoteblock =
		quote
		$(esc(model)).countsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).counts[$(esc(b))])
		end

	elseif expr.args[2] == :(:words)
		quoteblock =
		quote
		$(esc(model)).wordsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).words[$(esc(b))])
		end

	elseif expr.args[2] == :(:readers)
		quoteblock =
		quote
		$(esc(model)).readersbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).readers[$(esc(b))])
		end

	elseif expr.args[2] == :(:ratings)
		quoteblock =
		quote
		$(esc(model)).ratingsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).ratings[$(esc(b))])
		end

	elseif expr.args[2] == :(:views)
		quoteblock =
		quote
		$(esc(model)).viewsbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).views[$(esc(b))])
		end

	elseif expr.args[2] == :(:alpha)
		quoteblock =
		quote
		$(esc(model)).alphabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).alpha)
		end

	elseif expr.args[2] == :(:beta)
		quoteblock =
		quote
		$(esc(model)).betabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).beta)
		end

	elseif expr.args[2] == :(:newbeta)
		quoteblock =
		quote
		$(esc(model)).newbetabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=zeros(Float32, $(esc(model)).K, $(esc(model)).V))
		end

	elseif expr.args[2] == :(:gamma)
		quoteblock =
		quote
		$(esc(model)).gammabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=hcat($(esc(model)).gamma..., zeros(Float32, $(esc(model)).K, 64 - $(esc(model)).M % 64)))
		end

	elseif expr.args[2] == :(:phi)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).phibuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=ones(Float32, $(esc(model)).K, sum($(esc(model)).N[batch]) + 64 - sum($(esc(model)).N[batch]) % 64) / $(esc(model)).K)
		end

	elseif expr.args[2] == :(:Elogtheta)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).Elogthetabuf = cl.Buffer(Float32, $(esc(model)).context, :rw, $(esc(model)).K * (length(batch) + 64 - length(batch) % 64))
		end

	elseif expr.args[2] == :(:Elogthetasum)
		quoteblock =
		quote
		$(esc(model)).Elogthetasumbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).Elogthetasum)
		end

	elseif expr.args[2] == :(:C)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).Cbuf = cl.Buffer(Int, $(esc(model)).context, (:r, :copy), hostbuf=$(esc(model)).C[batch])
		end

	elseif expr.args[2] == :(:newtontemp)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).newtontempbuf = cl.Buffer(Float32, $(esc(model)).context, :rw, $(esc(model)).K^2 * (length(batch) + 64 - length(batch) % 64))
		end

	elseif expr.args[2] == :(:newtongrad)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).newtongradbuf = cl.Buffer(Float32, $(esc(model)).context, :rw, $(esc(model)).K * (length(batch) + 64 - length(batch) % 64))
		end

	elseif expr.args[2] == :(:newtoninvhess)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).newtoninvhessbuf = cl.Buffer(Float32, $(esc(model)).context, :rw, $(esc(model)).K^2 * (length(batch) + 64 - length(batch) % 64))
		end

	elseif expr.args[2] == :(:mu)
		quoteblock =
		quote
		$(esc(model)).mubuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).mu)
		end

	elseif expr.args[2] == :(:sigma)
		quoteblock =
		quote
		$(esc(model)).sigmabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).sigma)
		end

	elseif expr.args[2] == :(:invsigma)
		quoteblock =
		quote
		$(esc(model)).invsigmabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).invsigma)
		end

	elseif expr.args[2] == :(:lambda)
		quoteblock =
		quote
		$(esc(model)).lambdabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=hcat($(esc(model)).lambda..., zeros(Float32, $(esc(model)).K, 64 - $(esc(model)).M % 64)))
		end

	elseif expr.args[2] == :(:vsq)
		quoteblock =
		quote
		$(esc(model)).vsqbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=hcat($(esc(model)).vsq..., zeros(Float32, $(esc(model)).K, 64 - $(esc(model)).M % 64)))
		end

	elseif expr.args[2] == :(:lzeta)
		quoteblock =
		quote
		$(esc(model)).lzetabuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).lzeta)
		end

	elseif expr.args[2] == :(:alef)
		quoteblock =
		quote
		$(esc(model)).alefbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).alef)
		end

	elseif expr.args[2] == :(:newalef)
		quoteblock =
		quote
		$(esc(model)).newalefbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=fill($(esc(model)).a, $(esc(model)).K, $(esc(model)).V))
		end

	elseif expr.args[2] == :(:bet)
		quoteblock =
		quote
		$(esc(model)).betbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).bet)
		end

	elseif expr.args[2] == :(:gimel)
		quoteblock =
		quote
		$(esc(model)).gimelbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=hcat($(esc(model)).gimel..., zeros(Float32, $(esc(model)).K, 64 - $(esc(model)).M % 64)))
		end

	elseif expr.args[2] == :(:dalet)
		quoteblock =
		quote
		$(esc(model)).daletbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).dalet)
		end

	elseif expr.args[2] == :(:he)
		quoteblock =
		quote
		$(esc(model)).hebuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).he)
		end

	elseif expr.args[2] == :(:newhe)
		quoteblock =
		quote
		$(esc(model)).newhebuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=fill($(esc(model)).e, $(esc(model)).K, $(esc(model)).U))
		end

	elseif expr.args[2] == :(:vav)
		quoteblock =
		quote
		$(esc(model)).vavbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).vav)
		end

	elseif expr.args[2] == :(:zayin)
		quoteblock =
		quote
		$(esc(model)).zayinbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=hcat($(esc(model)).zayin..., zeros(Float32, $(esc(model)).K, 64 - $(esc(model)).M % 64)))
		end

	elseif expr.args[2] == :(:het)
		quoteblock =
		quote
		$(esc(model)).hetbuf = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).het)
		end

	elseif expr.args[2] == :(:xi)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		$(esc(model)).xibuf = cl.Buffer(Float32, $(esc(model)).context, :rw, 2 * $(esc(model)).K * (sum($(esc(model)).R[batch]) + 64 - sum($(esc(model)).R[batch]) % 64))
		end
	end
	return quoteblock
end

macro host(args...)
	if isa(args[1], Expr)
		expr = args[1]
	else
		b = args[1]
		expr = args[2]
	end

	model = expr.args[1]

	if expr.args[2] == :(:alphabuf)
		quoteblock =
		quote
		$(esc(model)).alpha = cl.read($(esc(model)).queue, $(esc(model)).alphabuf)
		end

	elseif expr.args[2] == :(:betabuf)
		quoteblock =
		quote
		$(esc(model)).beta = reshape(cl.read($(esc(model)).queue, $(esc(model)).betabuf), $(esc(model)).K, $(esc(model)).V)
		end

	elseif expr.args[2] == :(:gammabuf)
		quoteblock =
		quote
		hostgamma = reshape(cl.read($(esc(model)).queue, $(esc(model)).gammabuf), $(esc(model)).K, $(esc(model)).M + 64 - $(esc(model)).M % 64)
		$(esc(model)).gamma = [hostgamma[:,d] for d in 1:$(esc(model)).M]
		end

	elseif expr.args[2] == :(:phibuf)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		Npsums = $(esc(model)).Npsums[$(esc(b))]
		hostphi = reshape(cl.read($(esc(model)).queue, $(esc(model)).phibuf), $(esc(model)).K, sum($(esc(model)).N[batch]) + 64 - sum($(esc(model)).N[batch]) % 64)
		$(esc(model)).phi = [hostphi[:,Npsums[d]+1:Npsums[d+1]] for d in 1:length(batch)]
		end

	elseif expr.args[2] == :(:Elogthetabuf)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		hostElogtheta = reshape(cl.read($(esc(model)).queue, $(esc(model)).Elogthetabuf), $(esc(model)).K, length(batch) + 64 - length(batch) % 64)
		$(esc(model)).Elogtheta = [hostElogtheta[:,d] for d in 1:length(batch)]
		end

	elseif expr.args[2] == :(:Elogthetasumbuf)
		quoteblock =
		quote
		$(esc(model)).Elogthetasum = cl.read($(esc(model)).queue, $(esc(model)).Elogthetasumbuf)
		end

	elseif expr.args[2] == :(:mubuf)
		quoteblock =
		quote
		$(esc(model)).mu = cl.read($(esc(model)).queue, $(esc(model)).mubuf)
		end

	elseif expr.args[2] == :(:sigmabuf)
		quoteblock =
		quote
		$(esc(model)).sigma = reshape(cl.read($(esc(model)).queue, $(esc(model)).sigmabuf), $(esc(model)).K, $(esc(model)).K)
		end

	elseif expr.args[2] == :(:invsigmabuf)
		quoteblock =
		quote
		$(esc(model)).invsigma = reshape(cl.read($(esc(model)).queue, $(esc(model)).invsigmabuf), $(esc(model)).K, $(esc(model)).K)
		end

	elseif expr.args[2] == :(:lambdabuf)
		quoteblock =
		quote
		hostlambda = reshape(cl.read($(esc(model)).queue, $(esc(model)).lambdabuf), $(esc(model)).K, $(esc(model)).M + 64 - $(esc(model)).M % 64)
		$(esc(model)).lambda = [hostlambda[:,d] for d in 1:$(esc(model)).M]
		end

	elseif expr.args[2] == :(:vsqbuf)
		quoteblock =
		quote
		hostvsq = reshape(cl.read($(esc(model)).queue, $(esc(model)).vsqbuf), $(esc(model)).K, $(esc(model)).M + 64 - $(esc(model)).M % 64)
		$(esc(model)).vsq = [hostvsq[:,d] for d in 1:$(esc(model)).M]
		end

	elseif expr.args[2] == :(:lzetabuf)
		quoteblock =
		quote
		$(esc(model)).lzeta = cl.read($(esc(model)).queue, $(esc(model)).lzetabuf)
		end

	elseif expr.args[2] == :(:alefbuf)
		quoteblock =
		quote
		$(esc(model)).alef = reshape(cl.read($(esc(model)).queue, $(esc(model)).alefbuf), $(esc(model)).K, $(esc(model)).V)
		end

	elseif expr.args[2] == :(:betbuf)
		quoteblock =
		quote
		$(esc(model)).bet = cl.read($(esc(model)).queue, $(esc(model)).betbuf)
		end

	elseif expr.args[2] == :(:gimelbuf)
		quoteblock =
		quote
		hostgimel = reshape(cl.read($(esc(model)).queue, $(esc(model)).gimelbuf), $(esc(model)).K, $(esc(model)).M + 64 - $(esc(model)).M % 64)
		$(esc(model)).gimel = [hostgimel[:,d] for d in 1:$(esc(model)).M]
		end

	elseif expr.args[2] == :(:daletbuf)
		quoteblock =
		quote
		$(esc(model)).dalet = cl.read($(esc(model)).queue, $(esc(model)).daletbuf)
		end

	elseif expr.args[2] == :(:hebuf)
		quoteblock =
		quote
		$(esc(model)).he = reshape(cl.read($(esc(model)).queue, $(esc(model)).hebuf), $(esc(model)).K, $(esc(model)).U)
		end

	elseif expr.args[2] == :(:vavbuf)
		quoteblock =
		quote
		$(esc(model)).vav = cl.read($(esc(model)).queue, $(esc(model)).vavbuf)
		end

	elseif expr.args[2] == :(:zayinbuf)
		quoteblock =
		quote
		hostzayin = reshape(cl.read($(esc(model)).queue, $(esc(model)).zayinbuf), $(esc(model)).K, $(esc(model)).M + 64 - $(esc(model)).M % 64)
		$(esc(model)).zayin = [hostzayin[:,d] for d in 1:$(esc(model)).M]
		end

	elseif expr.args[2] == :(:hetbuf)
		quoteblock =
		quote
		$(esc(model)).het = cl.read($(esc(model)).queue, $(esc(model)).hetbuf)
		end

	elseif expr.args[2] == :(:xibuf)
		quoteblock =
		quote
		batch = $(esc(model)).batches[$(esc(b))]
		Rpsums = $(esc(model)).Rpsums[$(esc(b))]
		hostxi = reshape(cl.read($(esc(model)).queue, $(esc(model)).xibuf), 2 * $(esc(model)).K, sum($(esc(model)).R[batch]) + 64 - sum($(esc(model)).R[batch]) % 64)
		$(esc(model)).xi = [hostxi[:,Rpsums[m]+1:Rpsums[m+1]] for m in 1:length(batch)]
		end
	end
	return quoteblock
end

macro gpu(args...)
	if length(args) == 1
		@assert isa(args[1], Expr)
		batchsize = Inf
		expr = args[1]
	elseif length(args) == 2
		@assert isa(args[1], Integer)
		@assert isa(args[2], Expr)
		batchsize = args[1]
		expr = args[2]
	else
		throw(ArgumentError("Wrong number of arguments."))
	end

	@assert ispositive(batchsize)
	@assert expr.args[1] == :train! "GPU acceleration only applies to the train! function."

	quote
	local model = $(esc(expr.args[2]))
	local batchsize = Int(min($(esc(batchsize)), length(model.corp)))
	local kwargs = [(kw.args[1], kw.args[2]) for kw in $(esc(expr.args[3:end]))]

	if isa(model, LDA)
		fakecorp = Corpus(docs=[Document([1])], lex=["1"])
		gpumodel = gpuLDA(fakecorp, 1)

		gpumodel.corp = model.corp

		gpumodel.K = model.K
		gpumodel.M = model.M
		gpumodel.V = model.V
		gpumodel.N = model.N
		gpumodel.C = model.C

		gpumodel.batches = partition(1:model.M, batchsize)
		gpumodel.B = length(gpumodel.batches)

		gpumodel.topics = model.topics

		gpumodel.alpha = model.alpha
		gpumodel.beta = model.beta
		gpumodel.gamma = model.gamma
		gpumodel.phi = unshift!([ones(model.K, model.N[d]) / model.K for d in gpumodel.batches[1][2:end]], model.phi)
		gpumodel.elbo = model.elbo

		fixmodel!(gpumodel)
		train!(gpumodel; kwargs...)

		model.topics = gpumodel.topics
		model.alpha = gpumodel.alpha
		model.beta = gpumodel.beta
		model.gamma = gpumodel.gamma
		model.phi = gpumodel.phi[1]
		model.Elogtheta = gpumodel.Elogtheta[1]
		model.Elogthetasum = zeros(gpumodel.K)
		model.elbo = gpumodel.elbo

		model.beta ./= sum(model.beta, 2)
		model.phi ./= sum(model.phi, 1)
		nothing

	elseif isa(model, fLDA)
		nothing

	elseif isa(model, CTM)
		fakecorp = Corpus(docs=[Document([1])], lex=["1"])
		gpumodel = gpuCTM(fakecorp, 1)

		gpumodel.corp = model.corp

		gpumodel.K = model.K
		gpumodel.M = model.M
		gpumodel.V = model.V
		gpumodel.N = model.N
		gpumodel.C = model.C

		gpumodel.batches = partition(1:model.M, batchsize)
		gpumodel.B = length(gpumodel.batches)

		gpumodel.topics = model.topics

		gpumodel.mu = model.mu
		gpumodel.sigma = model.sigma
		gpumodel.beta = model.beta
		gpumodel.lambda = model.lambda
		gpumodel.vsq = model.vsq
		gpumodel.lzeta = fill(model.lzeta, model.M)
		gpumodel.phi = unshift!([ones(model.K, model.N[d]) / model.K for d in gpumodel.batches[1][2:end]], model.phi)
		gpumodel.elbo = model.elbo

		fixmodel!(gpumodel)
		train!(gpumodel; kwargs...)

		model.topics = gpumodel.topics
		model.mu = gpumodel.mu
		model.sigma = gpumodel.sigma
		model.invsigma = gpumodel.invsigma
		model.beta = gpumodel.beta
		model.lambda = gpumodel.lambda
		model.vsq = gpumodel.vsq
		model.lzeta = gpumodel.lzeta[1]
		model.phi = gpumodel.phi[1]
		model.elbo = gpumodel.elbo

		model.beta ./= sum(model.beta, 2)
		model.phi ./= sum(model.phi, 1)
		nothing

	elseif isa(model, fCTM)
		nothing

	elseif isa(model, DTM)
		nothing

	elseif isa(model, CTPF)
		fakecorp = Corpus(docs=[Document([1], readers=[1])], lex=["1"], users=["1"])
		gpumodel = gpuCTPF(fakecorp, 1)

		gpumodel.corp = model.corp

		gpumodel.batches = partition(1:model.M, batchsize)
		gpumodel.B = length(gpumodel.batches)

		gpumodel.topics = model.topics
		gpumodel.scores = model.scores
		gpumodel.libs = model.libs
		gpumodel.drecs = model.drecs
		model.urecs = model.urecs

		gpumodel.K = model.K
		gpumodel.M = model.M
		gpumodel.V = model.V
		gpumodel.U = model.U
		gpumodel.N = model.N
		gpumodel.C = model.C
		gpumodel.R = model.R

		gpumodel.a = model.a
		gpumodel.b = model.b
		gpumodel.c = model.c
		gpumodel.d = model.d
		gpumodel.e = model.e
		gpumodel.f = model.f
		gpumodel.g = model.g
		gpumodel.h = model.h
		gpumodel.alef = model.alef
		gpumodel.bet = model.bet
		gpumodel.gimel = model.gimel
		gpumodel.dalet = model.dalet
		gpumodel.he = model.he
		gpumodel.vav = model.vav
		gpumodel.zayin = model.zayin
		gpumodel.het = model.het
		gpumodel.phi = unshift!([ones(model.K, model.N[d]) / model.K for d in gpumodel.batches[1][2:end]], model.phi)
		gpumodel.xi = unshift!([ones(2model.K, model.R[d]) / 2model.K for d in gpumodel.batches[1][2:end]], model.xi)
		gpumodel.elbo = model.elbo

		fixmodel!(gpumodel)
		train!(gpumodel; kwargs...)

		model.topics = gpumodel.topics
		model.scores = gpumodel.scores
		model.drecs = gpumodel.drecs
		model.urecs = gpumodel.urecs

		model.alef = gpumodel.alef
		model.bet = gpumodel.bet
		model.gimel = gpumodel.gimel
		model.dalet = gpumodel.dalet
		model.he = gpumodel.he
		model.vav = gpumodel.vav
		model.zayin = gpumodel.zayin
		model.het = gpumodel.het
		model.phi = gpumodel.phi[1]
		model.xi = gpumodel.xi[1]
		model.elbo = gpumodel.elbo

		model.phi ./= sum(model.phi, 1)
		model.xi ./= sum(model.xi, 1)
		nothing

	else
		nothing
	end
	end
end


text:
module test_suite

using Test,
    BasisFunctions, BasisFunctions.Test

using DoubleFloats,
    FastTransforms,
    FFTW,
    LinearAlgebra,
    Random,
    StaticArrays,
    DomainSets,
    DomainIntegrals

import BasisFunctions.Test:
    supports_approximation,
    supports_interpolation,
    suitable_function,
    suitable_interpolation_grid

⊕ = BasisFunctions.:⊕

linspace(a,b,c) = range(a, stop=b, length=c)

my_rand = rand
Random.seed!(1234)

BF = BasisFunctions
const show_timings = false

# const LargeFloat = Double64
const LargeFloat = BigFloat

##########
# Testing
##########
ENV["JULIA_DEBUG"]="all"
include("test_scenariolist.jl")

Delimit("Utilities")
include("util/test_utilities.jl")
include("util/test_pgfplots.jl")

Delimit("Operators")
include("test_operators.jl")

Delimit("Dictionaries")
include("dictionaries/test_dictionaries.jl")
include("dictionaries/test_pwconstants.jl")

Delimit("Check evaluations, interpolations, extensions, expansions")
include("test_fourier.jl")
include("test_chebyshev.jl")
include("test_ops.jl")
include("test_gram.jl")

Delimit("Generic OPS")
include("test_generic_OPS.jl")

Delimit("Plots")
include("test_plots.jl")

println(" All tests passed!")

end # module


text:
using Plots
t = collect(-4:0.2:4)
function huber(t; λ = 1)
    flag = (abs.(t) .< λ)
    y1 = t.^2 / 2 .* flag .+ (λ * abs.(t) .- λ^2 / 2) .* (1 .- flag)
    y2 = λ * abs.(t) .- λ^2 / 2
    return y1, y2
end
y1, y2 = huber(t, λ = 2)
plot(t, y1, label = "Huber", legend = :bottomright)
plot!(t, y2, label = "Linear")


text:
# Autogenerated wrapper script for Poppler_jll for x86_64-unknown-freebsd
export libpoppler, libpoppler_cpp, libpoppler_glib, pdfattach, pdfdetach, pdffonts, pdfimages, pdfinfo, pdfseparate, pdftocairo, pdftohtml, pdftoppm, pdftops, pdftotext, pdfunite

using Cairo_jll
using Fontconfig_jll
using Glib_jll
using JpegTurbo_jll
using Libtiff_jll
using OpenJpeg_jll
using libpng_jll
JLLWrappers.@generate_wrapper_header("Poppler")
JLLWrappers.@declare_library_product(libpoppler, "libpoppler.so.113")
JLLWrappers.@declare_library_product(libpoppler_cpp, "libpoppler-cpp.so.0")
JLLWrappers.@declare_library_product(libpoppler_glib, "libpoppler-glib.so.8")
JLLWrappers.@declare_executable_product(pdfattach)
JLLWrappers.@declare_executable_product(pdfdetach)
JLLWrappers.@declare_executable_product(pdffonts)
JLLWrappers.@declare_executable_product(pdfimages)
JLLWrappers.@declare_executable_product(pdfinfo)
JLLWrappers.@declare_executable_product(pdfseparate)
JLLWrappers.@declare_executable_product(pdftocairo)
JLLWrappers.@declare_executable_product(pdftohtml)
JLLWrappers.@declare_executable_product(pdftoppm)
JLLWrappers.@declare_executable_product(pdftops)
JLLWrappers.@declare_executable_product(pdftotext)
JLLWrappers.@declare_executable_product(pdfunite)
function __init__()
    JLLWrappers.@generate_init_header(Cairo_jll, Fontconfig_jll, Glib_jll, JpegTurbo_jll, Libtiff_jll, OpenJpeg_jll, libpng_jll)
    JLLWrappers.@init_library_product(
        libpoppler,
        "lib/libpoppler.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpoppler_cpp,
        "lib/libpoppler-cpp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpoppler_glib,
        "lib/libpoppler-glib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        pdfattach,
        "bin/pdfattach",
    )

    JLLWrappers.@init_executable_product(
        pdfdetach,
        "bin/pdfdetach",
    )

    JLLWrappers.@init_executable_product(
        pdffonts,
        "bin/pdffonts",
    )

    JLLWrappers.@init_executable_product(
        pdfimages,
        "bin/pdfimages",
    )

    JLLWrappers.@init_executable_product(
        pdfinfo,
        "bin/pdfinfo",
    )

    JLLWrappers.@init_executable_product(
        pdfseparate,
        "bin/pdfseparate",
    )

    JLLWrappers.@init_executable_product(
        pdftocairo,
        "bin/pdftocairo",
    )

    JLLWrappers.@init_executable_product(
        pdftohtml,
        "bin/pdftohtml",
    )

    JLLWrappers.@init_executable_product(
        pdftoppm,
        "bin/pdftoppm",
    )

    JLLWrappers.@init_executable_product(
        pdftops,
        "bin/pdftops",
    )

    JLLWrappers.@init_executable_product(
        pdftotext,
        "bin/pdftotext",
    )

    JLLWrappers.@init_executable_product(
        pdfunite,
        "bin/pdfunite",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestEvent
  using Event
  using FactCheck

  facts("Observer creation should work") do
    @fact length(Event.create().events) --> 0
  end

  facts("Simple event binding should work") do
    flag = false
    obs  = Event.create()
    function handler() flag = true end

    Event.on(obs, "event", handler)
    Event.fire(obs, "event")

    @fact flag --> true
  end

  facts("Multiple event binding should work") do
    flag1 = false
    flag2 = false
    obs   = Event.create()
    function handler1() flag1 = true end
    function handler2() flag2 = true end

    Event.on(obs, "event", handler1)
    Event.on(obs, "event", handler2)
    Event.fire(obs, "event")

    @fact flag1 --> true
    @fact flag2 --> true
  end

  facts("Unbinding events should work") do
    flag1 = false
    flag2 = false
    obs   = Event.create()
    function handler1() flag1 = true end
    function handler2() flag2 = true end

    Event.on(obs, "event", handler1)
    Event.on(obs, "event", handler2)
    Event.off(obs, "event", handler1)
    Event.off(obs, "event", handler2)
    Event.fire(obs, "event")

    @fact flag1 --> false
    @fact flag2 --> false
  end

  facts("clear() should clear all handlers") do
    flag1 = false
    flag2 = false
    obs   = Event.create()
    function handler1() flag1 = true end
    function handler2() flag2 = true end

    Event.on(obs, "event", handler1)
    Event.on(obs, "event", handler2)
    Event.clear(obs)
    Event.fire(obs, "event")

    @fact flag1 --> false
    @fact flag2 --> false
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
using Term

import Pkg; Pkg.add("Documenter")
using Documenter

DocMeta.setdocmeta!(Term, :DocTestSetup, :(using Term); recursive = true)

makedocs(;
    modules = [Term],
    authors = "FedeClaudi <federicoclaudi@protonmail.com> and contributors",
    repo = "https://github.com/FedeClaudi/Term.jl/blob/{commit}{path}#{line}",
    sitename = "Term.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://FedeClaudi.github.io/Term.jl",
        assets = String[],
    ),
    pages = [
        "Home"=> "index.md",
        "Manual" => [
            "Styled text"           => "styled_text.md",
            "Renderables"           => "renderables.md",
            "Layout"                => "content_layout.md",
            "Logging"               => "logging.md",
            "Errors tracebacks"     => "errors_tracebacks.md",
            "Code introspection"    => "introspection.md",
        ],
        "API" => [
            "Term"          => "api_term.md",
            "Box"           => "api_box.md",
            "Color"         => "api_color.md",
            "Console"       => "api_console.md",
            "Inspect"       => "api_inspect.md",
            "Layout"        => "api_layout.md",
            "Logging"       => "api_logging.md",
            "Markup"        => "api_markup.md",
            "Measure"       => "api_measure.md",
            "Panel"         => "api_panel.md",
            "Renderables"   => "api_renderables.md",
            "Segment"       => "api_segment.md",
            "Style"         => "api_style.md",
        ],
    ],
)

deploydocs(; repo = "github.com/FedeClaudi/Term.jl", devbranch = "gh-pages")


text:
# This file contains functions to get the events data

Base.convert(::Type{JSON3.Array}, x::Vector{JSON3.Object}) = x

"""
    WyscoutEvent
Definition of a struct containing event stream data for a particular event.
"""
Base.@kwdef mutable struct WyscoutEvent
    event_id::Int
    game_id::Int
    milliseconds::Float64
    period_id::String
    player_id::Int
    positions::Vector{Dict}
    subtype_id::Int
    subtype_name::String
    tags::Vector{Dict}
    team_id::Int
    type_id::Int
    type_name::String
end

"""
    function events()
Return a Dataframe with all events in the game
"""
function events(events_data::PublicWyscoutLoader, game_id::Int)::Vector{WyscoutEvent}
    subset_test = subset_indexes_lazy(events_data.events, game_id)

    subset_dic = convert.(Dict, subset_test)
    vector_events = Vector{WyscoutEvent}()

    for i in eachindex(subset_dic)

        event = WyscoutEvent(
            event_id = subset_dic[i]["id"],
            game_id = subset_dic[i]["matchId"],
            milliseconds = subset_dic[i]["eventSec"],
            period_id =  string(subset_dic[i]["matchPeriod"]),
            player_id = subset_dic[i]["playerId"],
            positions = convert.(Dict,subset_dic[i]["positions"]),
            subtype_id = if typeof(subset_dic[i]["subEventId"]) == LazyJSON.String{String} 0 else subset_dic[i]["subEventId"] end,
            subtype_name = string(subset_dic[i]["subEventName"]),
            tags = convert.(Dict,subset_dic[i]["tags"]),
            team_id = convert(Int,subset_dic[i]["teamId"]),
            type_id = subset_dic[i]["eventId"],
            type_name = string(subset_dic[i]["eventName"])
        )
        push!(vector_events, event)
    end

    return vector_events
end

"""
    function subset_index_lazy

Use the vector Any v which is the JSON parsed lazily and the game_id to request data
with matchId correspoding to the game id.
"""
function subset_indexes_lazy(v::Vector{Any}, game_id::Int)
    json_data = Vector{Any}()
    @inbounds for i in eachindex(v)
        if v[i]["matchId"] == game_id
            push!(json_data, v[i])
        end
    end
    return json_data
end


text:
using Optimize
using NLPModels

using Plots
# pyplot()
plotly()

n_min = 2
n_max = 200

TestCUTEst = false

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

# Select solvers
using LineSearch
using Stopping
using LSDescentMethods

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

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

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

options = [Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => Newarmijo_wolfe)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => _hagerzhang2!, :τ₀ => 0.01, :mayterminate => true, :τ₁ => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => _morethuente2!, :τ₀ => 0.01, :τ₁ => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => TR_SecA_ls, :τ₀ => 0.000000010000000, :τ₁ => 0.092999999998146)
           Dict{Symbol,Any}(:verbose=>false, :stp => stp, :linesearch => TR_Nwt_ls, :τ₀ => 0.0001, :τ₁ => 0.001)
           ]

include("compare_solvers.jl")

P4 = compare_solvers_with_options2(solvers, options, labels, test_probs, n_min, n_max, printskip = false)[2]


text:
module package
using Cxx
export getAll, getPath, command

"""
    getAll()

Gets a vector of all packages.
"""
getAll() = convert(Vector{String}, icxx"std::vector<std::string> s; ros::package::getAll(s); return s;")

"""
    getPath(s::String)

Returns the fully-qualified path to a package, or an empty string if the package is not found.
"""
getPath(s::String) = unsafe_string(@cxx ros::package::getPath(pointer(s)))

"""
    command(s::String)

Runs a rospack command of the form 'rospack `s`', returning the output as a single string.
"""
command(s::String) = unsafe_string(@cxx ros::package::command(pointer(s)))

end


text:
function normalizedcut(A::AbstractMatrix{T}, b::AbstractVector{T}, k::Int,
                        isfun::Bool, tol::Float64) where {T}
    a = @view A[k, :]
    β = b[k]
    na = norm(a, 2)
    if isfun || na < tol
        a, β
    else
        a / na, β / na
    end
end

"""
    checkredundancy{T}(A::AbstractMatrix{T}, b::AbstractVector{T},
                       Anew::AbstractMatrix{T}, bnew::AbstractVector{T},
                       isfun::Bool, islb::Bool, tol::Float64, ident::Bool=false)

Check redundant cuts between the Polyhedra (A,b) and (Anew, bnew).
Return index of redundant cuts in `Anew`.

# Arguments
* `isfun::Bool`
    States if the Polyhedra defines a function
* `islb::Bool`
    States if the Polyhedra is a lower-bound or an upper-bound
* `tol::Float64`
    Tolerance of redundancy check
* `ident::Bool`
    States whether (A,b)==(Anew,bnew) if we want to remove redundant
    lines in a single Polyhedra

"""
function checkredundancy(A::AbstractMatrix{T}, b::AbstractVector{T},
                            Anew::AbstractMatrix{T}, bnew::AbstractVector{T},
                            isfun::Bool, islb::Bool, tol::Float64,
                            ident::Bool=false) where {T}
    # index of redundants cuts
    redundants = Int[]
    # number of new lines
    nnew = size(Anew, 1)

    for kk in 1:nnew
        a, β = normalizedcut(Anew, bnew, kk, isfun, tol)
        chk, indk = isinside(A, b, a, isfun, tol)
        # if ident is true, we need to take care of the possible match
        # of the cut a with itself in the matrix Anew
        if chk && (~ident || indk!=kk)
            ared, βred = normalizedcut(A, b, indk, isfun, tol)
            if islb ? β <= βred+tol : β+tol >= βred
                push!(redundants, kk)
            end
        end
    end

    redundants
end

"""Check if `λ` is a line of matrix `A`. `λ` might not have the same `eltype`
as `A` and `b` as it might have been scaled by `normalizecut`."""
function isinside(A::AbstractMatrix{T}, b::AbstractVector{T},
                     λ::AbstractVector, isfun::Bool, tol::Float64) where {T}
    nlines = size(A, 1)

    check = false
    k = 0
    while ~check && k < nlines
        k += 1
        a, β = normalizedcut(A, b, k, isfun, tol)
        check = (norm(a - λ, Inf) < tol)
    end
    check, k
end

checkredundancy(A, b, isfun, islb, tol)=checkredundancy(A, b, A, b, isfun, islb, tol, true)


text:
module PlusPlus

export ++

"""
++(a, b)

Concatenates strings much better than *
"""
++(a::AbstractString, b::AbstractString) = a*b

"""
++(a, b...)
Concatenates Vectors much better than vcat
"""
++(a::AbstractVector, b::AbstractVector...) = vcat(a, b...)

end # module


text:
using FFTW

struct CKKSEncoding{ScaleT} <: AbstractVector{Complex{Float64}}
    # TODO: Should this just be looked up by type?
    PlainT::Type
    # We permute the elements of the encoding in order to make the galois group
    # act cyclically on the plaintext elements.
    data::OffsetVector{Complex{Float64}}
end
Base.size(e::CKKSEncoding) = Base.size(e.data)
Base.axes(e::CKKSEncoding) = Base.axes(e.data)
Base.getindex(e::CKKSEncoding, i::Integer) = Base.getindex(e.data, i)
Base.setindex!(e::CKKSEncoding, v, i::Integer) = Base.setindex!(e.data, v, i)
scale(::Type{CKKSEncoding{ScaleT}}) where {ScaleT} = ScaleT
drop_last(::Type{CKKSEncoding{ScaleT}}) where {ScaleT} = CKKSEncoding{drop_last(ScaleT)}

"""
    ℤmstarPermutation

Recall that in a packed encoding the individual plaintext slots `c_i` correspond
to evaluation of the polynomial `p` at the i-th power of some m-th root of unity
``ζ_m``` for ``i ∈ ℤ_m^\\star`` where ``ℤ_m^\\star`` is the multiplicative group
of integers modulo ``m`` [1]. Alternatively, for m a power of 2, recall that
the ``c_i`` are the length ϕ(n) negacyclic fft of the coefficents of `p`
(the correspondence between these two is explained in the manual). This
identification is useful, because `ℤ_m^\\star` generally decomposes as a product
of cyclic groups. For example for `m=2^k` a power of two, we have:

```math
    \\mathbb{Z}_m^\\star \\simeq C_2 \\times C_{2^{k-2}}
```

Another useful representation of this cyclic group is a `2x(k-2)` matrix of
plaintext slots where the group action acts as as [`circshift`](@ref) on the
rows and columns respectively (corresponding to the action in the ``C_2``
and ``C_{2^{k-2}}`` factors respectively).

This type implements that representation, arranging the elements of
``ℤ_m^\\star`` (as integers) in the mentioned matrix representation.

[1] https://en.wikipedia.org/wiki/Multiplicative_group_of_integers_modulo_n
"""
struct ℤmstarPermutation <: AbstractArray{Int, 2}
    M::Int
    ℤmstarPermutation(m::Integer) = (@assert ispow2(m); new(convert(Int,m)))
end

Base.size(p::ℤmstarPermutation) = (2, div(p.M, 4))
function Base.getindex(p::ℤmstarPermutation, row::Integer, col::Integer)
    ℤmstar_generator = 3
    rs, cls = size(p)
    (row in 1:rs && col in 1:cls) || Base.throw_boundserror(p, (row, col))
    mod(row == 2 ? p.M-ℤmstar_generator^col : ℤmstar_generator^col, p.M)
end

function CKKSEncoding{ScaleT}(plain::PlainT) where {ScaleT, PlainT}
    # Decode
    scaled = map(x->reinterpret(ScaleT, x), NTT.coeffs_primal(plain))
    scaled = Float64[convert(Float64, x) for x in scaled]
    # Undo the root of unity premul
    multed = map(x->convert(Float64, x), scaled) .* [Complex{Float64}(exp(big(-2*k/(2*length(scaled))*pi*im))) for k in eachindex(scaled)]
    # FFT it and take only the non-conjugated coefficients
    plain = OffsetArray(fft(collect(multed)), 0:length(multed)-1)[ℤmstarPermutation(2length(multed))[1,:] .>> 1]
    plain = OffsetArray(plain, 0:(length(multed)÷2-1))
    CKKSEncoding{ScaleT}(PlainT, plain)
end

function Base.convert(::Type{NTT.RingElement}, s::CKKSEncoding{<:Any})
    convert(s.PlainT, s)
end

function Base.convert(T::Type{<:NTT.RingElement}, s::CKKSEncoding{ScaleT}) where {ScaleT}
    data = collect(s.data)

    cmplx = OffsetArray(zeros(Complex{Float64}, 2*length(s.data)), 0:2*length(s.data)-1)
    idxs = ℤmstarPermutation(4length(s.data)) .>> 1
    for i in axes(s.data, 1)
        cmplx[idxs[1, i+1]] = s.data[i]
        cmplx[idxs[2, i+1]] = conj(s.data[i])
    end

    ipoints = OffsetArray(ifft(cmplx.parent), 0:2*length(s.data)-1)
    # Make it negacyclic
    nipoints = ipoints .* [Complex{Float64}(exp(big(2*k/(2*length(ipoints))*pi*im))) for k in eachindex(ipoints)]

    # Project to real coefficents
    nipointsr = map(nipoints) do p
        @assert isapprox(imag(p), 0, atol=10^-10)
        real(p)
    end

    # Encode into the fixed fraction representation
    encoded = map(ScaleT{eltype(s.PlainT)}, nipointsr)

    plaintext = map(x->x.x, encoded)
    s.PlainT(plaintext, nothing)
end

# Multiplication with plaintext scalars
function Base.:*(a::CipherText{CKKSEncoding{Tscale}, P, T, N}, b::AbstractFloat) where {Tscale, P, T, N}
    scaled = convert(Integer, Tscale{coefftype(T)}(b).x)
    CipherText{CKKSEncoding{Tscale^2}, P, T, N}(a.params, map(c->c*scaled, a.cs))
end

function Base.broadcasted(::typeof(*), a::Array{<:AbstractFloat, 1}, c::CipherText{CKKSEncoding{Tscale}, P, T, N}) where {Tscale, P, T, N}
    plain = CKKSEncoding{Tscale}(zero(c.cs[1]))
    plain .= OffsetArray(a, 0:length(a)-1)
    re = convert(NTT.RingElement, plain)
    CipherText{CKKSEncoding{Tscale^2}, P, T, N}(c.params, map(c->c*re, c.cs))
end

function Base.broadcasted(::typeof(+), a::CipherText{CKKSEncoding{Tscale}, P, T, N}, b::AbstractFloat) where {Tscale, P, T, N}
    plain = CKKSEncoding{Tscale}(zero(T))
    plain .= b
    CipherText{CKKSEncoding{Tscale}, P, T, N}(a.params, (a.cs[1] + convert(T, plain), a.cs[2:end]...))
end

function Base.broadcasted(::typeof(+), a::CipherText{CKKSEncoding{Tscale}, P, T, N}, b::AbstractArray{<:AbstractFloat, 1}) where {Tscale, P, T, N}
    plain = CKKSEncoding{Tscale}(zero(T))
    plain .= b
    @show typeof(convert(T, plain))
    CipherText{CKKSEncoding{Tscale}, P, T, N}(a.params, (a.cs[1] + convert(T, plain), a.cs[2:end]...))
end

# Modswitching a ciphertext
function modswitch(c::CipherText{CKKSEncoding{Tscale}, P, T, N}) where {scale, Tscale<:FixedRational{scale}, CC<:CRTEncoded, P, T<:NTT.RingElement{<:Any, CC}, N}
    CipherText{CKKSEncoding{drop_last(FixedRational{scale/modulus(moduli(CC).parameters[end])})}}(modswitch_drop(c.params),
        map(modswitch, c.cs))
end

# Multiplication of ciphertexts
function *(c1::CipherText{C1, Enc, T}, c2::CipherText{C2, Enc, T}) where {C1 <: CKKSEncoding, C2 <: CKKSEncoding, Enc, T}
    CipherText{CKKSEncoding{scale(C1)*scale(C2)}}(c1.params, enc_mul(c1, c2))
end


text:
using LinearAlgebra
using NamedDims
using NamedDims: names
using Test

@testset "rename" begin
    nda = NamedDimsArray{(:a, :b, :c, :d)}(ones(10,1,1,20))
    new_nda = rename(nda, (:w, :x, :y, :z))

    @test names(new_nda) === (:w, :x, :y, :z)
    @test parent(new_nda) === parent(nda)
end

@testset "dropdims" begin
    nda = NamedDimsArray{(:a, :b, :c, :d)}(ones(10, 1, 1, 20))

    @test_throws ArgumentError dropdims(nda; dims=1)
    @test_throws ArgumentError dropdims(nda; dims=:a)

    @test dropdims(nda; dims=:b) == ones(10, 1, 20) == dropdims(nda; dims=2)
    @test names(dropdims(nda; dims=:b)) == (:a, :c, :d) == names(dropdims(nda; dims=2))

    @test dropdims(nda; dims=(:b,:c)) == ones(10, 20) == dropdims(nda; dims=(2, 3))
    @test names(dropdims(nda; dims=(:b, :c))) == (:a, :d) == names(dropdims(nda; dims=(2, 3)))
end

@testset "$f" for f in (adjoint, transpose, permutedims)
    @testset "Vector $f" begin
        ndv = NamedDimsArray{(:foo,)}([10,20,30])
        @test f(ndv) == [10 20 30]
        @test names(f(ndv)) == (:_, :foo)

        if f === permutedims
            # unlike adjoint and tranpose, permutedims should not be its own inverse
            # The new dimension should stick around
            @test f(f(ndv)) == reshape([10, 20, 30], Val(2))
            @test names(f(f(ndv))) == (:foo, :_)
        else
            # Make sure vector double adjoint gets you back to the start.
            @test f(f(ndv)) == [10, 20, 30]
            @test names(f(f(ndv))) == (:foo,)
        end
    end

    @testset "Matrix $f" begin
        ndm = NamedDimsArray{(:foo,:bar)}([10 20 30; 11 22 33])
        @test f(ndm) == [10 11; 20 22; 30 33]
        @test names(f(ndm)) == (:bar, :foo)

        # Make sure implementation of matrix double adjoint is correct
        # since it is easy for the implementation of vector double adjoint broke it
        @test f(f(ndm)) == [10 20 30; 11 22 33]
        @test names(f(f(ndm))) == (:foo, :bar)
    end
end

@testset "permutedims" begin
    nda = NamedDimsArray{(:w,:x,:y,:z)}(ones(10,20,30,40))
    @test names(permutedims(nda, (:w,:x,:y,:z))) == (:w,:x,:y,:z) == names(permutedims(nda, 1:4))
    @test size(permutedims(nda, (:w,:x,:y,:z))) == (10,20,30,40) == size(permutedims(nda, 1:4))

    @test names(permutedims(nda, (:w,:y,:x,:z))) == (:w,:y,:x,:z) == names(permutedims(nda, (1,3,2,4)))
    @test size(permutedims(nda, (:w,:y,:x,:z))) == (10,30,20,40) == size(permutedims(nda, (1,3,2,4)))

    @test_throws Exception permutedims(nda, (:foo,:x,:y,:z))
    @test_throws Exception permutedims(nda, (:x,:y,:z))
    @test_throws Exception permutedims(nda, (:x,:x,:y,:z))

    @test_throws Exception permutedims(nda, (0,1,2,3))
    @test_throws Exception permutedims(nda, (2,3,4))
    @test_throws Exception permutedims(nda, (2,2,3,4))
end

# We test pinv here as it is defined in src/function_dims.jl
# using the same logic as permutedims, transpose etc
@testset "pinv" begin
    @testset "Matrix" begin
        nda = NamedDimsArray{(:a, :b)}([1.0 2 3; 4 5 6])
        @test names(pinv(nda)) == (:b, :a)
        @test nda * pinv(nda) ≈ NamedDimsArray{(:a, :a)}([1.0 0; 0 1])
    end
    @testset "Vector" begin
        nda = NamedDimsArray{(:a,)}([1.0, 2, 3])
        @test names(pinv(nda)) == (:_, :a)

        @test names(pinv(pinv(nda))) == (:a,)
        @test pinv(pinv(nda)) ≈ nda

        # See issue https://github.com/invenia/NamedDims.jl/issues/36
        @test_broken nda * pinv(nda) ≈ 1.0
    end
end