text:
#does an initial filter to get down the file size
function filterToEndowments!(data::NCCSData)

  assetCutoff::Float64 = ASSET_CUTOFF
  assetLowerBound::Float64 = ASSET_LOWER_BOUND
  largeDifferenceCutoff::Float64 = LARGE_DIFFERENCE_CUTOFF
  #try to selection the desired endowments
  #remove missing data
  data.df = data.df[completecases(data.df[:,
    [:category, :netassets, :totrev, :totexp, :invrealizeddirect, :invrealizedindirect]]), :]
  data.df = data.df[data.df[:category] .== :B, :]
  data.df = data.df[data.df[:charitytype] .== :pc, :]
  N::Int = size(data.df,1)

  println("N1: $N")
  #filter out organizations which are too small
  data.df[:tokeep] = trues(size(data.df,1))
  sort!(data.df, [:ein, :fisyr])
  by(data.df, :ein) do subdf::SubDataFrame

    subN::Int = size(subdf,1)
    maxAssets::Float64 = maximum(subdf[:netassets])
    if maxAssets < assetCutoff || subN ≠ size(unique(subdf[:,[:ein, :fisyr]]),1)
      subdf[:tokeep] = falses(subN)
    end
  end

  #delete the invalid rows
  data.df = data.df[data.df[:tokeep] .== true, :]
  data.df = data.df[data.df[:netassets] .≥ assetLowerBound, :]
  N = size(data.df,1)
  println("N2: $N")

  data.df[:returnamt] = Vector{Union{Float64,Missing}}(missings(N))
  data.df[:return] = Vector{Union{Float64,Missing}}(missings(N))

  data.df[:returncheckamt] = Vector{Union{Float64,Missing}}(missings(N))
  data.df[:returncheck] = Vector{Union{Float64,Missing}}(missings(N))

  data.df[:largeDifference] = Vector{Union{Bool,Missing}}(missings(N))

  by(data.df, :ein) do subdf::SubDataFrame
    subN = size(subdf,1)

    @inbounds @simd for i ∈ 2:subN
      if subdf[i-1, :fisyr] == subdf[i, :fisyr] - 1

        #NOTE: returns = unrealized + inv income s.t. unrealized = Δassets - netincome
        returnCandidateAmt::Float64 = (subdf[i, :netassets] .- subdf[i-1, :netassets]
          ) .+ subdf[i, :invrealizeddirect] .- subdf[i, :netincome]
        returnCandidate::Float64 = returnCandidateAmt / subdf[i-1, :netassets]

        returnCheckCandidateAmt::Float64 = (subdf[i, :netassets] .- subdf[i-1, :netassets]
          ) .+ subdf[i, :invrealizedindirect] .- subdf[i, :netincome]
        returnCheckCandidate::Float64 = returnCheckCandidateAmt / subdf[i-1, :netassets]

        if (returnCandidate < 2. && returnCandidate > -0.9 &&
            abs(returnCandidate - returnCheckCandidate) ≤ largeDifferenceCutoff)# && abs(subdf[i, :totinv]) ≥ 1. #make sure we have valid numbers here
          subdf[i, :returnamt] = returnCandidateAmt
          subdf[i, :return] = returnCandidate

          subdf[i, :returncheckamt] = returnCheckCandidateAmt
          subdf[i, :returncheck] = returnCheckCandidate
        end
      end
    end
  end

  sort!(data.df, [:ein, :fisyr])

  return data
end

function makeWideEndowmentReturns(data::NCCSData)::DataFrame
  #form a wideDF
  widedf::DataFrame = DataFrame(ein = unique(data.df[:ein]))
  allYears::Dict = Dict(i=>Symbol(i) for i::Int ∈ collect(1989:2016))

  wideN::Int = size(widedf,1)
  N::Int = size(data.df,1)

  #make an index for the new dataframe
  einIndex::Dict = Dict(widedf[i, :ein]=>i for i::Int ∈ 1:wideN)

  #fill the rows for the wideDf
  for y::Int ∈ keys(allYears)
    widedf[allYears[y]] = Vector{Union{Missing, Float64}}(missings(wideN))
  end

  #this will hold the latest names
  widedf[:name] = Vector{Union{Missing, String}}(missings(wideN))

  sort!(data.df, [:ein, :fisyr])

  by(data.df, :ein) do subdf::SubDataFrame
    wideRow::Int = einIndex[subdf[1,:ein]]
    subN::Int = size(subdf,1)

    for i::Int ∈ 2:subN
      if !ismissing(subdf[i,:return])
        widedf[wideRow, allYears[subdf[i,:fisyr]]] = subdf[i,:return]
      end
    end

    #record the most recent name
    nameCandidates::Vector{String} = subdf[(!ismissing).(subdf[:,:name]), :name]
    if length(nameCandidates) > 0
      widedf[wideRow, :name] = nameCandidates[end]
    end
  end

  println(summary(widedf))

  return widedf
end

function computeEndowmentReturns(; refreshFilter::Bool = true,
    endowmentReturnsName::String = ENDOWMENT_RETURNS_NAME,
    writeEndowmentReturnsCSV::Bool = true,
    dataPath::String = DATA_PATH)::Nothing

  if refreshFilter
    data::NCCSData = deserializeNCCSData()
    data = filterToEndowments!(data)
    serializeNCCSData(data, dataName = endowmentReturnsName)
  else
    data = deserializeNCCSData(dataName = endowmentReturnsName)
  end

  widedf::DataFrame = makeWideEndowmentReturns(data)
  println("widedf size: $(size(widedf))")

  if writeEndowmentReturnsCSV
    outPath::String = "$dataPath\\$(endowmentReturnsName).csv"
    uCSV.write(outPath, widedf)

    outPath = "$dataPath\\$(endowmentReturnsName)_long.csv"
    uCSV.write(outPath, data.df)
  end

  return nothing

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
export DGSEM

struct DGSEM{L, C, G, A1, A2, A3, NF}
  law::L
  cell::C
  grid::G
  MJ::A1
  MJI::A2
  faceMJ::A3
  numericalflux::NF

  function DGSEM(; law, cell, grid, numericalflux)
    M = mass(cell)
    _, J = components(metrics(grid))
    MJ = M * J
    MJI = 1 ./ MJ

    faceM = facemass(cell)
    _, faceJ = components(facemetrics(grid))

    faceMJ = faceM * faceJ

    args = (law, cell, grid, MJ, MJI, faceMJ, numericalflux)
    new{typeof.(args)...}(args...)
  end
end

function (dg::DGSEM)(dq, q, time)
  cell = dg.cell
  grid = dg.grid
  device = Bennu.device(arraytype(cell))
  dim = ndims(cell)
  Nq = size(cell)[1]

  @assert all(size(cell) .== Nq)
  @assert(length(eltype(q)) == numberofstates(dg.law))

  comp_stream = Event(device)

  workgroup = ntuple(i -> i <= dim ? Nq : 1, 3)
  ndrange = (length(grid) * workgroup[1], Base.tail(workgroup)...)
  comp_stream = volumeterm!(device, workgroup)(
    dg.law,
    dq,
    q,
    derivatives_1d(cell)[1],
    metrics(grid),
    dg.MJ,
    dg.MJI,
    points(grid),
    Val(dim),
    Val(Nq),
    Val(numberofstates(dg.law));
    ndrange,
    dependencies = comp_stream
  )

  Nfp = Nq ^ (dim  - 1)
  workgroup_face = (Nfp,)
  ndrange = (Nfp * length(grid),)
  faceix⁻, faceix⁺ = faceindices(grid)
  facenormal, _ = components(facemetrics(grid))
  comp_stream = surfaceterm!(device, workgroup_face)(
    dg.law,
    dq,
    q,
    Val(Bennu.connectivityoffsets(cell, Val(2))),
    dg.numericalflux,
    dg.MJI,
    faceix⁻,
    faceix⁺,
    dg.faceMJ,
    facenormal,
    boundaryfaces(grid),
    points(grid),
    Val(dim);
    ndrange,
    dependencies = comp_stream
  )

  wait(comp_stream)
end

@kernel function volumeterm!(law,
                             dq,
                             q,
                             D,
                             metrics,
                             MJ,
                             MJI,
                             points,
                             ::Val{dim},
                             ::Val{Nq},
                             ::Val{Ns}) where {dim, Nq, Ns}
  @uniform begin
    FT = eltype(law)
    Nq1 = Nq
    Nq2 = dim > 1 ? Nq : 1
    Nq3 = dim > 2 ? Nq : 1
  end

  l_F = @localmem FT (Nq1, Nq2, Nq3, dim, Ns)
  dqijk = @private FT (Ns,)
  p_MJ = @private FT (1,)

  e = @index(Group, Linear)
  i, j, k = @index(Local, NTuple)

  @inbounds begin
    ijk = i + Nq * (j - 1 + Nq * (k - 1))

    g = metrics[ijk, e].g

    qijk = q[ijk, e]
    x⃗ijk = points[ijk, e]
    MJijk = MJ[ijk, e]

    fijk = flux(law, qijk, x⃗ijk)

    @unroll for s in 1:Ns
      @unroll for d in 1:dim
        l_F[i, j, k, d, s] = g[d, 1] * fijk[1, s]
        if dim > 1
          l_F[i, j, k, d, s] += g[d, 2] * fijk[2, s]
        end
        if dim > 2
          l_F[i, j, k, d, s] += g[d, 3] * fijk[3, s]
        end
        l_F[i, j, k, d, s] *= MJijk
      end
    end

    fill!(dqijk, -zero(FT))
    source!(law, dqijk, qijk, x⃗ijk)
    nonconservative_term!(law, dqijk, qijk, x⃗ijk)

    @synchronize

    ijk = i + Nq * (j - 1 + Nq * (k - 1))
    MJIijk = MJI[ijk, e]

    @unroll for n in 1:Nq
      Dni = D[n, i] * MJIijk
      Dnj = D[n, j] * MJIijk
      Dnk = D[n, k] * MJIijk
      @unroll for s in 1:Ns
        dqijk[s] += Dni * l_F[n, j, k, 1, s]
        if dim > 1
          dqijk[s] += Dnj * l_F[i, n, k, 2, s]
        end
        if dim > 2
          dqijk[s] += Dnk * l_F[i, j, n, 3, s]
        end
      end
    end

    dq[ijk, e] += dqijk
  end
end

@kernel function surfaceterm!(law,
                              dq,
                              q,
                              ::Val{faceoffsets},
                              numericalflux,
                              MJI,
                              faceix⁻,
                              faceix⁺,
                              faceMJ,
                              facenormal,
                              boundaryfaces,
                              points,
                              ::Val{dim}) where {faceoffsets, dim}
  @uniform begin
    FT = eltype(q)
    nfaces = 2 * dim
  end

  e⁻ = @index(Group, Linear)
  i = @index(Local, Linear)

  @inbounds begin
    @unroll for face in 1:nfaces
        j = i + faceoffsets[face]
        id⁻ = faceix⁻[j, e⁻]

        n⃗ = facenormal[j, e⁻]
        fMJ = faceMJ[j, e⁻]

        x⃗⁻ = points[id⁻]
        q⁻ = q[id⁻]

        boundarytag = boundaryfaces[face, e⁻]
        if boundarytag == 0
          id⁺ = faceix⁺[j, e⁻]
          q⁺ = q[id⁺]
        else
          q⁺ = boundarystate(law, n⃗, x⃗⁻, q⁻, boundarytag)
        end

        nf = surfaceflux(numericalflux, law, n⃗, x⃗⁻, q⁻, q⁺)
        dq[id⁻] -= fMJ * nf * MJI[id⁻]

        @synchronize(mod(face, 2) == 0)
    end
  end
end

ext: jl


