text:
"""
    getindex!(cache,a::AbstractArray,i...)

Returns the item of the array `a` associated with index `i`
by (possibly) using the scratch data passed in the `cache` object.

It defaults to

    getindex!(cache,a::AbstractArray,i...) = a[i...]

As for standard Julia arrays, the user needs to implement only one of the following signatures
depending on the `IndexStyle` of the array.

    getindex!(cache,a::AbstractArray,i::Integer)
    getindex!(cache,a::AbstractArray{T,N},i::Vararg{Integer,N}) where {T,N}

# Examples

Iterating over an array using the `getindex!` function

```jldoctest
using Gridap.Arrays

a = collect(10:15)

cache = array_cache(a)
for i in eachindex(a)
  ai = getindex!(cache,a,i)
  println("\$i -> \$ai")
end

# output
1 -> 10
2 -> 11
3 -> 12
4 -> 13
5 -> 14
6 -> 15
```

"""
@inline getindex!(cache,a::AbstractArray,i...) = a[i...]
@inline getindex!(cache,a::AbstractArray,i::CartesianIndex) = getindex!(cache,a,Tuple(i)...)
@inline getindex!(cache,a::AbstractArray,i::Integer) = _getindex_1d!(IndexStyle(a),cache,a,i)
@inline function getindex!(cache,a::AbstractArray{T,N},i::Vararg{Integer,N}) where {T,N}
  _getindex_nd!(IndexStyle(a),cache,a,CartesianIndex(i))
end
@inline _getindex_1d!(s::IndexLinear,cache,a,i) = a[i]
@inline _getindex_1d!(s::IndexCartesian,cache,a,i) = _getindex_nd!(s,cache,a,CartesianIndices(a)[i])
@inline _getindex_nd!(s::IndexLinear,cache,a,i) = _getindex_1d!(s,cache,a,LinearIndices(a)[i])
@inline _getindex_nd!(s::IndexCartesian,cache,a,i) = a[i]

"""
    array_cache(a::AbstractArray)

Returns a cache object to be used in the [`getindex!`](@ref) function.
It defaults to

    array_cache(a::T) where T = nothing

for types `T` such that `uses_hash(T) == Val(false)`, and

    function array_cache(a::T) where T
      hash = Dict{UInt,Any}()
      array_cache(hash,a)
    end

for types `T` such that `uses_hash(T) == Val(true)`, see the [`uses_hash`](@ref) function. In the later case, the
type `T` should implement the following signature:

    array_cache(hash::Dict,a::AbstractArray)

where we pass a dictionary (i.e., a hash table) in the first argument. This hash table can be used to test
if the object `a` has already built a cache and re-use it as follows

    id = objectid(a)
    if haskey(hash,id)
      cache = hash[id] # Reuse cache
    else
      cache = ... # Build a new cache depending on your needs
      hash[id] = cache # Register the cache in the hash table
    end

This mechanism is needed, e.g., to re-use intermediate results in complex lazy operation trees.
In multi-threading computations, a different hash table per thread has to be used in order
to avoid race conditions.
"""
array_cache(a::AbstractArray) = _default_array_cache(a,uses_hash(a))
array_cache(hash::Dict,a::AbstractArray) = _default_array_cache(hash,a,uses_hash(a))
_default_array_cache(a,s::Val{true}) = array_cache(Dict{UInt,Any}(),a)
_default_array_cache(a,s::Val{false}) = nothing
_default_array_cache(hash::Dict,a,s::Val{false}) = array_cache(a)
_default_array_cache(hash::Dict,a,s::Val{true}) = @abstractmethod

"""
    uses_hash(::Type{<:AbstractArray})

This function is used to specify if the type `T` uses the
hash-based mechanism to reuse caches.  It should return
either `Val(true)` or `Val(false)`. It defaults to

    uses_hash(::Type{<:AbstractArray}) = Val(false)

Once this function is defined for the type `T` it can also
be called on instances of `T`.

"""
uses_hash(::Type{<:AbstractArray}) = Val(false)
uses_hash(::T) where T = uses_hash(T)

"""
$(TYPEDSIGNATURES)

Returns an arbitrary instance of `eltype(a)`. The default returned value is the first entry
in the array if `length(a)>0` and `testvalue(eltype(a))` if `length(a)==0`
See the [`testvalue`](@ref) function.

# Examples

```jldoctest
using Gridap.Arrays

a = collect(3:10)
ai = testitem(a)

b = Int[]
bi = testitem(b)

(ai, bi)

# output
(3, 0)

```

"""
function testitem(a::AbstractArray{T}) where T
  #@check isconcretetype(T) "This array is type-instable"
  if length(a) >0
    first(a)
  else
    testvalue(T)
  end::T
end

function testitem(a::Fill)
  a.value
end

function testitem(a::Number)
  a
end

"""
    testvalue(::Type{T}) where T

Returns an arbitrary instance of type `T`. It defaults to `zero(T)` for
non-array types and to an empty array for array types.
It can be overloaded for new types `T` if `zero(T)` does not makes sense.
This function is used to compute  [`testitem`](@ref) for 0-length arrays.
"""
function testvalue end

@inline testvalue(::Type{T}) where T = zero(T)
@inline testvalue(v) = testvalue(typeof(v))

function testvalue(::Type{T}) where T<:AbstractArray{E,N} where {E,N}
   similar(T,tfill(0,Val(N))...)
end

function testvalue(::Type{T}) where T<:Transpose{E,A} where {E,A}
  a = testvalue(A)
  Transpose(a)
end

testvalue(::Type{Base.OneTo{T}}) where T = Base.OneTo(zero(T))

testvalue(::Type{Base.UnitRange{T}}) where T = UnitRange(one(T),zero(T))

function testvalue(::Type{T}) where T<:Fill{E,N,A} where {E,N,A}
  Fill(zero(E),testvalue(A))
end

function testvalue(::Type{<:Tuple})
  @notimplemented "testvalue on Tuple type only implemented up to 8 tuple elements"
end

#@fverdugo: use meta-programming here
function testvalue(::Type{Tuple{T1,T2,T3,T4,T5,T6,T7,T8}}) where {T1,T2,T3,T4,T5,T6,T7,T8}
  (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5),testvalue(T6),testvalue(T7),testvalue(T8))
end

function testvalue(::Type{Tuple{T1,T2,T3,T4,T5,T6,T7}}) where {T1,T2,T3,T4,T5,T6,T7}
  (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5),testvalue(T6),testvalue(T7))
end

function testvalue(::Type{Tuple{T1,T2,T3,T4,T5,T6}}) where {T1,T2,T3,T4,T5,T6}
  (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5),testvalue(T6))
end

function testvalue(::Type{Tuple{T1,T2,T3,T4,T5}}) where {T1,T2,T3,T4,T5}
  (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4),testvalue(T5))
end

function testvalue(::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4}
  (testvalue(T1),testvalue(T2),testvalue(T3),testvalue(T4))
end

function testvalue(::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3}
  (testvalue(T1),testvalue(T2),testvalue(T3))
end

function testvalue(::Type{Tuple{T1,T2}}) where {T1,T2}
  (testvalue(T1),testvalue(T2))
end

function testvalue(::Type{Tuple{T1}}) where {T1}
  (testvalue(T1),)
end

function testvalue(::Type{Tuple{}})
  ()
end

"""
    get_array(a::AbstractArray)

Returns `a`.
"""
function get_array(a::AbstractArray)
  a
end

# Test the interface

"""
    test_array(
      a::AbstractArray{T,N}, b::AbstractArray{S,N},cmp=(==)) where {T,S,N}

Checks if the entries in `a` and `b` are equal using the comparison function `cmp`.
It also stresses the new methods added to the `AbstractArray` interface.
"""
function test_array(
  a::AbstractArray{T,N}, b::AbstractArray{S,N},cmp=(==)) where {T,S,N}

  function _test_loop(indices)
    cache = array_cache(a)
    t = true
    for i in indices
      bi = b[i]
      ai = getindex!(cache,a,i)
      t = t && cmp(bi,ai)
    end
    @test t
  end

  @test cmp(a,b)
  _test_loop(eachindex(a))
  _test_loop(LinearIndices(a))
  _test_loop(CartesianIndices(a))
  cache = array_cache(a)
  t = true
  for i in eachindex(a)
    ai = getindex!(cache,a,i)
    t = t && (typeof(ai) <: eltype(a))
    t = t && (typeof(ai) <: T)
  end
  @test t
  @test isa(testitem(a),eltype(a))
  if length(a) > 0
    @test testitem(a) == first(a)
  end
  true
end


text:
import Pkg
import UUIDs

@inline function _strip(s::AbstractString)::String
    result::String = convert(String, strip(s))::String
    return result
end

@inline function _strip(v::AbstractVector{<:AbstractString})::Vector{String}
    result::Vector{String} = Vector{String}(undef, 0)
    for i = 1:length(v)
        v_i_stripped = _strip(v[i])::String
        if length(v_i_stripped) > 0
            push!(result, v_i_stripped)
        end
    end
    return result
end

@inline function is_stdlib(name::AbstractString)::Bool
    _stdlibs = Pkg.Types.stdlibs()
    return strip(name) in values(_stdlibs)
end

@inline function is_jll(name::AbstractString)::Bool
    return endswith(lowercase(strip(name)), "_jll")
end

@inline function with_temp_dir(f::Function)
    original_dir = pwd()
    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    cd(tmp_dir)

    result = f(tmp_dir)

    cd(original_dir)
    rm(tmp_dir; force = true, recursive = true)
    return result
end

ext: jl


text:
@testset "PLS2 Prediction Tests (in sample)" begin

	@testset "Single Column Prediction Test" begin

		X        = MLJ.table([1; 2; 3.0][:,:])
		Y        = MLJ.table([1 1; 2 2; 3.0 3.0])

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=1)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

		train = range(1,stop=length(X))
		MLJ.fit!(pls_machine, rows=train,force=true)
		pred = MLJ.predict(pls_machine, rows=train);
		pred = MLJ.matrix(pred)

		@test isequal(round.(pred),[1 1; 2 2; 3 3.0])

	end

	@testset "Linear Prediction Tests " begin

		X        = MLJ.table([1 2; 2 4; 4 6.0])
		Y        = MLJ.table([4 2;6 4;8 6.0])

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=train);
		pred = MLJ.matrix(pred)

		@test isequal(round.(pred),[4 2;6 4;8 6.0])

		X           = MLJ.table([1 -2; 2 -4; 4 -6.0])
		Y           = MLJ.table([-4 -2;-6 -4;-8 -6.0])

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=train);
		pred = MLJ.matrix(pred)

		@test isequal(round.(pred),[-4 -2;-6 -4;-8 -6.0])

	end

end

@testset "PLS2 Pediction Tests (out of sample)" begin

	@testset "Linear Prediction Tests (Ax + b) | A>0" begin

		X        = MLJ.table([1 2;2 4;3 6;6 12;7 14.0;4 8;5 10.0])
		Y        = MLJ.table([2 2;4 4;6 6;12 12;14 14.0;8 8;10 10.0])

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

		train = range(1,stop=5)
		test  = range(6,stop=7)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);
		pred = MLJ.matrix(pred)

		@test isequal(round.(pred),[8 8;10 10.0])

		X        = MLJ.table([1 2;2 4;3 6;6 12;7 14.0; 4 8;5 10.0])
		Y        = MLJ.table([2 4;4 6;6 8;12 14;14 16.0; 8 10;10 12.0])

    	pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

		train = range(1,stop=5)
		test  = range(6,stop=7)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);
		pred = MLJ.matrix(pred)

		@test isequal(round.(pred),[8 10;10 12.0])

	end

	@testset "Linear Prediction Tests (Ax + b) | A<0" begin

		X        = MLJ.table([1 -2;2 -4;3 -6;6 -12;7 -14.0; 4 -8;5 -10.0])
		Y        = MLJ.table([2 -2;4 -4;6 -6;12 -12;14 -14.0; 8 -8;10 -10.0])

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

		train = range(1,stop=5)
		test  = range(6,stop=7)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);
		pred = MLJ.matrix(pred)

	    @test isequal(round.(pred),[8 -8;10 -10.0])

		X        = MLJ.table([1 -2;2 -4;3 -6;6 -12;7 -14.0; 4 -8;5 -10.0])
		Y        = MLJ.table([2 -4;4 -6;6 -8;12 -14;14 -16.0; 8 -10;10 -12.0])

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

		train = range(1,stop=5)
		test  = range(6,stop=7)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);
		pred = MLJ.matrix(pred)

		@test isequal(round.(pred),[8 -10;10 -12.0])

	end

