text:
struct AdjointDiffCache{UF,PF,G,TJ,PJT,uType,JC,GC,PJC,JNC,PJNC,rateType,DG,DI,AI,FM}
  uf::UF
  pf::PF
  g::G
  J::TJ
  pJ::PJT
  dg_val::uType
  jac_config::JC
  g_grad_config::GC
  paramjac_config::PJC
  jac_noise_config::JNC
  paramjac_noise_config::PJNC
  f_cache::rateType
  dg::DG
  diffvar_idxs::DI
  algevar_idxs::AI
  factorized_mass_matrix::FM
  issemiexplicitdae::Bool
end

"""
    adjointdiffcache(g,sensealg,discrete,sol,dg;quad=false)

return (AdjointDiffCache, y)
"""
function adjointdiffcache(g::G,sensealg,discrete,sol,dg::DG,f;quad=false,noiseterm=false,needs_jac=false) where {G,DG}
  prob = sol.prob
  if prob isa DiffEqBase.SteadyStateProblem
    @unpack u0, p = prob
    tspan = (nothing, nothing)
  #elseif prob isa SDEProblem
  #  @unpack tspan, u0, p = prob
  else
    @unpack u0, p, tspan = prob
  end
  numparams = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(p)
  numindvar = length(u0)
  isautojacvec = get_jacvec(sensealg)

  issemiexplicitdae = false
  mass_matrix = sol.prob.f.mass_matrix
  if mass_matrix isa UniformScaling
    factorized_mass_matrix = mass_matrix'
  elseif mass_matrix isa Tuple{UniformScaling,UniformScaling}
    factorized_mass_matrix = (I',I')
  else
    mass_matrix = mass_matrix'
    diffvar_idxs = findall(x->any(!iszero, @view(mass_matrix[:, x])), axes(mass_matrix, 2))
    algevar_idxs = setdiff(eachindex(u0), diffvar_idxs)
    # TODO: operator
    M̃ = @view mass_matrix[diffvar_idxs, diffvar_idxs]
    factorized_mass_matrix = lu(M̃, check=false)
    issuccess(factorized_mass_matrix) || error("The submatrix corresponding to the differential variables of the mass matrix must be nonsingular!")
    isempty(algevar_idxs) || (issemiexplicitdae = true)
  end
  if !issemiexplicitdae
    diffvar_idxs = eachindex(u0)
    algevar_idxs = 1:0
  end

  if !needs_jac
    J = (issemiexplicitdae || !isautojacvec) ? similar(u0, numindvar, numindvar) : nothing
  else
    # Force construction of the Jacobian
    J = similar(u0, numindvar, numindvar)
  end

  if !discrete
    if dg !== nothing
      pg = nothing
      pg_config = nothing
      if dg isa Tuple && length(dg) == 2
        dg_val = (similar(u0, numindvar),similar(u0, numparams))
        dg_val[1] .= false
        dg_val[2] .= false
      else
        dg_val = similar(u0, numindvar) # number of funcs size
        dg_val .= false
      end
    else
      if !(prob isa RODEProblem)
        pg = UGradientWrapper(g,tspan[2],p)
      else
        pg = RODEUGradientWrapper(g,tspan[2],p,last(sol.W))
      end
      pg_config = build_grad_config(sensealg,pg,u0,p)
      dg_val = similar(u0, numindvar) # number of funcs size
      dg_val .= false
    end
  else
    dg_val = nothing
    pg = nothing
    pg_config = nothing
  end

  if DiffEqBase.has_jac(f) || (J === nothing)
    jac_config = nothing
    uf = nothing
  else
    if DiffEqBase.isinplace(prob)
      if !(prob isa RODEProblem)
        uf = DiffEqBase.UJacobianWrapper(f,tspan[2],p)
      else
        uf = RODEUJacobianWrapper(f,tspan[2],p,last(sol.W))
      end
      jac_config = build_jac_config(sensealg,uf,u0)
    else
      if !(prob isa RODEProblem)
        uf = DiffEqBase.UDerivativeWrapper(f,tspan[2],p)
      else
        uf = RODEUDerivativeWrapper(f,tspan[2],p,last(sol.W))
      end
      jac_config = nothing
    end
  end

  if prob isa DiffEqBase.SteadyStateProblem
    y = copy(sol.u)
  else
    y = copy(sol.u[end])
  end

  if typeof(prob.p) <: DiffEqBase.NullParameters
      _p = similar(y,(0,))
  else
      _p = prob.p
  end

  if sensealg.autojacvec isa ReverseDiffVJP ||
    (sensealg.autojacvec isa Bool && sensealg.autojacvec && DiffEqBase.isinplace(prob))

    if prob isa DiffEqBase.SteadyStateProblem
       if DiffEqBase.isinplace(prob)
        tape = ReverseDiff.GradientTape((y, _p)) do u,p
          du1 = p !== nothing && p !== DiffEqBase.NullParameters() ? similar(p, size(u)) : similar(u)
          f(du1,u,p,nothing)
          return vec(du1)
        end
      else
        tape = ReverseDiff.GradientTape((y, _p)) do u,p
          vec(f(u,p,nothing))
        end
      end
    elseif  noiseterm && (!StochasticDiffEq.is_diagonal_noise(prob) || isnoisemixing(sensealg))
      tape = nothing
    else
      if DiffEqBase.isinplace(prob)
        if !(prob isa RODEProblem)
          tape = ReverseDiff.GradientTape((y, _p, [tspan[2]])) do u,p,t
            du1 = (p !== nothing && p !== DiffEqBase.NullParameters()) ? similar(p, size(u)) : similar(u)
            f(du1,u,p,first(t))
            return vec(du1)
          end
        else
          tape = ReverseDiff.GradientTape((y, _p, [tspan[2]],last(sol.W))) do u,p,t,W
            du1 = p !== nothing && p !== DiffEqBase.NullParameters() ? similar(p, size(u)) : similar(u)
            f(du1,u,p,first(t),W)
            return vec(du1)
          end
        end
      else
        if !(prob isa RODEProblem)
          tape = ReverseDiff.GradientTape((y, _p, [tspan[2]])) do u,p,t
            vec(f(u,p,first(t)))
          end
        else
          tape = ReverseDiff.GradientTape((y, _p, [tspan[2]],last(sol.W))) do u,p,t,W
            return f(u,p,first(t),W)
          end
        end
      end
    end

    if compile_tape(sensealg.autojacvec)
      paramjac_config = ReverseDiff.compile(tape)
    elseif tape !== nothing && sensealg.autojacvec isa Bool && sensealg.autojacvec && DiffEqBase.isinplace(prob)
      compile = try
          !hasbranching(prob.f,copy(u0),u0,p,prob.tspan[2])
      catch
          false
      end
      if compile
          paramjac_config = ReverseDiff.compile(tape)
      else
          paramjac_config = tape
      end
    else
      paramjac_config = tape
    end

    pf = nothing
  elseif sensealg.autojacvec isa EnzymeVJP
    if typeof(prob.p) <: DiffEqBase.NullParameters
      paramjac_config = zero(y),prob.p,zero(y),zero(y)
    else
      paramjac_config = zero(y),zero(_p),zero(y),zero(y)
    end
    pf = let f = f.f
        if DiffEqBase.isinplace(prob) && prob isa RODEProblem
            function (out,u,_p,t,W)
                f(out, u, _p, t, W)
                nothing
            end
        elseif DiffEqBase.isinplace(prob)
            function (out,u,_p,t)
                f(out, u, _p, t)
                nothing
            end
        elseif !DiffEqBase.isinplace(prob) && prob isa RODEProblem
            function (out,u,_p,t,W)
                out .= f(u, _p, t, W)
                nothing
            end
        else !DiffEqBase.isinplace(prob)
            function (out,u,_p,t)
                out .= f(u, _p, t)
                nothing
            end
        end
    end
  elseif (DiffEqBase.has_paramjac(f) || isautojacvec || quad)
    paramjac_config = nothing
    pf = nothing
  else
    if DiffEqBase.isinplace(prob)
      if !(prob isa RODEProblem)
        pf = DiffEqBase.ParamJacobianWrapper(f,tspan[1],y)
      else
        pf = RODEParamJacobianWrapper(f,tspan[1],y,last(sol.W))
      end
      paramjac_config = build_param_jac_config(sensealg,pf,y,p)
    else
      if !(prob isa RODEProblem)
        pf = ParamGradientWrapper(f,tspan[2],y)
      else
        pf = RODEParamGradientWrapper(f,tspan[2],y,last(sol.W))
      end
      paramjac_config = nothing
    end
  end

  pJ = (quad || isautojacvec) ? nothing : similar(u0, numindvar, numparams)

  f_cache = DiffEqBase.isinplace(prob) ? deepcopy(u0) : nothing

  if noiseterm
    if (sensealg.noise isa ReverseDiffNoise ||
      (sensealg.noise isa Bool && sensealg.noise && DiffEqBase.isinplace(prob)))

      jac_noise_config = nothing
      paramjac_noise_config = []

      if DiffEqBase.isinplace(prob)
        for i in 1:numindvar
          function noisetape(indx)
            if StochasticDiffEq.is_diagonal_noise(prob)
              ReverseDiff.GradientTape((y, _p, [tspan[2]])) do u,p,t
                du1 = p !== nothing && p !== DiffEqBase.NullParameters() ? similar(p, size(u)) : similar(u)
                f(du1,u,p,first(t))
                return du1[indx]
              end
            else
              ReverseDiff.GradientTape((y, _p, [tspan[2]])) do u,p,t
                du1 = similar(p, size(prob.noise_rate_prototype))
                du1 .= false
                f(du1,u,p,first(t))
                return du1[:,indx]
              end
            end
          end
          tapei = noisetape(i)
          if compile_tape(sensealg.noise)
            push!(paramjac_noise_config, ReverseDiff.compile(tapei))
          elseif tapei !== nothing && sensealg.noise isa Bool && sensealg.noise && DiffEqBase.isinplace(prob)
              compile = try
                  !hasbranching(prob.f,copy(u0),u0,p,prob.tspan[2])
              catch
                  false
              end
              if compile
                  push!(paramjac_noise_config, ReverseDiff.compile(tapei))
              else
                  push!(paramjac_noise_config, tapei)
              end
          else
            push!(paramjac_noise_config, tapei)
          end
        end
      else
        for i in 1:numindvar
          function noisetapeoop(indx)
            if StochasticDiffEq.is_diagonal_noise(prob)
              ReverseDiff.GradientTape((y, _p, [tspan[2]])) do u,p,t
                f(u,p,first(t))[indx]
              end
            else
              ReverseDiff.GradientTape((y, _p, [tspan[2]])) do u,p,t
                f(u,p,first(t))[:,indx]
              end
            end
          end
          tapei = noisetapeoop(i)
          if compile_tape(sensealg.noise)
            push!(paramjac_noise_config, ReverseDiff.compile(tapei))
          else
            push!(paramjac_noise_config, tapei)
          end
        end
      end
    elseif (sensealg.noise isa Bool && !sensealg.noise)
      if DiffEqBase.isinplace(prob)
        if StochasticDiffEq.is_diagonal_noise(prob)
          pf = DiffEqBase.ParamJacobianWrapper(f,tspan[1],y)
          if isnoisemixing(sensealg)
            uf = DiffEqBase.UJacobianWrapper(f,tspan[2],p)
            jac_noise_config = build_jac_config(sensealg,uf,u0)
          else
            jac_noise_config = nothing
          end
        else
          pf = ParamNonDiagNoiseJacobianWrapper(f,tspan[1],y,prob.noise_rate_prototype)
          uf = UNonDiagNoiseJacobianWrapper(f,tspan[2],p,prob.noise_rate_prototype)
          jac_noise_config = build_jac_config(sensealg,uf,u0)
        end
        paramjac_noise_config = build_param_jac_config(sensealg,pf,y,p)
      else
        if StochasticDiffEq.is_diagonal_noise(prob)
          pf = ParamGradientWrapper(f,tspan[2],y)
          if isnoisemixing(sensealg)
            uf = DiffEqBase.UDerivativeWrapper(f,tspan[2],p)
          end
        else
          pf = ParamNonDiagNoiseGradientWrapper(f,tspan[1],y)
          uf = UNonDiagNoiseGradientWrapper(f,tspan[2],p)
        end
        paramjac_noise_config = nothing
        jac_noise_config = nothing
      end
      if StochasticDiffEq.is_diagonal_noise(prob)
        pJ = similar(u0, numindvar, numparams)
        if isnoisemixing(sensealg)
          J = similar(u0, numindvar, numindvar)
        end
      else
        pJ = similar(u0, numindvar*numindvar, numparams)
        J = similar(u0, numindvar*numindvar, numindvar)
      end

    else
      paramjac_noise_config = nothing
      jac_noise_config = nothing
    end
  else
    paramjac_noise_config = nothing
    jac_noise_config = nothing
  end

  adjoint_cache = AdjointDiffCache(uf,pf,pg,J,pJ,dg_val,
                          jac_config,pg_config,paramjac_config,
                          jac_noise_config,paramjac_noise_config,
                          f_cache,dg,diffvar_idxs,algevar_idxs,
                          factorized_mass_matrix,issemiexplicitdae)

  return adjoint_cache, y