text:
#
# # implements the online flexible least squares algorithm... modeled on Montana et al (2009):
# #   "Flexible least squares for temporal data mining and statistical arbitrage"
#
# # Our cost function: Cₜ(βₜ; μ) = (yₜ - xₜ'βₜ)² + μ ‖Δβₜ‖²
# #   Below we use Vω = μ⁻¹Iₚ  along with a nicer to use relationship: μ = (1 - δ) / δ
# #   We accept 0 <= δ <= 1 as a constructor argument which controls the weighting of new observations
# #   δ close to 0 corresponds to large μ, which means the parameter vector β changes slowly
# #   δ close to 1 corresponds to small μ, which means the parameter vector β changes quickly
#
# # OnlineFLS: Thomas Breloff (Cointegrated Technologies)
#
#
#
# # TODO: allow for time-varying Vω???
# #  to accomplish... lets represent Vω as a vector of Variance's (i.e. the diagonal of Vω)
#
# #-------------------------------------------------------# Type and Constructors
#
# type OnlineFLS <: OnlineStat
# 	p::Int  		# number of independent vars
# 	Vω::MatF    # pxp (const) covariance matrix of Δβₜ
# 	# Vω::Vector{Variance}
# 	Vε::Variance     # variance of error term
# 	yvar::Variance   # used for normalization
# 	xvars::Variances  # used for normalization
# 	# weighting::W  # weighting scheme for Variances... algo convergence is dictated by Vω
#
# 	n::Int
# 	β::VecF 		# the current estimate in: yₜ = xₜβₜ + εₜ
#
# 	# these are needed to update β
# 	R::MatF     # pxp matrix
# 	q::Float64  # called Q in paper
# 	K::VecF     # px1 vector (equivalent to Kalman gain)
#
# 	yhat::Float64  #most recent estimate of y
#
# 	function OnlineFLS(p::Integer, δ::Real = 0.0001, wgt::Weighting = default(Weighting))
#
# 		# calculate the covariance matrix Vω from the smoothing parameter δ
# 		@assert δ > 0. && δ <= 1.
# 		μ = (1. - δ) / δ
# 		Vω = eye(p) / μ
#
# 		Vε = Variance(wgt)
#
# 		yvar = Variance(wgt)
# 		xvars = Variances(p, wgt)
# 		# xvars = [Variance(wgt) for i in 1:p]
#
# 		# create and init the object
# 		o = new(p, Vω, Vε, yvar, xvars)
# 		empty!(o)
# 		o
# 	end
# end
#
#
# function OnlineFLS(x::AVecF, y::Real, δ::Real = 0.0001, wgt::Weighting = default(Weighting))
# 	p = length(x)
# 	o = OnlineFLS(p, δ, wgt)
# 	update!(o, x, y)
# 	o
# end
#
# function OnlineFLS(x::AMatF, y::AVecF, δ::Real = 0.0001, wgt::Weighting = default(Weighting))
# 	p = size(x,2)
# 	o = OnlineFLS(p, δ, wgt)
# 	update!(o, x, y)
# 	o
# end
#
# #-----------------------------------------------------------------------# state
#
# # nonzerostd(x) = if0then1(std(x))
#
# # state vars: [normalizedBeta, rawBeta, Variance(y), Variance(x), std(ε), mostRecentEstimateOfY, nobs]
# statenames(o::OnlineFLS) = [:βn, :β, :yvar, :xvars, :σε, :yhat, :nobs]
# state(o::OnlineFLS) = Any[copy(o.β), o.β * std(o.yvar) ./ if0then1(std(o.xvars)), copy(o.yvar), copy(o.xvars), std(o.Vε), o.yhat, nobs(o)]
#
# βn(o::OnlineFLS) = o.β
# StatsBase.coef(o::OnlineFLS) = (o.β * std(o.yvar)) ./ if0then1(std(o.xvars))
#
# #---------------------------------------------------------------------# update!
#
#
# # NOTE: assumes x mat is (T x p), where T is the number of observations
# # TODO: optimize
# # function update!(o::OnlineFLS, x::AMatF, y::AVecF)
# # 	@assert length(y) == size(x,1)
# # 	for i in length(y)
# # 		update!(o, row(x,i), y[i])
# # 	end
# # end
#
# function update!(o::OnlineFLS, x::AVecF, y::Real)
#
# 	# standardize y and x
# 	y = standardize!(o.yvar, y)
# 	x = standardize!(o.xvars, x)
#
# 	# calc error and update error variance
# 	yhat = dot(x, o.β)
# 	ε = y - yhat
# 	update!(o.Vε, ε)
#
# 	# update sufficient stats to get the Kalman gain
# 	o.R += o.Vω - (o.q * o.K) * o.K'
# 	Rx = o.R * x
# 	o.q = dot(x, Rx) + var(o.Vε)
# 	o.K = Rx / if0then1(o.q)
#
# 	# @LOG ε var(o.Vε)
# 	# @LOG diag(o.R)
# 	# @LOG Rx
# 	# @LOG o.q
# 	# @LOG o.K
#
# 	# update β
# 	o.β += o.K * ε
#
# 	@DEBUG o.β
#
# 	# save the denormalized estimate of y
# 	o.yhat = unstandardize(o.yvar, yhat)
#
# 	o.n += 1
# 	return
# end
#
# # NOTE: keeps consistent p... just resets state
# function Base.empty!(o::OnlineFLS)
# 	p = o.p
# 	empty!(o.Vε)
# 	o.n = 0
# 	o.β = zeros(p)
#
# 	# since Rₜ = Pₜ₋₁ + Vω, initialize with Vω
# 	o.R = copy(o.Vω)
#
# 	o.q = 0.
# 	o.K = zeros(p)
# 	o.yhat = 0.
# end
# #
# # function Base.merge!(o1::OnlineFLS, o2::OnlineFLS)
# # 	error("not implemented")
# # end
# #
# #
# # function StatsBase.coeftable(o::OnlineFLS)
# # 	error("not implemented")
# # end
# #
# # function StatsBase.confint(o::OnlineFLS, level::Real = 0.95)
# # 	error("not implemented")
# # end
#
# # predicts yₜ for a given xₜ
# function StatsBase.predict(o::OnlineFLS, x::AVecF)
# 	unstandardize(o.yvar, dot(o.β, standardize(o.xvars, x)))
# end
#
# # NOTE: uses most recent estimate of βₜ to predict the whole matrix
# function StatsBase.predict(o::OnlineFLS, x::AMatF)
# 	n = size(x,1)
# 	pred = zeros(n)
# 	for i in 1:n
# 		pred[i] = predict(o, row(x,i))
# 	end
# 	pred
# end


text:
using DataStructures

@testset "Vocabulary" begin

    words = ["a", "c", "-", "d", "c", "a", "b", "r", "a", "c", "d"]
    vocab = Vocabulary(words, 2, "</s>")
    @test vocab isa Vocabulary
    @test vocab.vocab isa Dict
    @test vocab.unk_cutoff isa Int
    @test vocab.unk_label isa String
    @test vocab.allword isa Array{String,1}
    @test length(vocab.vocab) == 4 #only 4 differnt string over word freq 2
    @test isequal(vocab.unk_cutoff, 2)
    @test vocab.unk_label == "</s>"
    @test isequal(vocab.allword ,["a", "c", "-", "d", "c", "a", "b", "r", "a", "c", "d", "</s>"])
    @test isequal(vocab.vocab, Dict{String,Int}("</s>"=>1,"c"=>3,"a"=>3,"d"=>2))
    #to check lookup function
    @test lookup(vocab,["a", "b", "c", "alien"]) == ["a", "</s>", "c", "</s>"]
    word_set = ["<unk>", "is", "already", "there"]
    @test_throws ErrorException Vocabulary(word_set, 1, "<unk>")
end

@testset "preprocessing" begin
    @testset "ngramizenew" begin
        sample_text = ["this", "is", "some", "sample", "text"]
        ngrams = TextAnalysis.ngramizenew(sample_text, 1)

        @test isequal(ngrams, ["this", "is", "some", "sample", "text"])

        ngrams = TextAnalysis.ngramizenew(sample_text,2)
        @test isequal(ngrams, ["this is", "is some", "some sample", "sample text"])

        ngrams = TextAnalysis.ngramizenew(sample_text,1,2)
        @test isequal(ngrams, ["this", "is", "some", "sample", "text", "this is", "is some", "some sample", "sample text"])
    end

    @testset "Padding function" begin
        example = ["1", "2", "3", "4", "5"]
        padded=padding_ngram(example, 2, pad_left=true, pad_right=true)
        @test isequal(padded,["<s> 1", "1 2", "2 3", "3 4", "4 5", "5 </s>"])
        @test isequal(example, ["<s>", "1", "2", "3", "4", "5", "</s>"])

        example = ["1", "2", "3", "4", "5"] #if used
        padded=padding_ngram(example, 2, pad_right=true)
        @test isequal(padded,["1 2", "2 3", "3 4", "4 5", "5 </s>"])
    end
    @testset "everygram function" begin
        example = ["1", "2", "3", "4", "5"]
        everyngms = everygram(example, min_len=1, max_len=2)
        @test isequal(everyngms, ["1", "2", "3", "4", "5", "1 2", "2 3", "3 4", "4 5"])
    end
end

@testset "counter" begin
    exam = ["To", "be", "or", "not", "to", "be","To", "be", "or", "not", "to", "be"]
    fit = (TextAnalysis.counter2(exam, 2, 2))
    @test fit isa DataStructures.DefaultDict
    @test length(fit) == 5 #length of unique words
end

@testset "language model" begin

    @testset "MLE" begin
        voc =["my", "name", "is", "salman", "khan", "and", "he", "is", "shahrukh", "Khan"]
        train = ["khan", "is", "my", "good", "friend", "and", "He", "is", "my", "brother"]
        model = MLE(voc)
        fit = model(train, 2, 2) #considering only bigrams
        unmaskedscore = score(model, fit, "is" ,"<unk>")
        @test unmaskedscore == 0.3333333333333333
        @test score(model, fit, "is", "alien") == Inf #context not in vocabulary
        @test score(model, fit, "alien", "is") == 0 # word not in vocabulary
    end

    @testset "Lidstone" begin
        voc =["my","name","is","salman","khan","and","he","is","shahrukh","Khan"]
        train = ["khan", "is", "my", "good", "friend", "and", "He", "is", "my", "brother"]
        model2 = Lidstone(voc, 1.0)
        fit = model2(train, 2, 2)
        @test score(model2, fit,"is", "alien") == 0.1
        @test score(model2, fit, "alien", "is") >= 0
    end
    @testset "Laplace" begin
        voc =["my","name","is","salman","khan","and","he","is","shahrukh","Khan"]
        train = ["khan", "is", "my", "good", "friend", "and", "He", "is", "my", "brother"]
        model3 = Laplace(voc)
        fit2 = model3(train, 2, 2)
        @test score(model3, fit2,"is", "alien") == 0.1
    end
    @testset "WittenBellInterpolated" begin
        voc =["my", "name", "is", "salman", "khan", "and", "he", "is", "shahrukh", "Khan"]
        train = ["khan", "is", "my", "good", "friend", "and", "He", "is", "my", "brother"]
        model = WittenBellInterpolated(voc)
        fit = model(train, 2, 2)
        @test score(model, fit,"is", "alien") == 0.2
        @test score(model, fit, "alien", "is") == 0.4
        @test score(model, fit,"alien") == 0.2 #should be non-zero
    end
    @testset "KneserNeyInterpolated" begin
        voc =["my", "name", "is", "salman", "khan", "and", "he", "is", "shahrukh", "Khan"]
        train = ["khan", "is", "my", "good", "friend", "and", "He", "is", "my", "brother"]
        model = KneserNeyInterpolated(voc, 0.1)
        fit = model(train, 2, 2)
        @test score(model, fit, "is", "alie") == 0.2
        @test score(model, fit, "alien", "is") == 0.11000000000000001
    end
end


text:
# -*- coding: utf-8 -*-
include("gerda_data_processing/init.jl")
include("gerda_data_processing/fct.jl")
include("fct-sim.jl")

plots_base_path = "pulses/sim/raw_0vbb/plots/"
base_path_raw   = "pulses/sim/raw_0vbb/"
base_path       = "pulses/sim/0vbb/"

parameters = JSON.parsefile("res/parameters.json")
bool_AE = false
std_intersect = 2;

# Search for available channels
channels = []
for c in 0:1:36 # change to 0:1:36 for ALL channels
    files = glob(base_path_raw*string(c)*"-"*channel_to_name[c]*"/*-wf*.h5");
    if length(files) > 0
        push!(channels, c)
    end
end

for ch in channels
    Base.run(`clear`)
#     IJulia.clear_output(true)
    str_ch = lpad(ch, 2, "0");
    @info("Processing Ch" * str_ch * " | "*channel_to_name[ch])
    @info(">---------------------------------------------------------------<")
#     files = glob(base_path_raw*string(ch)*"-"*channel_to_name[ch]*"/"*channel_to_name[ch]*"-wf*.h5");
    files = glob(base_path_raw*string(ch)*"-"*channel_to_name[ch]*"/"*channel_to_name[ch]*"-cl-wf*.h5");