end;


text:
"""
    RandomVariable(m::Model, mean::Number, name::AbstractString)
Return a `RandomVariable` DataType for given flexibility model `m` given the `mean` and
the `name`. An anonymous JuMP variable is added directly to the flexibility model and its
is appended to `FlexibilityData.RVcols`.

**Arguments**
- `m::Model` The flexibility model.
- `mean::Number` The variable mean.
- `name::AbstractString` The variable name.
"""
function RandomVariable(m::Model, mean::Number, name::AbstractString)
    # Get the data and update it
    flex_data = getflexibilitydata(m)
    flex_data.numRVs += 1
    push!(flex_data.RVmeans, mean)
    push!(flex_data.RVnames, name)

    # Add an associated anonymous variable directly to the model and save its index
    @variable(m)
    index = length(m.colVal)
    push!(flex_data.RVcols, index)

    return RandomVariable(m, flex_data.numRVs)
end

"""
    getmean(variable::RandomVariable)
Return the mean corresponding to a particular `RandomVariable`. Currently this only accepts a single random
variable and vector variables are not accepted directly.

**Arguments**
- `variable::RandomVariable` The random variable name, must be a single variable.

```julia
julia> getmean(T[1])
620
```
"""
function getmean(variable::RandomVariable)
    flex_data = getflexibilitydata(variable.m)
    return flex_data.RVmeans[variable.idx]
end

"""
    RecourseVariable(m::Model, name::AbstractString)
Return a `RecourseVariable` DataType for given flexibility model `m` given the `name`.
An anonymous JuMP variable is added directly to the flexibility model and its is
appended to `FlexibilityData.recourse_cols`.

**Arguments**
- `m::Model` The flexibility model.
- `name::AbstractString` The variable name.
"""
function RecourseVariable(m::Model, name::AbstractString)
    # Get the data and update it
    flex_data = getflexibilitydata(m)
    flex_data.num_recourse_vars += 1
    push!(flex_data.recourse_names, name)

    # Add an associated anonymous variable directly to the model and save its index
    @variable(m)
    index = length(m.colVal)
    push!(flex_data.recourse_cols, index)

    return RecourseVariable(m, flex_data.num_recourse_vars)
end

"""
    JuMP.getvalue(v::FlexibilityVariable)
Return the value of the a flexibility variable this is an extension of `JuMP.getvalue`.

**Arguments**
- `v::FlexibilityVariable` The flexibility variable.

```julia
julia> getvalue(T)
4-element Array{Float64,1}:
 620.0
 388.0
 581.0
 319.0
```
"""
JuMP.getvalue(v::FlexibilityVariable) = getvalue(Variable(v.m, v.idx))

"""
    JuMP.linearindex(v::FlexibilityVariable)
Return the index of the a flexibility variable this is an extension of `JuMP.linearindex`.

**Arguments**
- `v::FlexibilityVariable` The flexibility variable, must be a single variable.

```julia
julia> linearindex(Qc)
1
```
"""
JuMP.linearindex(v::FlexibilityVariable) = v.idx


text:
ENV["GKSwstype"] = "nul"  # For StatsPlots

import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))
using CytofDensityEstimation
const CDE = CytofDensityEstimation
using LaTeXStrings
using Distributions
using StatsPlots
import Random
using RCall
using BSON
import CytofDensityEstimation.Model.default_ygrid
using StatsFuns
@rimport stats as rstats
VD = Vector{Dict{Symbol, Any}}

ygrid = collect(range(-10, 8, length=200))
default_plotsize = (400, 400)

# Plot true data density.
function plot_true_data_density(simdata, imgdir; lw=3,
                                plotsize=default_plotsize,
                                ygrid=ygrid, legendpos=:topleft)
  plot(ygrid, pdf.(simdata[:mmC], ygrid), lw=lw,
       label=L"y_C", color=:blue)
  plot!(ygrid, pdf.(simdata[:mmT], ygrid), lw=lw,
        label=L"y_T", color=:red)
  plot!(size=plotsize, legend=:topleft)
  p0C, p0T = simdata[:gammaC], simdata[:gammaT]
  title!("prop. -∞ in C: $(p0C) | prop. -∞ in T: $(p0T)", titlefont=font(10))
  savefig(joinpath(imgdir, "data-true-density.pdf"))
  closeall()
end

# Plot observed data histogram.
function plot_observed_hist(yC, yT, imgdir; bins, binsT=nothing,
                            alpha=0.3, digits=5, plotsize=default_plotsize,
                            legendpos=:topleft)
  binsC = bins
  binsT == nothing && (binsT = bins)
  histogram(yC[isfinite.(yC)], bins=binsC, label=L"y_C", alpha=alpha,
            color=:blue, linealpha=0, normalize=true)
  histogram!(yT[isfinite.(yT)], bins=binsT, label=L"y_T", alpha=alpha,
             color=:red, linealpha=0, normalize=true)
  p0C = round(mean(isinf.(yC)), digits=digits)
  p0T = round(mean(isinf.(yT)), digits=digits)
  title!("prop. -∞ in C: $(p0C) | prop. -∞ in T: $(p0T)", titlefont=font(10))
  plot!(size=plotsize, legend=:topleft)
  savefig(joinpath(imgdir, "data-hist.pdf"))
  closeall()
end

function plot_simdata_with_hist(simdata, imgdir; bins, binsT=nothing,
                                lw=1, alpha=0.3, digits=5,
                                plotsize=default_plotsize, legendpos=:topleft)
  plot(ygrid, pdf.(simdata[:mmC], ygrid), lw=lw,
       label=nothing, color=:blue)
  plot!(ygrid, pdf.(simdata[:mmT], ygrid), lw=lw,
        label=nothing, color=:red)

  binsC = bins
  binsT == nothing && (binsT = bins)
  yC = simdata[:yC]
  yT = simdata[:yT]
  histogram!(yC[isfinite.(yC)], bins=binsC, label=L"y_C", alpha=alpha,
             color=:blue, linealpha=0, normalize=true)
  histogram!(yT[isfinite.(yT)], bins=binsT, label=L"y_T", alpha=alpha,
             color=:red, linealpha=0, normalize=true)

  p0C = round(mean(isinf.(yC)), digits=digits)
  p0T = round(mean(isinf.(yT)), digits=digits)
  title!("prop. -∞ in C: $(p0C) | prop. -∞ in T: $(p0T)", titlefont=font(10))
  plot!(size=plotsize, legend=:topleft)

  savefig(joinpath(imgdir, "simdata.pdf"))
  closeall()
end

function postprocess(chain, laststate, summarystats, yC, yT, imgdir;
                     bw_postpred=0.2, density_legend_pos=:best,
                     simdata=nothing, ygrid=collect(range(-8, 8, length=1000)))
  # Print summary statistics.
  CDE.Model.printsummary(chain, summarystats)
  flush(stdout)

  # Print KS Statistic.
  ks_fit = rstats.ks_test(yC, yT)
  println("KS-test p-value: ", ks_fit["p.value"][1])
  flush(stdout)

  # Plot results.
  @time CDE.Model.plotpostsummary(chain, summarystats, yC, yT, imgdir,
                                  bw_postpred=bw_postpred, ygrid=ygrid,
                                  density_legend_pos=density_legend_pos,
                                  simdata=simdata)
  flush(stdout)
end

function postprocess(chain0, chain1, data, imgdir, awsbucket;
                     simdata=nothing, bw_postpred=0.20,
                     density_legend_pos=:best, binsC=nothing, binsT=nothing,
                     ygrid=collect(range(-8, 8, length=1000)), p=nothing)
  mkpath(imgdir)
  CDE.Util.redirect_stdout_to_file(joinpath(imgdir, "bf.txt")) do
    println("Start merging results...")

    # Compute Bayes factor in favor of M1.
    _p = CDE.Model.group(:p, chain0)[1]
    p == nothing || (_p = p)
    @time lbf = CDE.Model.log_bayes_factor(data, VD(chain0[1]), VD(chain1[1]))
    pm1 = CDE.Model.posterior_prob1(lbf, logpriorodds=logit(_p))
    println("Log Bayes Factor: $(lbf) | P(M1|y): $(pm1)")

    # Print KS Statistic.
    ks_fit = rstats.ks_test(data.yC, data.yT)
    println("KS-test p-value: ", ks_fit["p.value"][1])
    flush(stdout)

    # Plot marginal posterior density
    CDE.Model.plot_posterior_predictive(data.yC, data.yT, chain0, chain1, pm1,
                                        imgdir, bw_postpred=bw_postpred,
                                        density_legend_pos=density_legend_pos,
                                        simdata=simdata, lw=.5, ls=:solid,
                                        binsC=binsC, binsT=binsT,
                                        xlims_=(-6, 6), digits=5, fontsize=7)

    CDE.Model.plot_gamma(data.yC, data.yT, chain0, chain1, pm1, imgdir)

    # DIC
    dic0, dic1 = CDE.dic(chain0[1], data), CDE.dic(chain1[1], data)
    println("(DIC0, DIC1): ($(round(dic0, digits=3)), $(round(dic1, digits=3)))")

    dic_average = CDE.dic(chain0[1], chain1[1], pm1, data)
    println("DIC average: $(round(dic_average, digits=3))")

    flush(stdout)
  end

  # Send results
  if awsbucket != nothing
    CDE.Util.s3sync(from=imgdir, to=awsbucket, tags=`--exclude '*.nfs'`)
  end
end

function defaults(yC, yT, K; seed=nothing)
  seed == nothing || Random.seed!(seed)

  data = CDE.Model.Data(yC, yT)
  prior = CDE.Model.Prior(K, mu=Normal(0, 3),
                          nu=LogNormal(3, .5),
                          p=Beta(100, 100),
                          psi=Normal(-1, .5),
                          a_omega=2.5,
                          tau=Gamma(.5, 1.),
                          eta=Dirichlet(K, 1/K))
  state = CDE.Model.State(data, prior)
  tuners = CDE.Model.Tuners(K, 1.0)
  return (state=state, data=data, prior=prior, tuners=tuners)
end

function run(config)
  resultsdir = config[:resultsdir]
  mkpath(resultsdir)
  path_to_log = joinpath(resultsdir, "log.txt")
  CDE.Util.redirect_stdout_to_file(path_to_log) do
    _run(config)
  end
end