end

getprob(S::SensitivityFunction) = (S isa ODEBacksolveSensitivityFunction) ? S.prob : S.sol.prob
inplace_sensitivity(S::SensitivityFunction) = isinplace(getprob(S))

struct ReverseLossCallback{λType,timeType,yType,RefType,FMType,AlgType,gType,cacheType}
  isq::Bool
  λ::λType
  t::timeType
  y::yType
  cur_time::RefType
  idx::Int
  F::FMType
  sensealg::AlgType
  g::gType
  diffcache::cacheType
end

function ReverseLossCallback(sensefun, λ, t, g, cur_time)
  @unpack sensealg, y = sensefun
  isq = (sensealg isa QuadratureAdjoint)

  @unpack factorized_mass_matrix = sensefun.diffcache
  prob = getprob(sensefun)
  idx = length(prob.u0)

  return ReverseLossCallback(isq, λ, t, y, cur_time, idx, factorized_mass_matrix, sensealg, g, sensefun.diffcache)
end

function (f::ReverseLossCallback)(integrator)
  @unpack isq, λ, t, y, cur_time, idx, F, sensealg, g = f
  @unpack diffvar_idxs, algevar_idxs, issemiexplicitdae, J, uf, f_cache, jac_config = f.diffcache

  p, u = integrator.p, integrator.u

  if sensealg isa BacksolveAdjoint
    copyto!(y,integrator.u[end-idx+1:end])
  end

  # Warning: alias here! Be careful with λ
  gᵤ = isq ? λ : @view(λ[1:idx])
  g(gᵤ,y,p,t[cur_time[]],cur_time[])

  if issemiexplicitdae
    jacobian!(J, uf, y, f_cache, sensealg, jac_config)
    dhdd = J[algevar_idxs, diffvar_idxs]
    dhda = J[algevar_idxs, algevar_idxs]
    # TODO: maybe need a `conj`
    Δλa = -dhda'\gᵤ[algevar_idxs]
    Δλd = dhdd'Δλa + gᵤ[diffvar_idxs]
  else
    Δλd = gᵤ
  end

  if F !== nothing
    F !== I && F !== (I,I) && ldiv!(F, Δλd)
  end

  u[diffvar_idxs] .+= Δλd
  u_modified!(integrator,true)
  cur_time[] -= 1
  return nothing
end

function generate_callbacks(sensefun, g, λ, t, t0, callback, init_cb,terminated=false)
  if sensefun isa NILSASSensitivityFunction
    @unpack sensealg = sensefun.S
  else
    @unpack sensealg = sensefun
  end

  if !sensefun.discrete
    cur_time = Ref(1)
  else
    cur_time = Ref(length(t))
  end

  reverse_cbs = setup_reverse_callbacks(callback,sensealg,g,cur_time,terminated)
  sensefun.discrete || return reverse_cbs, nothing

  # callbacks can lead to non-unique time points
  _t, duplicate_iterator_times = separate_nonunique(t)

  rlcb = ReverseLossCallback(sensefun, λ, t, g, cur_time)

  if eltype(_t) !== typeof(t0)
    _t = convert.(typeof(t0),_t)
  end
  cb = PresetTimeCallback(_t,rlcb)

  # handle duplicates (currently only for double occurances)
  if duplicate_iterator_times!==nothing
    # use same ref for cur_time to cope with concrete_solve
    cbrev_dupl_affect = ReverseLossCallback(sensefun, λ, t, g, cur_time)
    cb_dupl = PresetTimeCallback(duplicate_iterator_times[1],cbrev_dupl_affect)
    return CallbackSet(cb,reverse_cbs,cb_dupl), duplicate_iterator_times
  else
    return CallbackSet(cb,reverse_cbs), duplicate_iterator_times
  end
end

function separate_nonunique(t)
  # t is already sorted
  _t = unique(t)
  ts_with_occurances = [(i, count(==(i), t)) for i in _t]

  # duplicates (only those values which occur > 1 times)
  dupl = filter(x->last(x)>1, ts_with_occurances)

  ts = first.(dupl)
  occurances = last.(dupl)

  if isempty(occurances)
    itrs = nothing
  else
    maxoc = maximum(occurances)
    maxoc > 2 && error("More than two occurances of the same time point. Please report this.")
    # handle also more than two occurances
    itrs = [ts[occurances .>= i] for i=2:maxoc]
  end

  return _t, itrs
end

function out_and_ts(_ts, duplicate_iterator_times, sol)
  if duplicate_iterator_times === nothing
    ts = _ts
    out = sol(ts)
  else
    # if callbacks are tracked, there is potentially an event_time that must be considered
    # in the loss function but doesn't occur in saveat/t. So we need to add it.
    # Note that if it doens't occur in saveat/t we even need to add it twice
    # However if the callbacks are not saving in the forward, we don't want to compute a loss
    # value for them. This information is given by sol.t/checkpoints.
    # Additionally we need to store the left and the right limit, respectively.
    duplicate_times = duplicate_iterator_times[1] # just treat two occurances at the moment (see separate_nonunique above)
    _ts = Array(_ts)
    for d in duplicate_times
      (d ∉ _ts) && push!(_ts, d)
    end

    u1 = sol(_ts).u
    u2 = sol(duplicate_times,continuity=:right).u
    saveat = vcat(_ts,  duplicate_times...)
    perm = sortperm(saveat)
    ts = saveat[perm]
    u = vcat(u1, u2)[perm]
    out = DiffEqArray(u,ts)
  end
  return out, ts
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using IterativeLQR
using LinearAlgebra

# ## horizon
T = 11

# ## acrobot
num_state = 2
num_action = 1

function particle(x, u, w)
   A = [1.0 1.0; 0.0 1.0]
   B = [0.0; 1.0]
   return A * x + B * u[1]
end

# ## model
dynamics = Dynamics(particle, num_state, num_action)
model = [dynamics for t = 1:T-1]

# ## initialization
x1 = [0.0; 0.0]
xT = [1.0; 0.0]
ū = [1.0e-1 * randn(num_action) for t = 1:T-1]
x̄ = rollout(model, x1, ū)

# ## objective
ot = (x, u, w) -> 0.1 * dot(x, x) + 0.1 * dot(u, u)
oT = (x, u, w) -> 0.1 * dot(x, x)
ct = Cost(ot, num_state, num_action)
cT = Cost(oT, num_state, 0)
objective = [[ct for t = 1:T-1]..., cT]

# ## constraints
goal(x, u, w) = x - xT

cont = Constraint()
conT = Constraint(goal, num_state, 0)
# conT = Constraint(
#       (c, x, u, w) -> c .= x - xT,
#       (cx, x, u, w) -> cx .= 1.0 * I(length(x)),
#       (cu, x, u, w) -> cu .= zeros(length(x), length(u)),
#       length(xT), num_state, num_action)
constraints = [[cont for t = 1:T-1]..., conT]

# ## problem
prob = Solver(model, objective, constraints)
initialize_controls!(prob, ū)
initialize_states!(prob, x̄)

# ## solve
solve!(prob)

# ## solution
x_sol, u_sol = get_trajectory(prob)

ext: jl


text:
# following tests require DiscreteValueIteration
POMDPs.add("DiscreteValueIteration")
using DiscreteValueIteration

mdp = GridWorld()

pomdp = FullyObservablePOMDP(mdp)

@test observations(pomdp) == states(pomdp)
@test n_observations(pomdp) == n_states(pomdp)
@test state_type(pomdp) == obs_type(pomdp)

s_po = initial_state(pomdp, MersenneTwister(1))
s_mdp = initial_state(mdp, MersenneTwister(1))
@test s_po == s_mdp

solver = ValueIterationSolver(max_iterations = 100)
mdp_policy = solve(solver, mdp, verbose=false)
pomdp_policy = solve(solver, pomdp, verbose=false)
@test mdp_policy.util == pomdp_policy.util