#     append!(files, files2)
#     return length(files)
    # prog = Progress(length(files), dt=0.5,
    #                         barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
    #                         barlen=10)

    for file in files
        file_nr = findall(x->x==file, files)[1]
        filename  = split(file, basename(file))[1]
        filename *= channel_to_name[ch] * "_w_filter" * split(basename(file), channel_to_name[ch])[2]
        if !isfile(filename)
            @info("Processing file "*string(file_nr)*"/"*string(length(files)) * " for Ch " * string(ch) * " | " * channel_to_name[ch])
            @info(">---------------------------------------------------------------<")
            data_raw = HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "generated_waveforms")
            end

            t  = Table( waveform    = [],
                        evtno       = [],
                        event_mult  = [],
                        detno       = [],
                        det_ch      = [],
                        hits_totnum = [],
                        edep        = [],
                        pos         = []);
            temp = data_raw |> @filter(sum(_.waveform.value) > 0) |> Table

            if length(unique(data_raw.evtno)) != length(temp.evtno)
                evtno = unique(data_raw.evtno)
                @info("Sorting positive and negative waveforms")
                pro = Progress(length(evtno), dt=0.5,
                                barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                                barlen=10)
                for i in 1:1:length(evtno)
                    temp = data_raw |> @filter(_.evtno == evtno[i]) |> Table
                    if sum(temp[1].waveform.value) > 0
                        a = 1
                        b = 2
                    else
                        a = 2
                        b = 1
                    end
                    wf     = add_baseline_and_extend_tail.(temp.waveform, 6000, 12000)
                    time   = (0:10:11990)u"ns"
                    value  = wf[a].value[1:10:end]
                    wf_cut = [RDWaveform(time, value)]

                    append!(t, Table(waveform   = wf_cut,
                                    evtno       = [evtno[i]],
                                    event_mult  = [temp[1].event_mult],
                                    detno       = [temp[1].detno],
                                    det_ch      = [ch],
                                    hits_totnum = [temp[1].hits_totnum],
                                    edep        = [temp[1].edep],
                                    pos         = [temp[1].pos]))
                    next!(pro)
                end
            else
                waveforms_cut = []
                wf     = add_baseline_and_extend_tail.(temp.waveform, 6000, 12000)
                time   = (0:10:11990)u"ns"
                for i in 1:1:size(temp,1)
                    value  = wf[i].value[1:10:end]
                    push!(waveforms_cut, RDWaveform(time, value))
                end
                t  = Table(waveform = waveforms_cut,
                        evtno       = temp.evtno,
                        event_mult  = temp.event_mult,
                        detno       = temp.detno,
                        det_ch      = zeros(length(temp.evtno)) .+ ch,
                        hits_totnum = temp.hits_totnum,
                        edep        = temp.edep,
                        pos         = temp.pos);
            end

            events = Table( waveform    = ArrayOfRDWaveforms( Array{typeof(t.waveform[1]), 1}(t.waveform) ),
                            evtno       = Int32.(t.evtno),
                            event_mult  = Int32.(t.event_mult),
                            detno       = VectorOfArrays(Array{Array{Int32,1},1}(t.detno)),
                            det_ch      = Int32.(t.det_ch),
                            hits_totnum = Int32.(t.hits_totnum),
                            edep        = VectorOfArrays(Array{typeof(t[1].edep),1}(t.edep)),
                            pos         = VectorOfArrays(Array{typeof(t[1].pos),1}(t.pos)));

            @info("Store intermediate file")
            filename  = split(file, basename(file))[1]
            filename *= channel_to_name[ch] * "_w_bl" * split(basename(file), channel_to_name[ch])[2]
            HDF5.h5open(filename, "w") do h5f
                LegendHDF5IO.writedata(h5f, "data", events)
            end
            @info("Apply filter and time alignment")
            GBP = parameters[string(ch)]["filter"]["parameters"]["GBP"]*1e6
            tau = parameters[string(ch)]["filter"]["parameters"]["tau"]*1e-6
            Cd  = parameters[string(ch)]["filter"]["parameters"]["Cd"]*1e-12
            Cf  = parameters[string(ch)]["filter"]["parameters"]["Cf"]*1e-12

            if tau != 0 && GBP != 0
                data = Table( energy     = [],
                            multiplicity = [],
                            timestamp    = [],
                            run          = [],
                            channel      = [],
                            waveform     = [])
                filename  = split(file, basename(file))[1]
                filename *= channel_to_name[ch] * "_w_bl" * split(basename(file), channel_to_name[ch])[2]
                tt = HDF5.h5open(filename, "r") do h5f
                    LegendHDF5IO.readdata(h5f, "data")
                end
                pro = Progress(size(tt,1), dt=0.5,
                            barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                            barlen=10)
                for i in 1:1:size(tt,1)
                    pulse = tt[i].waveform.value
                    if sum(pulse) > 0
                        filtered_pulse = applyElectronics(pulse; Ts = 10e-9, GBP = GBP, tau = tau, Kv = 150e3, Cd = Cd, Cf = Cf, Rf = 500e6)
                        inter = find_intersect(filtered_pulse, maximum(filtered_pulse)/2, 5)
                        cut  = filtered_pulse[(inter-399):1:(inter+400)]
                        time = (0:10:7990)u"ns"
                        cut_pulse = RDWaveform(time, cut)

                        append!(data, Table(energy       = [tt[i].edep],
                                            multiplicity = [tt[i].event_mult],
                                            timestamp    = [tt[i].evtno],
                                            run          = [0],
                                            channel      = [tt[i].det_ch],
                                            waveform     = [cut_pulse]))
                        next!(pro)
                    end
                end
                @info(length(data.waveform))
                @info("Storing waveforms...")
                events = Table(energy        = VectorOfArrays(Array{typeof(data[1].energy),1}(data.energy)),
                                multiplicity = Int32.(data.multiplicity),
                                timestamp    = Int32.(data.timestamp),
                                run          = Int32.(data.run),
                                channel      = Int32.(data.channel),
                                waveform     = ArrayOfRDWaveforms( Array{typeof(data[1].waveform), 1}(data.waveform) ))

                filename  = split(file, basename(file))[1]
                filename *= channel_to_name[ch] * "_w_filter" * split(basename(file), channel_to_name[ch])[2]
                HDF5.h5open(filename, "w") do h5f
                    LegendHDF5IO.writedata(h5f, "data", events)
                end
            end
        else
            @info("File "*string(file_nr)*"/"*string(length(files))*" has already been processed!")
        end
        Base.run(`clear`)
#         IJulia.clear_output(true)

        #next!(prog)
    end
end


text:
using Test

const TallPosition  = CartesianIndex
const WidePosition  = NTuple{2,CartesianIndex}
const BlockPosition = Union{TallPosition,WidePosition}

"Compare `WidePosition`s in an order-agnostic way"
function Base.:(==)(left::WidePosition, right::WidePosition)
    same_order = left[1] == right[1] && left[2] == right[2]
    reversed   = left[1] == right[2] && left[2] == right[1]
    return same_order || reversed
end

#=------------------------------------------------------------------------------
| LEVEL ONE
|
| Start by modeling a `Block` and moving it at least one space up, down, left,
| or right.
|
| - You'll be given the shell of a `Block` object (or data structure), you'll
|   need to decide what data to store in it to keep track of the block's
|   location.
| - You will need to write the logic for a function (or class method) that will
|   move the block up, down, left, or right (no diagonal movement).
| - The tests will give the block's starting position as a list of the
|   coordinates of the squares the block is sitting on in (`row`, `col`)
|   format. One coordinate will indicate a block standing upright, two
|   coordinates will indicate a block lying on its side.
| - The tests will give directions in the form of a single character: 'U', 'D',
|   'L', or 'R'.
| - The tests will expect a list of the coordinates the block is sitting on, in
|   (`row`, `col`) format. One coordinate for a block standing upright, two
|   coordinates
|   for a block lying on its side.
------------------------------------------------------------------------------=#

struct Block
    # Create fields here, as needed
end

function Block(pos::BlockPosition)
    # Replace with a real constructor
    Block()
end

function level_one(block::Block, direction::Char)::BlockPosition
    # This function will be called for tests. Feel free
    # to write other functions to handle the logic for this
    # level, but use this function as your 'main' function.
    missing
end

# LEVEL ONE TESTS --------------------------------------------------------------

@testset "Can move once when the block is tall" begin
    start = CartesianIndex(0, 0)
    left  = (CartesianIndex(0, -2), CartesianIndex(0, -1))
    right = (CartesianIndex(0, 2),  CartesianIndex(0, 1))
    up    = (CartesianIndex(-2, 0), CartesianIndex(-1, 0))
    down  = (CartesianIndex(2, 0),  CartesianIndex(1, 0))

    @test level_one(Block(start), 'L') == left
    @test level_one(Block(start), 'R') == right
    @test level_one(Block(start), 'U') == up
    @test level_one(Block(start), 'D') == down
end

@testset "Can move once when the block is lying horizontally" begin
    start = (CartesianIndex(0, -1), CartesianIndex(0, 0))
    left  = CartesianIndex(0, -2)
    right = CartesianIndex(0, 1)
    up    = (CartesianIndex(-1, -1), CartesianIndex(-1, 0))
    down  = (CartesianIndex(1, -1),  CartesianIndex(1, 0))

    @test level_one(Block(start), 'L') == left
    @test level_one(Block(start), 'R') == right
    @test level_one(Block(start), 'U') == up
    @test level_one(Block(start), 'D') == down

end

@testset "Can move once when the block is lying vertically" begin
    start = (CartesianIndex(-1, 0), CartesianIndex(0, 0))
    left  = (CartesianIndex(-1, -1), CartesianIndex(0, -1))
    right = (CartesianIndex(-1, 1), CartesianIndex(0, 1))
    up    = CartesianIndex(-2, 0)
    down  = CartesianIndex(1, 0)

    @test level_one(Block(start), 'L') == left
    @test level_one(Block(start), 'R') == right
    @test level_one(Block(start), 'U') == up
    @test level_one(Block(start), 'D') == down

end

#=------------------------------------------------------------------------------
| LEVEL TWO
|
| Time to take that block for a ride. Now, given a list of directions, move
| the block around on an infinite flat grid and return the coordinates
| indicating the final direction.
|
| - Keep using your Block object (or data structure) from level 1.
| - Now, the tests will give directions as a list of 'U', 'D', 'L', or 'R'
|   characters.
| - The tests will expect a list of the coordinates the block is sitting on
|   after following the list of directions, in (`row`, `col`) format. One
|   coordinate for a block standing upright, two coordinates for a block
|   lying on its side.
------------------------------------------------------------------------------=#

function level_two(block::Block, directions::Vector{Char})::BlockPosition
    # This function will be called for tests. Feel free
    # to write other functions to handle the logic for this
    # level, but use this function as your 'main' function.
    missing
end

# LEVEL TWO TESTS -------------------------------------------------------------

@testset "Can move the block in a small circle" begin
    start      = CartesianIndex(0, 0)
    directions = ['U', 'R', 'D', 'L']
    stop       = (CartesianIndex(-1, 0), CartesianIndex(0, 0))
    @test level_two(Block(start), directions) == stop
end

@testset "Can move the block in a bigger circle" begin
    start      = CartesianIndex(0, 0)
    directions = ['U', 'U', 'L', 'L', 'D', 'D', 'R', 'R']
    stop       = CartesianIndex(0, 0)
    @test level_two(Block(start), directions) == stop
end

@testset "Can move diagonal there and back" begin
    start      = CartesianIndex(0, 0)
    directions = ['D', 'R', 'D', 'R', 'U', 'L', 'U', 'L']
    stop       = (CartesianIndex(0, 0), CartesianIndex(0, 0))
    @test level_two(Block(start), directions) == stop
end

@testset "Can make the block meander around" begin
    start      = CartesianIndex(0, 0)
    directions = ['R', 'U', 'R', 'D', 'R', 'U', 'R', 'D']
    stop       = (CartesianIndex(0, 5), CartesianIndex(0, 6))
    @test level_two(Block(start), directions) == stop
end

#=------------------------------------------------------------------------------
| LEVEL THREE
|
| Now that you can steer your block, it's time to start making it go where
| you want it to go. You may continue to assume an infinite flat grid, but
| now you need to take a coordinate and return the list of directions you
| moved your block to get it there. You should be able to identify a valid
| path from start to target, not necessarily the shortest path.
|
| - The level 3 tests will give you an initial block location and a target
|   location, as coordinates in the same format used throughout this
|   challenge.
| - In keeping with the game that inspired this challenge, targets will
|   alway be a single coordinate, and your block will need to stand upright
|   on that grid space in order to pass the test.
| - There are infinitely many paths that may be taken, considering the grid
|   is an infinite flat expanse, but remember that the tests will need to
|   run in a  reasonable amount of time.
| - The tests will expect a list of characters representing the directions
|   you moved, 'U', 'D', 'L', or 'R'.
------------------------------------------------------------------------------=#

function level_three(block::Block, target::BlockPosition)::Vector{Char}
    # This function will be called for tests. Feel free
    # to write other functions to handle the logic for this
    # level, but use this function as your 'main' function.
    missing
end

# LEVEL THREE TESTS -----------------------------------------------------------

@testset "Can move the block two spaces" begin
    start  = CartesianIndex(0, 0)
    target = CartesianIndex(-3, 0)
    path   = ['U', 'U']
    @test level_three(Block(start), target) == path
end

@testset "Can move the block one space in a slightly tricky way" begin
    start  = CartesianIndex(0, 0)
    target = CartesianIndex(0, 1)
    path   = level_three(Block(start), target)
    @test level_two(Block(start), path) == target
end

@testset "Can move the block pretty far away" begin
    start  = CartesianIndex(0, 0)
    target = CartesianIndex(100, -123)
    path   = level_three(Block(start), target)
    @test level_two(Block(start), path) == target