function _run(config)
  println("pid: $(getpid())"); flush(stdout)

  resultsdir = config[:resultsdir]
  awsbucket = config[:awsbucket]
  simdata = config[:simdata]
  yC = simdata[:yC]
  yT = simdata[:yT]
  K = config[:K]
  nsamps = config[:nsamps]
  nburn = config[:nburn]
  thin = config[:thin]

  # Make image path id needed.
  imgdir = joinpath(resultsdir, "img")
  mkpath(imgdir)

  # Create initial state, data, prior, and tuners objects.
  state, data, prior, tuners = defaults(yC, yT, K, seed=0)
  state.p = config[:p]
  state.beta = config[:beta]
  println("Priors:")
  foreach(fn -> println(fn, " => ", getfield(prior, fn)),
          fieldnames(CDE.Model.Prior))

  # Plot data.
  plot_observed_hist(yC, yT, imgdir, bins=100, binsT=100,
                     alpha=0.3, digits=5, legendpos=:topleft)
  plot_true_data_density(simdata, imgdir, lw=1, legendpos=:topleft)
  plot_simdata_with_hist(simdata, imgdir, bins=100, binsT=100)

  # Run analysis.
  println("Run Chain ..."); flush(stdout)
  fix=[:p, :beta, :zeta, :psi]
  state.zetaC .= 0
  state.zetaT .= 0
  state.psi .= 0
  @time chain, laststate, summarystats = CDE.fit(
      state, data, prior, tuners, nsamps=[nsamps], fix=fix,
      nburn=nburn, thin=thin, rep_aux=1, flags=CDE.Model.Flag[])
  flush(stdout)

  # Save results
  BSON.bson("$(resultsdir)/results.bson",
            Dict(:chain=>chain, :laststate=>laststate, :data=>data,
                 :prior=>prior, :summarystats=>summarystats,
                 :simdata=>simdata))

  # Load data:
  out = BSON.load("$(resultsdir)/results.bson")

  # Post process
  postprocess(out[:chain], out[:laststate], out[:summarystats], out[:data].yC,
              out[:data].yT, imgdir, bw_postpred=0.3, simdata=out[:simdata],
              ygrid=collect(range(-8, 8, length=200)),
              density_legend_pos=:topleft)

  # Send results
  if awsbucket != nothing
    CDE.Util.s3sync(from=resultsdir, to=awsbucket, tags=`--exclude '*.nfs'`)
  end

  # Print done.
  println("Done!"); flush(stdout)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using SimpleRoots, Test

using SimpleRoots: QuadraticError

@testset "test bracketing" begin
    @testset "sin 0.0 is 0.0" begin
        @test findzero(sin, (-1.0, 1.3); atol=1e-30)[1] ≈ 0.0 atol = 1e-30
        @test findzero(sin, (-1.0, 1.3); atol=1e-30)[2] == true
        @test findzero(sin, Secant(-1.0, 1.3); atol=1e-30)[1] ≈ 0.0 atol = 1e-30
        @test findzero(sin, Secant(-1.0, 1.3); atol=1e-30)[2] == true
        @test findzero(sin, Brent(-1.0, 1.3); atol=1e-30)[1] ≈ 0.0 atol = 1e-30
        @test findzero(sin, Brent(-1.0, 1.3); atol=1e-30)[2] == true
        @test findzero(sin, Bisection(-1.0, 1.3); atol=1e-100, max_iter=500)[1] ≈ 0.0 atol = 1e-30
        @test findzero(sin, Bisection(-1.0, 1.3); atol=1e-100, max_iter=500)[2] == true
    end

    @testset "sin with 1 iteration isn't found" begin
        @test findzero(sin, (-1.0, 1.3); atol=1e-30, max_iter=1)[2] == false
        @test findzero(sin, Secant(-1.0, 1.3); atol=1e-30, max_iter=1)[2] == false
        @test findzero(sin, Brent(-1.0, 1.3); atol=1e-30, max_iter=1)[2] == false
        @test findzero(sin, Bisection(-1.0, 1.3); atol=1e-100, max_iter=1)[2] == false
        # Except bisection when the bracket is symmetrical
        @test findzero(sin, Bisection(-1.0, -1.0); atol=1e-100, max_iter=1)[2] == true
    end

    @testset "sin(x) - x / 2" begin
        func1(x) = sin(x) - x / 2
        @test findzero(func1, [0.5pi, pi]; atol=1e-30)[1] == 1.895494267033981
        @test findzero(func1, Bisection([0.5pi, pi]); atol=1e-30)[1] == 1.895494267033981
        @test findzero(func1, Secant([0.5pi, pi]); atol=1e-30)[1] == 1.895494267033981
        @test findzero(func1, Brent([0.5pi, pi]); atol=1e-30)[1] == 1.895494267033981
    end

    @testset "inference" begin
        @inferred findzero(sin, (-1.0, 1.0); atol=1e-100, max_iter=100)
        @inferred findzero(sin, Secant(-1.0, 1.0); atol=1e-100, max_iter=100)
        @inferred findzero(sin, Brent(-1.0, 1.0); atol=1e-100, max_iter=100)
        @inferred findzero(sin, Bisection(-1.0, 1.0); atol=1e-100, max_iter=100)
    end
end

@testset "test quadratic polynomials" begin
    @test quad(1.0, 3.0, -4.0) == (-4.0, 1.0)
    @test quad(Both(), 1.0, 3.0, -4.0) == (-4.0, 1.0)
    @test quad(Lower(), 1.0, 3.0, -4.0) == -4.0
    @test quad(Upper(), 1.0, 3.0, -4.0) == 1.0

    @test quad(Both(), 0.0, 1.0, 1.0) == (-1.0, -1.0)
    @test quad(Lower(), 0.0, 1.0, 1.0) == -1.0
    @test quad(Upper(), 0.0, 1.0, 1.0) == -1.0

    @test quad(Both(), 1.0, 0.0, 0.0) == (0.0, 0.0)
    @test quad(Lower(), 1.0, 0.0, 0.0) == 0.0
    @test quad(Upper(), 1.0, 0.0, 0.0) == 0.0

    @test quad(Both(), 0.0, 0.0, 0.0) == (0.0, 0.0)
    @test quad(Lower(), 0.0, 0.0, 0.0) == 0.0
    @test quad(Upper(), 0.0, 0.0, 0.0) == 0.0

    @test_throws QuadraticError quad(1.0, 1.0, 1.0)
    @test_throws QuadraticError quad(0.0, 0.0, 1.0)

    @inferred quad(1.0, 3.0, -4.0)
end


text:
module Mamba

  using Reexport
  @reexport using Distributions

  #################### Imports ####################

  using DelimitedFiles
  using SpecialFunctions
  using Serialization
  using Distributed
  using Printf: @sprintf
  using LinearAlgebra
  using Calculus: gradient
  using Showoff: showoff

  import Base: Matrix, names
  import Compose: Context, context, cm, gridstack, inch, MeasureOrNumber, mm, pt, px
  import LinearAlgebra: cholesky, dot
  import Statistics: cor
  import Distributions:
         ## Generic Types
         Continuous, ContinuousUnivariateDistribution, Distribution,
         MatrixDistribution, Multivariate, MultivariateDistribution, PDiagMat,
         PDMat, ScalMat, Truncated, Univariate, UnivariateDistribution,
         ValueSupport,
         ## ContinuousUnivariateDistribution Types
         Arcsine, Beta, BetaPrime, Biweight, Cauchy, Chi, Chisq, Cosine,
         Epanechnikov, Erlang, Exponential, FDist, Frechet, Gamma, Gumbel,
         InverseGamma, InverseGaussian, Kolmogorov, KSDist, KSOneSided, Laplace,
         Levy, Logistic, LogNormal, NoncentralBeta, NoncentralChisq,
         NoncentralF, NoncentralT, Normal, NormalCanon, Pareto, Rayleigh,
         SymTriangularDist, TDist, TriangularDist, Triweight, Uniform, VonMises,
         Weibull,
         ## DiscreteUnivariateDistribution Types
         Bernoulli, Binomial, Categorical, DiscreteUniform, Geometric,
         Hypergeometric, NegativeBinomial, NoncentralHypergeometric, Pareto,
         PoissonBinomial, Skellam,
         ## MultivariateDistribution Types
         Dirichlet, Multinomial, MvNormal, MvNormalCanon, MvTDist,
         VonMisesFisher,
         ## MatrixDistribution Types
         InverseWishart, Wishart,
         ## Methods
         cdf, dim, gradlogpdf, insupport, isprobvec, logpdf, logpdf!, maximum,
         minimum, pdf, quantile, rand, sample!, support
  import Gadfly: draw, Geom, Guide, Layer, layer, PDF, PGF, Plot, plot, PNG, PS,
         render, Scale, SVG, Theme
  using LightGraphs: DiGraph, add_edge!, outneighbors,
         topological_sort_by_dfs, vertices
  import StatsBase: autocor, autocov, countmap, counts, describe, predict,
         quantile, sample, sem, summarystats

  include("distributions/pdmats2.jl")
  using .PDMats2

  #################### Types ####################

  ElementOrVector{T} = Union{T, Vector{T}}

  #################### Variate Types ####################

  abstract type ScalarVariate <: Real end
  abstract type ArrayVariate{N} <: DenseArray{Float64, N} end

  const AbstractVariate = Union{ScalarVariate, ArrayVariate}
  const VectorVariate = ArrayVariate{1}
  const MatrixVariate = ArrayVariate{2}

  #################### Distribution Types ####################

  const DistributionStruct = Union{Distribution,
                                   Array{UnivariateDistribution},
                                   Array{MultivariateDistribution}}

  #################### Dependent Types ####################

  mutable struct ScalarLogical <: ScalarVariate
    value::Float64
    symbol::Symbol
    monitor::Vector{Int}
    eval::Function
    sources::Vector{Symbol}
    targets::Vector{Symbol}
  end

  mutable struct ArrayLogical{N} <: ArrayVariate{N}
    value::Array{Float64, N}
    symbol::Symbol
    monitor::Vector{Int}
    eval::Function
    sources::Vector{Symbol}
    targets::Vector{Symbol}
  end

  mutable struct ScalarStochastic <: ScalarVariate
    value::Float64
    symbol::Symbol
    monitor::Vector{Int}
    eval::Function
    sources::Vector{Symbol}
    targets::Vector{Symbol}
    distr::UnivariateDistribution
  end

  mutable struct ArrayStochastic{N} <: ArrayVariate{N}
    value::Array{Float64, N}
    symbol::Symbol
    monitor::Vector{Int}
    eval::Function
    sources::Vector{Symbol}
    targets::Vector{Symbol}
    distr::DistributionStruct
  end

  const AbstractLogical = Union{ScalarLogical, ArrayLogical}
  const AbstractStochastic = Union{ScalarStochastic, ArrayStochastic}
  const AbstractDependent = Union{AbstractLogical, AbstractStochastic}

  #################### Sampler Types ####################

  mutable struct Sampler{T}
    params::Vector{Symbol}
    eval::Function
    tune::T
    targets::Vector{Symbol}
  end

  abstract type SamplerTune end

  struct SamplerVariate{T<:SamplerTune} <: VectorVariate
    value::Vector{Float64}
    tune::T

    function SamplerVariate{T}(x::AbstractVector, tune::T) where T<:SamplerTune
      v = new{T}(x, tune)
      validate(v)
    end

    function SamplerVariate{T}(x::AbstractVector, pargs...; kargs...) where T<:SamplerTune
      value = convert(Vector{Float64}, x)
      new{T}(value, T(value, pargs...; kargs...))
    end
  end

  #################### Model Types ####################

  struct ModelGraph
    graph::DiGraph
    keys::Vector{Symbol}
  end

  struct ModelState
    value::Vector{Float64}
    tune::Vector{Any}
  end

  mutable struct Model
    nodes::Dict{Symbol, Any}
    samplers::Vector{Sampler}
    states::Vector{ModelState}
    iter::Int
    burnin::Int
    hasinputs::Bool
    hasinits::Bool
  end

  #################### Chains Type ####################

  abstract type AbstractChains end

  struct Chains <: AbstractChains
    value::Array{Float64, 3}
    range::StepRange{Int, Int}
    names::Vector{AbstractString}
    chains::Vector{Int}
  end

  struct ModelChains <: AbstractChains
    value::Array{Float64, 3}
    range::StepRange{Int, Int}
    names::Vector{AbstractString}
    chains::Vector{Int}
    model::Model
  end

  #################### Includes ####################

  include("progress.jl")
  include("utils.jl")
  include("variate.jl")

  include("distributions/constructors.jl")
  include("distributions/distributionstruct.jl")
  include("distributions/extensions.jl")
  include("distributions/pdmatdistribution.jl")
  include("distributions/transformdistribution.jl")

  include("model/dependent.jl")
  include("model/graph.jl")
  include("model/initialization.jl")
  include("model/mcmc.jl")
  include("model/model.jl")
  include("model/simulation.jl")

  include("output/chains.jl")
  include("output/chainsummary.jl")
  include("output/discretediag.jl")
  include("output/fileio.jl")
  include("output/gelmandiag.jl")
  include("output/gewekediag.jl")
  include("output/heideldiag.jl")
  include("output/mcse.jl")
  include("output/modelchains.jl")
  include("output/modelstats.jl")
  include("output/rafterydiag.jl")
  include("output/stats.jl")
  include("output/plot.jl")

  include("samplers/sampler.jl")

  include("samplers/abc.jl")
  include("samplers/amm.jl")
  include("samplers/amwg.jl")
  include("samplers/bhmc.jl")
  include("samplers/bia.jl")
  include("samplers/bmc3.jl")
  include("samplers/bmg.jl")
  include("samplers/dgs.jl")
  include("samplers/hmc.jl")
  include("samplers/mala.jl")
  include("samplers/miss.jl")
  include("samplers/nuts.jl")
  include("samplers/rwm.jl")
  include("samplers/slice.jl")
  include("samplers/slicesimplex.jl")

  #################### Exports ####################

  export
    AbstractChains,
    AbstractDependent,
    AbstractLogical,
    AbstractStochastic,
    AbstractVariate,
    ArrayLogical,
    ArrayStochastic,
    ArrayVariate,
    Chains,
    Logical,
    MatrixVariate,
    Model,
    ModelChains,
    Sampler,
    SamplerTune,
    SamplerVariate,
    ScalarLogical,
    ScalarStochastic,
    ScalarVariate,
    Stochastic,
    VectorVariate

  export
    BDiagNormal,
    Flat,
    SymUniform

  export
    autocor,
    changerate,
    cor,
    describe,
    dic,
    discretediag,
    draw,
    gelmandiag,
    gettune,
    gewekediag,
    gradlogpdf,
    gradlogpdf!,
    graph,
    graph2dot,
    heideldiag,
    hpd,
    insupport,
    invlink,
    invlogit,
    link,
    logit,
    logpdf,
    logpdf!,
    mcmc,
    mcse,
    plot,
    predict,
    quantile,
    rafterydiag,
    rand,
    readcoda,
    relist,
    relist!,
    sample!,
    setinits!,
    setinputs!,
    setmonitor!,
    setsamplers!,
    summarystats,
    unlist,
    update!

  @static if VERSION >= v"1.0.0"
    export showall
  end

  export
    ABC,
    AMM, AMMVariate,
    AMWG, AMWGVariate,
    BHMC, BHMCVariate,
    BMC3, BMC3Variate,
    BMG, BMGVariate,
    DiscreteVariate,
    DGS, DGSVariate,
    HMC, HMCVariate,
    BIA, BIAVariate,
    MALA, MALAVariate,
    MISS,
    NUTS, NUTSVariate,
    RWM, RWMVariate,
    Slice, SliceMultivariate, SliceUnivariate,
    SliceSimplex, SliceSimplexVariate

  export
    cm,
    inch,
    mm,
    pt,
    px

  #################### Deprecated ####################

  include("deprecated.jl")