text:
col_1_to_30 = [-2.5   0.5         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0         0.0        0.0         0.0
1.0  -3.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0         0.0        0.0         0.0
0.0   0.0        -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.728864
0.0   0.0         0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       2.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         1.05443
0.0   0.0         0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        2.0        0.0        0.0         0.0         0.0         0.0        0.0         0.564472
0.0   0.0         0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.0        0.0         0.0         0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        2.0         0.0         0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         2.0         0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         2.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        2.0         0.0
0.0   0.0         0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.91546  -1.1331     0.217645   0.0        0.0         0.0         0.0         0.0        0.0         0.5
0.0   0.0         0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.21772  -2.19645   -0.021269   0.0        0.0         0.0         0.0         0.0        0.0         0.0
-0.0  -0.0         0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   2.25789  -0.845166   0.587279   0.0        0.0         0.0         0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546    -1.1331     0.217645    0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772    -2.19645   -0.021269    0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0       0.0        0.0        0.0        0.0         0.0         2.25789    -0.845166   0.587279    0.0
0.0   0.0         1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -3.50473   0.0        0.0        3.96463   -1.84856    -0.611345    0.0         0.0        0.0         1.45773
0.0   0.0         0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0      -3.50473   -1.55879    8.07172   -3.76354    -1.24466     0.0         0.0        0.0         2.10886
0.0   0.0         0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       1.55879   -3.50473   -0.142449   0.0664188   0.0219656   0.0         0.0        0.0         1.12894
0.0   0.0         0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0       -6.51418    0.0         0.0        11.8939     -5.54569   -1.83404     0.0
0.0   0.0         0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0       -6.51418    -4.67637    24.2152    -11.2906    -3.73397     0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        4.67637    -6.51418    -0.427348    0.199257   0.0658969   0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0        -2.56427     0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0        -2.56427   -0.584547    0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.584547  -2.56427     0.0
0.0   2.25709     0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.91546  -1.1331     0.217645   0.0        0.0         0.0         0.0         0.0        0.0        -5.25709
0.0   4.59528     0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.21772  -2.19645   -0.021269   0.0        0.0         0.0         0.0         0.0        0.0         0.0
-0.0  -0.0810972   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   2.25789  -0.845166   0.587279   0.0        0.0         0.0         0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0         0.0        0.0        17.8409
0.0   0.0         0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0         0.0        0.0        36.3227
0.0   0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0         0.0        0.0        -0.641022
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546    -1.1331     0.217645    0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772    -2.19645   -0.021269    0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0       0.0        0.0        0.0        0.0         0.0         2.25789    -0.845166   0.587279    0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.0
0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.0]

col_31_to_40 = [   0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-0.283275   0.0544112    0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-0.549114  -0.00531726   0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-0.211292   0.14682      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.728864   -0.283275    0.0544112    0.0        0.0        0.0          0.0   0.0
0.0        0.0          1.05443    -0.549114   -0.00531726   0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.564472   -0.211292    0.14682      0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.0          0.728864  -0.283275   0.0544112    0.0   0.0
0.0        0.0          0.0         0.0         0.0          1.05443   -0.549114  -0.00531726   0.0   0.0
0.0        0.0          0.0         0.0         0.0          0.564472  -0.211292   0.14682     -0.0  -0.0
0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.5        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.5          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.5         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.5         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.5          0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.0          0.5        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.0          0.0        0.5        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.5          0.0   0.0
-0.56655    0.108822     0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-1.09823   -0.0106345    0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-0.422583   0.29364      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0        0.0          1.45773    -0.56655     0.108822     0.0        0.0        0.0          0.0   0.0
0.0        0.0          2.10886    -1.09823    -0.0106345    0.0        0.0        0.0          0.0   0.0
0.0        0.0          1.12894    -0.422583    0.29364      0.0        0.0        0.0          0.0   0.0
0.0        0.0          0.0         0.0         0.0          1.45773   -0.56655    0.108822     0.0   0.564272
0.0        0.0          0.0         0.0         0.0          2.10886   -1.09823   -0.0106345    0.0   1.14882
0.0        0.0          0.0         0.0         0.0          1.12894   -0.422583   0.29364     -0.0  -0.0202743
0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-5.25709   -2.33819      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
2.33819   -5.25709      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-8.31853   -2.75105     -9.77127     0.0         0.0          0.0        0.0        0.0          0.0   0.0
-16.9359    -5.60095      0.0        -9.77127    -7.01456      0.0        0.0        0.0          0.0   0.0
0.298885   0.0988454    0.0         7.01456    -9.77127      0.0        0.0        0.0          0.0   0.0
0.0        0.0          2.23011    -1.03982    -0.343882    -3.84641    0.0        0.0          0.0   0.0
0.0        0.0          4.54034    -2.11699    -0.700119     0.0       -3.84641   -0.87682      0.0   0.0
0.0        0.0         -0.0801277   0.0373606   0.0123557    0.0        0.87682   -3.84641      0.0   0.0
0.0        0.0          0.0         0.0         0.0          1.31739   -0.614252  -0.203142    -2.5   2.0
0.0        0.0          0.0         0.0         0.0          2.63479   -1.2285    -0.406283     1.0  -2.0]

B_data = [col_1_to_30 col_31_to_40]


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Triangle(p1, p2, p3)