end

#=------------------------------------------------------------------------------
| LEVEL FOUR
|
| It's time to strive for efficiency! To accomplish this level, you should
| be able to identify the _shortest_ path from the block's start to the
| target.
|
| - The inputs and outputs of your function should be the same as for
|   level 3.
| - There may be more than one shortest path to get to the goal, you only
|   need to identify and return _one_ shortest path.
| - The tests will expect a list of characters representing the directions
|   you moved, 'U', 'D', 'L', or 'R'.
------------------------------------------------------------------------------=#

function level_four(block::Block, target::BlockPosition)::Vector{Char}
    # This function will be called for tests. Feel free
    # to write other functions to handle the logic for this
    # level, but use this function as your 'main' function.
    missing
end

# LEVEL FOUR TESTS -------------------------------------------------------------

@testset "The block can find the shortest short path" begin
    start  = CartesianIndex(0, 0)
    target = CartesianIndex(0, -1)
    path   = level_four(Block(start), target)
    @test level_two(Block(start), path) == target
    @test length(path) == 3
end

@testset "The block can find the shortest path down to the right" begin
    start  = CartesianIndex(0, 0)
    target = CartesianIndex(9, 9)
    path   = level_four(Block(start), target)
    @test level_two(Block(start), path) == target
    @test length(path) == 12
end

@testset "The block can find the shortest path up to the left" begin
    start  = CartesianIndex(0, 0)
    target = CartesianIndex(-8, -13)
    path   = level_four(Block(start), target)
    @test level_two(Block(start), path) == target
    @test length(path) == 15
end

#=------------------------------------------------------------------------------
| LEVEL FIVE
|
| Ok, smart guy/gal, you can navigate over an infinite, flat, unobstructed
| expanse, but how well does your block handle when there are obstacles? For
| this (final) level, you'll be given a representation of a finite grid
| containing obstacles, and will need to steer your block around those
| obstacles to the target. You'll also need to keep your block within the
| bounds of the grid.
|
| - You will need to parse a string into a grid representation. The tests
|   will give a string in the format ".....\n..x..\n.x...\x....\n...x.",
|   where '.' represents an empty space, and 'x' represents an obstacle.
| - Since the grid is now finite, you may assume that every coordinate not
|   explicity included in your grid representation as 'empty' is impassable.
| - You will continue to get a coordinate(s) for your block start and target
|   space.
| - The tests will continue to expect a list of characters representing the
|   directions you moved to get to the target.
------------------------------------------------------------------------------=#

function level_five(block::Block, target::BlockPosition, grid::String)::Vector{Char}
    # This function will be called for tests. Feel free
    # to write other functions to handle the logic for this
    # level, but use this function as your 'main' function.
    missing
end

# LEVEL FIVE TESTS -------------------------------------------------------------

@testset "The block can find the shortest short path with obstacles" begin
    start  = CartesianIndex(4, 5)
    target = CartesianIndex(4, 4)
    grid   = join(["..........",
                   "..........",
                   "...xx.....",
                   "..........",
                   "...xx.....",
                   "..........",
                   "..........",
                   "..........",
                   "..........",
                   ".........."], "\n")
    path   = level_four(Block(start), target)
    @test level_two(Block(start), path) == target
    @test length(path) == 7
end

@testset "The block can find the shortest longer path with obstacles" begin
    start  = CartesianIndex(3, 2)
    target = CartesianIndex(3, 8)
    grid   = join(["..........",
                   "..x.......",
                   "..x...x...",
                   "..x....x..",
                   "..x.......",
                   "..x.......",
                   "..x.......",
                   "..x.......",
                   "..........",
                   ".........."], "\n")
    path   = level_four(Block(start), target)
    @test level_two(Block(start), path) == target
    @test length(path) == 15
end

@testset "The block can find the shortest path with obstacles on a large grid" begin
    start  = CartesianIndex(14, 16)
    target = CartesianIndex(26, 4)
    grid = join(["...............................",
                 "...............................",
                 "...............x...............",
                 "...............x...............",
                 "...xxxxxxxxxxxxx...............",
                 "...............x...............",
                 "...............x...............",
                 "xxxxxxxxxxxxx..x...............",
                 "...............x...............",
                 "...............x...............",
                 "......xxxxxxxxxxxxxxxxxxxx.....",
                 "......x..................x.....",
                 "......x..................x.....",
                 "......x..................x.....",
                 "......x..................x.....",
                 "......x..................x.....",
                 "......x..................x.....",
                 "......x.xxxxxxxxxxxxxxxxxx.....",
                 "......x........................",
                 "......x........................",
                 "....xxxxxxxxxxxxxxxxxxxxxxxxxxx",
                 "...............................",
                 "xxxxxxxxxxxxxxxxxxxxxxxxxxx....",
                 "...............................",
                 "...............................",
                 "...............................",
                 "..............................."], "\n")
    path   = level_four(Block(start), target)
    @test level_two(Block(start), path) == target
    @test length(path) == 118
end

max_stars_repo_name: EvilTables/CodingDojo

max_issues_repo_name: EvilTables/CodingDojo

max_forks_repo_name: EvilTables/CodingDojo


text:
	include("MyStructs.jl")
function checkallcollision(posx::Float64,posy::Float64,vx::Float64,vy::Float64,radparticle::Float64, OF::Obstacles)

    nobs  = OF.n
    iscollided = false
    distancetravelled = Inf
    poscollidex = 0.0
    poscollidey = 0.0
    newvelocityx = 0.0
    newvelocityy = 0.0
    collidedi = 0
    # check if a collision takes place with every obstacle out there
    for i=1:nobs
        ofx,ofy,ofr = OF.x[i],OF.y[i],OF.r[i]
        # if the obstacle is too far away, it does not matter
        # if the particle collides when there is a collision taking
        # place earlier on

        if !((ofx-posx)^2+(ofy-posy)^2>distancetravelled^2+ofr^2+radparticle^2)
            _iscollided,_distancetravelled,_poscollidex,_poscollidey,
            _newvelocityx,_newvelocityy =
                checksinglecollision(posx,posy,vx,vy,radparticle,ofx,ofy,ofr)
            if _iscollided
                if _distancetravelled< distancetravelled
                    iscollided = true
                    distancetravelled = _distancetravelled
                    poscollidex = _poscollidex
                    poscollidey = _poscollidey
                    newvelocityx = _newvelocityx
                    newvelocityy = _newvelocityy
                    collidedi = i
                end
            end
        end
    end

    return iscollided,distancetravelled,poscollidex,poscollidey,newvelocityx,newvelocityy
end

function checksinglecollision(px::Float64,py::Float64,vx::Float64,vy::Float64,
                  radparticle::Float64,
                obsx::Float64,obsy::Float64,
                  radobs::Float64)

    # Checks whether the particle at posparticle with velocity veloparticle
    # and raidus radparticle collides
    # with the fixed obstacle at posobs with radius radobs.
    # If a collision takes place, we compute the distance to collision,
    # the position at which the collision takes place as well as the new
    # velocity that results from specular reflection

        # dist(s_min) = minimal distance between
        # posparticle+s*veloparticle and posobs
    vdotdiffpos = vx*(obsx-px) + vy*(obsy-py)
        if vdotdiffpos<= 0#dot(veloparticle,posobs .- posparticle)<= 0
            #s_min = -dot(posparticle .- posobs,veloparticle) /
            #dot(veloparticle,veloparticle)
            #if s_min<=radparticle+radobs
            # return if we are "inside" the obstacle or it is behind us
            return false,0.0,0.0,0.0,0.0,0.0
        end
        s_min = vdotdiffpos / (vx^2+vy^2)
        d_min = sqrt((px+s_min*vx-obsx)^2+(py+s_min*vy-obsy)^2)
        if d_min>=radparticle+radobs # no collision taking place
            return false,0.0,0.0,0.0,0.0,0.0
        end
        if s_min<=0# no collision taking place since the obs is behind the particle
            return false,0.0,0.0,0.0,0.0,0.0
        end
        if d_min> radobs+radparticle
            error("something went wrong")
        end
        if s_min<=0
            error("something went wrong here")
        end
        # d_min is the minimum,
        # but we collide earlier since radobs,radobs,radparticle>=0
        correction = sqrt((radobs+radparticle)^2 - d_min^2)

        distancetravelled = s_min - correction
        if distancetravelled<=0# no collision taking place since the obs is behind the particle
            return false,0.0,0.0,0.0,0.0,0.0
        end
        colx = px + distancetravelled*vx
        coly = py + distancetravelled*vy
        # new velocity: See Martin's Kinetic Theory lecture
        # only difference: particle may have finite size
        newvx,newvy = reflect(colx,coly,obsx,obsy,radobs,vx,vy)
        return true,distancetravelled*sqrt(vx^2+vy^2), colx,coly, newvx,newvy

end

function reflect(px,py,ox,oy,r,vx,vy)
    omegax = px-ox
    omegay = py-oy

    normomega = sqrt(omegax^2+omegay^2)
    omegax /= normomega
    omegay /= normomega

    dotprod = omegax*vx+omegay*vy
    return vx- 2.0*omegax*dotprod, vy-2.0*omegay*dotprod
end

function reflect(poscollide,posobs,radobs,v)
    Ω = (poscollide .- posobs)
    Ω /= norm(Ω)
    return v  .- 2.0*Ω*dot(Ω,v)
end

function checksinglecollision(posparticle::Array{Float64,1},
                  veloparticle::Array{Float64,1},
                  radparticle::Float64,
                  posobs::Array{Float64,1},
                  radobs::Float64)

    # Checks whether the particle at posparticle with velocity veloparticle
    # and raidus radparticle collides
    # with the fixed obstacle at posobs with radius radobs.
    # If a collision takes place, we compute the distance to collision,
    # the position at which the collision takes place as well as the new
    # velocity that results from specular reflection
    dim = size(posparticle,1)
    if dim==2
        # dist(s_min) = minimal distance between
        # posparticle+s*veloparticle and posobs
        if dot(veloparticle,posobs .- posparticle)<= 0
            #s_min = -dot(posparticle .- posobs,veloparticle) /
            #dot(veloparticle,veloparticle)
            #if s_min<=radparticle+radobs
            # return if we are "inside" the obstacle or it is behind us
            return false,0.0,[0.0;0.0],[0.0;0.0]
        end
        s_min = -dot(posparticle .- posobs,veloparticle) / dot(veloparticle,veloparticle)
        d_min = norm(posparticle .+ s_min*veloparticle .- posobs)
        if d_min>=radparticle+radobs # no collision taking place
            return false,0.0,[0.0;0.0],[0.0;0.0]
        end
        if s_min<=0# no collision taking place since the obs is behind the particle
            return false,0.0,[0.0;0.0],[0.0;0.0]
        end
        if d_min> radobs+radparticle
            error("something went wrong")
        end
        if s_min<=0
            error("something went wrong here")
        end
        # d_min is the minimum,
        # but we collide earlier since radobs,radobs,radparticle>=0
        correction = sqrt((radobs+radparticle)^2 - d_min^2)

        distancetravelled = s_min - correction
        if distancetravelled<=0# no collision taking place since the obs is behind the particle
            return false,0.0,[0.0;0.0],[0.0;0.0]
        end
        poscollide = posparticle .+ distancetravelled*veloparticle
        # new velocity: See Martin's Kinetic Theory lecture
        # only difference: particle may have finite size
        newvelocity = reflect(poscollide,posobs,radobs,veloparticle)
        return true,distancetravelled*norm(veloparticle), poscollide, newvelocity
    elseif dim==3
        error("Not yet implemented")
    end
end


text:
module HTensors

abstract type AbstractTensor{T} end

export HTensorNode, HTensorEvaluator
export isleaf, nmode, ncoeffs, shape, fulldim, array, getelt, getelts
export collapse
export preorder, postorder

include("ops.jl")
include("indexing.jl")
include("htensor.jl")
include("collapse.jl")
include("evaluator.jl")
include("ordering.jl")

end


text:
#==============================================================================#
#               ____ ____ _  _ ____ ___ ____ _  _ ___ ____                     #
#               |    |  | |\ | [__   |  |__| |\ |  |  [__                      #
#               |___ |__| | \| ___]  |  |  | | \|  |  ___]                     #
#==============================================================================#

@enum( FIRResponse, LOWPASS, BANDPASS, HIGHPASS, BANDSTOP )

#==============================================================================#
#          _  _ ____ _ ____ ____ ____    ___  ____ ____ _ ____ _  _            #
#          |_/  |__| | [__  |___ |__/    |  \ |___ [__  | | __ |\ |            #
#          | \_ |  | | ___] |___ |  \    |__/ |___ ___] | |__] | \|            #
#==============================================================================#