end


text:
## GPU runtime library

# reset the runtime cache from global scope, so that any change triggers recompilation
GPUCompiler.reset_runtime()

function signal_exception()
    return
end

function report_exception(ex)
    @veprintf("""
        ERROR: a %s was thrown during kernel execution.
        Run Julia on debug level 2 for device stack traces.
        """, ex)
    return
end

function report_oom(sz)
    @veprintf("ERROR: Out of dynamic VE memory (trying to allocate %i bytes)\n", sz)
    return
end

function report_exception_name(ex)
    @veprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex)
    return
end

function report_exception_frame(idx, func, file, line)
    @veprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    return
end


text:
using BifurcationInference
using ForwardDiff

using LaTeXStrings,Plots
using Plots.Measures

####################################################################################################
####################################################################################################
######################################################################################## saddle-node
include("minimal/saddle-node.jl"); X.roots .= [[[-1.0]]]
plot(F,[5/2,-1],X)
savefig("docs/figures/saddle-node.pdf")

####################################################################################################
####################################################################################################
########################################################################################## pitchfork
include("minimal/pitchfork.jl")
plot(F,[1/2,-1],X)
savefig("docs/figures/pitchfork.pdf")

####################################################################################################
####################################################################################################
################################################################################ bifurcation measure
measure(x) = 1/(1+abs(f(x)/ForwardDiff.derivative(f,x)))

layout = @layout [a;b{1.0w,0.5h}]
default(); default(grid=false,label="",margin=0mm)
plot(layout = layout, link = :x, size=(400,500) )

hline!([0],subplot=2,linewidth=1,ylim=(-5,5),color=:black, xlabel=L"\mathrm{arclength,}s", xmirror=true, topmargin=-5mm,
	ylabel=L"\mathrm{determinant}\,\quad\left|\!\!\!\!\frac{\partial F_{\theta}}{\partial u}\right|")

plot!([-1.52,-1.52],[0,1],subplot=1,linewidth=2,color=:gray)
plot!([-1.52,-1.52],[0,5],subplot=2,linewidth=2,color=:gray)

plot!([1.33,1.33],[0,1],subplot=1,linewidth=2,color=:gray)
plot!([1.33,1.33],[0,5],subplot=2,linewidth=2,color=:gray)

plot!([1.67,1.67],[0,1],subplot=1,linewidth=2,color=:gray)
plot!([1.67,1.67],[0,5],subplot=2,linewidth=2,color=:gray)

p = -2:0.01:2
f(p) = - 2 - (p-1/2)^3 + 3.1(p-1/2)
plot!( p, measure, subplot=1, linewidth=2, color=:red, ylim=(0,1) )
plot!( p, f, subplot=2, linewidth=2, color=:red )

f(p) = - 2 - p^3 + p
plot!( p, measure, subplot=1, linewidth=2, color=:gold, ylabel=L"\mathrm{Measure}\,\,\varphi_{\theta}(s)", ylim=(0,1) )
plot!( p, f, subplot=2, linewidth=2, color=:gold )

yticks!([0,1],subplot=1)
yticks!([-2,0,2],subplot=2)
xticks!([NaN],subplot=1)

scatter!([1.67],[0],subplot=2,marker=:star,color=:black,markersize=7)
scatter!([1.33],[0],subplot=2,marker=:star,color=:black,markersize=7)
scatter!([-1.52],[0],subplot=2,marker=:circle,color=:black,markersize=4)

savefig("docs/figures/bifurcation-measure.pdf")

####################################################################################################
####################################################################################################
############################################################################################ scaling
include("scaling/scaling.jl")
using StatsBase: mean,std
using LaTeXStrings

continuation_mean = Float64[]
continuation_std = Float64[]

forward_mean = Float64[]
forward_std = Float64[]

backward_mean = Float64[]
backward_std = Float64[]

for n ∈	41:100

	ts = [scaling_continuation(n,3) for _ ∈ 1:10]
	push!(continuation_mean,mean(ts))
	push!(continuation_std,std(ts))

	ts = [scaling_forward(n,3) for _ ∈ 1:10]
	push!(forward_mean,mean(ts))
	push!(forward_std,std(ts))

	# ts = [scaling_backward(n,3) for _ ∈ 1:10]
	# push!(backward_mean,mean(ts))
	# push!(backward_std,std(ts))
end

default();default(grid=false,label="", yscale=:log10, xscale=:log10,ms=3, ylim=(8e-2,1e1), msc=:auto, lc=:black,
	xlabel=L"\mathrm{States}\,\,N",ylabel=L"\mathrm{Execution\,\,Time}\,\,/\,\,\mathrm{sec}", size=(400,400),
	linewidth=2, legend=:topleft)

plot( 10:40, x->x^2/180, label=L"N^2", lc=:red )
scatter!( 1:40, backward_mean,yerror=backward_std, color=:black, label="Cost Gradient" )

# plot!( 1:100, x->x^3/205800, lc=:red )
# scatter!( 1:83, forward_mean, yerror=forward_std, color=:gray, label="Cost Function" )

# plot!( 20:60, x->x^2/1580, lc=:red )
# scatter!( 1:83, continuation_mean, yerror=continuation_std, color=:gold, label="Continuation" )

xticks!([1,3,5,10,20,40],["1","3","5","10","20","40"])
yticks!([1e-1,1e0,1e1],["0.1","1","10"])
savefig("docs/figures/scaling.pdf")

####################################################################################################
####################################################################################################
################################################################################ two-state-optima
using BifurcationInference,Plots,StaticArrays
using Clustering,JLD
using StatsBase: median

trajectories, losses, bifurcations, K = Vector{Float64}[], Float64[], Float64[], 0
for file ∈ readdir(join=true,"trajectories")
	data = load(file)

	mask = @. (data["bifurcations"]==2)&(data["losses"]<0.1)
	if ~any(x->any(isnan.(x)),data["trajectory"])

		append!(trajectories,data["trajectory"][mask])
		append!(bifurcations,data["bifurcations"][mask])
		append!(losses,data["losses"][mask])
		K += 1
	end
end

K/length(readdir("trajectories"))
optima = hcat(trajectories...)

@. optima[end,:] = abs(optima[end,:])
M,N = size(optima)
clusters = dbscan(optima,0.5,min_cluster_size=40)

layout = @layout [ a{0.7w} [b{0.5h}; b{0.5h} ] ]
default(); default(msc=:auto,label="",size=(600,450),grid=false)

fig = scatter([0],[0],markersize=0,title=L"\mathrm{Parameter\,\ Estimates}\,\,\theta^{*}", layout=layout, subplot=1,
	ylabel=L"\mathrm{Sensitivity/Degredation\,\,Ratio}\quad k/\mu_1",xlabel=L"\mathrm{Basal\,\,Production\,\,of}\,\,u_1\,,\quad a_1")

color = [:lightblue,:pink]
include("applied/two-state.jl")
for (j,cluster) ∈ enumerate(clusters)
	x,y = 10 .^(getindex(optima,3,cluster.core_indices)), getindex(optima,5,cluster.core_indices) ./ 10 .^getindex(optima,1,cluster.core_indices)
	@assert all( i -> sign(getindex(optima,3,i)) == sign(getindex(optima,4,i)), cluster.core_indices )

	scatter!(x,y,color=color[j],subplot=1,xscale=:log10,xlim=(1e-3,1e3),yscale=:log10,ylim=(1e-1,1e2))
	annotate!(median(x),median(y),j,subplot=1)

	parameters = getindex(optima,:,cluster.core_indices)
	if j == 1 X = StateSpace( 2, 0:0.002:10, [4,5] ) else X = StateSpace( 2, 0:0.01:10, [4,5] ) end

	plot!(F,parameters[:,rand(1:length(cluster.core_indices))],X,
		determinant=false,xlim=(0,10),ylim=(0.01,100),yscale=:log10,
		subplot=j+1, framestyle = :box)
	annotate!(8.8,30,j,subplot=j+1)
end

xlabel!("",subplot=2)
ylabel!("",subplot=2)
ylabel!("",subplot=3)
yticks!([NaN],subplot=2)
yticks!([NaN],subplot=3)

savefig("docs/figures/two-state-optima.pdf")


text:
"""
    strip_arg(arg::Symbol)
Strip anything extra (type annotations, default values, etc) from an argument.
For now this cannot handle keyword arguments (it will throw an error).
"""
strip_arg(arg::Symbol) = arg # once we have a symbol, we have stripped everything, so we can just return it
function strip_arg(arg_expr::Expr)
    if arg_expr.head == :parameters # keyword argument
        error("strip_arg can't handle keyword args yet (parsing arg expression $(arg_expr))")
    elseif arg_expr.head == :(::) # argument is type annotated, remove the annotation
        return strip_arg(arg_expr.args[1])
    elseif arg_expr.head == :kw # argument has a default value, remove the default
        return strip_arg(arg_expr.args[1])
    else
        error("strip_arg encountered something unexpected. arg_expr was $(arg_expr)")
    end
