text:
struct Characteristic <: AbstractCharacteristic
  ws::WeightSpec
  Fξ::Symbol
  X::Dict{Symbol, Symbol}
end

#builds a characteristic object, which includes cross-sectional weights and
#weight coefficients
function Characteristic(ws::WeightSpec, Fξ::Symbol)

  Fw = Symbol(:w, ws.weightprefix, Fξ)

  #a couple of helper funtions for consistent renaming
  lagname(s::Symbol) = Symbol(:L, s)
  retname(s::Symbol) = Symbol(:R, s)
  rettname(s::Symbol) = Symbol(:Rt, s)

  #WARNING: HARDCODE ALERT!!!
  X = Dict{Symbol, Symbol}()
  X[:Fw] = Fw
  X[:FLw] = lagname(X[:Fw])
  X[:FLLw] = lagname(X[:FLw])
  X[:FLLLw] = lagname(X[:FLLw])
  #not used for now
  X[:FRLw] = retname(X[:FLw])
  X[:FRLLw] = retname(X[:FLLw])

  X[:FRtw] = rettname(X[:Fw])
  return Characteristic(ws, Fξ, X)
end

weightfields(ξ::Characteristic) = collect(values(ξ.X))
Characteristic(ξ::Characteristic,  Fξ::Symbol, X::Dict
  ) = Characteristic(ξ.ws, Fξ, X)

demeanξ(ξraw::AbstractVector{<:Union{Real, Missing}}) = ξraw .- mean(skipmissing(ξraw))
#this gives the quantile fo each value in an array
#would have preferred the ecdf function, but it doesn't play nice with
#containers that can hold missing values, even if none are actually held
#=function quantiles(ξraw::AbstractVector{<:Union{Real,Missing}})
  (length(ξraw) .- competerank(ξraw,rev=true) .+ 1)./length(ξraw)
end=#

function quantiles(ξraw::AbstractVector{<:Union{Real,Missing}})
  (length(ξraw) .- competerank(ξraw,rev=true) .+ 1)./length(ξraw)
end

#alternate procedure
#=function quantiles(ξraw::AbstractVector{<:Union{Real,Missing}})
  F = ecdf(ξraw |> Vector{Float64})
  return map(F, ξraw)
end=#

#conditions each slice of the panel for a particular characteristic
function conditionξ!(F::Function, panel::AbstractDataFrame, Fξraw::Symbol, Fξ::Symbol)::Nothing

  (sum((!ismissing).(panel[!, Fξ])) > 0) && error("$Fξ in panel already contains data")
  spanels::GroupedDataFrame = groupby(view(panel, (!ismissing).(panel[!,Fξraw]), :), :date)
  Threads.@threads for i ∈ 1:length(spanels)
      spanel::SubDataFrame = spanels[i]
      Fξsubcol::SubArray = spanel[!, Fξ]
      Fξsubcol .= F(spanel[!, Fξraw])
  end

  return nothing
end

#=standard conditioning of a characteristic column
function conditionξ!(panel::AbstractDataFrame, Fξraw::Symbol;
  prefix::Symbol = PARAM[:charprefix],
  Fξ::Symbol = Symbol(prefix, Fξraw))::Nothing

  #allocate the data
  ξtype::Type = Union{eltype(panel[!,Fξraw]), Missing}
  panel[!, Fξ] = missings(ξtype, size(panel,1))

  conditionξ!(demeanξ, panel, Fξraw, Fξ)
  return nothing
end

#apply a quantile transformation before demeaning
function conditionξquantile!(panel::AbstractDataFrame, Fξraw::Symbol;
  prefix::Symbol = PARAM[:charprefix],
  Fξ::Symbol = Symbol(prefix, Fξraw))::Nothing

  #allocate the data
  panel[!, Fξ] = Vector{Union{eltype(panel[!,Fξraw])}}(undef, size(panel,1))

  conditionξ!(demeanξ ∘ quantiles, panel, Fξraw, Fξ)
  return nothing
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Pkg
Pkg.activate(@__DIR__)

using Dates
using HTTP
using MySQL
using JSON3
using StructTypes

struct jsonMsgObj
    message
end

struct jsonObj
    id
    randomNumber
end

StructTypes.StructType(::Type{jsonMsgObj}) = StructTypes.Struct()
StructTypes.StructType(::Type{jsonObj}) = StructTypes.Struct()

    function plaintext(req::HTTP.Request)
        headers = [ "Content-Type" => "text/plain",
                    "Server" => "Julia-HTTP",
                    "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT" ]

        return HTTP.Response(200, headers, body = "Hello, World!")
    end

    function jsonSerialization(req::HTTP.Request)
        headers = [ "Content-Type" => "application/json",
                    "Server" => "Julia-HTTP",
                    "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT" ]

        return HTTP.Response(200, headers, body = JSON3.write(jsonMsgObj("Hello, World!")))
    end

    function singleQuery(req::HTTP.Request)
        headers = [ "Content-Type" => "application/json",
                    "Server" => "Julia-HTTP",
                    "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT" ]

        randNum = rand(1:10000)

        conn = DBInterface.connect(MySQL.Connection, "tfb-database", "benchmarkdbuser", "benchmarkdbpass", db="hello_world")
        sqlQuery = "SELECT * FROM World WHERE id = $randNum"
        results = DBInterface.execute(conn, sqlQuery)
        row = first(results)
        dbNumber = row[2]

        DBInterface.close!(conn)
        return HTTP.Response(200, headers, body = JSON3.write(jsonObj(randNum, dbNumber)))
    end

    function multipleQueries(req::HTTP.Request)
        headers = [ "Content-Type" => "application/json",
                    "Server" => "Julia-HTTP",
                    "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT" ]

        numQueries = -1

        try
            numQueries = parse(Int64, (split(req.target, "="))[2])

        catch ArgumentError
            numQueries = 1

        finally
            if numQueries > 500
                numQueries = 500
            end

            if numQueries < 1
                numQueries = 1
            end
        end

        conn = DBInterface.connect(MySQL.Connection, "tfb-database", "benchmarkdbuser", "benchmarkdbpass", db="hello_world")

        responseArray = Array{jsonObj}(undef, numQueries)
        for i in 1:numQueries
            # randNum = randNumList[i]
            randNum = rand(1:10000)
            sqlQuery = "SELECT * FROM World WHERE id = $randNum"
            results = DBInterface.execute(conn, sqlQuery)
            row = first(results)
            dbNumber = row[2]
            responseArray[i] = jsonObj(randNum, dbNumber)
        end

        DBInterface.close!(conn)

        return HTTP.Response(200, headers, body = JSON3.write(responseArray))
    end

    function updates(req::HTTP.Request)
        headers = [ "Content-Type" => "application/json",
                    "Server" => "Julia-HTTP",
                    "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT" ]

        numQueries = -1

        try
            numQueries = parse(Int64, (split(req.target, "="))[2])

        catch ArgumentError
            numQueries = 1

        finally
            if numQueries > 500
                numQueries = 500
            end

            if numQueries < 1
                numQueries = 1
            end
        end

        conn = DBInterface.connect(MySQL.Connection, "tfb-database", "benchmarkdbuser", "benchmarkdbpass", db="hello_world")

        responseArray = Array{jsonObj}(undef, numQueries)
        for i in 1:numQueries
            randId = rand(1:10000)
            randNum = rand(1:10000)
            sqlQuery = "SELECT * FROM World WHERE id = $randId"
            results = DBInterface.execute(conn, sqlQuery)
            row = first(results)
            dbNumber = row[2]

            sqlQuery = "UPDATE World SET randomnumber = $randNum WHERE id = $randId"
            results = DBInterface.execute(conn, sqlQuery)
            responseArray[i] = jsonObj(randId, randNum)
        end

        DBInterface.close!(conn)
        return HTTP.Response(200, headers, body = JSON3.write(responseArray))
    end

    function fortunes(req::HTTP.Request)
        headers = [ "Content-Type" => "text/html; charset=utf-8",
                    "Server" => "Julia-HTTP",
                    "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT" ]

        fortunesList = []
        sqlQuery = "SELECT * FROM fortune"
        output = "<!DOCTYPE html><title>Fortunes</title>