function kaiserlength( transition::Real, attenuation::Real = 60; samplerate = 1.0 )

    transition = transition./samplerate
    numtaps      = ceil( Int, ( attenuation - 7.95 )/( 2*π*2.285*transition ))

    if attenuation > 50
        β = 0.1102*( attenuation - 8.7 )
    elseif attenuation >= 21
        β = 0.5842*( attenuation - 21 )^( 0.4 ) + 0.07886*( attenuation - 21 )
    else
        β = 0.0
    end

    return numtaps, β
end

#==============================================================================#
#          ____ _ ____    ___  ____ ____ ___ ____ ___ _   _ ___  ____          #
#          |___ | |__/    |__] |__/ |  |  |  |  |  |   \_/  |__] |___          #
#          |    | |  \    |    |  \ |__|  |  |__|  |    |   |    |___          #
#==============================================================================#
# numtaps       = Desired number of filter taps
#             = If FIRType is HIGHPASS, may return  samples to make it
#               a type 1 filter.
# F           = Cutoff frequency. Real for high-pass & lowpass, Vector{Real} for
#               band-pass & band-reject
# FIRResponse = The response of the filter: LOWPASS, BANDPASS, HIGHPASS, BANDSTOP

function firprototype( numtaps::Integer, F::Union{Real, Vector{Real}}; response::FIRResponse = LOWPASS )
    M = numtaps-1
    if response == LOWPASS
        prototype = [ 2*F*sinc(2*F*(n-M/2)) for n = 0:M ]
    elseif response == BANDPASS
        prototype = [ 2*(F[1]*sinc(2*F[1]*(n-M/2)) - F[2]*sinc(2*F[2]*(n-M/2))) for n = 0:M ]
    elseif response == HIGHPASS
        M = isodd( M ) ? M+1 : M
        prototype = [ sinc(n-M/2) - 2*F*sinc(2*F*(n-M/2)) for n = 0:M ]
    elseif response == BANDSTOP
        prototype = [ 2*(F[2]*sinc(2*F[2]*(n-M/2)) - F[1]*sinc(2*F[1]*(n-M/2))) for n = 0:M ]
    else
        error("Not a valid FIR_TYPE")
    end

    return prototype
end

#==============================================================================#
#                        ____ _ ____ ___  ____ ____                            #
#                        |___ | |__/ |  \ |___ [__                             #
#                        |    | |  \ |__/ |___ ___]                            #
#==============================================================================#

function firdes( numtaps::Integer, cutoff::Union{Real, Vector{Real}}, windowfunction::Function; response::FIRResponse = LOWPASS, samplerate = 1.0, beta = 6.75 )

    cutoff    = cutoff ./ samplerate
    prototype = firprototype( numtaps, cutoff, response=response )
    numtaps   = length( prototype )

    if windowfunction == kaiser
        return prototype .* kaiser( numtaps, beta )
    else
        return prototype .*  windowfunction( numtaps )
    end

end

function firdes( cutoff::Union{Real, Vector{Real}}, transitionwidth::Real, stopbandAttenuation::Real = 60; response::FIRResponse = LOWPASS, samplerate = 1.0 )

    ( numtaps, β ) = kaiserlength( transitionwidth, stopbandAttenuation; samplerate = samplerate )
    firdes( numtaps, cutoff, kaiser, response = response, samplerate=samplerate, beta=β )

end


text:
# OLS-UCB, "ordinary least squares, upper confidence bound".
# Based on https://arxiv.org/abs/1612.01859.

abstract type OLSUCBOptimisationAlgorithm end
function optimise_linear_sqrtquadratic(instance::CombinatorialInstance{T}, ::OLSUCBOptimisationAlgorithm,
                                       linear::Dict{T, Float64}, sqrtquadratic::Dict{Tuple{T, T}, Float64},
                                       ε::Float64, d::Int, verbose::Int; with_trace::Bool=false) where T end

mutable struct OLSUCB <: Policy
  ε::Float64 # Only for approximation algorithm. TODO: Move, like ESCB-2?
  λ::Float64 # Numerical parameter.
  Γ::Matrix{Float64} # Approximation of the covariance matrix.
  algo::OLSUCBOptimisationAlgorithm
  verbose::Int # 0: nothing. 1: summary of function. 2: iteration counter too.

  function OLSUCB(ε::Float64, λ::Float64, Γ::Matrix{Float64}, algo::OLSUCBOptimisationAlgorithm, verbose::Int=0)
    if any(Γ .< 0)
      error("The approximate covariance matrix must have nonnegative entries")
    end

    if ! isposdef(Γ)
      error("The approximate covariance matrix must be positive definite")
    end

    # Elements outside the diagonal are bounded.
    for i in 1:size(Γ, 1)
      for j in 1:size(Γ, 2)
        if i != j
          if Γ[i, j] > sqrt(Γ[i, i] * Γ[j, j])
            error("The approximate covariance matrix does not respect the bounding property Γij ≤ √(Γii Γjj) for i = $i and j = $j")
          end
        end
      end
    end

    if λ <= 0.
      error("The lambda (λ) parameter must be strictly positive")
    end

    return new(ε, λ, Γ, algo, verbose)
  end
end

mutable struct OLSUCBDetails <: PolicyDetails
  nIterations::Int
  nBreakpoints::Int
  bestLambda::Float64
  bestLinearObjective::Float64
  bestNonlinearObjective::Float64
  solverTimes::Vector{Float64}

  function OLSUCBDetails()
    return new(0, 0, 0.0, 0.0, 0.0, Float64[])
  end
end

function choose_action(instance::CombinatorialInstance{T}, policy::OLSUCB, state::State{T}; with_trace::Bool=false) where T
  # For initialisation, if some arms have never been tried, force them to be tried.
  # Without this, the standard deviation term cannot exist.
  if any(v == 0 for v in values(state.arm_counts))
    weights = Dict(arm => (state.arm_counts[arm] == 0.0) ? 1.0 : 0.0 for arm in keys(state.arm_counts))

    t0 = time_ns()
    sol = solve_linear(instance, weights)
    t1 = time_ns()

    if with_trace
      run_details = OLSUCBDetails()
      push!(run_details.solverTimes, (t1 - t0) / 1_000_000_000)
      return sol, run_details
    else
      return sol
    end
  end

  # If all arms have been tried, can compute the two parts of the objective function: the average reward and a standard deviation.
  t_dict = state.arm_counts
  w_dict = state.arm_average_reward
  n = state.round

  arms = collect(keys(t_dict))
  t = Float64[t_dict[arm] for arm in arms]
  w = Float64[w_dict[arm] for arm in arms]

  d = length(w)
  fn = log(n) + (d + 2) * log(log(n)) + d / 2 * log(1 + e / policy.λ)
  S = 2 * fn * inv(Diagonal(t)) * (policy.λ * Diagonal(policy.Γ) * Diagonal(t) + Diagonal(t) * Diagonal(policy.Γ)) * inv(Diagonal(t))

  S_dict = Dict{Tuple{T, T}, Float64}((arms[i], arms[j]) => S[i, j] for (i, j) in product(1:length(arms), 1:length(arms)))

  # Solve the maximisation w^T x + \sqrt{x^T S x} based on the chosen algorithm.
  return optimise_linear_sqrtquadratic(instance, policy.algo, w_dict, S_dict, policy.ε, d, policy.verbose, with_trace=with_trace)
end


text:
using Minerva

env = RandomEnvironment()
_ = reset!(env)
while true
    a = rand(env.action_space)
    s′, r, done, info = step!(env, a)
    @info s′, r, done, info
    done && break
end


text:
params = ["T1" * foldr(*, [",T$i" for i in 2:n], init="") for n in 1:20]

for ps in params
    eval(Meta.parse("zero(::Type{Tuple{$ps}}) where {$ps} = map(zero, ($ps,))"))
    eval(Meta.parse("one( ::Type{Tuple{$ps}}) where {$ps} = map(one,  ($ps,))"))

    eval(Meta.parse(
        "zero(::Type{NamedTuple{names, Tuple{$ps}}}) where {names, $ps} = NamedTuple{names}(map(zero, ($ps,)))"
    ))
    eval(Meta.parse(
        "one( ::Type{NamedTuple{names, Tuple{$ps}}}) where {names, $ps} = NamedTuple{names}(map(one,  ($ps,)))"
    ))
end

+(a::T, b::T) where T <: Union{Tuple, NamedTuple} = map(+, a, b)
*(a::T, b::T) where T <: Union{Tuple, NamedTuple} = map(*, a, b)

#############################
### general score methods ###
#############################

function set_scores(g::Grammar, scores::NamedTuple)
    Grammar(
        g.categories,
        g.start_categories,
        g.depcomp,
        g.depcomps,
        g.catidx2dcompidx,
        g.binary_rules,
        g.terminal_rules,
        g.all_rules,
        g.binary_dict,
        g.terminal_dict,
        scores
    )
end

function add_score(g::Grammar, name::Symbol, score)
    set_scores(g, merge(g.scores, NamedTuple{(name,)}((score,))))
end

function tabulate_score(g::Grammar, score_name::Symbol)
    collect(
        getfield(score(g, catidx, ruleidx), score_name)
        for catidx  in findallin(unique(g.depcomp, g.categories), g.categories),
            ruleidx in eachindex(g.all_rules)
    )
end

###################
### Count Score ###
###################

count_score(category, rule) = rule(category) === nothing ? 0 : 1

#############################
### Expected Counts Score ###
#############################

function expected_count_score(g::Grammar, base_score_name::Symbol)
    @closure function get_expected_counts(catidx, ruleidx)
        dcompidx = g.catidx2dcompidx[catidx]
        p        = getfield(score(g, catidx, ruleidx), base_score_name)
        m, n     = length(g.depcomps), length(g.all_rules)
        c        = sparse([dcompidx], [ruleidx], p, m, n)
        ExpectedCounts(p, c)
    end
end

function add_expected_count_score(
        g::Grammar, base_score_name::Symbol, score_name::Symbol=:exp_counts
    )
    add_score(g, score_name, expected_count_score(g, base_score_name))
end

struct ExpectedCounts{m, n, S, A <: AbstractMatrix{S}}
    inside_score    :: S
    expected_counts :: A

    function ExpectedCounts(inside_score::S, expected_counts::A) where
            {S, A <: AbstractMatrix{S}}
        m, n = size(expected_counts)
        new{m, n, S, A}(inside_score, expected_counts)
    end
end

function expected_counts_dataframe(
        g::Grammar, ec::ExpectedCounts{m, n, S, SparseMatrixCSC{S, Int}}
    ) where {m ,n, S}

    M = ec.expected_counts / ec.inside_score
    dcompidxs, ruleidxs, counts = findnz(M)
    DataFrame(
        dcomp = g.depcomps[dcompidxs],
        rule  = g.all_rules[ruleidxs],
        count = counts
    )
end

function zero(::Type{ExpectedCounts{m, n, S, Matrix{S}}}) where {m ,n, S}
    ExpectedCounts(zero(S), zeros(S, m, n))
end

function one(::Type{ExpectedCounts{m, n, S, Matrix{S}}}) where {m ,n, S}
    ExpectedCounts(one(S), zeros(S, m, n))
end

function zero(::Type{ExpectedCounts{m, n, S, SparseMatrixCSC{S, Int}}}) where {m ,n, S}
    ExpectedCounts(zero(S), spzeros(S, m, n))
end

function one(::Type{ExpectedCounts{m, n, S, SparseMatrixCSC{S, Int}}}) where {m ,n, S}
    ExpectedCounts(one(S), spzeros(S, m, n))
end

function +(c::ExpectedCounts, d::ExpectedCounts)
    ExpectedCounts(c.inside_score    + d.inside_score,
                   c.expected_counts + d.expected_counts)
end

function *(c::ExpectedCounts, d::ExpectedCounts)
    p, q = c.inside_score, d.inside_score
    ExpectedCounts(p * q, p * d.expected_counts + q * c.expected_counts)
end

#########################
### Random Prob Score ###
#########################

function zero_robust_dirichlet_sample(a::AbstractVector)
    non_zero_idxs = findall(x -> x > 0, a)
    probs = rand(Dirichlet(a[non_zero_idxs]))
    b = fill(zero(LogProb), length(a))
    b[non_zero_idxs] = probs
    b
end

function sample_dirichlet_rows(M::AbstractMatrix)
    mapslices(zero_robust_dirichlet_sample, M, dims=2)
end

function random_prob_score(g::Grammar, base_score_name::Symbol)
    params = collect(
        getfield(score(g, catidx, ruleidx), base_score_name)
        for catidx  in findallin(unique(g.depcomp, g.categories), g.categories),
            ruleidx in eachindex(g.all_rules)
    )

    probs = sample_dirichlet_rows(params)

    @closure (catidx, ruleidx) -> probs[g.catidx2dcompidx[catidx], ruleidx]
end

function add_random_prob_score(
        g::Grammar, base_score_name::Symbol, score_name::Symbol=:prob
    )
    add_score(g, score_name, random_prob_score(g, base_score_name))