end

"""
    required_func(signature)
Provide a default function implementation that throws an error when called.
"""
macro required_func(signature)
    if signature.head == :(=) # in this case a default implementation has already been supplied
        return esc(signature)
    end
    @assert signature.head == :call

    # get the names of all the arguments
    args = [strip_arg(expr) for expr in signature.args[2:end]]

    # get the name of the function
    fname = signature.args[1]

    error_string = "BayesNets.jl: No implementation of $fname for "

    # add each of the arguments to error string
    for (i,a) in enumerate(args)
        error_string *= "$a::\$(typeof($a))"
        if i == length(args)-1
            error_string *= ", and "
        elseif i != length(args)
            error_string *= ", "
        else
            error_string *= "."
        end
    end

    # if you are modifying this and want to debug, it might be helpful to print
    # println(error_string)

    body = Expr(:call, :error, parse("\"$error_string\""))

    return Expr(:function, esc(signature), esc(body))
end

"""
    infer_number_of_instantiations{I<:Int}(arr::AbstractVector{I})
Infer the number of instantiations, N, for a data type, assuming that it takes on the values 1:N
"""
function infer_number_of_instantiations{I<:Int}(arr::AbstractVector{I})
    lo, hi = extrema(arr)
    lo ≥ 1 || error("infer_number_of_instantiations assumes values in 1:N, value $lo found!")
    lo == 1 || warn("infer_number_of_instantiations assumes values in 1:N, lowest value is $(lo)!")
    hi
end

paramcount(::Bool) = 1
paramcount(::Real) = 1
paramcount(::Symbol) = 1
paramcount{R<:Real}(arr::Vector{R}) = length(arr)

function paramcount(tup::Tuple)
    retval = 0
    for v in tup
        retval += paramcount(v)
    end
    retval
end


text:
module IonLinearAlgebra

using LinearAlgebra

end


text:
# https://gridap.github.io/Tutorials/stable/pages/t009_stokes/

const E = 3.0e10 # Pa
const ν = 0.3 # dim-less
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
@law σe(ε) = λ*tr(ε)*one(ε) + 2*μ*ε # Pa
τ(ε) = sqrt(ε ⊙ σe(ε)) # Pa^(1/2)

const σ_u = 4.0e5 # Pa
const r_0 = σ_u / sqrt(E) # Pa^(1/2)
const H = 0.5 # dim-less

function d(r)
  1 - q(r)/r
end

function q(r)
  r_0 + H*(r-r_0)
end

@law function update(ε_in,r_in,d_in)
  τ_in = τ(ε_in)
  if τ_in <= r_in
    r_out = r_in
    d_out = d_in
    damaged = false
  else
    r_out = τ_in
    d_out = d(r_out)
    damaged = true
  end
  damaged, r_out, d_out
end

@law function σ(ε_in,r_in,d_in)

  _, _, d_out = update(ε_in,r_in,d_in)

  (1-d_out)*σe(ε_in)

end

@law function dσ(dε_in,ε_in,state)

  damaged, r_out, d_out = state

  if ! damaged
    return (1-d_out)*σe(dε_in)

  else
    c_inc = ((q(r_out) - H*r_out)*(σe(ε_in) ⊙ dε_in))/(r_out^3)
    return (1-d_out)*σe(dε_in) - c_inc*σe(ε_in)

  end
end

const b_max = VectorValue(0.0,0.0,-(9.81*2.5e3))

function project(q,trian,quad,order)

  a(u,v) = u*v
  l(v) = v*q
  t_Ω = AffineFETerm(a,l,trian,quad)

  V = TestFESpace(
    reffe=:Lagrangian, valuetype=Float64, order=order,
    triangulation=trian, conformity=:L2)

  U = TrialFESpace(V)
  op = AffineFEOperator(U,V,t_Ω)
  qh = solve(op)
  qh
end

function main(;n,nsteps)

  r = 12
  domain = (0,r,0,1,0,1)
  partition = (r*n,n,n)
  model = CartesianDiscreteModel(domain,partition)

  labeling = get_face_labeling(model)
  add_tag_from_tags!(labeling,"supportA",[1,3,5,7,13,15,17,19,25])
  add_tag_from_tags!(labeling,"supportB",[2,4,6,8,14,16,18,20,26])
  add_tag_from_tags!(labeling,"supports",["supportA","supportB"])

  order = 1

  V = TestFESpace(
    reffe=:Lagrangian, valuetype=VectorValue{3,Float64}, order=order,
    model=model,
    labels = labeling,
    dirichlet_tags=["supports"])

  U = TrialFESpace(V)

  trian = Triangulation(model)
  degree = 2*order
  quad = CellQuadrature(trian,degree)

  r = QPointCellField(r_0,trian,quad)
  d = QPointCellField(0.0,trian,quad)

  nls = NLSolver(show_trace=true, method=:newton)
  solver = FESolver(nls)

  function step(uh_in,factor)

    b = factor*b_max
    res(u,v) = ( ε(v) ⊙ σ(ε(u),r,d) ) - v⋅b
    function jac(u,du,v)
      state = update(ε(u),r,d)
      ε(v) ⊙ dσ(ε(du),ε(u),state)
    end
    t_Ω = FETerm(res,jac,trian,quad)
    op = FEOperator(U,V,t_Ω)

    uh_out, = solve!(uh_in,solver,op)

    update_state_variables!(quad,update,ε(uh_out),r,d)

    uh_out
  end

  factors = collect(1:nsteps)*(1/nsteps)
  uh = zero(V)

  for (istep,factor) in enumerate(factors)

    println("\n+++ Solving for load factor $factor in step $istep of $nsteps +++\n")

    uh = step(uh,factor)
    dh = project(d,trian,quad,order)
    rh = project(r,trian,quad,order)

    writevtk(
      trian,"results_$(lpad(istep,3,'0'))",
      cellfields=["uh"=>uh,"epsi"=>ε(uh),"damage"=>dh,
                  "threshold"=>rh,"sigma_elast"=>σe(ε(uh))])

  end

end

main(n=6,nsteps=20)

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging


text:
module ArrayAllez

include("cache.jl")

include("threads.jl")

include("inplace.jl")

@init @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
	include("inplace-flux.jl")
	include("prod+cumprod.jl")
end
include("inplace-zygote.jl")

include("dropdims.jl")

end


text:
module LandingGear

using LinearAlgebra
using StaticArrays
using ComponentArrays
using UnPack

using Flight.Modeling
import Flight.Modeling: init_x0, init_y0, init_u0, init_d0, f_cont!, f_disc!

using Flight.Plotting
import Flight.Plotting: plots

using Flight.Attitude
using Flight.Geodesy
using Flight.Terrain
using Flight.Kinematics
using Flight.Dynamics

using Flight.Components
import Flight.Components: WrenchTrait, AngularMomentumTrait, get_wr_b

export NoSteering, DirectSteering, get_steering_angle
export NoBraking, DirectBraking, get_braking_coefficient
export SimpleDamper, Strut
export LandingGearUnit

#basis elements, for convenience
const e1 = SVector{3,Float64}(1,0,0)
const e3 = SVector{3,Float64}(0,0,1)

########################### Steering #############################

abstract type AbstractSteering <: SystemDescriptor end

############## NoSteering ##############

struct NoSteering <: AbstractSteering end

get_steering_angle(::System{NoSteering}) = 0.0
f_cont!(::System{NoSteering}, args...) = nothing
f_disc!(::System{NoSteering}, args...) = false

############## DirectSteering ##############

Base.@kwdef struct DirectSteering <: AbstractSteering
    ψ_max::Float64 = π/6
end

Base.@kwdef struct DirectSteeringY
    ψ::Float64 = 0.0
end
#we need to make the contents of u mutable
init_u0(::DirectSteering) = Ref(0.0) #steering input
init_y0(::DirectSteering) = DirectSteeringY(0.0) #steering angle

function f_cont!(sys::System{DirectSteering})
    u = sys.u[]
    @assert abs(u) <= 1 "Steering input must be within [-1,1]"
    sys.y = DirectSteeringY(u * sys.params.ψ_max)
end

f_disc!(::System{DirectSteering}, args...) = false

get_steering_angle(sys::System{DirectSteering}) = sys.y.ψ

# function plots(t, data::AbstractVector{<:DirectSteeringY}; mode, save_path, kwargs...)
#     println("To be implemented")
# end

########################### Braking #############################

abstract type AbstractBraking <: SystemDescriptor end

############## NoBraking ###############

struct NoBraking <: AbstractBraking end

get_braking_coefficient(::System{NoBraking}) = 0.0
f_cont!(::System{NoBraking}, args...) = nothing
f_disc!(::System{NoBraking}, args...) = false

########### DirectBraking #############

Base.@kwdef struct DirectBraking <: AbstractBraking
   η_br::Float64 = 1.0 #braking efficiency
end

Base.@kwdef struct DirectBrakingY
   α_br::Float64 = 0.0 #braking coefficient
end

init_u0(::DirectBraking) = Ref(0.0)
init_y0(::DirectBraking) = DirectBrakingY()

function f_cont!(sys::System{DirectBraking})
    u = sys.u[]
    @assert (u <= 1 && u >= 0) "Braking input must be within [0,1]"
    sys.y = DirectBrakingY(u * sys.params.η_br)
end

f_disc!(::System{DirectBraking}, args...) = false

get_braking_coefficient(sys::System{DirectBraking}) = sys.y.α_br

function plots(t, data::AbstractVector{<:DirectBrakingY}; mode, save_path, kwargs...)
    # println("To be implemented")
end

########################### Damper #############################

abstract type AbstractDamper end #not a System!
get_damper_force(args...) = throw(MethodError(get_damper_force, args))

Base.@kwdef struct SimpleDamper <: AbstractDamper
    k_s::Float64 = 25000 #spring constant
    k_d_ext::Float64 = 1000 #extension damping coefficient
    k_d_cmp::Float64 = 1000 #compression damping coefficient
    ξ_min::Float64 = -2 #compression below which the shock absorber is disabled
end

#Force exerted by the damper along zs. The deformation ξ is positive along z_s
#(elongation). The resulting force can be negative, meaning the damper pulls the
#piston rod assembly upwards along the negative z_s axis. This can happen when
#ξ_dot > 0 and sufficiently large
function get_damper_force(c::SimpleDamper, ξ::Real, ξ_dot::Real)
    k_d = (ξ_dot > 0 ? c.k_d_ext : c.k_d_cmp)
    F = -(c.k_s * ξ + k_d * ξ_dot)
    F = F * (ξ > c.ξ_min)
    return F
end

########################## Strut #########################

Base.@kwdef struct Strut{D<:AbstractDamper} <: SystemDescriptor
    t_bs::FrameTransform = FrameTransform()
    damper::D = SimpleDamper()
    l_0::Float64 = 0.0 #natural length
end

Base.@kwdef struct StrutY
    wow::Bool = false #weight-on-wheel flag
    ξ::Float64 = 0.0 #damper deformation
    ξ_dot::Float64 = 0.0 #damper deformation rate
    t_sc::FrameTransform = FrameTransform() #strut to contact frame transform
    t_bc::FrameTransform = FrameTransform() #body to contact frame transform
    F_dmp::Float64 = 0.0 #damper force on the ground along the strut z axis
    v_eOc_c::SVector{3,Float64} = zeros(SVector{3}) #contact point velocity
    srf::SurfaceType = Terrain.DryTarmac #surface type at the reference point
end

init_y0(::Strut) = StrutY()