A triangle with points `p1`, `p2`, `p3`.
"""
struct Triangle{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim,T}
  vertices::V
end

measure(t::Triangle{2}) = abs(signarea(t))

function Base.in(p::Point{2}, t::Triangle{2})
  a, b, c = t.vertices
  abp = signarea(a, b, p)
  bcp = signarea(b, c, p)
  cap = signarea(c, a, p)
  areas = (abp, bcp, cap)
  all(areas .≥ 0) || all(areas .≤ 0)
end


text:
module StatisticsWithJulia

using Pkg, Markdown

include_dependency(joinpath("..", "Project.toml"))
include("functionsgeneral.jl")   # project_relative_path()
include("functionsnotebook.jl")  # list_notebooks()
include("functionspage.jl")      # footer(), header()
include("functionswebsite.jl")   # start(), update(), and create_sysimage()
include("notebooklist.jl")       # notebooklist

const _VERSION = VersionNumber(Pkg.TOML.parsefile(project_relative_path("Project.toml"))["version"])
@info """\n
    欢迎来到《Julia语言与数据科学初步》课程 v$(_VERSION)!
    请在 REPL 输入下面的命令开始课程：
    ```julia
    julia> StatisticsWithJulia.start()
    ````
\n"""

end # module


text:
# ---
# title: 868. Binary Gap
# id: problem868
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/binary-gap/description/>
# hidden: true
# ---
#
# Given a positive integer `n`, find and return _the **longest distance**
# between any two **adjacent** _`1` _' s in the binary representation of _`n` _.
# If there are no two adjacent_`1` _' s, return _`0` _._
#
# Two `1`'s are **adjacent** if there are only `0`'s separating them (possibly
# no `0`'s). The **distance** between two `1`'s is the absolute difference
# between their bit positions. For example, the two `1`'s in `"1001"` have a
# distance of 3.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 22
#     Output: 2
#     Explanation: 22 in binary is "10110".
#     The first adjacent pair of 1's is " _1_ 0 _1_ 10" with a distance of 2.
#     The second adjacent pair of 1's is "10 _11_ 0" with a distance of 1.
#     The answer is the largest of these two distances, which is 2.
#     Note that " _1_ 01 _1_ 0" is not a valid pair since there is a 1 separating the two 1's underlined.
#
#
# **Example 2:**
#
#
#
#     Input: n = 5
#     Output: 2
#     Explanation: 5 in binary is "101".
#
#
# **Example 3:**
#
#
#
#     Input: n = 6
#     Output: 1
#     Explanation: 6 in binary is "110".
#
#
# **Example 4:**
#
#
#
#     Input: n = 8
#     Output: 0
#     Explanation: 8 in binary is "1000".
#     There aren't any adjacent pairs of 1's in the binary representation of 8, so we return 0.
#
#
# **Example 5:**
#
#
#
#     Input: n = 1
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# RUN: julia --startup-file=no %s %t && llvm-link -S %t/* -o %t/module.ll
# RUN: cat %t/module.ll | FileCheck %s

## Notes:
# This script uses the `emit` function (defined llvmpasses.jl) to emit either
# optimized or unoptimized LLVM IR. Each function is emitted individually and
# `llvm-link` is used to create a single module that can be passed to opt.
# The order in which files are emitted and linked is important since `lit` will
# process the test cases in order.

include(joinpath("..", "testhelpers", "llvmpasses.jl"))

import Base.Experimental: Const, @aliasscope

# CHECK-LABEL: @julia_simple
function simple(A, B)
    @aliasscope @inbounds for I in eachindex(A, B)
        A[I] = Const(B)[I]
# CHECK: load double, {{.*}} !alias.scope [[SCOPE:![0-9]+]]
# CHECK: store double {{.*}} !noalias [[SCOPE]]
    end
    return 0 # return nothing causes japi1
end

# CHECK-LABEL: @julia_constargs
function constargs(A, B::Const)
    @aliasscope @inbounds for I in eachindex(A, B)
        A[I] = B[I]
# CHECK: load double, {{.*}} !alias.scope [[SCOPE2:![0-9]+]]
# CHECK: store double {{.*}} !noalias [[SCOPE2]]
    end
    return 0
end

# CHECK-LABEL: @"julia_micro_ker!
function micro_ker!(AB, Ac, Bc, kc, offSetA, offSetB)
    MR = 8; NR = 6;
    @inbounds @aliasscope for k in 1:kc
        for j in 1:NR, i in 1:MR
            AB[i+(j-1)*MR] = muladd(Const(Ac)[offSetA+i], Const(Bc)[offSetB+j], Const(AB)[i+(j-1)*MR])
# CHECK: load double, {{.*}} !alias.scope [[SCOPE3:![0-9]+]]
# CHECK: load double, {{.*}} !alias.scope [[SCOPE3]]
# CHECK: load double, {{.*}} !alias.scope [[SCOPE3]]
# CHECK: store double {{.*}} !noalias [[SCOPE3]]
        end
        offSetA += MR
        offSetB += NR
    end
    return
end

# CHECK: [[SCOPE]] = !{[[ALIASSCOPE:![0-9]+]]}
# CHECK: [[ALIASSCOPE]] = !{!"aliasscope", [[MDNODE:![0-9]+]]}
# CHECK: [[MDNODE]] = !{!"simple"}

emit(simple, Vector{Float64}, Vector{Float64})
emit(constargs, Vector{Float64}, Const{Float64, 1})
emit(micro_ker!, Matrix{Float64}, Vector{Float64}, Vector{Float64}, Int64, Int64, Int64)

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using Sampling
using Base.Test

function test_categoricalsampler(s, p::Vector{Float64}, n::Int, tol::Float64)
    # s: the sampler to be tested
    # p: the reference probability vector
    # n: the number of samples to be generated
    # tol: maximal allowable deviation

    K = length(p)
    @assert numcategories(s) == K

    cnts = zeros(Int, K)

    for i = 1:n
        k = rand(s)
        @assert 1 <= k <= K
        cnts[k] += 1
    end

    q = cnts ./ n
    # println(Base.maxabs(p - q))
    @test_approx_eq_eps p q tol
end

@test_throws ErrorException AliasTable(Float64[])

s = AliasTable([1.0])
for i = 1:10
    @test rand(s) == 1
end

for p in {[0.3, 0.7], [0.2, 0.3, 0.4, 0.1]}
    test_categoricalsampler(AliasTable(p), p, 10^5, 0.015)
end


text:
#####
##### Gradient AD implementation using ForwardDiff
#####

import .ReverseDiff

struct ReverseDiffLogDensity{L, C} <: ADGradientWrapper
    ℓ::L
    gradientconfig::C
end

"""
$(SIGNATURES)

AD via ReverseDiff.
"""
function ADgradient(::Val{:ReverseDiff}, ℓ)
    z = zeros(dimension(ℓ))
    cfg = ReverseDiff.GradientConfig(z)
    ReverseDiffLogDensity(ℓ, cfg)
end

Base.show(io::IO, ℓ::ReverseDiffLogDensity) = print(io, "ReverseDiff AD wrapper for ", ℓ.ℓ)

function logdensity_and_gradient(fℓ::ReverseDiffLogDensity, x::AbstractVector)
    @unpack ℓ, gradientconfig = fℓ
    buffer = _diffresults_buffer(ℓ, x)
    result = ReverseDiff.gradient!(buffer, x -> logdensity(ℓ, x), x, gradientconfig)
    _diffresults_extract(result)
end


text:
using BinaryBuilder

# Collection of sources required to build Nettle
name = "Nettle"
version = v"3.7.2"

sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/nettle/nettle-$(version).tar.gz",
                  "8d2a604ef1cde4cd5fb77e422531ea25ad064679ff0adf956e78b3352e0ef162"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nettle-*/
atomic_patch -p1 ../patches/alloca.patch
# Include patch for finding definition of `AT_HWCAP2` within the Linux
# kernel headers, rather than the glibc headers, since our glibc is too old
atomic_patch -p1 ../patches/AT_HWCAP2-linux_headers.patch

# Force c99 mode
export CFLAGS="${CFLAGS} -std=c99"

update_configure_scripts
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-include-path=${includedir}
make -j${nproc} SUBDIRS="tools"
make install    SUBDIRS="tools"
install_license COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libnettle", :libnettle),
    LibraryProduct("libhogweed", :libhogweed),
    ExecutableProduct("nettle-hash", :nettle_hash)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.2.0"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
try
    julia_py = ENV["_PYJULIA_JULIA_PY"]

    if Base.julia_cmd().exec[1] == julia_py
        @debug "Already monkey-patched. Skipping..." julia_py getpid()
    else
        @debug "Monkey-patching..." julia_py getpid()

        # Monkey patch `Base.package_slug`
        #
        # This is used for generating the set of precompilation cache paths
        # for PyJulia different to the standard Julia runtime.
        #
        # See also:
        # * Suggestion: Use different precompilation cache path for different
        #   system image -- https://github.com/JuliaLang/julia/pull/29914
        #
        if VERSION < v"1.4.0-DEV.389"
            Base.eval(
                Base,
                quote
                    function package_slug(uuid::UUID, p::Int = 5)
                        crc = _crc32c(uuid)
                        crc = _crc32c(unsafe_string(JLOptions().image_file), crc)
                        crc = _crc32c(unsafe_string(JLOptions().julia_bin), crc)
                        crc = _crc32c($julia_py, crc)
                        return slug(crc, p)
                    end
                end,
            )
        else
            Base.eval(Base, quote
                function package_slug(uuid::UUID, p::Int = 5)
                    crc = _crc32c(uuid)
                    crc = _crc32c($julia_py, crc)
                    return slug(crc, p)
                end
            end)
        end

        # Monkey patch `Base.julia_exename`.
        #
        # This is required for propagating the monkey patches to subprocesses.
        # This is important especially for the subprocesses used for
        # precompilation.
        #
        # See also:
        # * Request: Add an API for configuring julia_cmd --
        #   https://github.com/JuliaLang/julia/pull/30065
        #
        Base.eval(Base, quote
            julia_exename() = $julia_py
        end)
        @assert Base.julia_cmd().exec[1] == julia_py

        @debug "Successfully monkey-patched" julia_py getpid()
    end
catch err
    @error "Failed to monkey-patch `julia`" exception = (err, catch_backtrace()) getpid()
    rethrow()
end

max_stars_repo_name: dpinol/pyjulia

max_issues_repo_name: dpinol/pyjulia

max_forks_repo_name: dpinol/pyjulia


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

using Markdown
using InteractiveUtils

# ╔═╡ 4f0c141c-8ff2-4db7-b580-5f5b45762e26
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ 6ae7ec6a-bb4f-4e81-9766-8d0c4c191e82
md"""
# Problem 2.4.6
"""

# ╔═╡ 5dbb3805-6829-46ef-9ea1-cfd84d47e738
# Model
f(u, p, t) = p * (1.0-u)

# ╔═╡ 10097777-e4de-4d43-9342-87beabd7cb22
p = 1.0

# ╔═╡ 8e6e75ec-dc50-4b3d-a2f4-0f038f199bcc
u0 = 0.0

# ╔═╡ 4d91851b-4f1c-4d1b-a764-552a037166c3
tspan = 10.0

# ╔═╡ a8387560-b8cf-4331-a282-bd9b9625a739
prob = ODEProblem(f, u0, tspan, p)

# ╔═╡ e851eab8-e596-4dc8-8ac4-65e5d91bc6de
sol = solve(prob)

# ╔═╡ 350ddcf9-7a56-4530-8c95-bc2666460e74
plot(sol, xlabel="Time (arbitrary units)", ylabel="Concentration (arbitrary units)", legend=:bottomright)

# ╔═╡ 117275d0-bde0-11eb-2527-87d1644a8e38
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 19afef49-0616-473e-a515-b495b6e435ce
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ Cell order:
# ╠═6ae7ec6a-bb4f-4e81-9766-8d0c4c191e82
# ╠═5dbb3805-6829-46ef-9ea1-cfd84d47e738
# ╠═10097777-e4de-4d43-9342-87beabd7cb22
# ╠═8e6e75ec-dc50-4b3d-a2f4-0f038f199bcc
# ╠═4d91851b-4f1c-4d1b-a764-552a037166c3
# ╠═a8387560-b8cf-4331-a282-bd9b9625a739
# ╠═e851eab8-e596-4dc8-8ac4-65e5d91bc6de
# ╠═350ddcf9-7a56-4530-8c95-bc2666460e74
# ╠═117275d0-bde0-11eb-2527-87d1644a8e38
# ╠═19afef49-0616-473e-a515-b495b6e435ce
# ╠═4f0c141c-8ff2-4db7-b580-5f5b45762e26

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
module TestKwonly
using Kwonly
using Setfield
using Test

@settable struct AKW
    x
    y
    @add_kwonly AKW(x; y=2) = new(x, y)
end

@testset "kwonly" begin
    x0 = AKW(AKW(AKW(5), AKW(6, 7)))
    x1 = @set x0.x.x.x = 10
    x2 = @set x1.x.y.y = 20

    @test x2.x.x.x == 10
    @test x2.x.y.y == 20
    @test x2.x.y.x == x0.x.y.x == 6
end
end


text:
# please instantiate ./test/Manifest.toml for this

using JLD2
import Statistics

# ========================================
# Symmetric Fuzz Test Fixtures
# ========================================

# format: [(values, window_size, expected_medians)]

fixtures = Vector{Tuple{Vector{Float64}, Integer, Vector{Float64}}}(undef, 0)

function naive_symmetric_median(arr, window)
    output = similar(arr)
    offset = round((window - 1) / 2, Base.Rounding.RoundNearestTiesAway) |> Int
    for j in eachindex(arr)
        temp_offset = min(offset, j - 1, length(arr) - j)
        beginning = j - temp_offset
        ending = j + temp_offset
        to_median = arr[beginning:ending] |> skipmissing
        output[j] = Statistics.median(to_median)
    end
    output
end

function push_fixture(N, window_size)
    values = rand(N)
    expected_medians = naive_symmetric_median(values, window_size)
    push!(fixtures, (values, window_size, expected_medians))
end

push_fixture(100_000, 101)

for i in 1:50
    push_fixture(rand(3:20), rand(19:2:100))
    push_fixture(rand(10:1_000), rand(3:2:11))
    push_fixture(rand(100:10_000), rand(11:2:1_000))
end

filename = dirname(@__FILE__) * "/symmetric.jld2"
@save filename fixtures

max_stars_repo_name: Firionus/FastRunningMedian

max_issues_repo_name: Firionus/FastRunningMedian


text:
export CURANDError

struct CURANDError <: Exception
    code::curandStatus_t
    msg::AbstractString
end
Base.show(io::IO, err::CURANDError) = print(io, "CURANDError(code $(err.code), $(err.msg))")

function CURANDError(code::curandStatus_t)
    msg = status_message(code)
    return CURANDError(code, msg)
end

function status_message(status)
    if status == CURAND_STATUS_SUCCESS
        return "generator was created successfully"
    elseif status == CURAND_STATUS_VERSION_MISMATCH
        return "Header file and linked library version do not match"
    elseif status == CURAND_STATUS_NOT_INITIALIZED
        return "Generator not initialized"
    elseif status == CURAND_STATUS_ALLOCATION_FAILED
        return "Memory allocation failed"
    elseif status == CURAND_STATUS_TYPE_ERROR
        return "Generator is wrong type"
    elseif status == CURAND_STATUS_OUT_OF_RANGE
        return "Argument out of range"
    elseif status == CURAND_STATUS_LENGTH_NOT_MULTIPLE
        return "Length requested is not a multple of dimension"
    elseif status == CURAND_STATUS_DOUBLE_PRECISION_REQUIRED
        return "GPU does not have double precision required by MRG32k3a"
    elseif status == CURAND_STATUS_LAUNCH_FAILURE
        return "Kernel launch failure"
    elseif status == CURAND_STATUS_PREEXISTING_FAILURE
        return "Preexisting failure on library entry"
    elseif status == CURAND_STATUS_INITIALIZATION_FAILED
        return "Initialization of CUDA failed"
    elseif status == CURAND_STATUS_ARCH_MISMATCH
        return "Architecture mismatch, GPU does not support requested feature"
    elseif status == CURAND_STATUS_INTERNAL_ERROR
        return "Internal library error"
    else
        return "unknown status"
    end
end

macro check(func)
    quote
        local err::curandStatus_t
        err = $(esc(func::Expr))
        if err != CURAND_STATUS_SUCCESS
            throw(CURANDError(err))
        end
        err
    end
end


text:
using Documenter, GPUArrays

makedocs(
    modules = [GPUArrays],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "GPUArrays.jl",
    pages = [
        "Home" => "index.md",
    ],
    doctest = true
)

deploydocs(
    repo   = "github.com/JuliaGPU/GPUArrays.jl.git"
)


text:
# using MultiScaleTreeGraph
# using PlantGeom
# using CairoMakie

# Import / pre-compute
file = joinpath(dirname(dirname(pathof(PlantGeom))), "test", "files", "simple_plant.opf")
opf = read_opf(file)
meshes = get_ref_meshes(opf)
# Caching the meshes as we are plotting the opf several times:
transform!(opf, refmesh_to_mesh!)

# Make reference for ref meshes:
f, ax, p = viz(meshes)
save("reference_images/refmesh_basic.png", f)

f, ax, p = viz(meshes, color = [:burlywood4, :springgreen4, :burlywood4])
save("reference_images/refmesh_allcolors.png", f)

# Or just changing the color of some:
f, ax, p = viz(meshes, color = Dict(1 => :burlywood4, 3 => :burlywood4))
save("reference_images/refmesh_somecolors.png", f)

# One color for each vertex of the refmesh 0:
f, ax, p = viz(
    meshes,
    color = Dict(
        1 => 1:nvertices(meshes)[1],
        2 => 1:nvertices(meshes)[2],
        3 => 1:nvertices(meshes)[3]
    )
)
save("reference_images/refmesh_vertex_colors.png", f)

# OPF recipe
# Regular plot:
f, ax, p = viz(opf)
save("reference_images/opf_basic.png", f)

# With one shared color:
f, ax, p = viz(opf, color = :red)
save("reference_images/opf_one_color.png", f)

# One color per reference mesh:
f, ax, p = viz(opf, color = Dict(1 => :burlywood4, 2 => :springgreen4, 3 => :burlywood4))
save("reference_images/opf_one_color_per_ref.png", f)

# Or just changing the color of some:
f, ax, p = viz(opf, color = Dict(1 => :burlywood4))
save("reference_images/opf_one_color_one_ref.png", f)

# One color for each vertex of the refmesh 1:
f, ax, p = viz(opf, color = Dict(1 => 1:nvertices(get_ref_meshes(opf))[1]))
save("reference_images/opf_color_ref_vertex.png", f)

# Or coloring by opf attribute, e.g. using the mesh max Z coordinates (NB: need to use
# `refmesh_to_mesh!` before, see above):
transform!(opf, :geometry => (x -> zmax(x.mesh)) => :z_max, ignore_nothing = true)
f, ax, p = viz(opf, color = :z_max)
save("reference_images/opf_color_attribute.png", f)

# Or even coloring by the value of the Z coordinates of each vertex:
transform!(opf, :geometry => (x -> [i.coords[3] for i in x.mesh.points]) => :z, ignore_nothing = true)
f, ax, p = viz(opf, color = :z, showfacets = true)
save("reference_images/opf_color_attribute_vertex.png", f)

f, ax, p = viz(opf, color = :z)
colorbar(f[1, 2], p)
f
save("reference_images/opf_color_attribute_colorbar.png", f)

nothing # to avoid returning anything.


text:
# Register sumProductRule: SPPoissonOutNP, SPPoissonLPN
@sumProductRule(:node_type     => Poisson,
                :outbound_type => Message{Poisson},
                :inbound_types => (Nothing, Message{PointMass}),
                :name          => SPPoissonOutNP)

@sumProductRule(:node_type     => Poisson,
                :outbound_type => Message{Gamma},
                :inbound_types => (Message{PointMass}, Nothing),
                :name          => SPPoissonLPN)

# Register naiveVariationalRule: VBPoissonOut, VBPoissonL

@naiveVariationalRule(:node_type     => Poisson,
                      :outbound_type => Message{Poisson},
                      :inbound_types => (Nothing, ProbabilityDistribution),
                      :name          => VBPoissonOut)

@naiveVariationalRule(:node_type     => Poisson,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, Nothing),
                      :name          => VBPoissonL)


text:
# Autogenerated wrapper script for gmsh_jll for i686-w64-mingw32-cxx11
export gmsh_api, libgmsh

JLLWrappers.@generate_wrapper_header("gmsh")
JLLWrappers.@declare_file_product(gmsh_api)
JLLWrappers.@declare_library_product(libgmsh, "gmsh.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        gmsh_api,
        "lib\\gmsh.jl",
    )

    JLLWrappers.@init_library_product(
        libgmsh,
        "bin\\gmsh.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
##########################################
#         Importing
##########################################

include("FieldArithmetic.jl")

using .FieldArithmetic
import Base: +,-,*,/,^,==,inv,show,zeros,gcd,gcdx,iszero,isone,zero,one

##########################################
#
#Elliptic Curve in Weierstrass Normal Form
#
##########################################

struct Weierstrass_EC{T <: CRing, S <: CRingElem}
    parent::T
    A::S
    B::S
end

struct Weierstrass_ECElem{T <: CRingElem, S <: CRing}
    X::T
    Y::T
    Z::T
    curve::Weierstrass_EC{S, T}
end

# Curve-Constructor
function Weierstrass_EC(R::T, A::S, B::S) where {T <: CRing, S <: CRingElem}
    if isSingular(R,A,B)
        error("Curve is singular!")
    end
    return Weierstrass_EC{T,S}(R, A, B)
end

# Element-Constructor
function (E::Weierstrass_EC{T,S})(X::S, Y::S, Z::S) where {T <: CRing, S <: CRingElem}
    if !validPoint(X, Y, Z, E)
        error("Point is not valid!!!")
    end
    return Weierstrass_ECElem{S, T}(X,Y,Z,E)
end

###### Printing ######
function show(io::IO, R::Weierstrass_EC{T,S}) where {T <: CRing, S <: CRingElem}
    print(io, "Elliptic curve y^2 = x^3 + $(R.A)*x + $(R.B) defined over $(R.parent)")
end

function show(io::IO, P::Weierstrass_ECElem{T,S}) where {T <: CRingElem, S <: CRing}
    print(io, "($(P.X): $(P.Y): $(P.Z))")
end

###### Helper-functions ######
function isSingular(E::T, A::S, B::S) where {T <: CRing, S <: CRingElem}
    return iszero(4*A^3 + 27*B^2)
end

function validPoint(X::S, Y::S, Z::S, E::Weierstrass_EC{T,S}) where {T <: CRing, S <: CRingElem}
    return X^3 + X*Z^2*E.A + Z^3*E.B == Z*Y^2 && !isZeroCoordinate(X, Y, Z)
end

function isZeroCoordinate(X::S, Y::S, Z::S) where {S <: CRingElem}
    F = X.ring
    return (X,Y,Z)==(zero(F),zero(F),zero(F))
end

function identity(E::Weierstrass_EC{T,S}) where {T <: CRing, S <: CRingElem}
    F = E.parent
    return E(zero(F),one(F),zero(F))
end

function isIdentity(P::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    return P == identity(P.curve)
end

function sameCurve(P::Weierstrass_ECElem{S, T}, Q::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    return P.curve == Q.curve
end

function normalized(P::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    X, Y, Z = P.X, P.Y, P.Z
    E = P.curve
    F = E.parent
    if !iszero(Z)
        Zinv = inv(Z)
        return E(X*Zinv, Y*Zinv, one(F))
    else
        return identity(E)
    end
end

###### Basic operators #######

function -(P::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    E = P.curve
    return E(P.X, -P.Y, P.Z)
end

function +(P::Weierstrass_ECElem{S, T}, Q::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    if !(sameCurve(P,Q))
        error("Cannot add points from different curves!")
    end
    if isIdentity(P)
        return Q
    end
    if isIdentity(Q)
        return P
    end
    E = P.curve
    P = normalized(P)
    Q = normalized(Q)
    A = E.A
    X1, Y1, Z1 = P.X, P.Y, P.Z
    X2, Y2, Z2 = Q.X, Q.Y, Q.Z
    if X1 == X2 && Y1 == -Y2
        return identity(E)
    else
        if P == Q
            l = (3*X1^2 + A)*inv(2*Y1)
        else
            l = (Y2 - Y1)*inv(X2 - X1)
        end
        X3 = l^2 - X1 - X2
        Y3 = l*(X1 - X3) - Y1
        return E(X3, Y3, one(F))
    end
end

function -(P::Weierstrass_ECElem{S, T}, Q::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    return P+(-Q)
end

function *(n::U, P::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem, U <: Integer}
    if isIdentity(P) || n == 1
        return P
    end
    E = P.curve
    if n == 0
        return identity(E)
    end
    if n < 0
        P = -P
        n = -n
    end
    bit = U(1) << (ndigits(n, base = 2) - 1)
    Q = P
    bit >>= 1
    while bit != 0
        Q = Q+Q
        if (bit & n) != 0
            Q += P
        end
        bit >>=1
    end
    return Q
end

function *(P::Weierstrass_ECElem{S, T}, n::U) where {T <: CRing, S <: CRingElem, U <: Integer}
    return n*P
end

function ==(P::Weierstrass_ECElem{S, T}, Q::Weierstrass_ECElem{S, T}) where {T <: CRing, S <: CRingElem}
    P1 = normalized(P)
    Q1 = normalized(Q)
    return sameCurve(P,Q) && P1.X == Q1.X && P1.Y == Q1.Y && P1.Z == Q1.Z
end


text:
function inner_prod_metric(x)
	n = size(x, 2)
	D = zeros(eltype(x), n, n)
	@inbounds for i in 1:n
		D[i,i] = dot(view(x, :, i), view(x, :, i))
		for j in i+1:n
			D[i,j] = D[j,i] = dot(view(x, :, i), view(x, :, j))
		end
	end
	D
end

function inner_prod_metric_back(x, sym_dD)
	dD = 0.5*(sym_dD + sym_dD')
	m, n = size(x)
	dx = zeros(eltype(x), m, n)
	@inbounds for i in 1:m
		for j in 1:n
			for k in 1:n
				dx[i,j] += 2*dD[j,k]*x[i,k]
			end
		end
	end
	dx
end

function inner_prod_metric(x, xo)
	N1 = size(x, 2)
	N2 = size(xo, 2)
	D = zeros(eltype(x), N1, N2)
	@inbounds for i in 1:N1
		for j in 1:N2
			D[i,j] = dot(view(x, :, i), view(xo, :, j))
		end
	end
	D
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
module Experiment

export predict, loss, eval_model, train!, report_num_params, plot_task

using ..NeuralProcesses

using BSON
using CUDA
using Flux
using PyPlot
using Printf
using ProgressMeter
using Stheno
using Tracker

import StatsBase: std

include("checkpoint.jl")

"""
    eval_model(model, loss, data_gen, epoch::Integer; num_batches::Integer=256)

Evaluate model.

# Arguments
- `model`: Model to evaluate.
- `loss`: Loss function.
- `data_gen`: Data generator.
- `epoch::Integer`: Current epoch.

# Keywords
- `num_batches::Integer=256`: Number of batches to use.

# Returns
- `Tuple{Float32, Float32}`: Loss value and error.
"""
function eval_model(model, loss, data_gen, epoch::Integer; num_batches::Integer=256)
    model = NeuralProcesses.untrack(model)
    @time tuples = map(x -> loss(model, epoch, gpu.(x)...), data_gen(num_batches))
    values = map(x -> x[1], tuples)
    sizes = map(x -> x[2], tuples)

    # Compute and print loss.
    loss_value, loss_error = _mean_error(values)
    println("Losses:")
    @printf(
        "    %8.3f +- %7.3f (%d batches)\n",
        loss_value,
        loss_error,
        num_batches
    )

    # Normalise by average size of target set.
    @printf(
        "    %8.3f +- %7.3f (%d batches; normalised)\n",
        _mean_error(values ./ mean(sizes))...,
        num_batches
    )

    # Normalise by the target set size.
    @printf(
        "    %8.3f +- %7.3f (%d batches; global mean)\n",
        _mean_error(values ./ sizes)...,
        num_batches
    )

    return loss_value, loss_error
end

_mean_error(xs) = (mean(xs), 2std(xs) / sqrt(length(xs)))

_nanreport = Flux.throttle(() -> println("Encountered NaN loss! Returning zero."), 1)

function _nansafe(loss, xs...)
    value, value_size = loss(xs...)
    if isnan(value)
        _nanreport()
        return Tracker.track(identity, 0f0), value_size
    else
        return value, value_size
    end
end

"""
    train!(
        model,
        loss,
        data_gen,
        opt;
        bson=nothing,
        starting_epoch::Integer=1,
        epochs::Integer=100,
        tasks_per_epoch::Integer=1000,
        path="output"
    )

Train a model.

# Arguments
- `model`: Model to train.
- `loss`: Loss function.
- `data_gen`: Data generator.
- `opt`: Optimiser. See `Flux.Optimiser`.

# Keywords
- `bson`: Name of file to save model to. Set to `nothing` to not save the model.
- `starting_epoch::Integer=1`: Epoch to start training at.
- `epochs::Integer=100`: Number of epochs to train for.
- `tasks_per_epoch::Integer=1000`: Number of tasks to draw in each poch.
- `path="output"`: The model will be tested after every epochs, which produces plots. This
    specifies the directory to write plots to.
"""
function train!(
    model,
    loss,
    data_gen,
    opt;
    bson=nothing,
    starting_epoch::Integer=1,
    epochs::Integer=100,
    tasks_per_epoch::Integer=1000,
    path="output"
)
    CUDA.GPUArrays.allowscalar(false)

    # Divide out batch size to get the number of batches per epoch.
    batches_per_epoch = div(tasks_per_epoch, data_gen.batch_size)

    # Display the settings of the training run.
    @printf("Epochs:               %-6d\n", epochs)
    @printf("Starting epoch:       %-6d\n", starting_epoch)
    @printf("Tasks per epoch:      %-6d\n", batches_per_epoch * data_gen.batch_size)
    @printf("Batch size:           %-6d\n", data_gen.batch_size)
    @printf("Batches per epoch:    %-6d\n", batches_per_epoch)

    # Track the parameters of the model for training.
    model = NeuralProcesses.track(model)

    for epoch in starting_epoch:(starting_epoch + epochs - 1)
        # Perform epoch.
        CUDA.reclaim()
        @time begin
            ps = Flux.Params(Flux.params(model))
            @showprogress "Epoch $epoch: " for d in data_gen(batches_per_epoch)
                gs = Tracker.gradient(ps) do
                    first(_nansafe(loss, model, epoch, gpu.(d)...))
                end
                for p in ps
                    Tracker.update!(p, -Flux.Optimise.apply!(opt, p, Tracker.data(gs[p])))
                end
            end
        end

        # Evalute model.
        CUDA.reclaim()
        loss_value, loss_error = eval_model(
            NeuralProcesses.untrack(model),
            loss,
            data_gen,
            epoch
        )

        # Plot model.
        CUDA.reclaim()
        _plot_task(
            NeuralProcesses.untrack(model),
            data_gen,
            epoch,
            _make_plot_true(data_gen.process),
            path=path
        )

        # Save result.
        if !isnothing(bson)
            checkpoint!(
                bson,
                NeuralProcesses.untrack(model),
                epoch,
                loss_value,
                loss_error
            )
        end
    end
end

"""
    report_num_params(model)

Report the number of parameters of a model.

# Arguments
- `model`: Model to report the number of parameters of.
"""
function report_num_params(model)
    @printf("Number of parameters: %-6d\n", sum(map(length, Flux.params(model))))
end

function _plot_task(
    model,
    data_gen,
    epoch,
    plot_true = (plt, xc, yc, xt, σ²) -> nothing;
    path = "output",
    num_tasks = 5
)
    for i = 1:num_tasks
        x = collect(range(-3, 3, length=400))

        # Predict on a task.
        xc, yc, xt, yt = map(x -> x[:, 1, 1], data_gen(1)[1])
        μ, lower, upper, samples = predict(model, xc, yc, x, epoch=epoch)

        figure(figsize=(8, 4))

        # Scatter target and context set.
        scatter(xt, yt, c="tab:red", label="Target set")
        scatter(xc, yc, c="black", label="Context set")

        # Plot prediction of true, underlying model.
        plot_true(xc, yc, x, data_gen.σ²)

        # Plot prediction.
        if !isnothing(μ)
            plot(x, μ, c="tab:green", label="Model output")
            fill_between(x, lower, upper, alpha=0.2, facecolor="tab:green")
        end

        # Plot samples.
        if !isnothing(samples)
            plot(x, samples, c="tab:green", lw=0.5)
        end

        legend()
        tight_layout()

        if !isnothing(path)
            savefig("$path/epoch$epoch-$i.png", dpi=200)
        end

        close()
    end
end

_make_plot_true(process) = (xc, yc, xt, σ²) -> nothing

function _make_plot_true(process::GP)
    function plot_true(xc, yc, xt, σ²)
        xc, yc, xt = Float64.(xc), Float64.(yc), Float64.(xt)
        posterior = process | Obs(process(xc, σ²) ← yc)
        margs = marginals(posterior(xt))
        plot(xt, mean.(margs), c="tab:blue", label="GP")
        error = 2 .* sqrt.(std.(margs).^2 .+ σ²)
        plot(xt, mean.(margs) .- error, c="tab:blue", ls="--")
        plot(xt, mean.(margs) .+ error, c="tab:blue", ls="--")
    end
    return plot_true
end

end


text:
import Pkg

Pkg.add("ArgParse")
Pkg.add("Clustering")
Pkg.add("ClusterManagers")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Distances")
Pkg.add("Distributed")
Pkg.add("StatsBase")
Pkg.add("JuMP")
Pkg.add("Ipopt")
Pkg.add("MathOptInterface")
Pkg.add("Distributions")
Pkg.add("LineSearches")
Pkg.add("NLsolve")
Pkg.add("Optim")

using Pkg: @pkg_str
pkg"precompile"

max_stars_repo_name: ypradat/sigprofilerjulia

max_issues_repo_name: ypradat/sigprofilerjulia

max_forks_repo_name: ypradat/sigprofilerjulia


text:
# Interface -------

function stargazers(owner, repo; auth = AnonymousAuth(), options...)
    stargazers(auth, owner, repo; options...)
end

function stargazers(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/stargazers");
            headers = headers,
            options...)

    handle_error(r)

    data = JSON.parse(r.data)
    map!( u -> User(u), data)
end

function starred(user; auth = AnonymousAuth(), options...)
    starred(auth, user; options...)
end

function starred(auth::Authorization, user; headers = Dict(),
                                            query = Dict(),
                                            sort = nothing,
                                            direction = nothing,
                                            options...)
    authenticate_headers(headers, auth)

    sort != nothing && (query["sort"] = sort)
    direction != nothing && (query["direction"] = direction)

    r = get(URI(API_ENDPOINT; path = "/users/$user/starred"); query = query,
                                                              headers = headers,
                                                              options...)
    handle_error(r)

    data = JSON.parse(r.data)
    map!( u -> Repo(u), data)
end

function star(owner, repo; auth = AnonymousAuth(), options...)
    star(auth, owner, repo; options...)
end

function star(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)

    r = put(URI(API_ENDPOINT; path = "/user/starred/$owner/$repo"); headers = headers,
                                                                    options...)
    handle_error(r)
end

function unstar(owner, repo; auth = AnonymousAuth(), options...)
    unstar(auth, owner, repo; options...)
end

function unstar(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)

    r = delete(URI(API_ENDPOINT; path = "/user/starred/$owner/$repo"); headers = headers,
                                                                       options...)
    handle_error(r)
end


text:
#******************************************************************************#
#                                                                              #
#                                                         :::      ::::::::    #
#    code.jl                                            :+:      :+:    :+:    #
#                                                     +:+ +:+         +:+      #
#    By: mngomane <mngomane@student.42.fr>          +#+  +:+       +#+         #
#                                                 +#+#+#+#+#+   +#+            #
#    Created: 2015/03/13 22:44:39 by mngomane          #+#    #+#              #
#    Updated: 2015/03/13 22:44:39 by mngomane         ###   ########.fr        #
#                                                                              #
#******************************************************************************#

args = length(ARGS)

if	args == 1
	len = length(ARGS[1])

	if len > 0
		print(int(ARGS[1][1]))
	end
	if len > 1
		index = 2;

		while index <= len
			print(string(".", int(ARGS[1][index])))
			index += 1;
		end
	end
	print("\n")
end

max_stars_repo_name: mngomane/projects

max_issues_repo_name: mngomane/projects

max_forks_repo_name: mngomane/projects


text:
@testset "Bracket root" begin
    u = uk(3)
    u.ξk .= sort!(randn(5))
    u.σk .= σscale(u.ξk,2.0)
    u.coeffk .= rand(6);
    κ = 5.0
    xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])

    zlim = bracket(u, xlim[1], xlim[2])

    @test u(zlim[1])*u(zlim[2])<0.0

end
@testset "inverse_uk p=-1" begin
    u = uk(-1)
    zopt = inverse_uk(u,2.0,500.0)
    @test zopt == 2.0
    @test u(zopt) == 2.0
end

@testset "inverse_uk p=0" begin
    u = uk(0)
    u.coeffk .= [1.0;2.0]
    zopt = inverse_uk(u,2.0,500.0)
    @test zopt ==(2.0-1.0)/2.0
    @test u(zopt) == 2.0
end

@testset "inverse_uk p>0" begin
    u = uk(5)
    u.ξk .= sort!(randn(7))
    u.σk .= σscale(u.ξk,2.0)
    u.coeffk .= rand(8);

    zopt = inverse_uk(u,1.0,500.0)
    @test norm(u(zopt) -1.0)<1e-9

    # u = uk(5, randn(7), rand(7), [-0.5; rand(7)])
    #
    # zopt = inverse_uk(u,1.0;z0 = 0.1)
    # @test norm(u(zopt) -1.0)<1e-10
end

@testset "Invert S" begin
    Nx = 100
    Ne = 200
    p = 2
    γ = 2.0
    λ = 0.1
    δ = 1e-5
    κ = 10.0

    X = randn(Nx,Ne) .* randn(Nx,Ne)
    Xpost = zero(X)

    S = RadialMap(Nx, p, γ=γ, λ=λ, δ=δ, κ=κ)
    optimize(S, X);

    F = S(X)
    # Truth is a weak pertrubation
    ystar = X[1:50,1] + 0.01*cos.(randn(50))
    for i=1:Ne
        col = view(Xpost,:,i)
        inverse(col, view(F,:,i), S, ystar)
    end

    @test norm(S(Xpost)[51:Nx,:] - F[51:Nx,:])<1e-8
end

@testset "Invert S with Multi-threading" begin
    Nx = 100
    Ne = 200
    p = 2
    γ = 2.0
    λ = 0.1
    δ = 1e-5
    κ = 10.0

    X = randn(Nx,Ne) .* randn(Nx,Ne)
    Xserial = zero(X)
    Xthread = zero(X)
    Sserial = RadialMap(Nx, p, γ=γ, λ=λ, δ=δ, κ=κ)
    Sthread = RadialMap(Nx, p, γ=γ, λ=λ, δ=δ, κ=κ)
    optimize(Sserial, X; P = serial);
    optimize(Sthread, X; P = thread);

    F = Sserial(X)
    @test norm(Sserial(X)-Sthread(X))<1e-12
    # Truth is a weak pertrubation
    ystar = X[1:50,1] + 0.01*cos.(randn(50))
    Threads.@threads for i=1:Ne
        colserial = view(Xserial,:,i)
        colthread = view(Xthread,:,i)
        inverse(colserial, view(F,:,i), Sserial, ystar)
        inverse(colthread, view(F,:,i), Sthread, ystar)

    end

    @test norm(Sserial(Xserial)[51:Nx,:] - F[51:Nx,:])<1e-8
    @test norm(Sthread(Xthread)[51:Nx,:] - F[51:Nx,:])<1e-8

end

#
# @testset "Benchmark Algorithm" begin
# function compare_algo()
# @btime begin
# for i=1:1000
# u = uk(3)
# u.ξk .= sort!(randn(5))
# u.σk .= σscale(u.ξk,2.0)
# u.coeffk .= rand(6);
# κ = 5.0
# xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# zlim = bracket(u, xlim[1], xlim[2])
# end
# end
#
# @btime begin
# for i=1:1000
# u = uk(3)
# u.ξk .= sort!(randn(5))
# u.σk .= σscale(u.ξk,2.0)
# u.coeffk .= rand(6);
# κ = 5.0
# xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# zlim = bracket(u, xlim[1], xlim[2])
# Roots.bisection(u, zlim[1], zlim[2])
# end
# end
#
# @btime begin
# for i=1:1000
# u = uk(3)
# u.ξk .= sort!(randn(5))
# u.σk .= σscale(u.ξk,2.0)
# u.coeffk .= rand(6);
# κ = 5.0
# xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# zlim = bracket(u, xlim[1], xlim[2])
# find_zero(u, zlim, Bisection())
# end
# end
#
# #
# # @btime begin
# # for i=1:1000
# # u = uk(3)
# # u.ξk .= sort!(randn(5))
# # u.σk .= σscale(u.ξk,2.0)
# # u.coeffk .= rand(6);
# # κ = 5.0
# # xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# # zlim = bracket(u, xlim[1], xlim[2])
# # find_zero(u, zlim, Roots.A42())
# # end
# # end
#
# @btime begin
# for i=1:1000
# u = uk(3)
# u.ξk .= sort!(randn(5))
# u.σk .= σscale(u.ξk,2.0)
# u.coeffk .= rand(6);
# κ = 5.0
# xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# zlim = bracket(u, xlim[1], xlim[2])
# find_zero(u, zlim, Roots.Brent())
# end
# end
# #
# # @btime begin
# # for i=1:1000
# # u = uk(3)
# # u.ξk .= sort!(randn(5))
# # u.σk .= σscale(u.ξk,2.0)
# # u.coeffk .= rand(6);
# # κ = 5.0
# # xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# # zlim = bracket(u, xlim[1], xlim[2])
# # Roots.secant_method(u, zlim)
# # end
# # end
# #
# # @btime begin
# # for i=1:1000
# # u = uk(3)
# # u.ξk .= sort!(randn(5))
# # u.σk .= σscale(u.ξk,2.0)
# # u.coeffk .= rand(6);
# # κ = 5.0
# # xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# # zlim = bracket(u, xlim[1], xlim[2])
# # find_zero(u, zlim, Roots.Secant())
# # end
# # end
#
# # @btime begin
# # for i=1:1000
# # u = uk(3)
# # u.ξk .= sort!(randn(5))
# # u.σk .= σscale(u.ξk,2.0)
# # u.coeffk .= rand(6);
# # κ = 5.0
# # xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# # zlim = bracket(u, xlim[1], xlim[2])
# # find_zero((z->u(z), D(u),H(u)), mean(zlim), Roots.Newton())
# # end
# # end
# #
# # @btime begin
# # for i=1:1000
# # u = uk(3)
# # u.ξk .= sort!(randn(5))
# # u.σk .= σscale(u.ξk,2.0)
# # u.coeffk .= rand(6);
# # κ = 5.0
# # xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# # zlim = bracket(u, xlim[1], xlim[2])
# # find_zero((z->u(z), D(u),H(u)), mean(zlim), Roots.Newton())
# # end
# # end
# #
# # @btime begin
# # for i=1:1000
# # u = uk(3)
# # u.ξk .= sort!(randn(5))
# # u.σk .= σscale(u.ξk,2.0)
# # u.coeffk .= rand(6);
# # κ = 5.0
# # xlim = (u.ξk[1]-κ*u.σk[1], u.ξk[end]+κ*u.σk[end])
# # zlim = bracket(u, xlim[1], xlim[2])
# # find_zero((z->u(z), D(u),H(u)), mean(zlim), Roots.Halley())
# # end
# # end
#
# end
# compare_algo()
#
# end


text:
# Type that holds neighborhoods of a point
type Neighborhood{T<:FloatingPoint}
    inds::Vector{Int} # Indices of neighbors
    ds::Vector{T}     # Cost of successful connection
end

# Always assume in box 0,1?
type FMTParam
    X_lim::Int64
    Y_lim::Int64
    num_pts::Int64
    radius::Float64
    goal_rad::Float64
end

# Type that holds precomputed information for the path planner
type FMTPlanner
    params::FMTParam                    # Parameter object
    points::Array{Float64,2}            # Vector of points TODO parametrize dimension
    neighborhoods::Vector{Neighborhood} # Vector of all neighborhoods of points
end

max_stars_repo_name: stefantj/FMT

max_issues_repo_name: stefantj/FMT

max_forks_repo_name: stefantj/FMT


text:
module Bar

using Hyperspecialize

@concretize AlsoNotAType Set{Type}([Int64])
@concretize Float64 Set{Type}([UInt16])
@concretize NotAType Set{Type}([UInt64])
@concretize Float32 Set{Type}([UInt128])

end #module


text:
using CoffeeBreak
using Documenter

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

makedocs(;
    modules=[CoffeeBreak],
    authors="Markus DeMartini",
    repo="https://github.com/markusdemartini/CoffeeBreak.jl/blob/{commit}{path}#{line}",
    sitename="CoffeeBreak.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://markusdemartini.github.io/CoffeeBreak.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/markusdemartini/CoffeeBreak.jl",
    devbranch="main",
)


text:
using LoggingFacilities
using Logging

# Add some colors to the log

preferred_colors = Dict(
    Logging.Debug => ColorSpec(:grey, false),
    Logging.Info  => ColorSpec(:cyan, true),
    Logging.Warn  => ColorSpec(:yellow, true),
    Logging.Error => ColorSpec(:red, true),
)

clogger = ColorMessageTransformerLogger(current_logger(), preferred_colors);

with_logger(clogger) do
    x = 1; y = rand(); z = "Ryan"
    @info "hello world" x y
    @warn "pay attention!" y z
    @error "what!" x z
end

# Compose existing transformer loggers

normal_colors = Dict(
    Logging.Debug => ColorSpec(:light_black, false),
    Logging.Info  => ColorSpec(:cyan, false),
    Logging.Warn  => ColorSpec(:yellow, false),
    Logging.Error => ColorSpec(:red, false),
)

combo_logger = compose(
    MessageOnlyLogger(stderr, Logging.Debug),
    TimestampTransformerLogger(BeginningMessageLocation()),
    OneLineTransformerLogger,
    ColorMessageTransformerLogger(normal_colors)
)

with_logger(combo_logger) do
    x = 1; y = rand(); z = "Ryan"
    @debug "hey there"
    @info "hello world" x y
    @warn "pay attention!" y z
    @error "what!" x z
end


text:
"""
    infer_ClaDS(tree::Tree, n_reccord::Int64; ...)

Infer ClaDS parameters on a tree

# Arguments
- `tree::Tree`: the phylogeny on which to perform the inference.
- `n_reccord::Int64`: number of iterations between two computations of the gelman statistics. Default to `1_000`.

# Keyword arguments
- `former_run::CladsOutput`: the result of a run of `infer_ClaDS` on the tree. The new mcmc chains will be added to these ones.
- `f::Float64`: The sampling probability. It can be a `Float` (in which case the whole clade has the same sampling probability), or a vector length `n_tips(tree)`(in which case sampling probabilities are subclades specific andmust be given in the same order as the tip tip_labels(tree))` Default to `1.`.
- `goal_gelman::Float64`: The gelman parameter value below which the run is stoped. Default to `1.05`.
- `end_it::Int64`: Maximum number of MCMC iterations. If the number of iterations reaches this value, the run stop even if the gelman statistic is still above `goal_gelamn`. Default to `Inf`.
- `thin::Int64`: The thinning parameter. Default to `1.`.
- `burn::Float64`: The proportion of the mcmc that will be discarded befor computing the gelman statistics and point estiamtes for the parameters. Default to `0.25`.
- `n_trees::Int64`: Number of samples from the posterior distribution of complete phylogenies to be outputed. Default to `10`.
- `ltt_steps::Int64`: Number of time points at which the rate through time and diversity through time should be computed. Default to `50`.
- `print_state::Int64`: If `> 0`, the state of the chains is printed every `print_state`iteration. Default to `0`.
- `prior_ε::String` : The prior to be used for ε. Default to "uniform" (a uniform prior in [0,1000]), but a "lognormal" prior can be defined as an alternative. It can also be set to a flat prior on [0,Inf], as used in the paper, with "unifromInf".
- `logε0::Float64`: If `prior_ε = "lognormal"`, mean of the ε prior on the log scale.
- `sdε::Float64`: If `prior_ε = "lognormal"`, standard deviation of the ε prior on the log scale.
"""
function infer_ClaDS(tree::Tree, n_reccord=1000::Int64; ini_par = [], initialize_rates = 0, goal_gelman = 1.05,
    thin = 1, burn = 1/4, f = 1., plot_tree = 0, print_state = 0, max_node_number = 100, plot_chain = false,
    max_try = 10_000, it_edge_tree = 30, print_all = false, it_rates = 3, former_run = CladsOutput(), plot_burn = NaN, ltt_steps = 50,
    end_it = Inf, n_chains = 3, n_trees = 10, prior_ε = "uniform", logε0 = 0., sdε = 0.5)

    ntips = Int64((tree.n_nodes + 1)/2)
    if isnan(plot_burn)
        plot_burn = burn
    end

    if length(f) == 1
        fs = fill(f,tree.n_nodes-1)
    elseif length(f) == ntips
        fs = sample_fractions(tree,f)[2:end]
    else
        fs = f
    end
    if length(former_run.chains) == 0
        sampler = initialize_ClaDS2(tree , ini_par = ini_par, initialize_rates = initialize_rates, ltt_steps = ltt_steps, n_chains = n_chains, n_trees = n_trees)
        if initialize_rates > 0.
            sampler = add_iter_ClaDS2(sampler, 1, thin = initialize_rates, fs = fs, plot_tree = plot_tree, print_state = 1,
                max_node_number = max_node_number, max_try = max_try, it_edge_tree = 0.,
                print_all = print_all, it_rates = 1, n_trees = n_trees, ini = true)
        end
    else
        sampler = (former_run.chains,
            former_run.current_state[1],
            former_run.current_state[2],
            former_run.current_state[3],
            former_run.current_state[4],
            former_run.time_points,
            former_run.current_state[5],
            former_run.rtt_chains,
            former_run.enhanced_trees)

            former_run = CladsOutput()
    end

    min_branch = minimum(extract_branch_lengths(tree)[2:end])

    if min_branch < 1e-10
        println("Warning : some of the branches of your tree are very small ($(min_branch)), this may lead to numerical issues.")
    end

    gelman = 10.
    MAPS=[]
    nit = 0

    while (gelman > goal_gelman) & (nit < end_it)

        sampler = add_iter_ClaDS2(sampler, n_reccord, thin = thin, fs = fs, plot_tree = plot_tree, print_state = print_state,
            max_node_number = max_node_number, max_try = max_try, it_edge_tree = it_edge_tree,
            print_all = print_all, it_rates = it_rates, n_trees = n_trees, prior_ε = prior_ε, logε0 = logε0, sdε = sdε)

        nit += n_reccord

        chains = sampler[1]
        npar = sampler[4][1].n_nodes + 3
        @rput chains
        @rput npar

        if goal_gelman == 0.
            id_gelman = 0
            gelman = 10
        else
            gelman = gelman_est(sampler[1], npar, burn = burn)
            id_gelman = gelman[1]
            gelman = gelman[2]
        end

        if plot_chain
            plot_coda(sampler, burn = plot_burn, id_par=[1,2,3,max(4,id_gelman)])
        else
            chains = sampler[1]
            @rput chains
            @rput thin
            @rput burn
            id_par = 1
            @rput id_par
            reval("""
                require(coda)
                n_row = length(chains[[1]][[1]])
                ini = floor(burn * n_row)
                if(ini == 0) ini = 1
                    it = seq(ini, n_row, 1)
            """)
        end

        @rput ntips
        @rput n_chains

        reval("""
            npar2 = length(chains[[1]])
            n_row = length(chains[[1]][[1]])
            ini = floor(burn * n_row)
            if(ini == 0) ini = 1
            it = c(unique(ceiling(seq(ini, n_row-1,length.out = 100))),n_row)

            map_chains = mcmc.list(lapply(1:n_chains, function(i){
                mcmc(sapply(1:npar2, function(j){
                    if(j<=3 || j>(npar+ntips)) {
                        chains[[i]][[j]][it]
                    }else{
                        log(chains[[i]][[j]][it])
                    }
            }))}))

            unlist_chains = (sapply(1:npar2, function(k){
                x = c()
                for (n in 1:n_chains){
                    x = c(x,map_chains[[n]][,k] )
                    }
                return(x)
                }))

            MAPS=sapply(1:npar2, function(i){D=density(unlist_chains[,i]);
                return(D[[1]][which.max(D[[2]])])})
            #MAPS[2] = exp(MAPS[2])
        """)
        @rget MAPS
        println("iteration $(length(sampler[1][1][1])) ; gelman = $gelman for param $id_gelman")
        println("   σ = $(MAPS[1]), α = $(MAPS[2]), ε = $(MAPS[3]), λ0 = $(MAPS[4])")
        println("")
    end

    gelm = (0,0)
    npar = sampler[4][1].n_nodes + 3

    if goal_gelman > 0.
        gelm = gelman_est(sampler[1], npar, burn = burn)
    end

    @rput ntips
    @rput n_chains
    @rput npar

    reval("""
        npar2 = length(chains[[1]])
        n_row = length(chains[[1]][[1]])
        ini = floor(burn * n_row)
        if(ini == 0) ini = 1
        it = seq(ini, n_row, 1)

        map_chains = mcmc.list(lapply(1:n_chains, function(i){
            mcmc(sapply(1:npar2, function(j){
                if(j==1 || j==3 || j>(npar+ntips)) {
                    chains[[i]][[j]][it]
                }else{
                    log(chains[[i]][[j]][it])
                }
        }))}))

        unlist_chains = (sapply(1:npar2, function(k){
            x = c()
            for (n in 1:n_chains){
                x = c(x,map_chains[[n]][,k] )
                }
            return(x)
            }))

        MAPS=sapply(1:npar2, function(i){D=density(unlist_chains[,i]);
            return(D[[1]][which.max(D[[2]])])})
        MAPS[2] = exp(MAPS[2])
        MAPS[4:(npar+ntips)] = exp(MAPS[4:(npar+ntips)])

        mean_ltt = sapply((npar + ntips + 3):npar2, function(i){mean(unlist_chains[,i])})
        unlist_chains = list()
        map_chains = list()
    """)
    @rget MAPS
    @rget mean_ltt

    mean_rates_chains = sampler[8]
    @rput mean_rates_chains
    reval("""
        nmr = length(mean_rates_chains[[1]][[1]])
        map_chains = mcmc.list(lapply(1:n_chains, function(i){
            mcmc(sapply(1:nmr, function(j){
                res = c()
                for(t in it){
                    res = c(res,log(mean_rates_chains[[i]][[t]][j]))
                    }
                return(res)

                }))}))

        unlist_chains = (sapply(1:nmr, function(k){
            x = c()
            for (n in 1:n_chains){
                x = c(x,map_chains[[n]][,k] )
                }
            return(x)
            }))
        RTT_map = sapply(1:nmr, function(i){D=density(unlist_chains[,i]);
            return(D[[1]][which.max(D[[2]])])})
        RTT_map = exp(RTT_map)
    """)
    @rget RTT_map

    output = CladsOutput(
        tree,
        sampler[1],                         #chains
        mean_rates_chains,                  #rtt_chains
        MAPS[1], MAPS[2], MAPS[3], MAPS[4], #σ_map,α_map,ε_map,λ0_map
        MAPS[5:npar],                       #λi_map
        MAPS[(npar+1):(npar+ntips)],      #λtip_map
        mean_ltt,                           #DTT_mean
        RTT_map,                            #RTT_map
        sampler[6],                         #time_points
        sampler[9],                         #enhanced_trees
        gelm,                               #gelman statistics
        (sampler[2], sampler[3], sampler[4], sampler[5], sampler[7])
    #=
    current_state=#
    )
    return output

end


text:
@testset "Set Algorithms" begin

col_a1 = @Column A =   [1,   2,   4,   1,   4,   1]
col_b1 = @Column B = [1.0, 2.0, 1.0, 2.0, 4.0, 1.0]

col_a2 = @Column A =   [2,   3,   4,   4]
col_b2 = @Column B = [0.0, 0.0, 0.0, 4.0]

table1 = hcat(col_a1, col_b1)
table2 = hcat(col_a2, col_b2)

@test Set(@inferred(unique(col_a1)).data) == Set([1,2,4])
@test (tmp = copy(col_a2); unique!(tmp); Set(tmp.data) == Set([2,3,4]))

@test length(@inferred(unique(table1))) == 5
@test (tmp = copy(table1); unique!(tmp); tmp == unique(table1))

@test length(@inferred(union(col_a1,col_a2))) == 4
@test length(@inferred(setdiff(col_a1,col_a2))) == 1
@test length(@inferred(intersect(unique(col_a1),col_a2))) == 2 # there may be a problem with Julia's intersect for non-unique inputs...

@test length(@inferred(union(table1,table2))) == 8
@test length(@inferred(setdiff(table1,table2))) == 4
@test length(@inferred(intersect(table1,table2))) == 1

table1_copy = copy(table1)
@test (union!(table1_copy,table2); length(table1_copy) == 8)

end