end

##########################
### Enumerated Forests ###
##########################

enum_forest_score(catidx, ruleidx) = EnumForest([[(catidx, ruleidx)]])

struct EnumForest
    trees :: Vector{Vector{Tuple{Int, Int}}}
end

zero(::Type{EnumForest}) = EnumForest(Vector{Tuple{Int, Int}}[])
one(::Type{EnumForest})  = EnumForest([Tuple{Int, Int}[]])

function +(f::EnumForest, g::EnumForest)
    EnumForest( [f.trees; g.trees] )
end
function *(f::EnumForest, g::EnumForest)
    EnumForest( [[tf; tg] for tf in f.trees for tg in g.trees] )
end

######################
### Scored Forests ###
######################

function add_forest_score(
        g::Grammar, base_score_name::Symbol, score_name::Symbol=:forest)
    add_score(g, score_name, forest_score(getfield(g.scores, base_score_name)))
end

forest_score(score) = function (catidx, ruleidx)
    [RuleApp(catidx, ruleidx, score(catidx, ruleidx))]
end

struct RuleApp{S}
    catidx   :: Int
    ruleidx  :: Int
    children :: Vector{Vector{RuleApp{S}}}
    score    :: S
end

function RuleApp(catidx, ruleidx, score::S) where S
    RuleApp(catidx, ruleidx, Vector{Vector{RuleApp{S}}}(), score)
end

ScoredForest{S} = Vector{RuleApp{S}}
score(f::ScoredForest) = sum(ra.score for ra in f)

zero(::Type{ScoredForest{S}}) where S = RuleApp{S}[]

# +(f::ScoredForest, g::ScoredForest) = append!(f, g)
# +(f::ScoredForest, ra::RuleApp) = push!(f, ra)

# function *(ra::RuleApp, f::ScoredForest)
#     ra.score *= score(f)
#     push!(ra.children, f)
#     ra
# end

# function *(f1::ScoredForest, f2::ScoredForest)
#     for ra1 in f1
# 		ra1 * f2
#     end
#     f1
# end

+(f::ScoredForest, g::ScoredForest) = [f; g]
+(f::ScoredForest, ra::RuleApp) = [f; [ra]]

function *(f1::ScoredForest, f2::ScoredForest)
    [ra * f2 for ra in f1]
end

function *(ra::RuleApp, f::ScoredForest)
    RuleApp(
        ra.catidx,
        ra.ruleidx,
        [ra.children; [f]],
        ra.score * score(f)
    )
end

#############################
### Random Tree Iterators ###
#############################

struct RandomParseTreeIterator{S}
    f :: ScoredForest{S}
end

rand(f::ScoredForest) = RandomParseTreeIterator(f)

IteratorSize(::Type{RandomParseTreeIterator{S}}) where S = Base.SizeUnknown()
eltype(::Type{RandomParseTreeIterator{S}}) where S = Tuple{Int, Int}

function iterate(
        iter::RandomParseTreeIterator,
        state = [categorical_sample(iter.f, getfield.(iter.f, :score))]
    )
    if isempty(state)
        nothing
    else
        ra, ras = state[1], state[2:end]
        new_ras = map(ra.children) do f
            categorical_sample(f, getfield.(f, :score))
        end
        (ra.catidx, ra.ruleidx), prepend!(ras, new_ras)
    end
end

"""
    tree_struct(grammar, categories_and_ruleids)

convert a parse tree represented by a list of tuples of categories and rules to
a tree data structure representation
"""
function tree_struct(grammar, catidxs_and_ruleidxs)
    function apply_rules!(nodes, rules)
        if isempty(rules)
            nothing
        else
            r, rs = rules[1], rules[2:end]
            n, ns = nodes[1], nodes[2:end]
            if n.value in grammar.categories && isapplicable(r, n.value)
                for c in r(n.value)
                    add_child!(n, c)
                end
                apply_rules!([children(n); ns], rs)
            else
                apply_rules!(ns, rules)
            end
        end
        nothing
    end

	category = grammar.categories[catidxs_and_ruleidxs[1][1]]
	root     = Tree(category, Any)

	ruleidxs = getindex.(catidxs_and_ruleidxs, 2)
	rules    = grammar.all_rules[ruleidxs]

    apply_rules!([root], rules)
    root
end

function tree_structs(g::Grammar, f::EnumForest)
    tree_struct.(Ref(g), f.trees)
end


text:
export pairwise

function pairwise(x1::Var, x2::Var)
    @assert length(x1.batchdims) == length(x2.batchdims)
    @assert ndims(x1) == ndims(x2) == 2
    y = pairwise(x1.data, x1.batchdims, x2.data, x2.batchdims)
    batchdims = [x1.batchdims[i]*x2.batchdims[i] for i=1:length(x1.batchdims)]
    Var(y, batchdims, pairwise, (x1,x2))
end

function pairwise{T}(x1::Matrix{T}, batchdims1::Vector{Int}, x2::Matrix{T}, batchdims2::Vector{Int})
    cumdim1 = 0
    cumdim2 = 0
    m1 = size(x1, 1)
    m2 = size(x2, 1)
    n = 0
    for i = 1:length(batchdims1)
        n += batchdims1[i] * batchdims2[i]
    end

    y = Array{T}(m1+m2, n)
    yi = 1
    for k = 1:length(batchdims1)
        n1 = batchdims1[k]
        n2 = batchdims2[k]
        for i = cumdim1+1:cumdim1+n1
            for j = cumdim2+1:cumdim2+n2
                copyto!(y, yi, x1, (i-1)*m1+1, m1)
                yi += m1
                copyto!(y, yi, x2, (j-1)*m2+1, m2)
                yi += m2
            end
        end
        cumdim1 += n1
        cumdim2 += n2
    end
    y
end

function addgrad!(y::Var, ::typeof(pairwise), x1::Var, x2::Var)
    T = eltype(y.data)
    gx1 = isvoid(x1.grad) ? Array{T}(0,0) : x1.grad
    gx2 = isvoid(x2.grad) ? Array{T}(0,0) : x2.grad
    ∇pairwise!(y.grad, gx1, x1.batchdims, gx2, x2.batchdims)
end

function ∇pairwise!{T}(gy::Matrix{T}, gx1::Matrix{T}, batchdims1::Vector{Int}, gx2::Matrix{T}, batchdims2::Vector{Int})
    cumdim1 = 0
    cumdim2 = 0
    m1 = size(gx1, 1)
    m2 = size(gx2, 1)
    yi = 1
    for k = 1:length(batchdims1)
        n1 = batchdims1[k]
        n2 = batchdims2[k]
        for i = cumdim1+1:cumdim1+n1
            for j = cumdim2+1:cumdim2+n2
                isempty(gx1) || BLAS.axpy!(m1, T(1), pointer(gy,yi), 1, pointer(gx1,(i-1)*m1+1), 1)
                yi += m1
                isempty(gx2) || BLAS.axpy!(m2, T(1), pointer(gy,yi), 1, pointer(gx2,(j-1)*m2+1), 1)
                yi += m2
            end
        end
        cumdim1 += n1
        cumdim2 += n2
    end
end


text:
"""
    RBEstimator <: Estimator

Unbiased estimator for response rate `p` see also:

> Kunzmann K, Kieser M. Point estimation and p‐values in phase II adaptive two‐stage designs with a binary endpoint. Statistics in medicine. 2017 Mar 15;36(6):971-84.
"""
struct RBEstimator{TD<:Design} <: Estimator

  design::TD

  RBEstimator{TD}(design::TD) where {TD<:Design} = new(design)

end # RBEstimator

"""
    RBEstimator(design::TD) where {TD<:Design}

Create unbiased estimator for response rate `p` see also:

> Kunzmann K, Kieser M. Point estimation and p‐values in phase II adaptive two‐stage designs with a binary endpoint. Statistics in medicine. 2017 Mar 15;36(6):971-84.

"""
RBEstimator(design::TD) where {TD<:Design} = RBEstimator{TD}(design)

Base.show(io::IO, estimator::RBEstimator) = print("RBEstimator")

function estimate(estimator::RBEstimator, x1::T, x2::T) where {T<:Integer}

  checkx1x2(x1, x2, design(estimator))
  n1          = interimsamplesize(design(estimator))
  nominator   = BigInt(0.0) # we need big ints, as binomials get HUGE!
  denominator = BigInt(0.0)
  x = x1 + x2
  n = samplesize(design(estimator), x1)
  for xx1 in 0:(min(x, n1))
    if samplesize(design(estimator), xx1) == n
      nominator   += BigInt(binomial(n1 - 1, xx1 - 1))*BigInt(binomial(n - n1, x - xx1))
      denominator += BigInt(binomial(n1, xx1))*BigInt(binomial(n - n1, x - xx1))
    end
  end
  return convert(Float64, nominator/denominator)

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
command(dev::Device, data::AbstractString) = in(data, dev.commands) ? dev.io.command(data) : error("command $(data) not supported by device")

function scaled_values(sensor::LegoSensor)
    vals = Vector{Float64}(sensor.io.num_values())
    dec = sensor.io.decimals()
    multiplier = 10.0 ^ (-dec)
    for j = 1:length(vals)
        vals[j] = multiplier * getfield(sensor.io, symbol(:value, j-1))()
    end
    vals
end

function run_continuous(dev::TachoMotor, speed::Integer)
    dev.io.speed_regulation("on")
    dev.io.speed_sp(speed)
    command(dev, "run-forever")
end

function run_continuous(dev::TachoMotor, speed::Integer, timeout_ms::Integer)
    dev.io.time_sp(timeout_ms)
    dev.io.speed_regulation("on")
    dev.io.speed_sp(speed)
    command(dev, "run-timed")
end

function stop(dev::TachoMotor, stop_command::AbstractString="brake")
    dev.io.stop_command(stop_command)
    command(dev, "stop")
end

abstract Angular{T}

immutable Degrees{T} <: Angular{T}
    value::T
end

immutable Radians{T} <: Angular{T}
    value::T
end

convert(::Type{Degrees}, rad::Radians) = radians.value * 180 / pi
convert(::Type{Radians}, deg::Degrees) = degrees.value * pi / 180

function servo_absolute(dev::TachoMotor, degrees::Degrees, speed=nothing)
    ticks = round(Int, degrees.value / 360. * dev.io.count_per_rot())
    dev.io.speed_regulation("on")
    (speed !== nothing) && dev.io.speed_sp(speed)
    dev.io.position_sp(ticks)
    command(dev, "run-to-abs-pos")
end
servo_absolute(dev::TachoMotor, radians::Radians) = servo_absolute(dev, convert(Degrees, radians))

function servo_relative(dev::TachoMotor, degrees::Degrees, speed=nothing)
    ticks = round(Int, degrees.value / 360. * dev.io.count_per_rot())
    dev.io.speed_regulation("on")
    (speed !== nothing) && dev.io.speed_sp(speed)
    dev.io.position_sp(ticks)
    command(dev, "run-to-rel-pos")
end
servo_relative(dev::TachoMotor, radians::Radians) = servo_relative(dev, convert(Degrees, radians))


text:
using CUDAnative
using CuArrays

function applyfun(x::CuDeviceArray, F::CuDeviceArray, f)
    i = ( blockIdx().x - 1) * blockDim().x + threadIdx().x

    F[i] = f(x[i])
    return nothing
end

nthreads = 512
nblocks = 10
F = CuArrays.zeros(nthreads * nblocks)
x = CuArray( LinRange(0, 10, nthreads*nblocks) )

f(x, c) = c*x^2

#This workds
g1(x) = f(x,1)
@cuda blocks=nblocks threads=nthreads applyfun(x, F, g1 )

println("hey")

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
using Documenter
using DistributedFactorGraphs

makedocs(
    modules = [DistributedFactorGraphs],
    format = Documenter.HTML(),
    sitename = "DistributedFactorGraphs.jl",
    pages = Any[
        "Home" => "index.md",
        "Getting Started" => [
            "Introduction" => "getting_started.md",
            "Variables and Factors" => "variables_and_factors.md",
            "Common API Interface" => "ref_api.md",
            "Example" => "example.md"
        ],
        "DistributedFactorGraph API's" => [
            "Graphs.jl" => "apis/graphs.md",
            "MetaGraph.jl" => "apis/graphs.md",
            "GraffSDK.jl" => "apis/graphs.md",
            "CloudGraphs.jl" => "apis/graphs.md",
        ],
        "Reference" => "func_ref.md"
    ]
    # html_prettyurls = !("local" in ARGS),
    )

deploydocs(
    repo   = "github.com/JuliaRobotics/DistributedFactorGraphs.jl.git",
    target = "build",
    deps   = Deps.pip("mkdocs", "python-markdown-math")
)