function f_cont!(sys::System{<:Strut}, steering::System{<:AbstractSteering},
    terrain::AbstractTerrain, kin::KinData)

    @unpack t_bs, damper, l_0 = sys.params
    @unpack n_e, h_e, q_eb, q_nb = kin.pos
    @unpack v_eOb_b, ω_eb_b = kin.vel

    #strut frame axes
    q_bs = t_bs.q
    r_ObOs_b = t_bs.r #strut frame origin

    #compute contact reference point P
    r_OsP_s = l_0 * e3
    r_ObP_b = r_ObOs_b + q_bs(r_OsP_s)
    r_ObP_e = q_eb(r_ObP_b)
    r_OeOb_e = CartECEF(Geographic(n_e, h_e))
    r_OeP_e = r_OeOb_e + r_ObP_e
    P = Geographic(r_OeP_e)

    #get terrain data at the contact reference point (close enough to the actual
    #contact frame origin, which is still unknown)
    trn = get_terrain_data(terrain, P.l2d)
    h_trn = trn.altitude
    srf = trn.surface

    #compute damper deformation. if the projection of k_s onto z_n is close to
    #zero (strut nearly horizontal) or negative (strut upside down), set it to 0
    q_ns = q_nb ∘ q_bs
    k_s_zn = q_ns(e3)[3]
    Δh = AltOrth(P) - h_trn
    if k_s_zn > 1e-3
        ξ = min(0.0, Δh / k_s_zn)
    else
        ξ = 0.0 #0 instead of 0.0 causes type instability
    end
    wow = ξ < 0

    if !wow #we're done, assign only relevant outputs and return
        sys.y = StrutY(; wow, ξ, srf)
        return
    end

    #wheel frame axes
    ψ_sw = get_steering_angle(steering)
    q_sw = Rz(ψ_sw) #rotate strut axes to get wheel axes
    q_nw = q_ns ∘ q_sw #NED to contact axes rotation

    #contact frame axes
    k_c_n = trn.normal #NED components of contact z-axis
    i_w_n = q_nw(e1) #NED components of wheel x-axis
    i_w_n_trn = i_w_n - (i_w_n ⋅ k_c_n) * k_c_n #projection of i_w_n onto the terrain tangent plane
    i_c_n = normalize(i_w_n_trn) #NED components of contact x-axis
    j_c_n = k_c_n × i_c_n #NED components of contact y-axis
    R_nc = RMatrix(SMatrix{3,3}([i_c_n j_c_n k_c_n]), normalization = false)
    q_sc = q_ns' ∘ R_nc
    q_bc = q_bs ∘ q_sc

    #contact frame origin
    r_OsOc_s = e3 * (l_0 + ξ)
    r_OsOc_b = q_bs(r_OsOc_s)
    r_ObOc_b = r_OsOc_b + r_ObOs_b

    #construct contact frame transforms
    t_sc = FrameTransform(r_OsOc_s, q_sc)
    t_bc = FrameTransform(r_ObOc_b, q_bc)

    #contact frame origin velocity due to airframe motion
    v_eOc_afr_b = v_eOb_b + ω_eb_b × r_ObOc_b
    v_eOc_afr_c = q_bc'(v_eOc_afr_b)

    #compute the damper elongation rate required to cancel the airframe
    #contribution to the contact point velocity along the contact frame z axis
    q_sc = q_bs' ∘ q_bc
    k_s_c = q_sc'(e3)
    ξ_dot = -v_eOc_afr_c[3] / k_s_c[3]

    #compute contact point velocity
    v_eOc_dmp_c = k_s_c * ξ_dot
    v_eOc_c = v_eOc_afr_c + v_eOc_dmp_c

    F_dmp = get_damper_force(damper, ξ, ξ_dot)

    sys.y = StrutY(; wow, ξ, ξ_dot, t_sc, t_bc, F_dmp, v_eOc_c, srf)

end

f_disc!(::System{<:Strut}, args...) = false