id	message
"
        conn = DBInterface.connect(MySQL.Connection, "tfb-database", "benchmarkdbuser", "benchmarkdbpass", db="hello_world")
        results = DBInterface.execute(conn, sqlQuery)

        for row in results
            push!(fortunesList, [string(row[1]), row[2]])
        end

        push!(fortunesList, [string(0), "Additional fortune added at request time."])

        sort!(fortunesList, by = x -> x[2])

        for fortune in fortunesList
            id = fortune[1]
            message = HTTP.Strings.escapehtml(fortune[2])
            output = string(output, "
$id	$message
")
        end

        output = string(output, "
")

        DBInterface.close!(conn)
        return HTTP.Response(200, headers, body = output)
    end

const ROUTER = HTTP.Router()
HTTP.@register(ROUTER, "GET", "/plaintext", plaintext)
HTTP.@register(ROUTER, "GET", "/json", jsonSerialization)
HTTP.@register(ROUTER, "GET", "/db", singleQuery)
HTTP.@register(ROUTER, "GET", "/queries", multipleQueries)
HTTP.@register(ROUTER, "GET", "/updates", updates)
HTTP.@register(ROUTER, "GET", "/fortunes", fortunes)

HTTP.serve(ROUTER, "0.0.0.0" , 8080, reuseaddr=true)

ext: jl

max_stars_repo_name: efectn/FrameworkBenchmarks

max_issues_repo_name: efectn/FrameworkBenchmarks

max_forks_repo_name: efectn/FrameworkBenchmarks


text:
export TransferOperator,
    transferoperator,
    TransferOperatorGenerator,
    transopergenerator,
    GridBasedTransferOperator

import CausalityToolsBase.BinningScheme

""" Supertype of all transfer operator estimators. """
abstract type TransferOperator end

""" Supertype of all grid-based transfer operator estimators """
abstract type GridBasedTransferOperator <: TransferOperator end

"""
    TransferOperatorGenerator{E <: TransferOperator}(method::E, pts, init)

A transfer operator generator. Contains the partition information,
the method for transition probability estimation, the input `pts` and
`init`, which are things that can be pre-computed to speed up estimation
process.
"""
struct TransferOperatorGenerator{E <: TransferOperator, X, A}
    method::E # estimator with its input parameters
    pts::X    # the phase space / reconstruted state space points
    init::A   # pre-initialized things that speed up estimation process
end

"""
    transopergenerator(pts, method::TransferOperator) → to::TransferOperatorGenerator

Initialize a generator that creates transfer operators on demand, based on the given `method`.
This is efficient, because some things can be initialized and reused.

To approximate a transfer operator, call `to` as a function with the relevant arguments.

```julia
to = transopergenerator(x, SingleGrid(RectangularBinning(5)))
for i in 1:1000
    s = to()
    # do stuff with s and or x
    result[i] = stuff
end
```
"""
function transopergenerator end

"""
    transferoperator(pts, method::TransferOperator) → to::TransferOperatorApproximation
    transferoperator(pts, method::SingleGrid) → to::TransferOperatorApproximation
    transferoperator(pts, method::SimplexPoint; tol = 1e-8, randomsampling::Bool = false, n::Int = 100) → to::TransferOperatorApproximation
    transferoperator(pts, method::SimplexExact; tol::Real = 1e-8) → to::TransferOperatorApproximation

Create a transfer operator `to` from the phase/state space points `x` based on the given `method`.

Using the `SingleGrid` estimator will partition the state space into rectangular bins, while
the `SimplexPoint` and `SimplexExact` estimators will triangulate the state space.

- If using any of the triangulation estimators, `Simplices` must be imported after `PerronFrobenius`.
- If using `SimplexPoint` and `randomsampling==true`, then each simplex is represented by a
    point cloud uniformly distributed within each simplex. If `randomsampling==false`, then
    an even simplex splitting routine is used to generate interior points. `n` is the
    minimum number of sampling points (more points may be used if `randomsampling==false`,
    because the number of points generated by the splitting routine is dictated by the dimension
    and splitting factor).

## Examples

### Grid estimators

Use the `SingleGrid` estimator for long time series.

```julia
n = 30000
x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
y = [cos(2π*(t + sin(t)/10)) for t in 1:n] .+ 0.5rand(n)

τs = (2, 0, -2)
js = (2, 2, 1)
D = genembed([x, y], τs, js)

method = SingleGrid(RectangularBinning(5))
transferoperator(D, method)
```

### Triangulation estimators

Use the `SimplexExact` (slow) estimator or `SimplexPoint` (also slow,
but faster than `SimplexExact`) estimator for very short time series.
For these estimators to be imported, you need to load `Simplices` *after*
loading `PerronFrobenius`.

```julia
n = 30
x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
y = [cos(2π*(t + sin(t)/10)) for t in 1:n] .+ 0.5rand(n)

τs = (2, 0, -2)
js = (2, 2, 1)
D = genembed([x, y], τs, js)

# Compute transition probabilities using approximate simplex intersections.
# Sample each simplex using at least `n = 300` points.
method = SimplexApprox()
transferoperator(D, method, n = 300)
```
"""
function transferoperator(pts, method::TransferOperator)
    to = transopergenerator(pts, method)
    to()
end


text:
using Statistics

Statistics.varm(A::AbstractGPUArray{<:Real},m::AbstractArray{<:Real}; dims, corrected::Bool=true) =
    sum((A .- m).^2, dims=dims)/(prod(size(A)[[dims...]])::Int-corrected)

Statistics.stdm(A::AbstractGPUArray{<:Real},m::AbstractArray{<:Real}, dim::Int; corrected::Bool=true) =
    sqrt.(varm(A,m;dims=dim,corrected=corrected))

Statistics._std(A::AbstractGPUArray, corrected::Bool, mean, dims) =
    sqrt.(Statistics.var(A; corrected=corrected, mean=mean, dims=dims))

Statistics._std(A::AbstractGPUArray, corrected::Bool, mean, ::Colon) =
    sqrt.(Statistics.var(A; corrected=corrected, mean=mean, dims=:))

# Revert https://github.com/JuliaLang/Statistics.jl/pull/25
Statistics._mean(A::AbstractGPUArray, ::Colon)    = sum(A) / length(A)
Statistics._mean(f, A::AbstractGPUArray, ::Colon) = sum(f, A) / length(A)
Statistics._mean(A::AbstractGPUArray, dims)    = mean!(Base.reducedim_init(t -> t/2, +, A, dims), A)
Statistics._mean(f, A::AbstractGPUArray, dims) = sum(f, A, dims=dims) / mapreduce(i -> size(A, i), *, unique(dims); init=1)

function Statistics.covzm(x::AbstractGPUMatrix, vardim::Int=1; corrected::Bool=true)
    C = Statistics.unscaled_covzm(x, vardim)
    T = promote_type(typeof(one(eltype(C)) / 1), eltype(C))
    A = convert(AbstractArray{T}, C)
    b = 1//(size(x, vardim) - corrected)
    A .*= b
    return A
end

function Statistics.cov2cor!(C::AbstractGPUMatrix{T}, xsd::AbstractGPUArray) where T
    nx = length(xsd)
    size(C) == (nx, nx) || throw(DimensionMismatch("inconsistent dimensions"))
    tril!(C, -1)
    C += adjoint(C)
    C = Statistics.clampcor.(C ./ (xsd * xsd'))
    C[diagind(C)] .= oneunit(T)
    return C
end

function Statistics.corzm(x::AbstractGPUMatrix, vardim::Int=1)
    c = Statistics.unscaled_covzm(x, vardim)
    return Statistics.cov2cor!(c, sqrt.(diag(c)))
end


text:
module FETerms

using Gridap
using Gridap.Helpers

export FETerm
export AffineFETerm
export FESource
export LinearFETerm
export NonLinearFETerm

export setup_cell_jacobian
export setup_cell_residual
export setup_cell_matrix
export setup_cell_vector
export setup_cell_ids

# Interfaces

abstract type FETerm end

"""
Returns an object (e.g. a CellMatrix for single field problems) representing
the contribution to the Jacobian of the given term. Returns nothing if the term
has not contribution to the Jacobian (typically for source terms)
"""
setup_cell_jacobian(t::FETerm,uh,v,du) = @abstractmethod

"""
Returns an object (e.g. a CellVector) representing the contribution to the
residual of the given term. Returns always something.
"""
setup_cell_residual(t::FETerm,uh,v) = @abstractmethod

setup_cell_ids(t::FETerm)::CellNumber = @abstractmethod

abstract type AffineFETerm <: FETerm end

"""
Returns an object (e.g. a CellMatrix) representing the contribution to the
system matrix of the given term. Returns nothing if the term has not
contribution (typically for source terms)
"""
setup_cell_matrix(t::AffineFETerm,v,u) = @abstractmethod

"""
Returns an object (e.g. a CellVector) representing the contribution to the
system rhs of the given term. Returns nothing if the term has not
contribution (typically for linear terms)
"""
setup_cell_vector(t::AffineFETerm,v,uhd) = @abstractmethod

function setup_cell_jacobian(t::AffineFETerm,uh,v,du)
  setup_cell_matrix(t,v,du)
end

# Pretty printing

import Base: show

function show(io::IO,self::FETerm)
  print(io,"$(nameof(typeof(self))) object")
end

function show(io::IO,::MIME"text/plain",self::FETerm)
  show(io,self)
end

# Concrete implementations

struct AffineFETermFromIntegration <: AffineFETerm
  biform::Function
  liform::Function
  trian::Triangulation
  quad::CellQuadrature
end

function AffineFETerm(
  biform::Function,liform::Function,trian::Triangulation,quad::CellQuadrature)
  AffineFETermFromIntegration(biform,liform,trian,quad)
end

function setup_cell_matrix(t::AffineFETermFromIntegration,v,u)
  _v = _restrict_if_needed(v,t.trian)
  _u = _restrict_if_needed(u,t.trian)
  integrate(t.biform(_v,_u),t.trian,t.quad)
end

function setup_cell_vector(t::AffineFETermFromIntegration,v,uhd)
  _v = _restrict_if_needed(v,t.trian)
  _uhd = _restrict_if_needed(uhd,t.trian)
  integrate(t.liform(_v)-t.biform(_v,_uhd),t.trian,t.quad)
end

function setup_cell_residual(t::AffineFETermFromIntegration,uh,v)
  _v = _restrict_if_needed(v,t.trian)
  _uh = _restrict_if_needed(uh,t.trian)
  integrate(t.biform(_v,_uh)-t.liform(_v),t.trian,t.quad)
end

setup_cell_ids(t::AffineFETermFromIntegration) = _setup_cell_ids(t.trian)

struct FESourceFromIntegration <: AffineFETerm
  liform::Function
  trian::Triangulation
  quad::CellQuadrature
end

function FESource(
  liform::Function,
  trian::Triangulation,
  quad::CellQuadrature)
  FESourceFromIntegration(liform,trian,quad)
end

setup_cell_matrix(t::FESourceFromIntegration,v,u) = nothing

function setup_cell_vector(t::FESourceFromIntegration,v,uhd)
  _v = _restrict_if_needed(v,t.trian)
  integrate(t.liform(_v),t.trian,t.quad)
end

function setup_cell_residual(t::FESourceFromIntegration,uh,v)
  _v = _restrict_if_needed(v,t.trian)
  integrate(-t.liform(_v),t.trian,t.quad)
end

setup_cell_ids(t::FESourceFromIntegration) = _setup_cell_ids(t.trian)

struct LinearFETermFromIntegration <: AffineFETerm
  biform::Function
  trian::Triangulation
  quad::CellQuadrature
end

function LinearFETerm(
  biform::Function,trian::Triangulation,quad::CellQuadrature)
  LinearFETermFromIntegration(biform,trian,quad)
end

function setup_cell_matrix(t::LinearFETermFromIntegration,v,u)
  _v = _restrict_if_needed(v,t.trian)
  _u = _restrict_if_needed(u,t.trian)
  integrate(t.biform(_v,_u),t.trian,t.quad)
end

function setup_cell_vector(t::LinearFETermFromIntegration,v,uhd)
  _v = _restrict_if_needed(v,t.trian)
  _uhd = _restrict_if_needed(uhd,t.trian)
  integrate(-t.biform(_v,_uhd),t.trian,t.quad)
end

function setup_cell_residual(t::LinearFETermFromIntegration,uh,v)
  _v = _restrict_if_needed(v,t.trian)
  _uh = _restrict_if_needed(uh,t.trian)
  integrate(t.biform(_v,_uh),t.trian,t.quad)
end

setup_cell_ids(t::LinearFETermFromIntegration) = _setup_cell_ids(t.trian)

struct NonLinearFETerm <: FETerm
  res::Function
  jac::Function
  trian::Triangulation
  quad::CellQuadrature
end

function setup_cell_jacobian(t::NonLinearFETerm,uh,v,du)
  _v = _restrict_if_needed(v,t.trian)
  _uh = _restrict_if_needed(uh,t.trian)
  _du = _restrict_if_needed(du,t.trian)
  integrate(t.jac(_uh,_v,_du),t.trian,t.quad)
end

function setup_cell_residual(t::NonLinearFETerm,uh,v)
  _v = _restrict_if_needed(v,t.trian)
  _uh = _restrict_if_needed(uh,t.trian)
  integrate(t.res(_uh,_v),t.trian,t.quad)
end

setup_cell_ids(t::NonLinearFETerm) = _setup_cell_ids(t.trian)

# Dealing with several FETerms

function setup_cell_jacobian(uh,v,du,terms::Vararg{<:FETerm})
  m = []
  for term in terms
    r = _jac(term,uh,v,du)
    if r != nothing
      c = _cellids(term)
      _append_matrix_contribution!(m,r,c)
    end
  end
  m
end

function setup_cell_residual(uh,v,terms::Vararg{<:FETerm})
  m = []
  for term in terms
    r = setup_cell_residual(term,uh,v)
    c = _cellids(term)
    _append_vector_contribution!(m,r,c)
  end
  m
end

function setup_cell_matrix(v,u,terms::Vararg{<:AffineFETerm})
  m = []
  for term in terms
    r = _mat(term,v,u)
    if r != nothing
      c = _cellids(term)
      _append_matrix_contribution!(m,r,c)
    end
  end
  m
end

function setup_cell_vector(v,uhd,terms::Vararg{<:AffineFETerm})
  m = []
  for term in terms
    r = _vec(term,v,uhd)
    if r != nothing
      c = _cellids(term)
      _append_vector_contribution!(m,r,c)
    end
  end
  m
end

# Helpers

function _append_vector_contribution!(m,r,c)
  push!(m,(r,c))
end

function _append_vector_contribution!(m,r::SkeletonCellVector,c)
  push!(m,(r.cellvector1,c[1]))
  push!(m,(r.cellvector2,c[2]))
end

function _append_matrix_contribution!(m,r,c)
  push!(m,(r,c,c))
end

function _append_matrix_contribution!(m,r::SkeletonCellMatrix,c)
  push!(m,(r.cellmatrix11,c[1],c[1]))
  push!(m,(r.cellmatrix12,c[1],c[2]))
  push!(m,(r.cellmatrix21,c[2],c[1]))
  push!(m,(r.cellmatrix22,c[2],c[2]))
end

const _jac = setup_cell_jacobian

const _mat = setup_cell_matrix

const _vec = setup_cell_vector

const _cellids = setup_cell_ids

_restrict_if_needed(u,trian) = u

_restrict_if_needed(u,trian::BoundaryTriangulation) = restrict(u,trian)

_restrict_if_needed(u,trian::SkeletonTriangulation) = restrict(u,trian)

_setup_cell_ids(trian) = IdentityCellNumber(Int,ncells(trian))

_setup_cell_ids(trian::BoundaryTriangulation) = trian.descriptor.facet_to_cell

function _setup_cell_ids(trian::SkeletonTriangulation)
  n1 = trian.descriptor1.facet_to_cell
  n2 = trian.descriptor2.facet_to_cell
  (n1, n2)
end

end # module


text:
import BenchmarkProfiles

# Customized performance_profile method, derived from BenchmarkProfiles.jl
#Copyright (c) 2016: Dominique Orban.
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
function performance_profile(T :: Array{Float64,2}, labels :: Vector{String};
                             logscale :: Bool=true,
                             title :: String="",
                             ymax = 1.1,
                             xmax = NaN,
                             linewidth=1,
                             kwargs...)

  (ratios, max_ratio) = BenchmarkProfiles.performance_ratios(T, logscale=logscale)
  (np, ns) = size(ratios)

  ratios = [ratios; 2.0 * max_ratio * ones(1, ns)]
  xs = [1:np+1;] / np
  length(labels) == 0 && (labels = [@sprintf("column %d", col) for col = 1 : ns])
  profile = Plots.plot(; kwargs...)  # initial empty plot
  for s = 1 : ns
    Plots.plot!(ratios[:, s], xs, t=:steppost, label=labels[s], linewidth=linewidth)
  end
  if isfinite(xmax)
      Plots.xlims!(logscale ? 0.0 : 1.0, xmax)
  else
      Plots.xlims!(logscale ? 0.0 : 1.0, 1.1 * max_ratio)
  end
  Plots.ylims!(0, ymax)
  Plots.xlabel!("\$F\$" * (logscale ? " (log scale)" : ""))
  Plots.ylabel!("\$P\$")
  Plots.title!(title)
  return profile
end

max_stars_repo_name: mlubin/PajaritoSupplement

max_issues_repo_name: mlubin/PajaritoSupplement

max_forks_repo_name: mlubin/PajaritoSupplement


text:
type IIDMeodInput <: AbstractInput end

################################################################################
################################################################################
## Untested

## Read inputs
function read_input(intput_type::T, filenames::Array{String}) where T<:Type{IIDMeodInput}
  data_branches = readdlm(filenames[1])
  data_buses = readdlm(filenames[2])
  data_generators = readdlm(filenames[3])
  data_limits = readdlm(filenames[4])
  data_loads = readdlm(filenames[5])
  data_substations = readdlm(filenames[6])

  bus = Dict{String, Dict{String, Any}}()

  for i=1:size(data_buses, 1)
      busid = Int(data_buses[i, 1])
      substationid = Int(data_buses[i, 2])
      V = data_buses[i, 4] # pu
      θ = data_buses[i, 5] # radians
      P, Q = data_buses[i, 6:7] # MW, MVar
      basekV, minV, maxV = data_substations[substationid, 4:6] # kV, pu, pu

      busname = "BUS_$busid"
      data_buses[i, 3] == 0 || warn("Buses - $busname : parameter cc!=0 ($(data_buses[i, 3])), not handled")
      data_buses[i, 8] == 0 || warn("Buses - $busname : parameter fault!=0 ($(data_buses[i, 8])), not handled")
      data_buses[i, 9] == 0 || warn("Buses - $busname : parameter curative!=0 ($(data_buses[i, 9])), not handled")

      bus[busname] = Dict{String, Any}()
      bus_elems = bus[busname]
      bus_elems["Volt"] = IIDMVolt(busname, busid, string(substationid), basekV, minV, maxV)
      bus_elems["VoltMeasure"] = IIDMVoltMeasure(V*exp(im*θ), P+im*Q)
  end

  for i=1:size(data_loads, 2)
      busid = data_loads[i, 2]
      P, Q = data_loads[i, 4:5] #MW, MVar
      busname = "BUS_$busid"

      data_loads[i, 6] == 0 || warn("Loads - $busname : parameter fault!=0 ($(data_loads[i, 6])), not handled")
      data_loads[i, 7] == 0 || warn("Loads - $busname : parameter curative!=0 ($(data_loads[i, 7])), not handled")

      bus_elems = bus["$busname"]
      nb = count([ismatch(r"Load_", x) for x in keys(bus_elems)])
      bus_elems["Load_$(nb+1)"] = IIDMLoad(string(nb+1), P+im*Q)
  end

  for i=1:size(data_generators, 1)
      busid = data_generators[i, 2]
      minP, maxP = data_generators[i, 5:6] # MW
      volt_regulator_on, targetV = data_generators[i, 11:12] # bool, pu
      minQmaxP, minQminP, maxQmaxP, maxQminP = data_generators[i, 7:10] # MVar
      S_bounds = [(minP+im*minQminP, minP+im*maxQminP), (maxP+im*minQmaxP, maxP+im*maxQmaxP)]
      busname = "BUS_$busid"

      data_generators[i, 2] == data_generators[i, 3] || warn("Generators - $busname : 'bus' is not 'con. bus', not handled")
      data_generators[i, 15] == 0 || warn("Loads - $busname : parameter fault!=0 ($(data_generators[i, 15])), not handled")
      data_generators[i, 16] == 0 || warn("Loads - $busname : parameter curative!=0 ($(data_generators[i, 16])), not handled")

      bus_elems = bus["$busname"]
      nb = count([ismatch(r"Gen_", x) for x in keys(bus_elems)])
      bus_elems["Gen_$(nb+1)"] = IIDMGenerator(string(nb+1), minP, maxP, targetV, volt_regulator_on, S_bounds)
  end

  link = Dict{Link, Dict{String, Any}}()
  for i=1:size(data_branches, 1)
      busid1, busid2 = data_branches[i, 2:3]
      r, x, g1, g2, b1, b2 = data_branches[i, 7:12] # pu
      p1, p2, q1, q2 = data_branches[i, 16:19] # MW, MW, MVar, MVar

      linkname = Link("BUS_$busid1", "BUS_$busid2")
      if !haskey(link, linkname)
          link[linkname] = Dict{String, Any}()
      end
      link_elems = link[linkname]

      basekV_dest = basekV_orig = lim = 0

      nb = count([ismatch(r"Branch_", x) for x in keys(link_elems)])
      link_elems["IIDMLine_π_$(nb+1)"] = IIDMLine_π(r, x, g1+im*b1, g2+im*b2, basekV_orig, basekV_dest, lim)
      link_elems["IIDMLine_πMeasure_$(nb+1)"] = IIDMLine_πMeasure(p1+im*q1, p2+im*q2)
  end

  ds = DataSource(bus, link)

  ## Building GridStructure
  node_linksin = node_linksout = Dict{String, Set{Link}}()
  node_vars = Dict{String, Dict{String, Variable}}()
  link_vars = Dict{Link, Dict{Type, Variable}}()
  gs = GridStructure("BaseCase", node_linksin, node_linksout)

  node_formulations = Dict{String, Dict{Tuple{Type, String}, Symbol}}()
  link_formulations = Dict{Link, Dict{Tuple{Type, String}, Symbol}}()
  mp = MathematicalProgramming(node_formulations, link_formulations, node_vars, link_vars)
  return OPFProblems("BaseCase"=>Scenario(ds, gs, mp))
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
### Following R. Nelsen 'An Introduction to Copulas', Springer Science & Business Media, 1999 - 216,
### for bivariate Archimedean copulas `C(u₁,u₂)` data can be generated as follow:
### * draw `u₁ = rand()`,
### * define `w = ∂C(u₁, u₂)\∂u₁` and inverse `u₂ = f(w, u₁)`,
### * draw  `w = rand()`
### * return a pair u₁, u₂.

### This method is applied here for Clayton, Frank and Ali-Mikhail-Haq copula. If we use
### this method recursively, we can get `n`-variate data with uniform marginals on
### `[0,1]`, where each neighbour pair
### of marginals `uᵢ uⱼ` for `j = i+1` are draw form a bivariate subcopula with
### parameter `θᵢ`, the only condition for `θᵢ`
### is such as for a corresponding bivariate copula.

"""
  rand2cop(u1::Real, θ::Real, copula::String, w::Real)

Returns Real  generated using copula::String given u1 and w from uniform distribution
and copula parameter θ.
"""
function rand2cop(u1::T, θ::T, copula::String, w::T) where T <: Real
  copula in ["clayton", "amh", "frank"] || throw(AssertionError("$(copula) copula is not supported"))
  if copula == "clayton"
    return (u1^(-θ)*(w^(-θ/(1+θ)) -1) +1)^(-1/θ)
  elseif copula == "frank"
    return -1/θ*log(1 +(w*(1-exp(-θ)))/(w*(exp(-θ*u1) -1) -exp(-θ*u1)))
  elseif copula == "amh"
    a = 1 -u1
    b = 1 -θ *(1 +2 *a *w)+2*θ^2*a^2 *w
    c = 1 -θ *(2 -4*w +4 *a *w)+θ^2 *(1 -4 *a *w+4 *a^2 *w)
    return 2*w *(a *θ -1)^2 /(b+sqrt(c))
  end
end
"""
    Chain_of_Archimedeans

Subsequent pairs of marginals are modeled by bi-variate copulas form the Archimedean
family, following copulas are supported: "Clayton", "Frank", "Ali-Mikhail-Haq",
"Gumbel" copula is not supported.

Fields:
  - n::Int - number of variables
  - θ::Vector{Real} - a vector of parameters
  - copulas::Vector{String} - a vector indicating bi-variate copulas.
  possible elements "clayton", "frank", "amh"

Constructor

    Chain_of_Archimedeans(θ::Vector{Real}, copulas::Vector{String})

requires length(θ) = length(copulas) and limitations on θ for particular bi-variate copulas

Constructor

    Chain_of_Archimedeans(θ::Vector{Real}, copulas::String)
one copula family for all subsequent pairs of marginals.

Constructors

    Chain_of_Archimedeans(θ::Vector{Real}, copulas::Vector{String}, cor::Type{<:CorrelationType})

    Chain_of_Archimedeans(θ::Vector{Real}, copulas::String, cor::Type{<:CorrelationType})

For computing copula parameters from expected correlation use empty type cor::Type{<:CorrelationType}.
cor ∈ {SpearmanCorrelation, KendallCorrelation} uses these correlations in the place of θ  in the constructor
to compute θ.

In all cases n = length(θ)+1.

```jldoctest

julia> c = Chain_of_Archimedeans([4., 11.], "frank")
Chain_of_Archimedeans(3, [4.0, 11.0], ["frank", "frank"])

julia> c = Chain_of_Archimedeans([.5, .7], ["frank", "clayton"], "Kendall")
Chain_of_Archimedeans(3, [5.736282707019972, 4.666666666666666], ["frank", "clayton"])

```
"""
struct Chain_of_Archimedeans{T}
  n::Int
  θ::Vector{T}
  copulas::Vector{String}
  function(::Type{Chain_of_Archimedeans})(θ::Vector{T}, copulas::Vector{String}) where T <: Real
      n = length(θ)+1
      n-1 == length(copulas) || throw(BoundsError("length(θ) ≠ length(copulas)"))
      map(i -> testbivθ(θ[i], copulas[i]), 1:n-1)
      for copula in copulas
        copula in ["clayton", "amh", "frank"] || throw(AssertionError("$(copula) copula is not supported"))
      end
      new{T}(n, θ, copulas)
  end
  function(::Type{Chain_of_Archimedeans})(θ::Vector{T}, copulas::Vector{String}, cor::Type{<:CorrelationType}) where T <: Real
      n = length(θ)+1
      n - 1 == length(copulas) || throw(BoundsError("length(θ) ≠ length(copulas)"))
      for copula in copulas
        copula in ["clayton", "amh", "frank"] || throw(AssertionError("$(copula) copula is not supported"))
      end
      θ = map(i -> usebivρ(θ[i], copulas[i], cor), 1:n-1)
      new{T}(n, θ, copulas)
  end
  function(::Type{Chain_of_Archimedeans})(θ::Vector{T}, copula::String) where T <: Real
      n = length(θ)+1
      map(i -> testbivθ(θ[i], copula), 1:n-1)
      copula in ["clayton", "amh", "frank"] || throw(AssertionError("$(copula) copula is not supported"))
      new{T}(n, θ, [copula for i in 1:n-1])
  end
  function(::Type{Chain_of_Archimedeans})(θ::Vector{T}, copula::String, cor::Type{<:CorrelationType}) where T <: Real
      n = length(θ)+1
      map(i -> testbivθ(θ[i], copula), 1:n-1)
      copula in ["clayton", "amh", "frank"] || throw(AssertionError("$(copula) copula is not supported"))
      θ = map(i -> usebivρ(θ[i], copula, cor), 1:n-1)
      new{T}(n, θ, [copula for i in 1:n-1])
  end
end

"""
    simulate_copula(t::Int, copula::Chain_of_Archimedeans; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations of multivariate data modeled by the chain of bivariate
Archimedean copulas, i.e.

    Chain_of_Archimedeans(θ::Vector{Flota64}, copulas::Union{String, Vector{String}})

```jldoctest
julia> Random.seed!(43);

julia> c = Chain_of_Archimedeans([4., 11.], "frank")
Chain_of_Archimedeans(3, [4.0, 11.0], ["frank", "frank"])

julia> simulate_copula(1, c)
1×3 Array{Float64,2}:
 0.180975  0.492923  0.679345

julia> c = Chain_of_Archimedeans([.5, .7], ["frank", "clayton"], KendallCorrelation)
Chain_of_Archimedeans(3, [5.736282707019972, 4.666666666666666], ["frank", "clayton"])

julia> Random.seed!(43);

julia> simulate_copula(1, c)
1×3 Array{Float64,2}:
 0.180975  0.408582  0.646887
```
"""
function simulate_copula(t::Int, copula::Chain_of_Archimedeans{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    U = zeros(T, t, copula.n)
    simulate_copula!(U, copula; rng = rng)
    U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Chain_of_Archimedeans; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, returns size(U,1) realizations of data modeled by the chain
of bivariate Archimedean copulas. N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> c = Chain_of_Archimedeans([4., 11.], "frank")
Chain_of_Archimedeans(3, [4.0, 11.0], ["frank", "frank"])

julia> u = zeros(1,3)
1×3 Array{Float64,2}:
 0.0  0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(u,c)

julia> u
1×3 Array{Float64,2}:
 0.180975  0.492923  0.679345

 julia> u = zeros(1,3)
1×3 Array{Float64,2}:
 0.0  0.0  0.0

julia> Random.seed!(43);

julia> c = Chain_of_Archimedeans([.5, .7], ["frank", "clayton"], KendallCorrelation)
Chain_of_Archimedeans(3, [5.736282707019972, 4.666666666666666], ["frank", "clayton"])

julia> simulate_copula!(u,c)

julia> u
1×3 Array{Float64,2}:
 0.180975  0.408582  0.646887
```
"""
function simulate_copula!(U::Matrix{T}, copula::Chain_of_Archimedeans{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    θ = copula.θ
    copulas = copula.copulas
    size(U, 2) == copula.n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    for j in 1:size(U,1)
      u = rand(rng, T)
      w = rand(rng, T, copula.n-1)
      for i in 1:copula.n-1
        u = hcat(u, rand2cop(u[i], θ[i], copulas[i], w[i]))
      end
      U[j,:] = u
    end
end

"""
  testbivθ(θ::Union{Real}, copula::String)

Tests bivariate copula parameter

clayton bivariate sub-copulas with parameters (θᵢ ≥ -1) ^ ∧ (θᵢ ≠ 0).
amh -- Ali-Mikhail-Haq bi-variate sub-copulas with parameters -1 ≥ θᵢ ≥ 1
Frank bi-variate sub-copulas with parameters (θᵢ ≠ 0)
"""
function testbivθ(θ::Union{T, Int}, copula::String) where T <: Real
  !(0. in θ)|(copula == "amh") || throw(DomainError("not supported for θ = 0"))
  if copula == "clayton"
    θ >= -1 || throw(DomainError("not supported for θ < -1"))
  elseif copula == "amh"
    -1 <= θ <= 1|| throw(DomainError("amh biv. copula supported only for -1 ≤ θ ≤ 1"))
  end
  Nothing
end

"""
    usebivρ(ρ::Real, copula::String, cor::Type{<:CorrelationType})

Returns Real, a copula parameter given the Spearman or the Kendall correlation

For Clayton or Frank copula the correlation must fulfill (-1 > ρᵢ > 1) ∧ (ρᵢ ≠ 0)
For the AMH copula Spearman must fulfill -0.2816 > ρᵢ >= .5, while Kendall -0.18 < τ < 1/3
"""
function usebivρ(ρ::T, copula::String, cor::Type{<:CorrelationType}) where T <: Real
  if copula == "amh"
      -0.2816 < ρ <= 0.5 || throw(DomainError("correlation coeficiant must fulfill -0.2816 < ρ <= 0.5"))
    if cor == KendallCorrelation
      -0.18 < ρ < 1/3 || throw(DomainError("correlation coeficiant must fulfill -0.2816 < ρ <= 1/3"))
    end
  else
    -1 < ρ < 1 || throw(DomainError("correlation coeficiant must fulfill -1 < ρ < 1"))
    !(0. in ρ) || throw(DomainError("not supported for ρ = 0"))
  end
  if cor == KendallCorrelation
      return τ2θ(ρ, copula)
  elseif cor == SpearmanCorrelation
      return ρ2θ(ρ, copula)
  end
end

# chain frechet copulas
"""
    Chain_of_Frechet

Chain of bi-variate Frechet copulas. Models each subsequent pair of marginals by the
bi-variate Frechet copula.

Fields:
    - n::Int - number of marginals
    - α::Vector{Real}  - vector of parameters for the maximal copula
    - β::Vector{Real} - vector of parameters for the minimal copula

Here α[i] and β[i] parameterized bi-variate Frechet copula between i th and i+1 th marginals.

Constructors

    Chain_of_Frechet(α::Vector{Real})
here β = zero(0)

    Chain_of_Frechet(α::Vector{Real}, β::Vector{Real})

```jldoctest
julia> Chain_of_Frechet([0.2, 0.3, 0.4])
Chain_of_Frechet(4, [0.2, 0.3, 0.4], [0.0, 0.0, 0.0])

julia> Chain_of_Frechet([0.2, 0.3, 0.4], [0.1, 0.1, 0.1])
Chain_of_Frechet(4, [0.2, 0.3, 0.4], [0.1, 0.1, 0.1])
```
"""
struct Chain_of_Frechet{T}
  n::Int
  α::Vector{T}
  β::Vector{T}
  function(::Type{Chain_of_Frechet})(α::Vector{T}) where T <: Real
      n = length(α)+1
      β = zero(α)
      minimum(α) >= 0 || throw(DomainError("negative α parameter"))
      maximum(α) <= 1 || throw(DomainError("α parameter greater than 1"))
      new{T}(n, α, β)
  end
  function(::Type{Chain_of_Frechet})(α::Vector{T}, β::Vector{T}) where T <: Real
      n = length(α)+1
      n == length(β) +1 || throw(AssertionError("length(α) ≠ length(β)"))
      minimum(α) >= 0 || throw(DomainError("negative α parameter"))
      minimum(β) >= 0 || throw(DomainError("negative β parameter"))
      maximum(α+β) <= 1 || throw(DomainError("α[i] + β[i] > 0"))
      new{T}(n, α, β)
  end
end

"""
    simulate_copula(t::Int, copula::Chain_of_Frechet; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations modeled by the chain of bivariate two parameter Frechet copulas

```jldoctest
julia> Random.seed!(43)

julia> simulate_copula(10, Chain_of_Frechet([0.6, 0.4], [0.3, 0.5]))
  10×3 Array{Float64,2}:
  0.996764  0.996764  0.996764
  0.204033  0.795967  0.204033
  0.979901  0.979901  0.0200985
  0.120669  0.879331  0.120669
  0.453027  0.453027  0.453027
  0.800909  0.199091  0.800909
  0.54892   0.54892   0.54892
  0.933832  0.933832  0.0661679
  0.396943  0.396943  0.396943
  0.804096  0.851275  0.955881
```
"""
function simulate_copula(t::Int, copula::Chain_of_Frechet{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
  α = copula.α
  β = copula.β
  n = copula.n
  fncopulagen(α, β, rand(rng, T, t, n))
end

"""
  fncopulagen(α::Vector{Real}, β::Vector{Real}, u::Matrix{Real})

```jldoctest

julia> fncopulagen(2, [0.2, 0.4], [0.1, 0.1], [0.2 0.4 0.6; 0.3 0.5 0.7])
2×3 Array{Float64,2}:
 0.6  0.4  0.2
 0.7  0.5  0.3

```
"""
function fncopulagen(α::Vector{T}, β::Vector{T}, u::Matrix{T}) where T <: Real
  p = invperm(sortperm(u[:,1]))
  u = u[:,end:-1:1]
  lx = floor.(Int, size(u,1).*α)
  li = floor.(Int, size(u,1).*β) + lx
  for j in 1:size(u, 2)-1
    u[p[1:lx[j]],j+1] = u[p[1:lx[j]], j]
    r = p[lx[j]+1:li[j]]
    u[r,j+1] = 1 .-u[r,j]
  end
  u
end


text:
# const codons = ["AAA", "AAC", "AAG", "AAT", "ACA", "ACC", "ACG", "ACT", "AGA", "AGC", "AGG", "AGT", "ATA", "ATC", "ATG", "ATT", "CAA", "CAC", "CAG", "CAT", "CCA", "CCC", "CCG", "CCT", "CGA", "CGC", "CGG", "CGT", "CTA", "CTC", "CTG", "CTT", "GAA", "GAC", "GAG", "GAT", "GCA", "GCC", "GCG", "GCT", "GGA", "GGC", "GGG", "GGT", "GTA", "GTC", "GTG", "GTT", "TAA", "TAC", "TAG", "TAT", "TCA", "TCC", "TCG", "TCT", "TGA", "TGC", "TGG", "TGT", "TTA", "TTC", "TTG", "TTT"]

const CodonQualityTriplet = Tuple{UInt8,UInt8,UInt8}
const CodonIndex = Tuple{Int,Int,Int}

codonqualitytriplet() = (0x00,0x00,0x00)::CodonQualityTriplet
codonindex() = (0,0,0)::CodonIndex

valid(ci::CodonIndex) = ci[1]>0 && ci[2]>0 && ci[3]>0
shiftin(t::Tuple{T,T,T},s::T) where {T} = (t[2],t[3],s)

function codon(ci::CodonIndex)
	n = ['A','C','G','T']
	string(ci[1]>0 ? n[ci[1]] : 'x',
		   ci[2]>0 ? n[ci[2]] : 'x',
		   ci[3]>0 ? n[ci[3]] : 'x')
end

const CodonQualityDict = Dict{CodonQualityTriplet,Int}

# dimensions: c, b, a, pos
# where abc is a codon
const CodonQualityMap = Array{CodonQualityDict,4}

# init Dict for each codon and position
codonqualitymap(len::Int) = CodonQualityDict[CodonQualityDict() for n1=1:4,n2=1:4,n3=1:4, i=1:len]

const CodonQualityAccumulator = Vector{CodonQualityMap}

codonqualityaccumulator(lengths::Vector{Int}) = CodonQualityMap[codonqualitymap(l) for l in lengths]
codonqualityaccumulator(len::Int) = codonqualityaccumulator([len])

# increase the count of a codon quality triplet
function inc!(m::CodonQualityMap,refPos::Int,codonInd::CodonIndex,
			  q::CodonQualityTriplet)
	d = m[codonInd[3],codonInd[2],codonInd[1],refPos] # the dict for this position and codon
	d[q] = get!(d,q,0)+1 # is there a better way to increase a value for a key?
end
function inc!(a::CodonQualityAccumulator,refInd::Int,refPos::Int,
			  codonInd::CodonIndex,q::CodonQualityTriplet)
	inc!(a[refInd],refPos,codonInd,q)
end

# qualityfilter! replaces all quality values below threshold with 0
#                and merges keys that are equal
function qualityfilter!(d::CodonQualityDict, threshold::Integer, scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[])
	isempty(d) && return

	resize!(scratch,0)

	# get all that should be filtered
	#for (k,v) in d
	for k in keys(d)
		if k[1]<threshold || k[2]<threshold || k[3]<threshold
			push!(scratch,k)
		end
	end

	# for each filtered entry
	for k in scratch
		v = pop!(d,k) # remove and get value

		# update quality triplet
		qt = (k[1]>=threshold ? k[1] : 0x00,
			  k[2]>=threshold ? k[2] : 0x00,
			  k[3]>=threshold ? k[3] : 0x00)

		if qt[1]>0 || qt[2]>0 || qt[3]>0
			d[qt] = get!(d,qt,0) + v
		end
	end
end
function qualityfilter!(m::CodonQualityMap, threshold::Integer, scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[])
	for d in m
		qualityfilter!(d,threshold,scratch)
	end
end
function qualityfilter!(a::CodonQualityAccumulator, threshold::Integer, scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[])
	for m in a
		qualityfilter!(m,threshold,scratch)
	end
end

iscomplete(qt::CodonQualityTriplet) = qt[1]>0 && qt[2]>0 && qt[3]>0
hascompletecodon(d::CodonQualityDict) = any(iscomplete,keys(d))
function nbrcompletecodons(d::CodonQualityDict)
	s = 0
	for (k,v) in d
		iscomplete(k) && (s+=v)
	end
	s
end

# function removeambiguous!(d::CodonQualityDict,complete::Array{Bool,3},scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[])
# end

# handles 4x4x4 matrix for a given position in the reference
function removeambiguous!(a::AbstractArray{T,3},scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[]) where {T}
	complete = map(hascompletecodon,a)

	for c1=1:4
		for c2=1:4
			for c3=1:4
				d = a[c3,c2,c1]

				resize!(scratch,0)
				for k in keys(d)

					# TODO: code is not type stable (Colon/Int for indexing), fix
					if !any(complete[ k[3]==0 ? (:) : c3,
						              k[2]==0 ? (:) : c2,
						              k[1]==0 ? (:) : c1 ])
						#println((c3,c2,c1))
						push!(scratch,k)
					end
				end

				for k in scratch
					delete!(d,k)
				end
			end
		end
	end
end

function removeambiguous!(m::CodonQualityMap,scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[])
	for p=1:size(m,4)
		#println(p)
		removeambiguous!(view(m,:,:,:,p),scratch)
	end
end
function removeambiguous!(a::CodonQualityAccumulator,scratch::Vector{CodonQualityTriplet}=CodonQualityTriplet[])
	for m in a
		removeambiguous!(m,scratch)
	end
end

function processread!(cqa::CodonQualityAccumulator, r::BamRead)
	# println(pos(r), " ", cigar(r))
	# println(seq(r))

	refInd = refID(r)+1 # 0 to 1-based index

	refPos = -1
	codonIndex = codonindex()
	qt = codonqualitytriplet()

	for c in r
		o = op(c)

		if o==CIGAR_M || o==CIGAR_X || o==CIGAR_Eq # the only chunks that are aligned to both read and reference

			# reset counter if the last chunk didn't end where the current start
			p = refpos(c)
			refPos+1 != p && (codonIndex=codonindex()) # invalidate codon
			refPos = p

			cSeq = seq(r,c)
			cQual = qual(r,c)

			for (i,n,q) in zip(1:len(c),cSeq,cQual)
				codonIndex = shiftin(codonIndex,base_index(n)) # anything but A=1,C=2,G=3,T=4 will be 0 and invalidate
				qt = shiftin(qt,q)

				if valid(codonIndex)
					# refPos-2 to get first nuc in codon
					inc!(cqa,refInd,refPos-2,codonIndex,qt)
				end

				refPos += 1
			end
		else
			# not a chunk that is aligned to both read and reference
			refPos = -1 # reset pos since we have to start over on a new codon
		end
	end
end

# strands is one of: :both, :forward, :reverse
function codonqualitycount(bamFile::BamFile; strands=:both, mappingQualityThreshold::Int=30, ignoreChimeric::Bool=true)
	@assert strands∈[:both, :forward, :reverse] "Strand must be one of :both, :forward, :reverse."

	seqs = sequences(bamFile)
	seqLengths = Int[s[2] for s in seqs]

	flagsUnset = BAM_FUNMAP | BAM_FSECONDARY | BAM_FDUP
	flagsSet = 0

	strands == :forward && (flagsUnset |= BAM_FREVERSE)
	strands == :reverse && (flagsSet   |= BAM_FREVERSE)

	cqa = codonqualityaccumulator(seqLengths)
	processbam!(bamFile,cqa, flagsSet=flagsSet, flagsUnset=flagsUnset,
		        mappingQualityThreshold=mappingQualityThreshold,
		        ignoreChimeric=ignoreChimeric)
	cqa
end

# test function to compare with older basecountextractor results
# TODO: remove
function savetxt(filename::String,cqa::CodonQualityAccumulator,refInd=1)
	os = open(filename,"w")
	m = cqa.list[refInd]

	for refPos=1:size(m.map,4)
		for c1=1:4
			for c2=1:4
				for c3=1:4
					# output dict
					d = m[refPos,(c1,c2,c3)]

					dKeys = collect(keys(d))
					dVals = collect(values(d))

					ind = sortperm(dKeys)
					dKeys = dKeys[ind]
					dVals = dVals[ind]

					for k in dKeys
						print(os, Int(k[1]), ' ', Int(k[2]), ' ', Int(k[3]), ' ')
					end
					println(os)
					for v in dVals
						print(os, v, ' ')
					end
					println(os)

					# wrong order
					# for k in keys(d)
					# 	print(os, Int(k[1]), ' ', Int(k[2]), ' ', Int(k[3]), ' ')
					# end
					# println(os)
					# for v in values(d)
					# 	print(os, v, ' ')
					# end
					# println(os)
				end
			end
		end
	end

	close(os)
end


text:
# using Revise, FLPDevelopment, BrowseTables

pathMac = "/Volumes/GoogleDrive/My Drive/Flipping/Datasets/Pharmacology/CaspaseDevelopment"
pathLin = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Pharmacology/CaspaseDevelopment"
if ispath(pathMac)
    path = pathMac
elseif ispath(pathLin)
    path = pathLin
else
    error("unknown path")
end
##
Cas_d, Cas_p, Cas_b, Cas_s = process_dataset(path)

for df in (Cas_p, Cas_b, Cas_s)
    df[!,:Virus] = [get(VirusDict,x,"Missing") for x in df.MouseID]
    df[!,:Sex] = [x in females ? "F" : "M" for x in df.MouseID]
    # df[!,:PreInterpoke] = [ismissing(x) ? 0.0 : x for x in df.PreInterpoke]
    transform!(df, :Virus => categorical, renamecols=false)
    gd = groupby(df,:Session)
    transform!(gd, :Streak => maximum => :Performance)
end

levels!(Cas_p.Virus,["tdTomato", "Caspase"])
levels!(Cas_b.Virus,["tdTomato", "Caspase"])
levels!(Cas_s.Virus,["tdTomato", "Caspase"])

for dd in [Cas_s, Cas_p]
    dd[!,:Gen] = [g == "HET" ? "Rbp4-cre" : "Wild Type" for g in dd.Gen]
    dd[!,:Combo] = dd.Gen .* "\n" .* string.(dd.Virus)
    dd[!,:Group] = [r.Gen == "Wild Type" ? "Wild Type" : r.Gen .* "\n" .* string(r.Virus) for r in eachrow(dd)]
end

Afreq = countmap(Cas_s.AfterLast)
Aprob = Dict()
for (a,f) in Afreq
    Aprob[a] = round(f/nrow(Cas_s),digits = 5)
end
Cas_s[!,:IncorrectStart] = [!x for x in Cas_s.CorrectStart]
Cas_s[!,:IncorrectLeave] = [!x for x in Cas_s.CorrectLeave]
Cas_s[!,:P_AfterLast] = [Aprob[a] for a in Cas_s.AfterLast]
gd = groupby(Cas_s,:Session)
transform!(gd, :AfterLast => frequency)
transform!(gd, :AfterLast => (x -> x .<= quantile(x,0.95)) => :Limit)
transform!(gd, :Num_Rewards => cumsum => :Cum_Rewards)
Cas_s[!,:RewRate] = Cas_s.Cum_Rewards ./ Cas_s.Stop
Cas_s[!,:TrialRewRate] = Cas_s.Cum_Rewards ./ Cas_s.Streak
gdp = groupby(Cas_p,[:Session, :Streak])
transform!(gdp,:PokeIn => (x -> x .- x[1]) => :In)
transform!(gdp,[:PokeIn, :PokeOut] => ((x,y) -> y .- x[1]) => :Out)
# for (df,name) in zip([Cas_p, Cas_b, Cas_s],["Pokes", "Bouts", "Streaks"])
#     filename = joinpath(path, "Results_" * string(today()),"FullInfo" * name * ".csv")
#     CSV.write(filename,df)
# end


text:
"""
A pointer to a `T` in some manually managed region of memory.
"""
# TODO do we want to also keep the page address/size in here? If we complicate the loading code a little we could avoid writing it to the page, so it would only exist on the stack.
struct Manual{T}
    ptr::Ptr{Void}

    function Manual{T}(ptr::Ptr{Void}) where {T}
        @assert isbits(T)
        new(ptr)
    end
end

function rewrite_address(expr)
    if !(expr isa Expr)
        esc(expr)
    elseif expr.head == :.
        (object, field) = expr.args
        if field isa QuoteNode
            fieldname = field.value
        elseif field isa Expr && field.head == :quote
            fieldname = field.args[1]
        else
            error("Impossible?")
        end
        :(get_address($(rewrite_address(object)), $(Val{fieldname})))
    elseif expr.head == :ref
        object = expr.args[1]
        :(get_address($(rewrite_address(object)), $(map(esc, expr.args[2:end])...)))
    elseif expr.head == :macrocall
        rewrite_address(macroexpand(expr))
    else
        error("Don't know how to compute address for $expr")
    end
end

"""
    @a man.x[2].y

Get a `Manual` pointing at the *address* of `man.x[2].y`.
"""
macro a(expr)
    rewrite_address(expr)
end

function rewrite_value(expr)
    if (expr isa Expr) && (expr.head == :(=))
        if length(expr.args) == 2
            :(unsafe_store!($(rewrite_address(expr.args[1])), $(esc(expr.args[2]))))
        else
            error("Don't know how to compute assignment $expr")
        end
    else
        :(unsafe_load($(rewrite_address(expr))))
    end
end

"""
    @v man.x[2].y

Get the *value* at `man.x[2].y`.

    @v man.x[2].y = 42

Set the *value* at `man.x[2].y`.

NOTE macros bind tightly, so:

    # invalid syntax
    @v man.x[2].y < 42

    # valid syntax
    (@v man.x[2].y) < 42
"""
macro v(expr)
    rewrite_value(expr)
end

@generated function get_address(man::Manual{T}, ::Type{Val{field}}) where {T, field}
    i = findfirst(fieldnames(T), field)
    @assert i != 0 "$T has no field $field"
    quote
        $(Expr(:meta, :inline))
        Manual{$(fieldtype(T, i))}(man.ptr + $(fieldoffset(T, i)))
    end
end

@generated function Base.unsafe_load(man::Manual{T}) where {T}
    if isempty(fieldnames(T))
        # is a primitive type
        quote
            $(Expr(:meta, :inline))
            unsafe_load(convert(Ptr{T}, man.ptr))
        end
    else
        # is a composite type - recursively load its fields
        # so that specializations of this method can hook in and alter loading
        quote
            $(Expr(:meta, :inline))
            $(Expr(:new, T, @splice (i, field) in enumerate(fieldnames(T)) quote
                unsafe_load(get_address(man, $(Val{field})))
            end))
        end
    end
end

@generated function Base.unsafe_store!(man::Manual{T}, value::T) where {T}
    if isempty(fieldnames(T))
        # is a primitive type
        quote
            $(Expr(:meta, :inline))
            unsafe_store!(convert(Ptr{T}, man.ptr), value)
            value
        end
    else
        # is a composite type - recursively store its fields
        # so that specializations of this method can hook in and alter storing
        quote
            $(Expr(:meta, :inline))
            $(@splice (i, field) in enumerate(fieldnames(T)) quote
                unsafe_store!(get_address(man, $(Val{field})), value.$field)
            end)
            value
        end
    end
end

# if the value is the wrong type, try to convert it (just like setting a field normally)
@inline function Base.unsafe_store!(man::Manual{T}, value) where {T}
    unsafe_store!(man, convert(T, value))
end

# pointers to other parts of the region need to be converted into offsets

@inline function Base.unsafe_load(man::Manual{Manual{T}}) where {T}
    offset = unsafe_load(Manual{UInt64}(man.ptr))
    Manual{T}(man.ptr + offset)
end

@inline function Base.unsafe_store!(man::Manual{Manual{T}}, value::Manual{T}) where {T}
    offset = value.ptr - man.ptr
    unsafe_store!(Manual{UInt64}(man.ptr), offset)
end


text:
struct RelationType <: AbstractRelationType
    label::Label
    is_root::Bool
end

# Remote functions

function create(x::RemoteConcept{RelationType})
    result = execute(x.transaction, RelationTypeRequestBuilder.create_req(x.concept.label))
    return Relation(result.relation_type_create_res.relation)
end


text:
function ImageSim(target,scope,n_img)
    # target  - 1xn array of SNRs for targets
    #            - SNR
    #            - velocity magnitude (pixels/image)
    # scope   - Telescope and Camera parameters
    #            - f (m)
    #            - D (m)
    #            - σ (pixels) sensor noise
    #            - pixel pitch (m)
    #            - λ (m)
    #            - 𝓂 (m1,m2)
    # n_img   - number of images to create in sequence

    # Define output dictionary
    ImgSeq = Dict()

    # Extract Target Parameters
    SNR = target.SNR
    vel = target.vel

    # Extract Telescope Parameters
    f  = scope.f
    D  = scope.D
    σ  = scope.σ
    p  = scope.p
    λ  = scope.λ
    m1 = scope.𝓂[1]  # Image size x
    m2 = scope.𝓂[2]  # Image size y
    σ_s = scope.σ_s

    # Derived Values
    n    = length(SNR)   # number of objects in image
    𝓂   = m1*m2         # total number of pixels
    μx   = rand(1:m1,n)  # random centroid of each object, x
    μy   = rand(1:m2,n)  # random centroid of each object, y
    W_m1 = 2*ceil(3*σ_s) # reasonable size of the tracking window
    W_m2 = W_m1          # reasonable size of the tracking window
    W𝓂  = W_m1*W_m2     # tracking window pixels
    s    = [snr*σ*√W𝓂 for snr in SNR]
    # TODO: Add loop for each μx and μy for case with multiple objects in single image
    vdir = ([m1/2 m2/2] - [μx μy])/norm([m1/2 m2/2] - [μx μy]) # Velocity direction towards center of frame
    V    = round.(vel.*vdir)

    # Noise - Normal Distribution
    n_pdf = Normal(0, σ)

    # Loop over Image Sequence
    for seq_i in 1:n_img

        # Image
        Img = Float64[] # Start with an empty array

        # Add Sensor Noise to Image
        for pixel in 1:𝓂
            push!(Img,rand(n_pdf))
        end

        # Reshape the Image into the correct size
        Img = reshape(Img, (m1,m2))

        # Create Signal Distributions - Multivariate Normal Distribution
        s_pdf = []
        if μx[1] > 0 || μx[1] <= m1
            if μy[1] > 0 || μy[1] <= m2
                for i in 1:n
                    push!(s_pdf, MvNormal([μx[i];μy[i]], σ_s^2*Matrix(I, 2, 2)))
                end
            end
        end

        # Photons - Estimate where each photon lands
        sᵢ = []
        if μx[1] > 0 || μx[1] <= m1
            if μy[1] > 0 || μy[1] <= m2
                for i in 1:n
                    push!(sᵢ,Int.(round.(clamp.(rand(s_pdf[i],Int(s[i])),1,m1))))
                    [Img[sᵢ[i][1,j],sᵢ[i][2,j]]+=1 for j in 1:size(sᵢ[i],2)]
                end
            end
        end

        #display(heatmap(Img, color = :greys, title = string("Velocity = ",V)))

        # Object truth information
        Objects = Dict(i => (μx[i],μy[i],s[i],s_pdf[i], V) for i in 1:n)

        # Update Object Location
        [μx[ii] += V[1] for ii in 1:n]
        [μy[ii] += V[2] for ii in 1:n]

        # Dictionary to return all images and truth locations
        push!(ImgSeq, seq_i => (Img , Objects))
    end

    return sort(OrderedDict(ImgSeq))
end


text:
"""
Example simulating the car racing environment.

Dependencies:
 - MPOPIS
 - Printf
 - Random
 - Plots
 - ProgressMeter
 - Dates
 - LinearAlgebra
 - Distributions
 - example_utils.jl
"""

"""
simulate_car_racing(; kwargs...)
    Simulate a car racing scenarion with 1 or multiple cars
kwargs:
 - num_trials = 1,                                   # Number of trials
 - num_steps = 200,                                  # Max number of steps per scenario
 - num_cars = 1,                                     # Number of cars
 - laps = 2,                                         # Max number of laps (if using curve.csv track)
 - policy_type = :cemppi,                            # Type of policy (see `get_policy` for options)
 - num_samples = 150,                                # Number of initial samples for the policy
 - horizon = 50,                                     # Time step horizon
 - λ = 10.0,                                         # Inverve temperature settings for IT weighting
 - α = 1.0,                                          # Control cost parameter
 - U₀ = zeros(Float64, num_cars*2),                  # Default initial contorl
 - cov_mat = block_diagm([0.0625, 0.1], num_cars),   # Control covariance matrix
 - ais_its = 10,                                     # Number of AIS algo iterations
 - λ_ais = 20.0,                                     # Inverse temperature for AIS algo (if applicable)
 - ce_elite_threshold = 0.8,                         # CE elite threshold (if applicable)
 - ce_Σ_est = :ss,                                   # CE Σ estimation methods (if applicable)
 - cma_σ = 0.75,                                     # CMA step factor (if applicable)
 - cma_elite_threshold = 0.8,                        # CMA elite threshold (if applicable)
 - state_x_sigma = 0.0,                              # Add normal noise std dev to x position at each time step (only for single car)
 - state_y_sigma = 0.0,                              # Add normal noise std dev to y position at each time step (only for single car)
 - state_ψ_sigma = 0.0,                              # Add normal noise std dev to heading at each time step (only for single car)
 - seed = Int(rand(1:10e10)),                        # Seed algorithm and envrionment (incremtented by trial number)
 - log_runs = true,                                  # Output results at each of each trial
 - plot_steps = false,                               # Plot each step (NOT RECOMMENDED FOR LARGE RUNS)
 - pol_log = false,                                  # Turn on policy logger (slows down the process)
 - plot_traj = false,                                # Plot the samples trajectories
 - plot_traj_perc = 1.0,                             # Percentage of samples trajectories to plot (if applicable)
 - text_with_plot = true,                            # Plot car 1 state info on plot
 - text_on_plot_xy = (80.0, -60.0)                   # XY position of output text (if applicable)
 - save_gif = false,                                 # Save gif
"""
function simulate_car_racing(;
    num_trials = 1,
    num_steps = 200,
    num_cars = 1,
    policy_type = :cemppi,
    laps = 2,
    num_samples = 150,
    horizon = 50,
    λ = 10.0,
    α = 1.0,
    U₀ = zeros(Float64, num_cars*2),
    cov_mat = block_diagm([0.0625, 0.1], num_cars),
    ais_its = 10,
    λ_ais = 20.0,
    ce_elite_threshold = 0.8,
    ce_Σ_est = :ss,
    cma_σ = 0.75,
    cma_elite_threshold = 0.8,
    state_x_sigma = 0.0,
    state_y_sigma = 0.0,
    state_ψ_sigma = 0.0,
    seed = Int(rand(1:10e10)),
    log_runs = true,
    plot_steps = false,
    pol_log = false,
    plot_traj = false,
    plot_traj_perc = 1.0,
    text_with_plot = true,
    text_on_plot_xy = (80.0, -60.0),
    save_gif = false,
)

    if num_cars > 1
        sim_type = :mcr
    else
        sim_type = :cr
    end

    @printf("\n")
    @printf("%-30s%s\n", "Sim Type:", sim_type)
    @printf("%-30s%d\n", "Num Cars:", num_cars)
    @printf("%-30s%d\n", "Num Trails:", num_trials)
    @printf("%-30s%d\n", "Num Steps:", num_steps)
    @printf("%-30s%d\n", "Max Num Laps:", laps)
    @printf("%-30s%s\n","Policy Type:", policy_type)
    @printf("%-30s%d\n", "Num samples", num_samples)
    @printf("%-30s%d\n", "Horizon", horizon)
    @printf("%-30s%.2f\n", "λ (inverse temp):", λ)
    @printf("%-30s%.2f\n", "α (control cost param):", α)
    if policy_type != :mppi && policy_type != :gmppi
        @printf("%-30s%d\n", "# AIS Iterations:", ais_its)
        if policy_type ∈ [:μΣaismppi, :μaismppi, :pmcmppi]
            @printf("%-30s%.2f\n", "λ_ais (ais inverse temp):", λ_ais)
        elseif policy_type == :cemppi
            @printf("%-30s%.2f\n", "CE Elite Threshold:", ce_elite_threshold)
            @printf("%-30s%s\n", "CE Σ Est Method:", ce_Σ_est)
        elseif policy_type == :cmamppi
            @printf("%-30s%.2f\n", "CMA Step Factor (σ):", cma_σ)
            @printf("%-30s%.2f\n", "CMA Elite Perc Thres:", cma_elite_threshold)
        end
    end
    @printf("%-30s[%.4f, ..., %.4f]\n", "U₀", U₀[1], U₀[end])
    @printf("%-30s%s([%.4f %.4f; %.4f %.4f], %d)\n", "Σ", "block_diagm",
        cov_mat[1,1], cov_mat[1,2], cov_mat[2,1], cov_mat[2,2], num_cars)
    if num_cars == 1
        @printf("%-30s%.4f\n", "Noise, State X σ:", state_x_sigma)
        @printf("%-30s%.4f\n", "Noise, State Y σ:", state_y_sigma)
        @printf("%-30s%.4f\n", "Noise, Heading σ:", state_ψ_sigma)
    end
    @printf("%-30s%d\n", "Seed:", seed)
    @printf("\n")

    # Must have policy log on if plotting trajectories
    if plot_traj
        pol_log = true
    end

    gif_name = "$sim_type-$num_cars-$policy_type-$num_samples-$horizon-$λ-$α-"
    if policy_type != :mppi && policy_type != :gmppi
        gif_name = gif_name * "$ais_its-"
    end
    if policy_type == :cemppi
        gif_name = gif_name * "$ce_elite_threshold-"
        gif_name = gif_name * "$ce_Σ_est-"
    elseif policy_type ∈ [:μΣaismppi, :μaismppi, :pmcmppi]
        gif_name = gif_name * "$λ_ais-"
    elseif policy_type == :cmamppi
        gif_name = gif_name * "$cma_σ-"
        gif_name = gif_name * "$cma_elite_threshold-"
    end
    gif_name = gif_name * "$num_trials-$laps.gif"
    anim = Animation()

    rews = zeros(Float64, num_trials)
    steps = zeros(Float64, num_trials)
    rews_per_step = zeros(Float64, num_trials)
    lap_ts = [zeros(Float64, num_trials) for _ in 1:laps]
    mean_vs = zeros(Float64, num_trials)
    max_vs = zeros(Float64, num_trials)
    mean_βs = zeros(Float64, num_trials)
    max_βs = zeros(Float64, num_trials)
    β_viols = zeros(Float64, num_trials)
    T_viols = zeros(Float64, num_trials)
    C_viols = zeros(Float64, num_trials)
    exec_times = zeros(Float64, num_trials)

    @printf("Trial    #: %12s : %7s: %12s", "Reward", "Steps", "Reward/Step")
    for ii ∈ 1:laps
        @printf(" : %6s%d", "lap ", ii)
    end
    @printf(" : %7s : %7s", "Mean V", "Max V")
    @printf(" : %7s : %7s", "Mean β", "Max β")
    @printf(" : %7s : %7s", "β Viol", "T Viol")
    if sim_type == :mcr
        @printf(" : %7s", "C Viol")
    end
    @printf(" : %7s", "Ex Time")
    @printf("\n")

    for k ∈ 1:num_trials

        if sim_type == :cr
            env = CarRacingEnv(rng=MersenneTwister())
        elseif sim_type == :mcr
            env = MultiCarRacingEnv(num_cars, rng=MersenneTwister())
        end

        pol = get_policy(
            policy_type,
            env,num_samples, horizon, λ, α, U₀, cov_mat, pol_log,
            ais_its,
            λ_ais,
            ce_elite_threshold, ce_Σ_est,
            cma_σ, cma_elite_threshold,
        )

        seed!(env, seed + k)
        seed!(pol, seed + k)

        pm = Progress(num_steps, 1, "Trial $k ....", 50)
        # Start timer
        time_start = Dates.now()

        lap_time = zeros(Int, laps)
        v_mean_log = Vector{Float64}()
        v_max_log = Vector{Float64}()
        β_mean_log = Vector{Float64}()
        β_max_log = Vector{Float64}()
        rew, cnt, lap, prev_y = 0, 0, 0, 0
        trk_viol, β_viol, crash_viol = 0, 0, 0

        # Main simulation loop
        while !env.done && cnt <= num_steps
            # Get action from policy
            act = pol(env)
            # Apply action to envrionment
            env(act)
            cnt += 1
            # Get reward at the step
            step_rew = reward(env)
            rew += step_rew

            # Plot or collect the plot for the animation
            if plot_steps || save_gif
                if plot_traj
                    p = plot(env, pol, plot_traj_perc, text_output=text_with_plot, text_xy=text_on_plot_xy)
                else
                    p = plot(env, text_output=text_with_plot, text_xy=text_on_plot_xy)
                end
                if save_gif frame(anim) end
                if plot_steps display(p) end
            end

            if sim_type == :cr
                env.state[1] += state_x_sigma * randn(env.rng)
                env.state[2] += state_y_sigma * randn(env.rng)

                δψ = state_ψ_sigma * randn(env.rng)
                env.state[3] += δψ

                # Passive rotation matrix
                rot_mat = [ cos(δψ) sin(δψ) ;
                           -sin(δψ) cos(δψ) ]
                V′ = rot_mat*[env.state[4]; env.state[5]]
                env.state[4:5] = V′
            end

            next!(pm)

            # Get logging information
            curr_y = env.state[2]
            if sim_type == :mcr
                curr_y = minimum([en.state[2] for en ∈ env.envs])
                vs = [norm(en.state[4:5]) for en ∈ env.envs]
                βs = [abs(calculate_β(en)) for en ∈ env.envs]
            else
                vs = norm(env.state[4:5])
                βs = abs(calculate_β(env))
            end
            push!(v_mean_log, mean(vs))
            push!(v_max_log, maximum(vs))
            push!(β_mean_log, mean(βs))
            push!(β_max_log, maximum(βs))

            # Determine if violations occurred
            if step_rew < -4000
                ex_β = exceed_β(env)
                within_t = sim_type == :cr ? within_track(env).within : within_track(env)
                if ex_β β_viol += 1 end
                if !within_t trk_viol += 1 end
                temp_rew = step_rew + ex_β*5000 + !within_t*1000000
                if temp_rew < -10500 crash_viol += 1 end
            end

            if sim_type == :mcr
                # Not exact, but works
                d = minimum([norm(en.state[1:2]) for en ∈ env.envs])
            else
                d = norm(env.state[1:2])
            end

            # Estimate to increment lap count on curve.csv
            if prev_y < 0.0 && curr_y >= 0.0 && d <= 15.0
                lap += 1
                lap_time[lap] = cnt
            end
            if lap >= laps || trk_viol > 10 || β_viol > 50
                env.done = true
            end
            prev_y = curr_y
        end

        # Stop timer
        time_end = Dates.now()
        seconds_ran = Dates.value(time_end - time_start) / 1000

        rews[k] = rew
        steps[k] = cnt-1
        rews_per_step[k] = rews[k]/steps[k]
        exec_times[k] = seconds_ran
        if sim_type ∈ (:cr, :mcr)
            for ii ∈ 1:laps
                lap_ts[ii][k] = lap_time[ii]
            end
            mean_vs[k] = mean(v_mean_log)
            max_vs[k] = maximum(v_max_log)
            mean_βs[k] = mean(β_mean_log)
            max_βs[k] = maximum(β_max_log)
            β_viols[k] = β_viol
            T_viols[k] = trk_viol
            C_viols[k] = crash_viol
        end

        # For clearing the progress bar
        if cnt > num_steps
            print("\u1b[1F") # Moves cursor to beginning of the line n lines up
            print("\u1b[0K") # Clears  part of the line. n=0: clear from cursor to end
        else
            print("\e[2K") # clear whole line
            print("\e[1G") # move cursor to column 1
        end
        if log_runs
            @printf("Trial %4d: %12.2f : %7d: %12.2f", k, rew, cnt-1, rew/(cnt-1))
            for ii ∈ 1:laps
                @printf(" : %7d", lap_time[ii])
            end
            @printf(" : %7.2f : %7.2f", mean(v_mean_log), maximum(v_max_log))
            @printf(" : %7.2f : %7.2f",  mean(β_mean_log), maximum(β_max_log))
            @printf(" : %7d : %7d", β_viol, trk_viol)
            if sim_type == :mcr
                @printf(" : %7d", crash_viol)
            end
            @printf(" : %7.2f", seconds_ran)
            @printf("\n")
        end
    end

    # Output summary results
    @printf("-----------------------------------\n")
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "AVE", mean(rews), mean(steps), mean(rews_per_step))
    for ii ∈ 1:laps
        @printf(" : %7.2f", mean(lap_ts[ii]))
    end
    @printf(" : %7.2f : %7.2f", mean(mean_vs), mean(max_vs))
    @printf(" : %7.2f : %7.2f",  mean(mean_βs), mean(max_βs))
    @printf(" : %7.2f : %7.2f", mean(β_viols), mean(T_viols))
    if sim_type == :mcr
        @printf(" : %7.2f", mean(C_viols))
    end
    @printf(" : %7.2f\n", mean(exec_times))
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "STD", std(rews), std(steps), std(rews_per_step))
    for ii ∈ 1:laps
        @printf(" : %7.2f", std(lap_ts[ii]))
    end
    @printf(" : %7.2f : %7.2f", std(mean_vs), std(max_vs))
    @printf(" : %7.2f : %7.2f",  std(mean_βs), std(max_βs))
    @printf(" : %7.2f : %7.2f", std(β_viols), std(T_viols))
    if sim_type == :mcr
        @printf(" : %7.2f", std(C_viols))
    end
    @printf(" : %7.2f\n", std(exec_times))
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "MED",
        quantile_ci(rews)[2], quantile_ci(steps)[2], quantile_ci(rews_per_step)[2])
    for ii ∈ 1:laps
        @printf(" : %7.2f", quantile_ci(lap_ts[ii])[2])
    end
    @printf(" : %7.2f : %7.2f", quantile_ci(mean_vs)[2], quantile_ci(max_vs)[2])
    @printf(" : %7.2f : %7.2f",  quantile_ci(mean_βs)[2], quantile_ci(max_βs)[2])
    @printf(" : %7.2f : %7.2f", quantile_ci(β_viols)[2], quantile_ci(T_viols)[2])
    if sim_type == :mcr
        @printf(" : %7.2f", quantile_ci(C_viols)[2])
    end
    @printf(" : %7.2f\n", quantile_ci(exec_times)[2])
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "L95",
        quantile_ci(rews)[1], quantile_ci(steps)[1], quantile_ci(rews_per_step)[1])
    for ii ∈ 1:laps
        @printf(" : %7.2f", quantile_ci(lap_ts[ii])[1])
    end
    @printf(" : %7.2f : %7.2f", quantile_ci(mean_vs)[1], quantile_ci(max_vs)[1])
    @printf(" : %7.2f : %7.2f",  quantile_ci(mean_βs)[1], quantile_ci(max_βs)[1])
    @printf(" : %7.2f : %7.2f", quantile_ci(β_viols)[1], quantile_ci(T_viols)[1])
    if sim_type == :mcr
        @printf(" : %7.2f", quantile_ci(C_viols)[1])
    end
    @printf(" : %7.2f\n", quantile_ci(exec_times)[1])
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "U95",
        quantile_ci(rews)[3], quantile_ci(steps)[3], quantile_ci(rews_per_step)[3])
    for ii ∈ 1:laps
        @printf(" : %7.2f", quantile_ci(lap_ts[ii])[3])
    end
    @printf(" : %7.2f : %7.2f", quantile_ci(mean_vs)[3], quantile_ci(max_vs)[3])
    @printf(" : %7.2f : %7.2f",  quantile_ci(mean_βs)[3], quantile_ci(max_βs)[3])
    @printf(" : %7.2f : %7.2f", quantile_ci(β_viols)[3], quantile_ci(T_viols)[3])
    if sim_type == :mcr
        @printf(" : %7.2f", quantile_ci(C_viols)[3])
    end
    @printf(" : %7.2f\n", quantile_ci(exec_times)[3])
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "MIN",
        minimum(rews), minimum(steps), minimum(rews_per_step))
    for ii ∈ 1:laps
        @printf(" : %7.2f", minimum(lap_ts[ii]))
    end
    @printf(" : %7.2f : %7.2f", minimum(mean_vs), minimum(max_vs))
    @printf(" : %7.2f : %7.2f",  minimum(mean_βs), minimum(max_βs))
    @printf(" : %7.2f : %7.2f", minimum(β_viols), minimum(T_viols))
    if sim_type == :mcr
        @printf(" : %7.2f", minimum(C_viols))
    end
    @printf(" : %7.2f\n", minimum(exec_times))
    @printf("Trials %3s: %12.2f : %7.2f: %12.2f", "MAX", maximum(rews), maximum(steps), maximum(rews_per_step))
    for ii ∈ 1:laps
        @printf(" : %7.2f", maximum(lap_ts[ii]))
    end
    @printf(" : %7.2f : %7.2f", maximum(mean_vs), maximum(max_vs))
    @printf(" : %7.2f : %7.2f",  maximum(mean_βs), maximum(max_βs))
    @printf(" : %7.2f : %7.2f", maximum(β_viols), maximum(T_viols))
    if sim_type == :mcr
        @printf(" : %7.2f", maximum(C_viols))
    end
    @printf(" : %7.2f\n", maximum(exec_times))

    if save_gif
        println("Saving gif...$gif_name")
        gif(anim, gif_name, fps=10)
    end
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
#used to generate the seat graphs.
using Gadfly
using Colors: RGB
using DataFrames

#sql connections to come later
function graph()
    Gadfly.with_theme(:dark) do
        #data = data()
        names = ["Libs", "NDP", "UCP", "RSP", "Ind"]
        seats = [5, 13, 5, 2, 1]
        colours = [RGB(184/255,44/255,44/255), RGB(255/255,129/255,0),
        RGB(65/255,88/255,169/255), RGB(236/255,45/255,74/255), RGB(102/255,102/255,102/255)]
        df = DataFrames.DataFrame(Parties = names, Seats = seats, Colours = colours)
        p = Gadfly.plot(df ,x=:Parties, y=:Seats, color=:Colours, Geom.bar,
        Guide.xlabel("Parties"), Guide.ylabel("Seat Counts"))
        img = SVG("seatchart.svg", 10cm, 14cm)
        draw(img, p)
    end
end

graph()

max_stars_repo_name: cmhoc/backbone

max_issues_repo_name: cmhoc/backbone

max_forks_repo_name: cmhoc/backbone


text:
"""
Generates test cases, which are sets of arguments to use for testing functions.
"""
module UnitTestDesign

export IPOG
export GND
export Excursion
export generate_tuples
export all_tuples
export all_values
export all_pairs
export all_triples
export values_excursion
export pairs_excursion
export triples_excursion
export full_factorial

include("combinations.jl")
include("coverage_matrix.jl")
include("coverage_set.jl")
include("greedy_tuples.jl")
include("parameter_order.jl")
include("full_factorial.jl")
include("excursions.jl")
include("factorial_interface.jl")

end


text:
struct GradientDescentIterable{T, F, S}
    w0::T
    f::F
    stepsize::S
end

GradientDescentIterable(w0, f; stepsize) = GradientDescentIterable(w0, f, to_iterator(stepsize))

Base.IteratorSize(::Type{<:GradientDescentIterable}) = Base.IsInfinite()

mutable struct GradientDescentState{T, R, G, S}
    w::T
    f_w::R
    grad_f_w::G
    stepsize_iterator::S
end

function Base.iterate(iter::GradientDescentIterable)
    w = copy(iter.w0)
    grad_f_w, f_w = gradient(iter.f, w)
    state = GradientDescentState(w, f_w, grad_f_w, Iterators.Stateful(iter.stepsize))
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

function Base.iterate(iter::GradientDescentIterable, state::GradientDescentState)
    stepsize = popfirst!(state.stepsize_iterator)
    state.w .-= stepsize .* state.grad_f_w
    state.grad_f_w, state.f_w = gradient(iter.f, state.w)
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

struct GradientDescent
    kwargs
    GradientDescent(; kwargs...) = new(kwargs)
end

(alg::GradientDescent)(args...) = GradientDescentIterable(args...; alg.kwargs...)


text:
using TrajectoryOptimization
using Random
using BenchmarkTools
Random.seed!(7)
### Solver options ###
dt = 0.1
opts = TrajectoryOptimization.SolverOptions()
opts.square_root = false
opts.verbose = false
opts.cache = false
opts.c1 = 1e-8
opts.c2 = 10.0
opts.cost_tolerance_intermediate = 1e-4
opts.constraint_tolerance = 1e-4
opts.cost_tolerance = 1e-4
opts.iterations_outerloop = 50
opts.iterations = 250
opts.iterations_linesearch = 25
opts.constraint_decrease_ratio = 0.25
opts.penalty_scaling = 10.0
opts.bp_reg_initial = 1.0
opts.outer_loop_update_type = :individual
opts.use_static = false
opts.resolve_feasible = false
opts.λ_second_order_update = false
opts.bp_reg_type = :state
######################

### Set up model, objective, solver ###
# Model, objective (unconstrained)
model_pendulum, obj_uncon_pendulum = TrajectoryOptimization.Dynamics.pendulum!
model_dubins, obj_uncon_dubins = TrajectoryOptimization.Dynamics.dubinscar!
model_cartpole, obj_uncon_cartpole = TrajectoryOptimization.Dynamics.cartpole_udp

## Constraints

# pendulum
u_min_pendulum = -2
u_max_pendulum = 2
x_min_pendulum = [-20;-20]
x_max_pendulum = [20; 20]

# dubins car
u_min_dubins = [-1; -1]
u_max_dubins = [1; 1]
x_min_dubins = [0; -100; -100]
x_max_dubins = [1.0; 100; 100]

# cartpole
u_min_cartpole = -20
u_max_cartpole = 40
x_min_cartpole = [-10; -1000; -1000; -1000]
x_max_cartpole = [10; 1000; 1000; 1000]

# -Constrained objective
obj_con_pendulum = ConstrainedObjective(obj_uncon_pendulum, u_min=u_min_pendulum, u_max=u_max_pendulum, x_min=x_min_pendulum, x_max=x_max_pendulum)
obj_con_dubins = ConstrainedObjective(obj_uncon_dubins, u_min=u_min_dubins, u_max=u_max_dubins, x_min=x_min_dubins, x_max=x_max_dubins)
obj_con_cartpole = ConstrainedObjective(obj_uncon_cartpole, u_min=u_min_cartpole, u_max=u_max_cartpole, x_min=x_min_cartpole, x_max=x_max_cartpole)

# Solver
intergrator = :rk3
solver_pendulum = Solver(model_pendulum,obj_con_pendulum,integration=intergrator,dt=dt,opts=opts)
solver_dubins = Solver(model_dubins,obj_con_dubins,integration=intergrator,dt=dt,opts=opts)
solver_cartpole = Solver(model_cartpole,obj_con_cartpole,integration=intergrator,dt=dt,opts=opts)

# -Initial state and control trajectories
X_interp_pendulum = line_trajectory(solver_pendulum.obj.x0,solver_pendulum.obj.xf,solver_pendulum.N)
X_interp_dubins = line_trajectory(solver_dubins.obj.x0,solver_dubins.obj.xf,solver_dubins.N)
X_interp_cartpole = line_trajectory(solver_cartpole.obj.x0,solver_cartpole.obj.xf,solver_cartpole.N)

U_pendulum = rand(solver_pendulum.model.m,solver_pendulum.N)
U_dubins = rand(solver_dubins.model.m,solver_dubins.N)
U_cartpole = rand(solver_cartpole.model.m,solver_cartpole.N)

#######################################

### Solve ###
@time results_pendulum, stats_pendulum = solve(solver_pendulum,U_pendulum)
@time results_dubins, stats_dubins = solve(solver_dubins,U_dubins)
@time results_cartpole, stats_cartpole = solve(solver_cartpole,U_cartpole)

function cost_constraints_thread(solver::Solver, res::ConstrainedIterResults)
    N = solver.N
    J = zeros(N)
    Threads.@threads for k = 1:N-1
        J[k] = 0.5*res.C[k]'*res.Iμ[k]*res.C[k] + res.LAMBDA[k]'*res.C[k]
    end

    if solver.control_integration == :foh
        J[N] = 0.5*res.C[N]'*res.Iμ[N]*res.C[N] + res.LAMBDA[N]'*res.C[N]
    end

    J[N] += 0.5*res.CN'*res.IμN*res.CN + res.λN'*res.CN

    return sum(J)
end

function _cost_thread(solver::Solver,res::SolverVectorResults,X=res.X,U=res.U)
    # pull out solver/objective values
    N = solver.N; Q = solver.obj.Q; xf::Vector{Float64} = solver.obj.xf; Qf::Matrix{Float64} = solver.obj.Qf; m = solver.model.m; n = solver.model.n
    obj = solver.obj
    dt = solver.dt

    R = getR(solver)

    J = zeros(N)
    Threads.@threads for k = 1:N-1
        if solver.control_integration == :foh
            Xm = res.xmid[k]
            Um = (U[k] + U[k+1])/2

            J[k] = solver.dt/6*(stage_cost(X[k],U[k],Q,R,xf) + 4*stage_cost(Xm,Um,Q,R,xf) + stage_cost(X[k+1],U[k+1],Q,R,xf)) # Simpson quadrature (integral approximation) for foh stage cost
        else
            J[k] = solver.dt*stage_cost(X[k],U[k],Q,R,xf)
        end
    end

    J[N] = 0.5*(X[N] - xf)'*Qf*(X[N] - xf)

    return sum(J)
end

# threaded functions
function calculate_midpoints_thread!(results::SolverVectorResults, solver::Solver, X::Vector, U::Vector)
    n,m,N = get_sizes(solver)
    Threads.@threads for k = 1:N-1
        results.xmid[k] = cubic_midpoint(results.X[k],results.xdot[k],results.X[k+1],results.xdot[k+1],solver.dt)
    end
end

function calculate_derivatives_thread!(results::SolverVectorResults, solver::Solver, X::Vector, U::Vector)
    n,m,N = get_sizes(solver)
    Threads.@threads for k = 1:N
        solver.fc(results.xdot[k],X[k],U[k][1:m])
    end
end

function calculate_jacobians_thread!(res::ConstrainedIterResults, solver::Solver)::Nothing #TODO change to inplace '!' notation throughout the code
    N = solver.N
    Threads.@threads for k = 1:N-1
        if solver.control_integration == :foh
            res.fx[k], res.fu[k], res.fv[k] = solver.Fd(res.X[k], res.U[k], res.U[k+1])
            res.Ac[k], res.Bc[k] = solver.Fc(res.X[k], res.U[k])
        else
            res.fx[k], res.fu[k] = solver.Fd(res.X[k], res.U[k])
        end
        solver.c_jacobian(res.Cx[k], res.Cu[k], res.X[k],res.U[k])
    end

    if solver.control_integration == :foh
        res.Ac[N], res.Bc[N] = solver.Fc(res.X[N], res.U[N])
        solver.c_jacobian(res.Cx[N], res.Cu[N], res.X[N],res.U[N])
    end

    solver.c_jacobian(res.Cx_N, res.X[N])
    return nothing
end

r1 = copy(results_pendulum)
r2 = copy(results_pendulum)
solver = solver_pendulum

println("Jacobians")
@btime update_jacobians!(r1,solver)
@btime calculate_jacobians_thread!(r2,solver)

println("Derivatives")
@btime calculate_derivatives!(r1,solver,r1.X,r1.U)
@btime calculate_derivatives_thread!(r2,solver,r2.X,r2.U)

println("Midpoints")
@btime calculate_midpoints!(r1,solver,r1.X,r1.U)
@btime calculate_midpoints_thread!(r2,solver,r2.X,r2.U)

println("Cost")
@btime _cost(solver,r1,r1.X,r1.U)
@btime _cost_thread(solver,r1,r1.X,r1.U)

println("Cost (constraints)")
@btime cost_constraints(solver,r1)
@btime cost_constraints_thread(solver,r1)

# profiler
using Profile
using Juno
@profiler results_pendulum, stats_pendulum = solve(solver_pendulum,U_pendulum)
#
Juno.profiletree()
Juno.profiler()

# @time results_pendulum, stats_pendulum = solve(solver_pendulum,U_pendulum)


text:
struct ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}
	H::Int
	h::T
	κ::Vector{T}
	q::Vector{Vector{T}}   # trajectory of q's   length=H+2
	u::Vector{Vector{T}}   # trajectory of u's   length=H
	w::Vector{Vector{T}}   # trajectory of w's   length=H
	γ::Vector{Vector{T}}   # trajectory of γ's   length=H
	b::Vector{Vector{T}}   # trajectory of b's   length=H
	z::Vector{Vector{T}}   # trajectory of z's   length=H
	θ::Vector{Vector{T}}   # trajectory of θ's   length=H
	iq0::SizedArray{Tuple{nq},Int,1,1,Vector{Int}}
	iq1::SizedArray{Tuple{nq},Int,1,1,Vector{Int}}
	iu1::SizedArray{Tuple{nu},Int,1,1,Vector{Int}}
	iw1::SizedArray{Tuple{nw},Int,1,1,Vector{Int}}
	iq2::SizedArray{Tuple{nq},Int,1,1,Vector{Int}}
	iγ1::SizedArray{Tuple{nc},Int,1,1,Vector{Int}}
	ib1::SizedArray{Tuple{nb},Int,1,1,Vector{Int}}
end

function contact_trajectory(model::Model, env::Environment, H::Int, h::T; κ::T=0.0) where {T}
	nq = model.nq
    nu = model.nu
    nw = model.nw
    nc = model.nc
    nb = nc * friction_dim(env)
	nz = num_var(model, env)
	nθ = num_data(model)

	q = [zeros(nq) for k=1:H+2]
	u = [zeros(nu) for k=1:H]
	w = [zeros(nw) for k=1:H]
	γ = [zeros(nc) for k=1:H]
	b = [zeros(nb) for k=1:H]
	z = [zeros(nz) for k=1:H]
	θ = [[zeros(nθ-1); copy(h)] for k=1:H]
	κ = [κ]
	off = 0
	iq0 = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q0
	iq1 = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q1
	iu1 = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iw1 = SizedVector{nw}(off .+ (1:nw)); off += nw # index of the disturbance w1
	off = 0
	iq2 = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iγ1 = SizedVector{nc}(off .+ (1:nc)); off += nc # index of the impact γ1
    ib1 = SizedVector{nb}(off .+ (1:nb)); off += nb # index of the linear friction b1

	return ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}(H,h,κ,q,u,w,γ,b,z,θ,iq0,iq1,iu1,iw1,iq2,iγ1,ib1)
end

function update_z!(traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}, t::Int) where {T,nq,nu,nw,nc,nb,nz,nθ}
	if t in (1:traj.H)
		traj.z[t][traj.iq2] = traj.q[t+2]
		traj.z[t][traj.iγ1] = traj.γ[t]
		traj.z[t][traj.ib1] = traj.b[t]
	end
	return nothing
end

function update_z!(traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}) where {T,nq,nu,nw,nc,nb,nz,nθ}
	for t in eachindex((1:traj.H))
		update_z!(traj, t)
	end
	return nothing
end

function update_θ!(traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}, t::Int) where {T,nq,nu,nw,nc,nb,nz,nθ}
	if t in (1:traj.H)
		traj.θ[t][traj.iq0] = traj.q[t]
		traj.θ[t][traj.iq1] = traj.q[t+1]
		traj.θ[t][traj.iu1] = traj.u[t]
		traj.θ[t][traj.iw1] = traj.w[t]
	end
	return nothing
end

function update_θ!(traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}) where {T,nq,nu,nw,nc,nb,nz,nθ}
	for t in eachindex((1:traj.H))
		update_θ!(traj, t)
	end
	return nothing
end

function repeat_ref_traj(traj::ContactTraj, N::Int;
       idx_shift = (1:0))

    shift = (traj.q[end] - traj.q[2])[idx_shift]

    q = [deepcopy(traj.q)...]
    u = [deepcopy(traj.u)...]
    w = [deepcopy(traj.w)...]
    γ = [deepcopy(traj.γ)...]
    b = [deepcopy(traj.b)...]
    z = [deepcopy(traj.z)...]
    θ = [deepcopy(traj.θ)...]
    κ = deepcopy(traj.κ)

    for i = 1:N-1
		for t = 1:traj.H
			qt = copy(traj.q[t + 2])
			qt[idx_shift] .+= i * shift
			push!(q, qt)
			push!(u, copy(traj.u[t]))
			push!(w, copy(traj.w[t]))
			push!(γ, copy(traj.γ[t]))
			push!(b, copy(traj.b[t]))
			push!(z, copy(traj.z[t]))
			push!(θ, copy(traj.θ[t]))
		end
    end

    return ContactTraj{typeof(traj).parameters...}(traj.H * N, traj.h, traj.κ,
        q, u, w, γ, b, z, θ,
		traj.iq0, traj.iq1, traj.iu1, traj.iw1, traj.iq2, traj.iγ1, traj.ib1)
end

function sub_traj(traj::ContactTraj, idx::AbstractVector{Int})

    q = deepcopy(traj.q[collect([[idx..., idx[end] + 1, idx[end] + 2]]...)])
    u = deepcopy(traj.u[idx])
    w = deepcopy(traj.w[idx])
    γ = deepcopy(traj.γ[idx])
    b = deepcopy(traj.b[idx])
    z = deepcopy(traj.z[idx])
    θ = deepcopy(traj.θ[idx])
    κ = deepcopy(traj.κ)

    return ContactTraj{typeof(traj).parameters...}(length(idx), traj.h, traj.κ,
        q, u, w, γ, b, z, θ,
		traj.iq0, traj.iq1, traj.iu1, traj.iw1, traj.iq2, traj.iγ1, traj.ib1)
end

function update_friction_coefficient!(traj::ContactTraj, model::Model, env::Environment)
	for t = 1:traj.H
		q2, γ1, b1, ψ1, s1, η1, __ = unpack_z(model, env, traj.z[t])
		q0, q1, u1, w1, _, h = unpack_θ(model, traj.θ[t])
		traj.z[t] .= pack_z(model, env, q2, γ1, b1, ψ1, η1)
		traj.θ[t] .= pack_θ(model, q0, q1, u1, w1, copy(model.μ_world), h)
	end
	nothing
end

function get_trajectory(model::Model, env::Environment, gait_path::String;
		load_type::Symbol = :split_traj, update_friction::Bool = false)
	#TODO: assert model exists

	nq = model.nq
	nu = model.nu
	nw = model.nw
	nc = model.nc
	nb = nc * friction_dim(env)
	res = JLD2.jldopen(gait_path)

	if load_type == :split_traj
		q, u, γ, b, h = res["q"], res["u"], res["γ"], res["b"], mean(res["h̄"])
		ū = res["ū"]
		ψ = [ut[nu + nc + nb .+ (1:nc)] for ut in ū]
		η = [ut[nu + nc + nb + nc .+ (1:nb)] for ut in ū]

		T = length(u)

		traj = contact_trajectory(T, h, model, env)
		traj.q .= deepcopy(q)
		traj.u .= deepcopy(u)
		traj.γ .= deepcopy(γ)
		traj.b .= deepcopy(b)
		traj.z .= [pack_z(model, env, q[t+2], γ[t], b[t], ψ[t], η[t]) for t = 1:T]
		traj.θ .= [pack_θ(model, q[t], q[t+1], u[t], zeros(nw), model.μ_world, h) for t = 1:T]
	elseif load_type == :split_traj_alt
		q, u, γ, b, ψ, η, μ, h = res["qm"], res["um"], res["γm"], res["bm"], res["ψm"], res["ηm"], res["μm"], res["hm"]
		T = length(u)

		traj = contact_trajectory(model, env, T, h)
		traj.q .= deepcopy(q)
		traj.u .= deepcopy(u)
		traj.γ .= deepcopy(γ)
		traj.b .= deepcopy(b)
		traj.z .= [pack_z(model, env, q[t+2], γ[t], b[t], ψ[t], η[t]) for t = 1:T]
		traj.θ .= [pack_θ(model, q[t], q[t+1], u[t], zeros(nw), μ, h) for t = 1:T]
	elseif load_type == :joint_traj
		traj = res["traj"]
	end
	update_friction && update_friction_coefficient!(traj, model, env)
	return traj
end

# Check tracking performance
function tracking_error(ref_traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ},
		sim_traj, N_sample::Int; idx_shift=(1:0)) where {T,nq,nu,nw,nc,nb,nz,nθ}

	# Horizons
	H_ref = ref_traj.H
	H_sim = length(sim_traj.u)
	H̄_sim = H_sim / N_sample
	dupl_traj = repeat_ref_traj(ref_traj, Int(ceil(H̄_sim / H_ref)), idx_shift=idx_shift)
	H_dupl = dupl_traj.H

	q_error = 0.0
	u_error = 0.0
	γ_error = 0.0
	b_error = 0.0
	cnt = 0
    for t = 1:H_dupl
		cnt += 1
		(t-1)*N_sample + 1 > H_sim && break
		q_error += norm(dupl_traj.q[t+2] - sim_traj.q[(t-1)*N_sample + 3], 1) / nq
		u_error += norm(dupl_traj.u[t] - sim_traj.u[(t-1)*N_sample + 1], 1) / nu
		γ_error += norm(dupl_traj.γ[t] - sim_traj.γ[(t-1)*N_sample + 1], 1) / nc
		b_error += norm(dupl_traj.b[t] - sim_traj.b[(t-1)*N_sample + 1], 1) / nb
    end

	q_error /= cnt
	u_error /= cnt
	γ_error /= cnt
	b_error /= cnt
    return q_error, u_error, γ_error, b_error
end

function initial_conditions(traj::ContactTraj)
	q1 = traj.q[2]
	v1 = (traj.q[2] - traj.q[1]) ./ traj.h

	return q1, v1
end


text:
@testset "load_gmsh_file" begin
    @test_throws ErrorException("File a does not exist.") load_gmsh_file("a")
    @test_throws ErrorException("Gmsh mesh file format saved in binary mode. Please convert to ASCII mode.") load_gmsh_file("data/gmsh/gmsh41_triangles_coarse_only_mesh_bin.msh")

    include("data/gmsh/gmsh41_triangles_coarse_only_mesh.jl")
    @test load_gmsh_file("data/gmsh/gmsh41_triangles_coarse_only_mesh.msh") == test_msh

    include("data/gmsh/gmsh41_triangles_coarse.jl")
    @test load_gmsh_file("data/gmsh/gmsh41_triangles_coarse.msh") == test_msh

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    @test load_gmsh_file("data/gmsh/gmsh41_tetrahedrons_coarse.msh") == test_msh
end

@testset "gmsh_do_physicalnames" begin
    include("data/gmsh/gmsh41_triangles_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_triangles_coarse.msh")
    @test gmsh_do_physicalnames(msh.raw_PhysicalNames) == test_physicals

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_tetrahedrons_coarse.msh")
    @test gmsh_do_physicalnames(msh.raw_PhysicalNames) == test_physicals
end

@testset "gmsh_do_entities" begin
    include("data/gmsh/gmsh41_triangles_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_triangles_coarse.msh")
    @test gmsh_do_entities(msh.raw_Entities) == test_entities

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_tetrahedrons_coarse.msh")
    @test gmsh_do_entities(msh.raw_Entities) == test_entities
end

@testset "gmsh_do_nodes" begin
    include("data/gmsh/gmsh41_triangles_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_triangles_coarse.msh")
    @test gmsh_do_nodes(msh.raw_Nodes) == test_nodes

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_tetrahedrons_coarse.msh")
    @test gmsh_do_nodes(msh.raw_Nodes) == test_nodes
end

@testset "gmsh_do_elements" begin
    include("data/gmsh/gmsh41_triangles_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_triangles_coarse.msh")
    @test gmsh_do_elements(msh.raw_Elements) == test_elements

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    msh = load_gmsh_file("data/gmsh/gmsh41_tetrahedrons_coarse.msh")
    @test gmsh_do_elements(msh.raw_Elements) == test_elements
end

@testset "load_gmsh" begin
    include("data/gmsh/gmsh41_triangles_coarse.jl")
    @test load_gmsh("data/gmsh/gmsh41_triangles_coarse.msh") == (physicalTagtoName=test_physicals,
         tagToPhysicalTag=test_entities,
         nodes=test_nodes,
         elements=test_elements)

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    @test load_gmsh("data/gmsh/gmsh41_tetrahedrons_coarse.msh") == (physicalTagtoName=test_physicals,
         tagToPhysicalTag=test_entities,
         nodes=test_nodes,
         elements=test_elements)
end

@testset "import_gmsh" begin
    include("data/gmsh/gmsh41_triangles_coarse.jl")
    ir20 = import_gmsh("data/gmsh/gmsh41_triangles_coarse.msh")

    @test ir20._v == test_ir20._v
    @test ir20.name == test_ir20.name

    @test ir20.left.name == test_ir20.left.name
    @test ir20.left.attributes == test_ir20.left.attributes
    @test ir20.left.nshapes == test_ir20.left.nshapes
    @test ir20.left.shapedesc == test_ir20.left.shapedesc

    @test ir20.right.name == test_ir20.right.name
    @test ir20.right.attributes == test_ir20.right.attributes
    @test ir20.right.nshapes == test_ir20.right.nshapes
    @test ir20.right.shapedesc == test_ir20.right.shapedesc

    include("data/gmsh/gmsh41_tetrahedrons_coarse.jl")
    ir30 = import_gmsh("data/gmsh/gmsh41_tetrahedrons_coarse.msh")

    @test ir30._v == test_ir30._v
    @test ir30.name == test_ir30.name

    @test ir30.left.name == test_ir30.left.name
    @test ir30.left.attributes == test_ir30.left.attributes
    @test ir30.left.nshapes == test_ir30.left.nshapes
    @test ir30.left.shapedesc == test_ir30.left.shapedesc

    @test ir30.right.name == test_ir30.right.name
    @test ir30.right.attributes == test_ir30.right.attributes
    @test ir30.right.nshapes == test_ir30.right.nshapes
    @test ir30.right.shapedesc == test_ir30.right.shapedesc
end