text:
function DampedSpring(a::Body, b::Body, anchorA::Vect, anchorB::Vect, restlength::Real, stiffness::Real, damping::Real)
    ccall(dlsym(libchipmunk, :cpDampedSpringNew), Ptr{Void}, (Ptr{Void}, Ptr{Void}, Vect, Vect, Cdouble, Cdouble, Cdouble,), a.ptr, b.ptr, anchorA, anchorB, restlength, stiffness, damping)
end

function get_anchorA(spring::DampedSpring)
    ccall(dlsym(libchipmunk, :cpDampedSpringGetAnchorA), Vect, (Ptr{Void},), spring.ptr)
end

function set_anchorA(spring::DampedSpring, anchor::Vect)
    ccall(dlsym(libchipmunk, :cpDampedSpringSetAnchorA), Void, (Ptr{Void}, Vect,), spring.ptr, anchor)
end

function get_anchorB(spring::DampedSpring)
    ccall(dlsym(libchipmunk, :cpDampedSpringGetAnchorB), Vect, (Ptr{Void},), spring.ptr)
end

function set_anchorB(spring::DampedSpring, anchor::Vect)
    ccall(dlsym(libchipmunk, :cpDampedSpringSetAnchorB), Void, (Ptr{Void}, Vect,), spring.ptr, anchor)
end

function get_restlength(spring::DampedSpring)
    ccall(dlsym(libchipmunk, :cpDampedSpringGetRestLength), Cdouble, (Ptr{Void},), spring.ptr)
end

function set_restlength(spring::DampedSpring, restlength::Real)
    ccall(dlsym(libchipmunk, :cpDampedSpringSetRestLength), Void, (Ptr{Void}, Cdouble,), spring.ptr, restlength)
end

function get_stiffness(spring::DampedSpring)
    ccall(dlsym(libchipmunk, :cpDampedSpringGetStiffness), Cdouble, (Ptr{Void},), spring.ptr)
end

function set_stiffness(spring::DampedSpring, stiffness::Real)
    ccall(dlsym(libchipmunk, :cpDampedSpringSetStiffness), Void, (Ptr{Void}, Cdouble,), spring.ptr, stiffness)
end

function get_damping(spring::DampedSpring)
    ccall(dlsym(libchipmunk, :cpDampedSpringGetDamping), Cdouble, (Ptr{Void},), spring.ptr)
end

function set_damping(spring::DampedSpring, damping::Real)
    ccall(dlsym(libchipmunk, :cpDampedSpringSetDamping), Void, (Ptr{Void}, Cdouble,), spring.ptr, damping)
end

export DampedSpring, get_anchorA, set_anchorA, get_anchorB, set_anchorB, get_restlength, set_restlength,
get_stiffness, set_stiffness, get_damping, set_damping


text:
function TurnAroundRangeTroposphericDelayModifier(arg0::DiscreteTroposphericModel)
    return TurnAroundRangeTroposphericDelayModifier((DiscreteTroposphericModel,), arg0)
end

function get_parameters_drivers(obj::TurnAroundRangeTroposphericDelayModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::TurnAroundRangeTroposphericDelayModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
#!/usr/bin/env julia
using CairoMakie

CairoMakie.activate!(type = "svg")
scene = scatter(rand(10,2),figure=(;resolution=(300,300)))
save("plot2.svg",scene)

# 55s
# julia --compile=min hello_svg.jl  27.43s

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
utc = FixedTimeZone("UTC")
warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)

# Converting a ZonedDateTime into a DateTime
dt = DateTime(2015, 1, 1, 0)
zdt = ZonedDateTime(dt, warsaw)
@test DateTime(zdt) == dt

# Converting from ZonedDateTime to DateTime isn't possible as it is always inexact.
@test_throws MethodError convert(DateTime, zdt)

# Vectorized accessors
arr = repmat([zdt], 10)
@test Dates.DateTime(arr) == repmat([dt], 10)

# now function
dt = Dates.unix2datetime(time())  # Base.now in UTC
zdt = now(warsaw)
@test zdt.timezone == warsaw
@test isapprox(map(Dates.datetime2unix, [dt, TimeZones.utc(zdt)])...)

# Changing time zones
dt = DateTime(2015, 1, 1, 0)
zdt_utc = ZonedDateTime(dt, utc; from_utc=true)
zdt_warsaw = ZonedDateTime(dt, warsaw; from_utc=true)

# Identical since ZonedDateTime is immutable
@test astimezone(zdt_utc, warsaw) === zdt_warsaw
@test astimezone(zdt_warsaw, utc) === zdt_utc

# ZonedDateTime to Unix timestamp (and vice versa)
@test TimeZones.zdt2unix(ZonedDateTime(1970, utc)) == 0
@test TimeZones.unix2zdt(0) == ZonedDateTime(1970, utc)

for dt in (DateTime(2013, 2, 13), DateTime(2016, 8, 11))
    zdt = ZonedDateTime(dt, warsaw)
    offset = TimeZones.value(zdt.zone.offset)   # Total offset in seconds
    @test TimeZones.zdt2unix(zdt) == datetime2unix(dt) - offset
end

@test isa(TimeZones.zdt2unix(ZonedDateTime(1970, utc)), Float64)
@test isa(TimeZones.zdt2unix(Float32, ZonedDateTime(1970, utc)), Float32)
@test isa(TimeZones.zdt2unix(Int64, ZonedDateTime(1970, utc)), Int64)
@test isa(TimeZones.zdt2unix(Int32, ZonedDateTime(1970, utc)), Int32)