function plots(t, data::AbstractVector{<:StrutY}; mode, save_path, kwargs...)

    @unpack ξ, ξ_dot, F_dmp = StructArray(data)

    pd = Dict{String, Plots.Plot}()

    splt_ξ = thplot(t, ξ;
        title = "Elongation", ylabel = L"$\xi \ (m)$",
        label = "", kwargs...)

    splt_ξ_dot = thplot(t, ξ_dot;
        title = "Elongation Rate", ylabel = L"$\dot{\xi} \ (m/s)$",
        label = "", kwargs...)

    splt_F_dmp = thplot(t, F_dmp;
        title = "Force", ylabel = L"$F_{dmp} \ (N)$",
        label = "", kwargs...)

    pd["01_dmp"] = plot(splt_ξ, splt_ξ_dot, splt_F_dmp;
        plot_title = "Damper",
        layout = (1,3), link = :none,
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    save_plots(pd; save_path)

end

########################### Contact #############################

#static / dynamic friction coefficient set
struct StaticDynamic
    static::Float64
    dynamic::Float64
end

#static / dynamic friction interpolation from breakout factor
get_μ(f::StaticDynamic, α_bo::Real)::Float64 = α_bo * f.dynamic + (1-α_bo) * f.static

#pure rolling friction coefficient set
Base.@kwdef struct Rolling
    sd::StaticDynamic= StaticDynamic(0.03, 0.02)
end

get_μ(f::Rolling, ::SurfaceType)::StaticDynamic = f.sd

#pure skidding friction coefficient set
Base.@kwdef struct Skidding
    dry::StaticDynamic = StaticDynamic(0.75, 0.25)
    wet::StaticDynamic = StaticDynamic(0.25, 0.15)
    icy::StaticDynamic = StaticDynamic(0.075, 0.025)
end

function get_μ(f::Skidding, srf::SurfaceType)::StaticDynamic
    if srf === Terrain.DryTarmac
        return f.dry
    elseif srf === Terrain.WetTarmac
        return f.wet
    elseif srf === Terrain.IcyTarmac
        return f.icy
    else
        error("Invalid surface type")
    end
end

get_μ(f::Union{Rolling,Skidding}, srf::SurfaceType, α_bo::Real)::Float64 =
    get_μ(get_μ(f, srf), α_bo)

#this allocates. why??
# get_sd(f::Skidding, ::Val{Terrain.DryTarmac}) = f.dry
# get_sd(f::Skidding, ::Val{Terrain.WetTarmac}) = f.wet
# get_sd(f::Skidding, ::Val{Terrain.IcyTarmac}) = f.icy

#note: setting k_l > 0 asymptotically removes residual force fighting between
#contacts in equilibrium at the expense of some position creep under persistent
#external loads. k_l = 0.2 gives reasonable practical results

Base.@kwdef struct Contact <: SystemDescriptor
    rolling::Rolling = Rolling() #rolling friction coefficients
    skidding::Skidding = Skidding() #skidding friction coefficients
    v_bo::SVector{2,Float64} = [0.005, 0.01] #breakout velocity interval
    ψ_skid::Float64 = deg2rad(10) #skidding slip angle thresold
    k_p::Float64 = 5.0 #proportional gain for contact velocity regulator
    k_i::Float64 = 400.0 #integral gain for contact velocity regulator
    k_l::Float64 = 0.0 #integrator leak factor for contact velocity regulator
end

WrenchTrait(::System{<:Contact}) = HasNoWrench()

Base.@kwdef struct ContactY
    v::SVector{2,Float64} = zeros(SVector{2}) #contact plane velocity
    s::SVector{2,Float64} = zeros(SVector{2}) #contact plane velocity integral
    α_p::SVector{2,Float64} = zeros(SVector{2}) #proportional μ scale factor
    α_i::SVector{2,Float64} = zeros(SVector{2}) #integral μ scale factor
    α_raw::SVector{2,Float64} = zeros(SVector{2}) #total scale factor, raw
    α::SVector{2,Float64} = zeros(SVector{2}) #total scale factor, clipped
    sat::SVector{2,Bool} = zeros(SVector{2,Bool}) #scale factor saturation flag
    α_bo::Float64 = 0.0 #breakout factor
    μ_roll::Float64 = 0.0 #rolling friction coefficient
    μ_skid::Float64 = 0.0 #skidding friction coefficient
    α_br::Float64 = 0.0 #braking factor
    ψ_cv::Float64 = 0.0 #tire slip angle
    μ_max::SVector{2,Float64} = zeros(SVector{2}) #maximum friction coefficient
    μ::SVector{2,Float64} = zeros(SVector{2}) #scaled friction coefficient
    f_c::SVector{3,Float64} = zeros(SVector{3}) #normalized contact force
    F_c::SVector{3,Float64} = zeros(SVector{3}) #contact force
    wr_b::Wrench = Wrench() #resulting airframe Wrench
end

init_x0(::Contact) = ComponentVector(x = 0.0, y = 0.0) #v regulator integrator states
init_y0(::Contact) = ContactY()
init_d0(::Contact) = Ref(true) #contact active?

function f_cont!(sys::System{Contact}, strut::System{<:Strut},
                braking::System{<:AbstractBraking})

    @unpack rolling, skidding, v_bo, ψ_skid, k_p, k_i, k_l = sys.params
    @unpack wow, t_sc, t_bc, F_dmp, v_eOc_c, srf = strut.y

    if !wow #better than !sys.d[], which is set at f_disc! so has a one-step delay
        sys.ẋ .= 0
        sys.y = ContactY()
        return
    end

    v = SVector{2,Float64}(v_eOc_c[1], v_eOc_c[2]) #contact point velocity
    s = SVector{2,Float64}(sys.x) #velocity integrator state

    #breakout factor
    α_bo = clamp((norm(v) - v_bo[1]) / (v_bo[2] - v_bo[1]), 0, 1)

    μ_roll = get_μ(rolling, srf, α_bo)
    μ_skid = get_μ(skidding, srf, α_bo)

    #longitudinal friction coefficient
    α_br = get_braking_coefficient(braking)
    @assert (α_br >= 0 && α_br <= 1)
    μ_x = μ_roll + (μ_skid - μ_roll) * α_br

    #tire slip angle
    if α_bo < 1 #prevents chattering in μ_y for near-zero contact velocity
        ψ_cv = ψ_skid
    else
        ψ_cv = atan(v[2], v[1])
    end
    ψ_abs = abs(ψ_cv)
    @assert (ψ_abs <= π)

    #lateral friction coefficient
    if ψ_abs < ψ_skid
        μ_y = μ_skid * ψ_abs / ψ_skid
    elseif ψ_abs > π - ψ_skid
        μ_y = μ_skid * (1 - (ψ_skid + ψ_abs - π)/ ψ_skid)
    else
        μ_y = μ_skid
    end

    #maximum friction coefficient vector
    μ_max = @SVector [μ_x, μ_y]

    #the magnitude of μ_max cannot exceed μ_skid; if it does, scale down its components
    μ_max *= min(1, μ_skid / norm(μ_max))

    α_p = -k_p * v
    α_i = -k_i * s
    α_raw = α_p + α_i #raw μ scaling
    α = clamp.(α_raw, -1, 1) #clipped μ scaling
    μ = α .* μ_max

    #if not saturated, integrator accumulates
    sat = abs.(α_raw) .> abs.(α) #saturated?
    sys.ẋ .= (v - k_l * s) .* .!sat

    #normalized contact force projected on the contact frame
    f_c = SVector{3,Float64}(μ[1], μ[2], -1)
    f_s = t_sc.q(f_c) #project normalized force onto the strut frame
    @assert f_s[3] < 0

    N = -F_dmp / f_s[3]
    N = max(0, N) #cannot be negative (could happen with large xi_dot >0)
    F_c = f_c * N

    wr_c = Wrench(F = F_c)
    wr_b = t_bc(wr_c)

    sys.y = ContactY(; v, s, α_p, α_i, α_raw, α, sat, α_bo, μ_roll, μ_skid,
        α_br, ψ_cv, μ_max, μ, f_c, F_c, wr_b)

end

function f_disc!(sys::System{Contact}, strut::System{<:Strut})

    if strut.y.wow
        sys.d[] = true #declare contact as active
    else #no wow
        if sys.d[] #contact was active, declare it inactive and reset integrator states
            sys.d[] = false
            sys.x .= 0 #reset integrator states
            return true #continuous state was modified
        end
    end
    return false
end

function plots(t, data::AbstractVector{<:ContactY}; mode, save_path, kwargs...)

    @unpack v, s, α_p, α_i, α_raw, α, sat, α_bo, μ_roll, μ_skid,
        α_br, ψ_cv, μ_max, μ, f_c, F_c, wr_b = StructArray(data)

    extract_xy = (input) -> (input |> StructArray |> StructArrays.components)

    (v_x, v_y), (α_p_x, α_p_y), (α_i_x, α_i_y), (sat_x, sat_y),
        (α_raw_x, α_raw_y), (α_x, α_y), (μ_max_x, μ_max_y), (μ_x, μ_y) =
         map(extract_xy, (v, α_p, α_i, sat, α, α, μ_max, μ))

    pd = Dict{String, Plots.Plot}()

    splt_v_mag = thplot(t, norm.(v); title = "Contact Velocity Magnitude",
        ylabel = L"$\alpha_{bo}$", label = "", kwargs...)
    splt_α_bo = thplot(t, α_bo; title = "Breakout Factor",
        ylabel = L"$\alpha_{bo}$", label = "", kwargs...)
    splt_μ_roll = thplot(t, μ_roll; title = "Rolling Friction Coefficient",
        ylabel = L"$\mu_{roll}$", label = "", kwargs...)
    splt_μ_skid = thplot(t, μ_skid; title = "Skidding Friction Coefficient",
        ylabel = L"$\mu_{skid}$", label = "", kwargs...)

    pd["01_srf"] = plot(splt_v_mag, splt_α_bo, splt_μ_roll, splt_μ_skid;
        plot_title = "Rolling and Skidding Friction",
        layout = (2,2), link = :none,
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_v_x = thplot(t, v_x; title = "Contact Point Velocity",
        ylabel = L"$v^x \ (m/s)$", label = "", kwargs...)
    splt_α_p_x = thplot(t, α_p_x; title = "Proportional Term",
        ylabel = L"$\alpha_p^x$", label = "", kwargs...)
    splt_α_i_x = thplot(t, α_i_x; title = "Integral Term",
        ylabel = L"$\alpha_i^x$", label = "", kwargs...)
    splt_sat_x = thplot(t, sat_x; title = "Saturation",
        ylabel = L"$S^x$", label = "", kwargs...)
    splt_α_raw_x = thplot(t, α_raw_x; title = "Raw Output",
        ylabel = L"$\alpha_{raw}^{x}$", label = "", kwargs...)
    splt_α_x = thplot(t, α_x; title = "Clipped Output",
        ylabel = L"$\alpha^{x}$", label = "", kwargs...)

    pd["02_reg_x"] = plot(splt_v_x, splt_α_p_x, splt_α_i_x,
                        splt_sat_x, splt_α_raw_x, splt_α_x;
        plot_title = "Longitudinal Friction Regulator",
        layout = (2,3),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_α_br = thplot(t, α_br; title = "Braking Coefficient",
        ylabel = L"$\alpha_{br}$", label = "", kwargs...)
    splt_μ_max_x = thplot(t, μ_max_x; title = "Maximum Friction Coefficient",
        ylabel = L"$\mu_{max}^{x}$", label = "", kwargs...)
    splt_μ_x = thplot(t, μ_x; title = "Effective Friction Coefficient",
        ylabel = L"$\mu^{x}$", label = "", kwargs...)

    pd["03_mu_x"] = plot(splt_α_br, splt_μ_max_x, splt_μ_x;
        plot_title = "Longitudinal Friction",
        layout = (1,3),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_v_y = thplot(t, v_y; title = "Contact Point Velocity",
        ylabel = L"$v^y \ (m/s)$", label = "", kwargs...)
    splt_α_p_y = thplot(t, α_p_y; title = "Proportional Term",
        ylabel = L"$\alpha_p^y$", label = "", kwargs...)
    splt_α_i_y = thplot(t, α_i_y; title = "Integral Term",
        ylabel = L"$\alpha_i^y$", label = "", kwargs...)
    splt_sat_y = thplot(t, sat_y; title = "Saturation",
        ylabel = L"$S^y$", label = "", kwargs...)
    splt_α_raw_y = thplot(t, α_raw_y; title = "Raw Output",
        ylabel = L"$\alpha_{raw}^{y}$", label = "", kwargs...)
    splt_α_y = thplot(t, α_y; title = "Clipped Output",
        ylabel = L"$\alpha^{y}$", label = "", kwargs...)

    pd["04_reg_y"] = plot(splt_v_y, splt_α_p_y, splt_α_i_y,
                        splt_sat_y, splt_α_raw_y, splt_α_y;
        plot_title = "Lateral Contact Regulator",
        layout = (2,3),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    splt_ψ_cv = thplot(t, rad2deg.(ψ_cv); title = "Tire Slip Angle",
        ylabel = L"$\psi_{cv} \ (deg)$", label = "", kwargs...)
    splt_μ_max_y = thplot(t, μ_max_y; title = "Maximum Friction Coefficient",
        ylabel = L"$\mu_{max}^{y}$", label = "", kwargs...)
    splt_μ_y = thplot(t, μ_y; title = "Effective Friction Coefficient",
        ylabel = L"$\mu^{y}$", label = "", kwargs...)

    pd["05_mu_y"] = plot(splt_ψ_cv, splt_μ_max_y, splt_μ_y;
        plot_title = "Lateral Friction",
        layout = (1,3),
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    pd["06_Fn_c"] = thplot(t, f_c;
        plot_title = "Normalized Contact Force",
        ylabel = [L"$f_{Oc \ (trn)}^{c}$" L"$f_{Oc \ (trn)}^{c}$" L"$f_{Oc \ (trn)}^{c}$"],
        th_split = :h, link = :none,
        kwargs...)

    pd["07_F_c"] = thplot(t, F_c;
        plot_title = "Contact Force",
        ylabel = [L"$F_{Oc \ (trn)}^{c} \ (N)$" L"$F_{Oc \ (trn)}^{c} \ (N)$" L"$F_{Oc \ (trn)}^{c} \ (N)$"],
        th_split = :h, link = :none,
        kwargs...)

    pd["08_wr_b"] = thplot(t, wr_b;
        plot_title = "Wrench [Airframe]",
        wr_source = "trn", wr_frame = "b",
        kwargs...)

    save_plots(pd; save_path)

end

########################## LandingGearUnit #########################

Base.@kwdef struct LandingGearUnit{L<:Strut, S <: AbstractSteering,
                            B <: AbstractBraking} <: SystemGroupDescriptor
    strut::L = Strut()
    contact::Contact = Contact()
    steering::S = NoSteering()
    braking::B = NoBraking()
end

#if we avoid the generic fallback for SystemGroup, we don't need to define
#traits for Steering, Braking, Contact or Strut
WrenchTrait(::System{<:LandingGearUnit}) = HasWrench()
AngularMomentumTrait(::System{<:LandingGearUnit}) = HasNoAngularMomentum()

get_wr_b(sys::System{<:LandingGearUnit}) = sys.y.contact.wr_b

function f_cont!(sys::System{<:LandingGearUnit}, kinematics::KinData,
                terrain::AbstractTerrain)

    @unpack strut, contact, steering, braking = sys.subsystems

    f_cont!(steering)
    f_cont!(braking)
    f_cont!(strut, steering, terrain, kinematics)
    f_cont!(contact, strut, braking)

    sys.y = (strut = strut.y, contact = contact.y, steering = steering.y,
            braking = braking.y)

end

function f_disc!(sys::System{<:LandingGearUnit})

    @unpack strut, contact = sys.subsystems

    # println("Called LdgUnit f_disc!")
    contact_modified = f_disc!(contact, strut)
    return contact_modified

end

#need to override the default AirframeNode implementation, since the underlying
#subcomponents are not AirframeComponents and only Contact contributes wr_b.
#could also define trivial get_wr_b and get_hr_b methods for steering, braking
#and strut, and let the default AirframeNode implementation do its thing
# get_wr_b(sys::System{<:LandingGearUnit}) = sys.y.contact.wr_b

end #module


text:
function no_processing(key, val, global_p, current_p)
    k = Symbol(key)
    current_p = push!!(current_p, k=>val)
    global_p, current_p
end

process_value(::Val{:boundaries}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:fields}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:particles_from_file}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:dist_fn}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:probe}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:collisions}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:qed}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:bremsstrahlung}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:injector}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)

process_value(::Val{:subset}, key, val, global_p, current_p) =
    no_processing(key, val, global_p, current_p)


text:
using GPUAcceleratedTracking, DrWatson, DataFrames
@quickactivate "GPUAcceleratedTracking"

raw_data_df = collect_results(datadir("benchmarks/codereplica"))

fig = plot_replica_benchmark(raw_data_df)

save( plotsdir("benchmark_textmem.pdf"), fig)

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
using Base.Test
using DataFrames

let
	# TODO:
	# - grouped variables in formulas with interactions
	# - is it fast?  Can expand() handle DataFrames?
	# - deal with intercepts
	# - implement ^2 for datavector's
	# - support more transformations with I()?

	#test_context("Formula")
	#test_group("Formula")

	d = DataFrame()
	d["y"] = [1:4]
	d["x1"] = PooledDataArray([5:8])
	d["x2"] = [9:12]
	d["x3"] = [11:14]
	d["x4"] = [12:15]
	f = Formula(:(y ~ x1 * (log(x2) + x3)))
	mf = model_frame(f, d)
	mm = model_matrix(mf)
	@test mm.model_colnames == [
	 "(Intercept)"
	 "x1:6"
	 "x1:7"
	 "x1:8"
	 "log(x2)"
	 "x3"
	 "x1:6&log(x2)"
	 "x1:6&x3"
	 "x1:7&log(x2)"
	 "x1:7&x3"
	 "x1:8&log(x2)"
	 "x1:8&x3" ]

	tmp = d["x2"]

	#test_group("Basic tests")

	d = DataFrame()
	d["y"] = [1:4]
	d["x1"] = [5:8]
	d["x2"] = [9:12]
	d["x3"] = [13:16]
	d["x4"] = [17:20]

	x1 = [5.:8]
	x2 = [9.:12]
	x3 = [13.:16]
	x4 = [17.:20]
	f = Formula(:(y ~ x1 + x2))
	mf = model_frame(f, d)
	mm = model_matrix(mf)
	@test mm.response_colnames == ["y"]
	@test mm.model_colnames == ["(Intercept)","x1","x2"]
	@test mm.response == transpose([1. 2 3 4])
	@test mm.model[:,1] == ones(4)
	@test mm.model[:,2:3] == [x1 x2]

	#test_group("expanding a PooledVec into a design matrix of indicators for each dummy variable")

	a = expand(PooledDataArray(x1), "x1", DataFrame())
	@test a[:,1] == DataVector([0, 1., 0, 0])
	@test a[:,2] == DataVector([0, 0, 1., 0])
	@test a[:,3] == DataVector([0, 0, 0, 1.])
	@test colnames(a) == ["x1:6.0", "x1:7.0", "x1:8.0"]

	#test_group("create a design matrix from interactions from two DataFrames")

	b = DataFrame()
	b["x2"] = DataVector(x2)
	df = interaction_design_matrix(a,b)
	@test df[:,1] == DataVector([0, 10., 0, 0])
	@test df[:,2] == DataVector([0, 0, 11., 0])
	@test df[:,3] == DataVector([0, 0, 0, 12.])

	#test_group("expanding an singleton expression/symbol into a DataFrame")

	df = deepcopy(d)
	r = expand(:x2, df)
	@test isa(r, DataFrame)
	@test r[:,1] == DataVector([9,10,11,12])  # TODO: test float vs int return

	df = deepcopy(d)
	ex = :(log(x2))
	r = expand(ex, df)
	@test isa(r, DataFrame)
	@test r[:,1] == DataVector(log([9,10,11,12]))

	# ex = :(x1 & x2)
	# r = expand(ex, df)
	# @test isa(r, DataFrame)
	# @test ncol(r) == 1
	# @test r[:,1] == DataArray([45, 60, 77, 96])

	r = expand(:(x1 + x2), df)
	@test isa(r, DataFrame)
	@test ncol(r) == 2
	@test r[:,1] == DataVector(df["x1"])
	@test r[:,2] == DataVector(df["x2"])

	df["x1"] = PooledDataArray(x1)
	r = expand(:x1, df)
	@test isa(r, DataFrame)
	@test ncol(r) == 3
	@test r == expand(PooledDataArray(x1), "x1", DataFrame())

	r = expand(:(x1 + x2), df)
	@test isa(r, DataFrame)
	@test ncol(r) == 4
	@test r[:,1:3] == expand(PooledDataArray(x1), "x1", DataFrame())
	@test r[:,4] == DataVector(df["x2"])

	df["x2"] = PooledDataArray(x2)
	r = expand(:(x1 + x2), df)
	@test isa(r, DataFrame)
	@test ncol(r) == 6
	@test r[:,1:3] == expand(PooledDataArray(x1), "x1", DataFrame())
	@test r[:,4:6] == expand(PooledDataArray(x2), "x2", DataFrame())

	#test_group("Creating a model matrix using full formulas: y ~ x1 + x2, etc")

	df = deepcopy(d)
	f = Formula(:(y ~ x1 & x2))
	mf = model_frame(f, df)
	mm = model_matrix(mf)
	@test mm.model == [ones(4) x1.*x2]

	f = Formula(:(y ~ x1 * x2))
	mf = model_frame(f, df)
	mm = model_matrix(mf)
	@test mm.model == [ones(4) x1 x2 x1.*x2]

	df["x1"] = PooledDataArray(x1)
	x1e = [[0, 1, 0, 0] [0, 0, 1, 0] [0, 0, 0, 1]]
	f = Formula(:(y ~ x1 * x2))
	mf = model_frame(f, df)
	mm = model_matrix(mf)
	@test mm.model == [ones(4) x1e x2 [0, 10, 0, 0] [0, 0, 11, 0] [0, 0, 0, 12]]

	#test_group("Basic transformations")

	df = deepcopy(d)
	f = Formula(:(y ~ x1 + log(x2)))
	mf = model_frame(f, df)
	mm = model_matrix(mf)
	@test mm.model == [ones(4) x1 log(x2)]

	df = deepcopy(d)
	df["x1"] = PooledDataArray([5:8])
	f = Formula(:(y ~ x1 * (log(x2) + x3)))
	mf = model_frame(f, df)
	mm = model_matrix(mf)
	@test mm.model_colnames == [
	 "(Intercept)"
	 "x1:6"
	 "x1:7"
	 "x1:8"
	 "log(x2)"
	 "x3"
	 "x1:6&log(x2)"
	 "x1:6&x3"
	 "x1:7&log(x2)"
	 "x1:7&x3"
	 "x1:8&log(x2)"
	 "x1:8&x3" ]

	#test_group("Model frame response variables")

	f = Formula(:(x1 + x2 ~ y + x3))
	mf = model_frame(f, d)
	@test mf.y_indexes == [1, 2]
	@test isequal(mf.formula.lhs, [:(x1 + x2)])
	@test isequal(mf.formula.rhs, [:(y + x3)])

	f = Formula(:(x1 + x2 ~ y + x3))
	mf = model_frame(f, d)
	@test mf.y_indexes == [1, 2]
	@test isequal(mf.formula.lhs, [:(x1 + x2)])
	@test isequal(mf.formula.rhs, [:(y + x3)])

	# unique_symbol tests
	#@test unique_symbols(:(x1 + x2)) ==  {"x1"=>:x1, "x2"=>:x2}
	#unique_symbols(:(y ~ x1 + x2 + x3))

	# additional tests from Tom
	y = [1., 2, 3, 4]
	mm = model_matrix(model_frame(Formula(:(y ~ x2)), d))
	@test mm.model == [ones(4) x2]
	@test mm.response == y''

	df = deepcopy(d)
	df["x1"] = PooledDataArray(df["x1"])

	mm = model_matrix(model_frame(Formula(:(y ~ x2 + x3 + x3*x2)), df))
	@test mm.model == [ones(4) x2 x3 x2.*x3]
	mm = model_matrix(model_frame(Formula(:(y ~ x3*x2 + x2 + x3)), df))
	@test mm.model == [ones(4) x3 x2 x2.*x3]
	mm = model_matrix(model_frame(Formula(:(y ~ x1 + x2 + x3 + x4)), df))
	@test mm.model[:,2] == [0, 1., 0, 0]
	@test mm.model[:,3] == [0, 0, 1., 0]
	@test mm.model[:,4] == [0, 0, 0, 1.]
	@test mm.model[:,5] == x2
	@test mm.model[:,6] == x3
	@test mm.model[:,7] == x4

	mm = model_matrix(model_frame(Formula(:(y ~ x2 + x3 + x4)), df))
	@test mm.model == [ones(4) x2 x3 x4]
	mm = model_matrix(model_frame(Formula(:(y ~ x2 + x2)), df))
	@test mm.model == [ones(4) x2]
	mm = model_matrix(model_frame(Formula(:(y ~ x2*x3 + x2&x3)), df))
	@test mm.model == [ones(4) x2 x3 x2.*x3]
	mm = model_matrix(model_frame(Formula(:(y ~ x2*x3*x4)), df))
	@test mm.model == [ones(4) x2 x3 x4 x2.*x3 x2.*x4 x3.*x4 x2.*x3.*x4]
	mm = model_matrix(model_frame(Formula(:(y ~ x2&x3 + x2*x3)), df))
	@test mm.model == [ones(4) x2.*x3 x2 x3]  # TODO This disagrees with R
	mm = model_matrix(model_frame(Formula(:(y ~ x2 & x3 & x4)), df))
	@test mm.model == [ones(4) x2.*x3.*x4]

	#test_group("Column groups in formulas")
	set_group(d, "odd_predictors", ["x1","x3"])
	@test expand(:odd_predictors, d) == d["odd_predictors"]
	mf = model_frame(Formula(:(y ~ odd_predictors)), d)
	@test mf.df[:,1] == d["y"]
	@test mf.df[:,2] == d["x1"]
	@test mf.df[:,3] == d["x3"]
	@test ncol(mf.df) == 3
	mf = model_frame(Formula(:(y ~ odd_predictors * x2)), d)
	mm = model_matrix(mf)
	@test mm.model == [ones(4) x1 x3 x2 x1.*x2 x3.*x2]
end


text:
using PokerHands
using Test

@testset "PokerHands.jl" begin
    @test chances("AhAs") > chances("KcKd")
end


text:
export PointMassInSphere

# Simple double integrator
mutable struct PointMassInSphere{T<:AbstractFloat} <: Robot
  mass::T
  r::T
  hard_limit_vel::T
  hard_limit_accel::T
  btCollisionObject
end
function PointMassInSphere{T}() where T
  # Default values
  mass = 1.0
  r = 0.1
  hard_limit_vel = 1.0
  hard_limit_accel = 1.0

  btCollisionObject = BulletCollision.convex_hull_cylinder(SVector{3}(zeros(3)), SVector{3}(0.,0.,1.), r)

  return PointMassInSphere{T}(mass, r, hard_limit_vel, hard_limit_accel, btCollisionObject)
end
PointMassInSphere(::Type{T} = Float64; kwargs...) where {T} = PointMassInSphere{T}(; kwargs...)


text:
"""
    ParticleRepresentation{S} <: AbstractFilterRepresentation{S}

Abstract type for the representation of the conditional distribution over the hidden state by particles (samples), weighted or unweighted.
"""
abstract type ParticleRepresentation{S} <: AbstractFilterRepresentation{S} end

"""
    get_pos(ensemble::ParticleRepresentation, i)

Return a view of the position of the particle with index `i` from `ensemble`.
"""
function get_pos(ensemble::ParticleRepresentation, i) end

"""
    list_of_particles(ensemble::ParticleRepresentation)

Return a list of particle position views from `ensemble`.
"""
list_of_pos(ens::ParticleRepresentation) = [get_particle(ens, i) for i in 1:no_of_particles(ens)]

"""
    no_of_particles(ensemble::ParticleRepresentation)

Return the number of particles in `ensemble`.
"""
function no_of_particles(ensemble::ParticleRepresentation) end

"""
    eff_no_of_particles(ensemble::ParticleRepresentation)

Return the effective number of particles in `ensemble`. For unweighted ensembles, this is equal to the number of particles.
For weighted ensembles, this returns the inverse of the sum of the squared weights
"""
eff_no_of_particles(ensemble::ParticleRepresentation) = no_of_particles(ensemble)

"""
    particle_dim(ensemble::ParticleRepresentation)

Return the dimension of individual particles in `ensemble`.
"""
function particle_dim(ensemble::ParticleRepresentation) end

dim(ensemble::ParticleRepresentation) = particle_dim(ensemble) * no_of_particles(ensemble)

"""
    eltype(ensemble::ParticleRepresentation)

Return the type of individual particles in `ensemble`.
"""
Base.eltype(ensemble::ParticleRepresentation{S}) where S = S

#function propagate!(ens::ParticleRepresentation{S}, model::HiddenStateModel{S}) where S end
#function propagate!(ens::ParticleRepresentation{S}, model::ContinuousTimeHiddenStateModel{S}, dt) where S end

# How to add new particle representations:
# * Add a struct which is a subtype of ParticleRepresentation{S}, where S is the type to be represented
# * Implement a method for dim which returns the dimensionality of the representation.
# * Implement methods for
#   - get_pos
#   - no_of_particles
#   - particle_dim

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using Test
using DiBitVectors

@testset "DiBitVectors" begin
    d0 = DiBitVector()
    d1 = DiBitVector(10)
    d2 = DiBitVector(10, 0)

    @test length(d0) == 0
    @test isempty(d0)
    @test_throws ArgumentError pop!(d0)
    push!(d0, 1)
    @test length(d0) == 1
    @test pop!(d0) == 1
    @test length(d0) == 0
    @test_throws ArgumentError pop!(d0)

    @test length(d1) == length(d2) == 10
    @test d1 == d2
    @test all(d1 .== 0)
    @test all(d2 .== 0)

    @test size(d1) == size(d2) == (10,)

    d3 = DiBitVector(30, 3)
    @test all(d3 .== 3)
    @test d3[1] == d3[end] == 3

    push!(d3, 0)
    @test length(d3) == 31 && length(d3.data) == 1
    push!(d3, 1)
    @test length(d3) == 32 && length(d3.data) == 1
    push!(d3, 2)
    @test length(d3) == 33 && length(d3.data) == 2
    push!(d3, 3)
    @test length(d3) == 34 && length(d3.data) == 2

    @test pop!(d3) == 3
    @test length(d3) == 33 && length(d3.data) == 2
    @test pop!(d3) == 2
    @test length(d3) == 32 && length(d3.data) == 1
    @test pop!(d3) == 1
    @test length(d3) == 31 && length(d3.data) == 1
    @test pop!(d3) == 0
    @test length(d3) == 30 && length(d3.data) == 1
    @test pop!(d3) == 3
    @test length(d3) == 29 && length(d3.data) == 1

    @test zero(d3) == DiBitVector(length(d3))

    @test_throws BoundsError d3[0]
    @test_throws BoundsError d3[-1]
    @test_throws BoundsError d3[99991]
end