@test TimeZones.zdt2unix(ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0.75
@test TimeZones.zdt2unix(Float32, ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0.75
@test TimeZones.zdt2unix(Int64, ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0
@test TimeZones.zdt2unix(Int32, ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0

@testset "round-trip" begin
    zdt = ZonedDateTime(2010, 1, 2, 3, 4, 5, utc)
    round_trip = TimeZones.unix2zdt(TimeZones.zdt2unix(Int64, zdt))
    @test round_trip == zdt
end

@testset "millisecond loss" begin
    zdt = ZonedDateTime(2010, 1, 2, 3, 4, 5, 999, utc)
    round_trip = TimeZones.unix2zdt(TimeZones.zdt2unix(Int64, zdt))
    @test round_trip != zdt
    @test round_trip == floor(zdt, Dates.Second)
end

@testset "timezone loss" begin
    zdt = ZonedDateTime(2010, 1, 2, 3, 4, 5, warsaw)
    round_trip = TimeZones.unix2zdt(TimeZones.zdt2unix(Int64, zdt))
    @test round_trip == zdt
    @test !isequal(round_trip, zdt)
end


text:
# features.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2014 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
export ngrams, count, tfnorm, sparse_count, norm, znorm, ngram_iterator, ngrams!

struct NgramStringIterator
  string :: AbstractString
  order :: Int32
  truncated_start :: Bool
end
Base.iteratorsize(::Type{NgramStringIterator}) = Base.SizeUnknown()

mutable struct StringPosition
  start  :: Int32
  fin    :: Int32
  nth    :: Int32
end

function start(ngi :: NgramStringIterator)
  if ngi.truncated_start
    idx = 1
    for i = 1:(ngi.order-1)
      idx = nextind(ngi.string, idx)
    end
    return StringPosition(1, idx, ngi.order)
  else
    return StringPosition(1, 1, 1)
  end
end

done(ngi :: NgramStringIterator, position) = position.fin > endof(ngi.string)
function next(ngi :: NgramStringIterator, position)
  str = make_string(ngi.string, position.start, position.fin)
  if position.nth >= ngi.order
    position.start = nextind(ngi.string, position.start)
  end
  position.nth += 1
  position.fin  = nextind(ngi.string, position.fin)
  return str, position
end

# -------------------------------------------------------------------------------------------------------------------------
# feature extractors
# -------------------------------------------------------------------------------------------------------------------------
make_string(words :: AbstractString, b, e) = SubString(words, b, e)
make_string(words :: Array, b, e) = join(words[b:e], " ")

function ngrams(words::Array; order = 2, truncated_start = false)
  ret = AbstractString[]

  if !truncated_start
    for wi = 1:min(order - 1, length(words))
      push!(ret, make_string(words, 1, wi))
    end
  end

  for wi = order:length(words)
    push!(ret, make_string(words, wi - order + 1, wi))
  end
  return ret
end

function ngrams(words :: AbstractString; order = 2, truncated_start = false)
  ret = AbstractString[]
  return ngrams!(ret, words, order = order, truncated_start = truncated_start)
end

function ngrams!(ret :: Array, words :: AbstractString; order = 2, truncated_start = false)
  for x in ngram_iterator(words, order = order, truncated_start = truncated_start)
    push!(ret, x)
  end
  return ret
end

ngram_iterator(words :: AbstractString; order = 2, truncated_start = false) = NgramStringIterator(words, order, truncated_start)

# -------------------------------------------------------------------------------------------------------------------------
# feature vector operations
# -------------------------------------------------------------------------------------------------------------------------
function sparse_count(text, bkg)
  vec = spzeros(vocab_size(bkg), 1)
  for token in text
    vec[bkg[token]] += 1
  end
  return vec
end

function dict_count(tokens)
  map = DefaultDict{AbstractString,Int32}()
  for w in tokens
    map[w] += 1
  end
end

function znorm(vector, mean, std)
  @devec ret = (vector - mean) ./ std
  return ret
end

function norm(counts, bkg; cutoff = 0.0)
  for r in counts.rowval
    counts[r] *= bkg.stats[r]
  end

  # prune
  nn = nonzeros(counts)
  for i = 1:length(nn)
    if nn[i] < cutoff
      nn[i] = 0.0
    end
  end

  return counts
end


text:
module Tests
using PenaltyFunctions, LinearAlgebra, Test
P = PenaltyFunctions

element_penalties = [NoPenalty(), L1Penalty(), L2Penalty(),
                     ElasticNetPenalty(.7), ElasticNetPenalty(.7f0),
                     MCPPenalty(1.), MCPPenalty(1f0), MCPPenalty(1),
                     LogPenalty(1.), LogPenalty(1f0),
                     SCADPenalty(3.7), SCADPenalty(3.7f0)]
array_penalties = [NuclearNormPenalty(), GroupLassoPenalty(),
                   MahalanobisPenalty(rand(3,2))]

penalty_list = vcat(element_penalties, array_penalties)

@info "Show methods:"
for p in penalty_list
    printstyled("  > $(string(p))\n";             color = :green)
    printstyled("  > $(string(scaled(p, .1)))\n"; color = :blue)
end

#---------------------------------------------------------------------------# Begin Tests
println("\n")
@info "Begin Actual Tests"
@testset "Common" begin
    @test P.soft_thresh(1.0, 0.5) == 0.5
    @test P.soft_thresh!(ones(5), .5) == .5 * ones(5)
    @test P.name(L1Penalty()) == "L1Penalty"
end

@testset "ElementPenalty" begin
    function test_element_penalty(p::Penalty, θ::Number, s::Number, v1, v2, v3)
        r(x) = Float32(round(x; digits=5)) # avoids floating point issues for different parameter types
        @testset "$(P.name(p))" begin
            @test r(@inferred(value(p, θ))) ≈ r(v1)
            @test r(@inferred(deriv(p, θ))) ≈ r(v2)
            @test r.(value.(Ref(p), fill(θ, 5))) ≈ r.(fill(v1, 5))
            @test r.(deriv.(Ref(p), fill(θ, 5))) ≈ r.(fill(v2, 5))
            @test value(p, θ) == p(θ)
            @test value(p, θ, s) == p(θ, s)
            if isa(p, P.ProxableElementPenalty)
                @test r(@inferred(prox(p, θ, s))) ≈ r(v3)
                @test r.(prox.(Ref(p), fill(θ, 5), Ref(s))) ≈ r.(fill(v3, 5))
            end
        end
    end

    θ, s = rand(), rand()

    @testset "Constructors" begin
        @inferred NoPenalty()
        @inferred L1Penalty()
        @inferred L2Penalty()

        @inferred ElasticNetPenalty()
        @inferred ElasticNetPenalty(1.0)

        @inferred SCADPenalty()
        @inferred SCADPenalty(4.0)
        @inferred SCADPenalty(4)
        @inferred SCADPenalty(4, 1.)
        @inferred SCADPenalty(4.0, 1.0)
        @test_throws Exception SCADPenalty(1.0)
        @test_throws Exception SCADPenalty(3.7, -1.)
        @test typeof(SCADPenalty())         <: SCADPenalty{Float64}
        @test typeof(SCADPenalty(4))        <: SCADPenalty{Int}
        @test typeof(SCADPenalty(4, 1))     <: SCADPenalty{Int}
        @test typeof(SCADPenalty(4f0))      <: SCADPenalty{Float32}
        @test typeof(SCADPenalty(4f0, 1))   <: SCADPenalty{Float32}
        @test typeof(SCADPenalty(4f0, 1.0)) <: SCADPenalty{Float64}
        @test typeof(SCADPenalty(4, 1.0))   <: SCADPenalty{Float64}
        @test typeof(SCADPenalty(4.0, 1.0)) <: SCADPenalty{Float64}

        @inferred LogPenalty()
        @inferred LogPenalty(1.0)
        @test_throws Exception LogPenalty(-1.0)

        @inferred MCPPenalty()
        @inferred MCPPenalty(1.0)
        @test_throws Exception MCPPenalty(-1.0)
    end

    for T in (Float32, Float64), S in (Float32, Float64)
        test_element_penalty(NoPenalty(), T(θ), S(s), 0.0, 0.0, θ)
        test_element_penalty(L1Penalty(), T(θ), S(s), abs(θ), sign(θ), max(0.0, θ - sign(θ) * s))
        test_element_penalty(L2Penalty(), T(θ), S(s), .5 * θ^2, θ, θ / (1 + s))
        test_element_penalty(ElasticNetPenalty(.4), T(θ), S(s),
            .4 * value(L1Penalty(), θ) + .6 * value(L2Penalty(), θ),
            .4 * deriv(L1Penalty(), θ) + .6 * deriv(L2Penalty(), θ),
            prox(L2Penalty(), prox(L1Penalty(), θ, .4s), .6s)
        )

        test_element_penalty(LogPenalty(1.0), T(θ), S(s), log(1 + θ), 1 / (1 + θ), nothing)
        test_element_penalty(MCPPenalty(2.0), T(1.0), S(s), 1 - 1/4, 1 - 1/2, nothing)
        test_element_penalty(MCPPenalty(1.0), T(2.0), S(s), 1/2, 0.0, nothing)
    end

    @testset "SCADPenalty" begin
        for T in (Float32, Float64), S in (Float32, Float64)
            r(x) = round(x, digits=8)
            @test r(@inferred(value(SCADPenalty(T(3.8), T(.2)), S(.1)))) ≈ r(.02)
            @test r(@inferred(value(SCADPenalty(T(3.8), T(.1)), S(.2)))) ≈ r(-.5 * (.2^2 - .2^2 * 3.8 + .01) / (2.8))
            @test r(@inferred(value(SCADPenalty(T(3.8), T(.1)), S(9.)))) ≈ r(.5 * 4.8 * .01)
            @test r(@inferred(value(SCADPenalty(T(3.8), T(.1)), 9)))     ≈ r(.5 * 4.8 * .01)

            @test r(@inferred(deriv(SCADPenalty(T(3.8), T(.2)), S(.1)))) ≈ r(.2)
            @test r(@inferred(deriv(SCADPenalty(T(3.8), T(.1)), S(.2)))) ≈ r(.1 * (3.8 * .1 - .2) / (2.8 * .1))
            @test r(@inferred(deriv(SCADPenalty(T(3.8), T(.1)), S(9.)))) ≈ r(0.)
            @test r(@inferred(deriv(SCADPenalty(T(3.8), T(.1)), 9)))     ≈ r(0.)
        end
    end

    @testset "ElementPenalty methods" begin
        p = L1Penalty()
        for T in (Float32, Float64), S in (Float32, Float64)
            θ, s = rand(T, 10), rand(S, 10)
            @testset "value: $T, $S" begin
                @test round(@inferred(value(p, θ, s[1])), digits=4) ≈ round(s[1] * sum(abs, θ), digits=4)
                @test @inferred(value(p, θ))       ≈ sum(abs, θ)
                @test @inferred(value(p, θ, s))    ≈ sum(s .* abs.(θ))
            end
            @testset "deriv/grad: $T, $S" begin
                @test @inferred(deriv(p, θ[1], s[1]))  ≈ s[1] * sign(θ[1])
                @test @inferred(grad(p, θ))            ≈ sign.(θ)
                @test @inferred(grad(p, θ, s[1]))      ≈ s[1] * sign.(θ)
                @test @inferred(grad(p, θ, s))         ≈ s .* sign.(θ)
                @test deriv.(Ref(p), θ)       == grad(p, θ)
                @test deriv.(Ref(p), θ, s[1]) == grad(p, θ, s[1])
                @test deriv.(Ref(p), θ, s)    == grad(p, θ, s)

                buffer = rand(10)
                grad!(buffer, p, θ); @test buffer       ≈ sign.(θ)
                grad!(buffer, p, θ, s[1]); @test buffer ≈ s[1] * sign.(θ)
                grad!(buffer, p, θ, s); @test buffer    ≈ s .* sign.(θ)

                addgrad(buffer[1], p, θ[1]) ≈ buffer[1] + sign(θ[1])
                addgrad(buffer[1], p, θ[1], s[1]) ≈ buffer[1] + s[1] * sign(θ[1])
                ∇ = rand(10)
                ∇2 = copy(∇)
                addgrad!(∇, p, θ); @test ∇ ≈ ∇2 + sign.(θ)
                ∇ = rand(10)
                ∇2 = copy(∇)
                addgrad!(∇, p, θ, s[1]); @test ∇ ≈ ∇2 + s[1] * sign.(θ)
                ∇ = rand(10)
                ∇2 = copy(∇)
                addgrad!(∇, p, θ, s); @test ∇ ≈ ∇2 + s .* sign.(θ)
            end
            @testset "prox: $T, $S" begin
                p = L1Penalty()
                θ2 = copy(θ)
                prox!(p, θ, s[1]); @test θ ≈ P.soft_thresh.(θ2, s[1])
                θ = copy(θ2)
                prox!(p, θ, s); @test θ ≈ P.soft_thresh.(θ2, s)

                θ = rand(5)
                s = rand(5)
                @test @inferred(prox(p, θ, s)) == map((x,y) -> prox(p, x, y), θ, s)
            end
        end
    end
end

@testset "ScaledElementPenalty" begin
    for p in element_penalties
        s = scaled(p, .1)
        x = randn(5)
        # FIXME: @inference broken during tests (not in REPL though)
        #        only for ElasticNetPenalty{Float64} for some reason
        @test value(s, x) ≈ value(p, x, .1)
        @test value(s, x) == value(.1 * p, x)
        @test @inferred(deriv(s, x[1])) ≈ deriv(p, x[1], .1)
        @test @inferred(grad(s, x))     ≈ grad(p, x, .1)
        if typeof(p) <: P.ProxableElementPenalty
            @test @inferred(prox(s, x)) ≈ prox(p, x, .1)
        end
    end

    p = ElasticNetPenalty(.7)
    s = scaled(p, .2)
    x = randn(5)
    # FIXME: @inference broken during tests (not in REPL though)
    @test value(s, x) ≈ .2 * value(p, x)
    @test @inferred(deriv(s, x[1])) ≈ .2 * deriv(p, x[1])
    @test @inferred(grad(s, x)) ≈ .2 * grad(p, x)
    @test deriv.(Ref(s), x) ≈ .2 * deriv.(Ref(p), x)
    @test @inferred(prox(s, x)) ≈ prox(p, x, .2)
    @test @inferred(prox(s, x[1])) ≈ prox(p, x[1], .2)

    @test_throws ArgumentError scaled(p, -1.)
end

@testset "ArrayPenalty" begin
    @testset "NuclearNormPenalty" begin
        p = NuclearNormPenalty()
        Θ = randn(10, 5)
        s = .05
        # FIXME: @inference broken. seems like a type instability
        @test value(p, Θ) ≈ sum(svd(Θ).S)
        @test value(p, Θ) == p(Θ)
        @test value(p, Θ, s) ≈ s * sum(svd(Θ).S)
        prox!(p, Θ, s)
    end
    @testset "GroupLassoPenalty" begin
        p = GroupLassoPenalty()
        Θ = randn(10)
        s = .05
        @test @inferred(value(p, Θ)) ≈ norm(Θ)
        @test value(p, Θ) == p(Θ)
        prox!(p, Θ, s)

        Θ = .01 * ones(10)
        prox!(p, Θ, 10.) == zeros(10)
    end
    @testset "MahalanobisPenalty" begin
        C = randn(5, 10)
        p = MahalanobisPenalty(C)
        θ = rand(10)
        @test value(p, θ) == p(θ)
        s = .05
        @test @inferred(value(p, θ)) ≈ 0.5 * dot(C * θ, C * θ)
        prox!(p, θ, s)
    end
    @testset "ScaledArrayPenalty" begin
        p = GroupLassoPenalty()
        s = scaled(p, .1)
        Θ = randn(10)
        @test value(p, Θ) == p(Θ)
        @test @inferred(value(p, Θ, .1)) ≈ value(s, Θ)

        Θ2 = copy(Θ)
        prox!(p, Θ, .1); prox!(s, Θ2); @test Θ ≈ Θ2
    end
end

end  # module


text:
import Base: show

"""
# SampleModel

Create a SampleModel.

### Required arguments
```julia
* `name::AbstractString`               : Name for the model
* `model::AbstractString`              : Stan model source
```

### Optional arguments
```julia
* `n_chains::Vector{Int64}=[4]`        : Optionally updated in stan_sample()
* `seed::RandomSeed`                   : Random seed settings
* `output::Output`                     : File output options
* `init::Init`                         : Default interval bound for parameters
* `tmpdir::AbstractString`             : Directory where output files are stored
* `output_base::AbstractString`        : Base name for output files
* `exec_path::AbstractString`          : Path to cmdstan executable
* `data_file::vector{AbstractString}`  : Path to per chain data file
* `init_file::Vector{AbstractString}`  : Path to per chain init file
* `cmds::Vector{Cmd}`                  : Path to per chain init file
* `sample_file::Vector{String}         : Path to per chain samples file
* `log_file::Vector{String}            : Path to per chain log file
* `diagnostic_file::Vector{String}     : Path to per chain diagnostic file
* `summary=true`                       : Create computed stan summary
* `printsummary=true`                  : Show computed stan summary
* `sm::StanRun.StanModel`              : StanRun.StanModel
* `method::Sample                      : Will be Sample()
```

"""
mutable struct SampleModel <: CmdStanModels
  @shared_fields_stanmodels
  method::Sample
end

function SampleModel(
  name::AbstractString,
  model::AbstractString,
  n_chains=[4];
  method = Sample(),
  seed = StanBase.RandomSeed(),
  init = StanBase.Init(),
  output = StanBase.Output(),
  tmpdir = mktempdir())

  !isdir(tmpdir) && mkdir(tmpdir)

  StanBase.update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))
  sm = StanModel(joinpath(tmpdir, "$(name).stan"))

  output_base = StanRun.default_output_base(sm)
  exec_path = StanRun.ensure_executable(sm)

  stan_compile(sm)

  SampleModel(name, model, n_chains, seed, init, output,
    tmpdir, output_base, exec_path, String[], String[],
    Cmd[], String[], String[], String[], false, false, sm, method)
end

function Base.show(io::IO, ::MIME"text/plain", m::SampleModel)
  println(io, "  name =                    \"$(m.name)\"")
  println(io, "  n_chains =                $(StanBase.get_n_chains(m))")
  println(io, "  output =                  Output()")
  println(io, "    file =                    \"$(split(m.output.file, "/")[end])\"")
  println(io, "    diagnostics_file =        \"$(split(m.output.diagnostic_file, "/")[end])\"")
  println(io, "    refresh =                 $(m.output.refresh)")
  println(io, "  tmpdir =                  \"$(m.tmpdir)\"")
  println(io, "  method =                  Sample()")
  println(io, "    num_samples =             ", m.method.num_samples)
  println(io, "    num_warmup =              ", m.method.num_warmup)
  println(io, "    save_warmup =             ", m.method.save_warmup)
  println(io, "    thin =                    ", m.method.thin)
  if isa(m.method.algorithm, Hmc)
    println(io, "    algorithm =               HMC()")
    if isa(m.method.algorithm.engine, Nuts)
      println(io, "      engine =                  NUTS()")
      println(io, "        max_depth =               ", m.method.algorithm.engine.max_depth)
    elseif isa(m.method.algorithm.engine, Static)
      println(io, "      engine =                  Static()")
      println(io, "        int_time =                ", m.method.algorithm.engine.int_time)
    end
    println(io, "      metric =                  ", typeof(m.method.algorithm.metric))
    println(io, "      stepsize =                ", m.method.algorithm.stepsize)
    println(io, "      stepsize_jitter =         ", m.method.algorithm.stepsize_jitter)
  else
    if isa(m.method.algorithm, Fixed_param)
      println(io, "    algorithm =               Fixed_param()")
    else
      println(io, "    algorithm =               Unknown")
    end
  end
  println(io, "    adapt =                   Adapt()")
  println(io, "      gamma =                   ", m.method.adapt.gamma)
  println(io, "      delta =                   ", m.method.adapt.delta)
  println(io, "      kappa =                   ", m.method.adapt.kappa)
  println(io, "      t0 =                      ", m.method.adapt.t0)
  println(io, "      init_buffer =             ", m.method.adapt.init_buffer)
  println(io, "      term_buffer =             ", m.method.adapt.term_buffer)
  println(io, "      window =                  ", m.method.adapt.window)
end