text:
@testset "Tracker" begin
    @testset "tracking - AD : $AD" for AD = 0:3
        @var x a y b
        F = System([x^2 - a, x * y - a + b], [x, y], [a, b])

        tracker = Tracker(
            ParameterHomotopy(F, [1, 0], [2, 4]),
            options = TrackerOptions(automatic_differentiation = AD),
        )
        s = [1, 1]
        res = track(tracker, s, 1, 0)
        @test is_success(res)
        @test steps(res) ≤ 5
        @test isa(solution(res), Vector{ComplexF64})
        @test solution(res) ≈ [sqrt(2), -sqrt(2)]
        @test length(solution(res)) == 2
        @test is_success(track(tracker, res, 0, 1))

        code = track!(tracker, s, 1, 0)
        @test is_success(code)
        s0 = copy(tracker.state.x)
        @unpack μ, ω = tracker.state
        @test is_success(track(tracker, s0, 0, 1))
        @test is_success(track(tracker, s0, 0, 1, μ = μ, ω = ω))

        s = @SVector [1, 1]
        @test is_success(track(tracker, s, 1, 0))
    end

    @testset "projective tracking- AD: $AD" for AD = 0:3
        @var x a y b z
        F = System([x^2 - a * z^2, x * y + (b - a) * z^2], [x, y, z], [a, b])
        H = ParameterHomotopy(F, [1, 0], [2, 4])
        tracker = Tracker(
            on_affine_chart(H, (2,)),
            options = TrackerOptions(automatic_differentiation = AD),
        )

        s = [1, 1, 1]
        res = track(tracker, s, 1, 0)
        @test is_success(res)
        @test isa(solution(res), Vector{ComplexF64})
        x₀ = abs(solution(res)[end])
        @test affine_chart(PVector(solution(res))) ≈ [sqrt(2), -sqrt(2)] rtol = 1e-12 / x₀
    end

    @testset "iterator" begin
        @var x a y b
        F = System([x^2 - a, x * y - a + b], [x, y], [a, b])
        tracker = Tracker(ParameterHomotopy(F, [1, 0], [2, 4]))
        s = [1, 1]

        # path iterator
        typeof(first(iterator(tracker, s, 1.0, 0.0))) == Tuple{Vector{ComplexF64},Float64}

        tracker.options.max_step_size = 0.01
        @test length(collect(iterator(tracker, s, 1.0, 0.0))) ≥ 101

        F = System([x - a], [x], [a])
        ct = Tracker(
            ParameterHomotopy(F, [1], [2]),
            options = TrackerOptions(max_step_size = 0.015625),
        )
        Xs = Vector{ComplexF64}[]
        for (x, t) in iterator(ct, [1.0], 1.0, 0.0)
            push!(Xs, x)
        end

        @test length(Xs) ≥ length(1:0.015625:2)
    end

    @testset "path info" begin
        @var x a y b
        F = System([x^2 - a, x * y - a + b], [x, y], [a, b])
        tracker = Tracker(ParameterHomotopy(F, [1, 0], [2, 4]))
        info = path_info(tracker, [1, 1], 1, 0)
        @test !isempty(sprint(show, info))
    end

    @testset "Change parameters" begin
        @var x a y b
        F = System([x^2 - a, x * y - a + b]; parameters = [a, b])
        s = [1.0, 1.0 + 0im]
        tracker = Tracker(ParameterHomotopy(F, [2.2, 3.2], [2.2, 3.2]))
        start_parameters!(tracker, [1, 0])
        target_parameters!(tracker, [2, 4])
        res = track(tracker, s, 1.0, 0.0)
        @test is_success(res)
    end

    @testset "Straight Line Homotopy" begin
        @var x y
        F = System([x^2 + y^2 - 2.3, 2x + 3y - 4], [x, y])
        G = System((0.2 + 0.4im) .* [x^2 - 1, y - 1], [x, y])
        H = StraightLineHomotopy(G, F)
        S = [[1, 1], [-1, 1]]
        tracker = Tracker(H, options = (automatic_differentiation = 3,))

        @test is_success(track(tracker, S[1], 1, 0))
        @test is_success(track(tracker, S[2], 1, 0))
    end

    include("test_cases/steiner_higher_prec.jl")
    include("test_cases/four_bar.jl")
end


text:
##############################################################
# Hurdle Distributed Multinomial Regression (HDMR)
##############################################################

# default model used when not specified explicitly, as in fit(HDMR,...)
const HDMR_DEFAULT_MODEL = InclusionRepetition

"Abstract HDMR returned object"
abstract type HDMR{M<:Union{Missing, <:TwoPartModel}} <: DCR end

"""
Relatively heavy object used to return HDMR results when we care about the regulatrization paths.
"""
struct HDMRPaths{M<:Union{Missing, <:TwoPartModel}, X} <: HDMR{M}
  nlpaths::Vector{M} # independent TwoPartModel{GammaLassoPath} for each phrase
  intercept::Bool               # whether to include an intercept in each Poisson regression
                                # (only kept with remote cluster, not with local cluster)
  n::Int                        # number of observations. May be lower than provided after removing all zero obs.
  d::Int                        # number of categories (terms/words/phrases)
  inpos::X                      # indices of covars columns included in positives model
  inzero::X                     # indices of covars columns included in zeros model
end

"""
Relatively light object used to return HDMR results when we only care about estimated coefficients.
"""
struct HDMRCoefs{M<:TwoPartModel, T<:AbstractMatrix, S<:SegSelect, X} <: HDMR{M}
  coefspos::T                   # positives model coefficients
  coefszero::T                  # zeros model coefficients
  intercept::Bool               # whether to include an intercept in each Poisson regression
  n::Int                        # number of observations. May be lower than provided after removing all zero obs.
  d::Int                        # number of categories (terms/words/phrases)
  inpos::X                      # indices of covars columns included in positives model
  inzero::X                     # indices of covars columns included in zeros model
  select::S                     # path segment selector
end

HDMRCoefs{M}(coefspos::T, coefszero::T, intercept::Bool, n::Int, d::Int,
  inpos::X, inzero::X, select::S) where {M<:TwoPartModel, T<:AbstractMatrix, S<:SegSelect, X} =
  HDMRCoefs{M,T,S,X}(coefspos, coefszero, intercept, n, d, inpos, inzero, select)

function HDMRCoefs(m::HDMRPaths{M}, select::SegSelect=defsegselect) where {M<:TwoPartModel}
  coefspos, coefszero = coef(m, select)
  HDMRCoefs{M}(coefspos, coefszero, m.intercept, m.n, m.d, m.inpos, m.inzero, select)
end

"""
    fit(HDMR,covars,counts; <keyword arguments>)
    hdmr(covars,counts; <keyword arguments>)

Fit a Hurdle Distributed Multinomial Regression (HDMR) of counts on covars.

HDMR fits independent hurdle lasso regressions to each column of counts to
approximate a multinomial, picks a segement of each path, and
returns a coefficient matrix (wrapped in HDMRCoefs) representing point estimates
for the entire multinomial (includes the intercept if one was included).

# Example:
```julia
  m = fit(HDMR,covars,counts)
```
# Arguments
- `covars` n-by-p matrix of covariates
- `counts` n-by-d matrix of counts (usually sparse)

# Keywords
- `inpos=1:p` indices of covars columns included in model for positive counts
- `inzero=1:p` indices of covars columns included in model for zero counts
- `intercept::Bool=false` include a intercepts in each hurdle regression
- `parallel::Bool=true` parallelize the poisson fits
- `local_cluster::Bool=true` use local_cluster mode that shares memory across
    parallel workers that is appropriate on a single multicore machine, or
    remote cluster mode that is more appropriate when distributing across machines
    for which sharing memory is costly.
- `verbose::Bool=true`
- `showwarnings::Bool=false`
- `select::SegSelect=MinAICc()` path segment selection criterion
- `kwargs...` additional keyword arguments passed along to fit(Hurdle,...)
"""
function StatsBase.fit(::Type{<:HDMR{M}}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
  kwargs...) where {T<:AbstractFloat, M<:TwoPartModel}

  hdmr(covars, counts, M; kwargs...)
end

"""
    fit(HDMRPaths,covars,counts; <keyword arguments>)
    hdmrpaths(covars,counts; <keyword arguments>)

Fit a Hurdle Distributed Multinomial Regression (HDMR) of counts on covars, and returns
the entire regulatrization paths, which may be useful for plotting or picking
coefficients other than the AICc optimal ones. Same arguments as
[`fit(::HDMR)`](@ref).
"""
function StatsBase.fit(::Type{HDMRPaths{M}}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
  local_cluster=false, # ignored. will always assume remote_cluster
  kwargs...) where {T<:AbstractFloat, M<:TwoPartModel}

  hdmrpaths(covars, counts, M; kwargs...)
end

# default is to use HDMR_DEFAULT_MODEL
StatsBase.fit(::Type{HDMR}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
  kwargs...) where {T<:AbstractFloat} = fit(HDMRCoefs{HDMR_DEFAULT_MODEL}, covars, counts; kwargs...)
StatsBase.fit(::Type{HDMRCoefs}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
    kwargs...) where {T<:AbstractFloat} = fit(HDMRCoefs{HDMR_DEFAULT_MODEL}, covars, counts; kwargs...)
StatsBase.fit(::Type{HDMRPaths}, covars::AbstractMatrix{T}, counts::AbstractMatrix;
  kwargs...) where {T<:AbstractFloat} = fit(HDMRPaths{HDMR_DEFAULT_MODEL}, covars, counts; kwargs...)

# fit wrapper that takes a model (two formulas) and dataframe instead of the covars matrix
# e.g. @model(h ~ x1 + x2, c ~ x1)
# h and c on the lhs indicate model for zeros and positives, respectively.
"""
    fit(HDMR,@model(h ~ x1 + x2, c ~ x1),df,counts; <keyword arguments>)

Fits a HDMR but takes a model formula and dataframe instead of the covars matrix.
See also [`fit(::HDMR)`](@ref).

`h` and `c` on the lhs indicate the model for zeros and positives, respectively.
"""
function StatsBase.fit(::Type{T}, m::Model, df, counts::AbstractMatrix, args...;
  contrasts::Dict{Symbol,<:Any} = Dict{Symbol,Any}(), kwargs...) where {T<:HDMR}

  # parse and merge rhs terms
  trmszero = getrhsterms(m, :h)
  trmspos = getrhsterms(m, :c)
  trms, inzero, inpos = mergerhsterms(trmszero,trmspos)

  # create model matrix
  covars, counts, as = modelcols(trms, df, counts; model=T, contrasts=contrasts)

  # inzero and inpos may be different in the applied schema with factor variables
  inzero, inpos = mapins(inzero, inpos, as)

  # fit and wrap in TableCountsRegressionModel
  TableCountsRegressionModel(fit(T, covars, counts, args...;
    inzero=inzero, inpos=inpos, kwargs...), df, counts, m, as)
end

"""
    coef(m::HDMRCoefs)

Returns the coefficient matrices fitted with HDMR using
the segment selected during fit (MinAICc by default).

# Example:
```julia
  m = fit(HDMR,covars,counts)
  coefspos, coefszero = coef(m)
```
"""
StatsBase.coef(m::HDMRCoefs) = m.coefspos, m.coefszero

coeffill!(coefszero, coefspos, path::Missing, pzero, ppos, j, select::SegSelect) = nothing

function coeffill!(coefszero, coefspos, path::TwoPartModel, pzero, ppos, j, select::SegSelect)
  cjpos, cjzero = coef(path, select)
  coeffill!(coefspos, cjpos, ppos, j, select)
  coeffill!(coefszero, cjzero, pzero, j, select)
end

"""
    coef(m::HDMRPaths, select::SegSelect=MinAICc())

Returns all or selected coefficient matrices fitted with HDMR.

# Example:
```julia
  m = fit(HDMRPaths,covars,counts)
  coefspos, coefszero = coef(m, MinCVKfold{MinCVmse}(5))
```
"""
function StatsBase.coef(m::HDMRPaths, select::SegSelect=defsegselect)
  # get dims
  d = length(m.nlpaths)
  d < 1 && return nothing, nothing

  # drop missing paths
  nonmsngpaths = skipmissing(m.nlpaths)

  # get number of variables from paths object
  ppos = ncoefspos(m)
  pzero = ncoefszero(m)

  # establish maximum path lengths
  nλzero = nλpos = 0
  if !isempty(nonmsngpaths)
    nλzero = maximum(size(nlpath.mzero)[2] for nlpath in nonmsngpaths)
    nλpos = maximum(size(nlpath.mpos)[2] for nlpath in nonmsngpaths)
  end

  nλ = max(nλzero,nλpos)
  # p = pzero + ppos

  # allocate space
  coefszero = coefspace(pzero, d, nλ, select)
  coefspos = coefspace(ppos, d, nλ, select)

  # iterate over paths
  for j=1:d
    path = m.nlpaths[j]
    coeffill!(coefszero, coefspos, path, pzero, ppos, j, select::SegSelect)
  end

  coefspos, coefszero
end

"""
    predict(m,newcovars; <keyword arguments>)

Predict counts using a fitted HDMRPaths object and given newcovars.

# Example:
```julia
  m = fit(HDMRPaths,covars,counts)
  newcovars = covars[1:10,:]
  countshat = predict(m, newcovars; select=MinAICc())
```

# Arguments
- `m::HDMRPaths` fitted DMRPaths model (HDMRCoefs currently not supported)
- `newcovars` n-by-p matrix of covariates of same dimensions used to fit m.

# Keywords
- `select=MinAICc()` See [`coef(::RegularizationPath)`](@ref).
- `kwargs...` additional keyword arguments passed along to predict() for each
  category j=1..size(counts,2)
"""
function StatsBase.predict(m::HDMRPaths, newcovars::AbstractMatrix{T};
  select=defsegselect, kwargs...) where {T<:AbstractFloat}

  covarspos, covarszero = incovars(newcovars,m.inpos,m.inzero)

  _predict(m,covarszero;Xpos=covarspos,select=select,kwargs...)
end

function StatsBase.predict(m::M, newcovars::AbstractMatrix{T};
  select=defsegselect, kwargs...) where {T<:AbstractFloat, M<:HDMR}

  error("predict(m::HDMR,...) can currently only be evaluated for HDMRPaths structs returned from fit(HDMRPaths,...)")
end

"Number of covariates used for HDMR estimation of zeros model"
ncovarszero(m::HDMR) = length(m.inzero)

"Number of covariates used for HDMR estimation of positives model"
ncovarspos(m::HDMR) = length(m.inpos)

"Number of coefficient potentially including intercept used by model for zeros"
ncoefszero(m::HDMR) = ncovarszero(m) + (hasintercept(m) ? 1 : 0)

"Number of coefficient potentially including intercept used by model for positives"
ncoefspos(m::HDMR) = ncovarspos(m) + (hasintercept(m) ? 1 : 0)

"Destandardize two part model estimated using a potentially standardized covars matrix"
function destandardize!(tpm::TwoPartModel, covarsnorm::AbstractVector{T},
    inpos, inzero, standardize) where T

  if standardize
    destandardize!(tpm.mzero, covarsnorm[inzero], standardize)
    destandardize!(tpm.mpos, covarsnorm[inpos], standardize)
  end

  tpm
end

"Shorthand for fit(HDMRPaths,covars,counts). See also [`fit(::HDMRPaths)`](@ref)"
function hdmrpaths(covars::AbstractMatrix{T},counts::AbstractMatrix,::Type{M}=HDMR_DEFAULT_MODEL;
      inpos=1:size(covars,2), inzero=1:size(covars,2),
      intercept=true,
      parallel=true,
      verbose=true, showwarnings=false,
      standardize=true,
      m=nothing, l=nothing, D=nothing,
      kwargs...) where {T<:AbstractFloat,M<:TwoPartModel}
  # get dimensions
  n, d = size(counts)
  n1,p = size(covars)
  @assert n==n1 "counts and covars should have the same number of observations"

  ppos = length(inpos)
  pzero = length(inzero)

  verbose && @info("fitting $n observations on $d categories \n$ppos covariates for positive and $pzero for zero counts")

  # add one coef for intercept
  ncoefpos = ppos + (intercept ? 1 : 0)
  ncoefzero = pzero + (intercept ? 1 : 0)

  covars, counts, μpos, μzero, n = shifters(M, covars, counts, showwarnings, m, l, D)

  # standardize covars only once if needed
  covars, covarsnorm = Lasso.standardizeX(covars, standardize)

  covarspos, covarszero = incovars(covars,inpos,inzero)

  function tryfith(countsj::AbstractVector)
    try
      # we make it dense remotely to reduce communication costs
      # we use the same offsets for pos and zeros
      tpm = fit(M,GammaLassoPath,covarszero,Vector(countsj); Xpos=covarspos, offsetpos=μpos, offsetzero=μzero, verbose=false, standardize=false,
        # showwarnings=showwarnings, # uncomment do help with debugging
        kwargs...)
      destandardize!(tpm, covarsnorm, inpos, inzero, standardize)
    catch e
      showwarnings && @warn("fit($M...) failed for countsj with frequencies $(sort(countmap(countsj))) and will return missing path ($e)")
      missing
    end
  end

  # counts generator
  countscols = (counts[:,j] for j=1:d)

  if parallel
    verbose && @info("distributed $M run on remote cluster with $(nworkers()) nodes")
    mapfn = pmap
  else
    verbose && @info("serial $M run on a single node")
    mapfn = map
  end

  nlpaths = mapfn(tryfith,countscols)

  HDMRPaths(nlpaths, intercept, n, d, inpos, inzero)
end

"Returns a (covarspos, covarszero) tuple from covars"
function incovars(covars,inpos,inzero)
  #NOTE: we don't use views for typesafety
  inall = 1:size(covars,2)

  if inzero == inall
    covarszero = covars
  else
    covarszero = covars[:,inzero]
  end

  if inzero == inpos
    covarspos = covarszero
  elseif inpos == inall
    covarspos = covars
  else
    covarspos = covars[:,inpos]
  end

  covarspos, covarszero
end

function shifters(::Type{Hurdle}, covars::AbstractMatrix, counts::AbstractMatrix, showwarnings::Bool,
  prespecm::Union{Nothing, AbstractVector}, prespecl, prespecd)

  showwarnings && @warn("Old Hurdle-DMR model does not use prespecified l")
  shifters(Hurdle, covars, counts, showwarnings, prespecm, nothing, nothing)
end

# same as DMR, and uses same shifters for both model parts
function shifters(::Type{Hurdle}, covars::AbstractMatrix, counts::AbstractMatrix, showwarnings::Bool,
  prespecm::Union{Nothing, AbstractVector}, prespecl::Nothing, prespecd)

  covars, counts, μ, n = shifters(DMR, covars, counts, showwarnings, prespecm)
  covars, counts, μ, μ, n
end

vocab(counts, prespecd::Nothing) = size(counts, 2)
vocab(counts, prespecd::Int) = prespecd

function shifters(::Type{InclusionRepetition}, covars::AbstractMatrix, counts::AbstractMatrix{C}, showwarnings::Bool,
  prespecm::Union{Nothing, AbstractVector},
  prespecl::Union{Nothing, AbstractVector},
  prespecd::Union{Nothing, Int}) where C

  # standardize counts matrix to conform to GLM.FP
  counts = fpcounts(counts)

  # mi = total word count per observation
  m = totalcounts(counts, prespecm)

  # li = vocabulary per observation
  l = totalcounts(posindic(counts), prespecl)

  if any(iszero,m)
      # omit observations with no counts
      ixposm = findall(x->x!=zero(C), m)
      showwarnings && @warn("omitting $(length(m)-length(ixposm)) observations with no counts")
      m = m[ixposm]
      l = l[ixposm]
      counts = counts[ixposm,:]
      covars = covars[ixposm,:]
  end

  # get new dimensions
  # NOTE: d here is total vocabulary as opposed to obs-specific lexicon l
  d = vocab(counts, prespecd)
  n = length(m)

  # μipos = log(mi-li)
  μpos = broadcast((mi,li)->log(mi-li), m, l)

  # μizero = log(li/(D-li))
  μzero = broadcast(li->log(li/(d-li)), l)

  covars, counts, μpos, μzero, n
end

"Fits a regularized hurdle regression counts[:,j] ~ covars saving the coefficients in coefs[:,j]"
function hurdle_regression!(::Type{M}, coefspos::AbstractMatrix{T}, coefszero::AbstractMatrix{T}, j::Int, covars::AbstractMatrix{T},counts::AbstractMatrix{V},
            inpos, inzero, offsetpos, offsetzero;
            select=defsegselect,
            kwargs...) where {T<:AbstractFloat,V,M<:TwoPartModel}
  cj = Vector(counts[:,j])
  covarspos, covarszero = incovars(covars,inpos,inzero)
  # we use the same offsets for pos and zeros
  path = fit(M,GammaLassoPath,covarszero,cj;
    Xpos=covarspos, offsetpos=offsetpos, offsetzero=offsetzero, kwargs...)
  (coefspos[:,j], coefszero[:,j]) = coef(path, select)
  nothing
end

"Wrapper for hurdle_regression! that catches exceptions in which case it sets coefs to zero"
function tryfith!(::Type{M}, coefspos::AbstractMatrix{T}, coefszero::AbstractMatrix{T}, j::Int, covars::AbstractMatrix{T},counts::AbstractMatrix{V},
            inpos, inzero, offsetpos, offsetzero;
            showwarnings = false,
            kwargs...) where {T<:AbstractFloat,V,M<:TwoPartModel}
  try
    hurdle_regression!(M, coefspos, coefszero, j, covars, counts, inpos, inzero,
      offsetpos, offsetzero;
      # showwarnings=showwarnings, # uncomment do help with debugging
      kwargs...)
  catch e
    showwarnings && @warn("hurdle_regression! failed on count dimension $j with frequencies $(sort(countmap(counts[:,j]))) and will return zero coefs ($e)")
    # redudant ASSUMING COEFS ARRAY INTIAILLY FILLED WITH ZEROS, but can happen in serial mode
    for i=1:size(coefszero,1)
      coefszero[i,j] = zero(T)
    end
    for i=1:size(coefspos,1)
      coefspos[i,j] = zero(T)
    end
  end
end

"Shorthand for fit(HDMR,covars,counts). See also [`fit(::HDMR)`](@ref)"
function hdmr(covars::AbstractMatrix{T},counts::AbstractMatrix{V},::Type{M}=HDMR_DEFAULT_MODEL;
          inpos=1:size(covars,2), inzero=1:size(covars,2),
          intercept=true,
          parallel=true, local_cluster=true,
          verbose=true, showwarnings=false,
          kwargs...) where {T<:AbstractFloat,V,M<:TwoPartModel}
  if local_cluster || !parallel
    hdmr_local_cluster(M, covars,counts,inpos,inzero,intercept,parallel,verbose,showwarnings; kwargs...)
  else
    hdmr_remote_cluster(M, covars,counts,inpos,inzero,intercept,parallel,verbose,showwarnings; kwargs...)
  end
end

"""
This version is built for local clusters and shares memory used by both inputs
and outputs if run in parallel mode.
"""
function hdmr_local_cluster(::Type{M}, covars::AbstractMatrix{T},counts::AbstractMatrix{V},
          inpos,inzero,intercept,parallel,verbose,showwarnings;
          select=defsegselect,
          m=nothing, l=nothing, D=nothing,
          standardize=true, kwargs...) where {T<:AbstractFloat,V,M<:TwoPartModel}
  # get dimensions
  n, d = size(counts)
  n1,p = size(covars)
  @assert n==n1 "counts and covars should have the same number of observations"

  ppos = length(inpos)
  pzero = length(inzero)

  verbose && @info("fitting $n observations on $d categories \n$ppos covariates for positive and $pzero for zero counts")

  # add one coef for intercept
  ncoefpos = ppos + (intercept ? 1 : 0)
  ncoefzero = pzero + (intercept ? 1 : 0)

  covars, counts, μpos, μzero, n = shifters(M, covars, counts, showwarnings, m, l, D)

  # standardize covars only once if needed
  covars, covarsnorm = Lasso.standardizeX(covars, standardize)

  # fit separate GammaLassoPath's to each dimension of counts j=1:d and pick its min AICc segment
  if parallel
    verbose && @info("distributed $M run on local cluster with $(nworkers()) nodes")
    scounts = convert(SharedArray,counts)
    scoefszero = SharedMatrix{T}(ncoefzero,d)
    scoefspos = SharedMatrix{T}(ncoefpos,d)
    scovars = convert(SharedArray,covars)
    # μ = convert(SharedArray,μ) incompatible with GLM

    @sync @distributed for j=1:d
      tryfith!(M, scoefspos, scoefszero, j, scovars, scounts, inpos, inzero, μpos, μzero;
        verbose=false, showwarnings=showwarnings, intercept=intercept,
        standardize=false, select=select, kwargs...)
    end

    coefszero = convert(Matrix{T}, scoefszero)
    coefspos = convert(Matrix{T}, scoefspos)
  else
    verbose && @info("serial $M run on a single node")
    coefszero = Matrix{T}(undef,ncoefzero,d)
    coefspos = Matrix{T}(undef,ncoefpos,d)
    for j=1:d
      tryfith!(M, coefspos, coefszero, j, covars, counts, inpos, inzero, μpos, μzero;
        verbose=false, showwarnings=showwarnings, intercept=intercept,
        standardize=false, select=select, kwargs...)
    end
  end

  if standardize
    # destandardize coefs only once if needed
    destandardize!(coefspos, covarsnorm[inpos], standardize, intercept)
    destandardize!(coefszero, covarsnorm[inzero], standardize, intercept)
  end

  HDMRCoefs{M}(coefspos, coefszero, intercept, n, d, inpos, inzero, select)
end

"""
This version does not share memory across workers, so may be more efficient for
 small problems, or on remote clusters.
"""
function hdmr_remote_cluster(::Type{M}, covars::AbstractMatrix{T},counts::AbstractMatrix{V},
          inpos,inzero,intercept,parallel,verbose,showwarnings;
          select=defsegselect, kwargs...) where {T<:AbstractFloat,V,M<:TwoPartModel}
  paths = hdmrpaths(covars, counts, M; inpos=inpos, inzero=inzero, parallel=parallel, verbose=verbose, showwarnings=showwarnings, kwargs...)
  HDMRCoefs(paths, select)
end

"""
    posindic(A)

Returns an array of the same dimensions of indicators for positive entries in A.
"""
function posindic(A::AbstractArray{T}) where T
  # find positive y entries
  ixpos = (LinearIndices(A))[findall(x->x!=zero(T), A)]

  # build positive indicators matrix
  Ia = deepcopy(A)
  Ia[ixpos] .= one(T)
  Ia
end

function dropexplicitzeros(I, J, V::F) where {T, F<:AbstractVector{T}}
  if any(iszero, V)
    ixnonzero = broadcast(x->!iszero(x), V)
    @inbounds I = I[ixnonzero]
    @inbounds J = J[ixnonzero]
    nnonzero = sum(ixnonzero)
  else
    nnonzero = length(V)
  end

  I, J, nnonzero
end

"Sparse version simply replaces all the non-zero values with ones."
function posindic(A::SparseMatrixCSC{T}) where T
  m,n = size(A)
  I,J,V = findnz(A)
  I, J, nnonzero = dropexplicitzeros(I, J, V)
  sparse(I, J, fill(one(T),nnonzero), m, n)
end

ext: jl


text:
"""
    padded_tuple(v::AbstractVector, s)

Internal function used to compute shifts. Return a `Tuple` with as many element
as the dimensions of `v`. The first `length(s)` entries are filled with values
from `s`, the remaining entries are `0`. `s` should be an integer, in which case
`length(s) == 1`, or a container of integers with keys `1:length(s)`.

# Examples

```jldoctest padded_tuple
julia> ShiftedArrays.padded_tuple(rand(10, 10), 3)
(3, 0)

julia> ShiftedArrays.padded_tuple(rand(10, 10), (4,))
(4, 0)

julia> ShiftedArrays.padded_tuple(rand(10, 10), (1, 5))
(1, 5)
```
"""
padded_tuple(v::AbstractArray, s) = ntuple(i -> i ≤ length(s) ? s[i] : 0, ndims(v))

"""
    ShiftedArray(parent::AbstractArray, shifts, default)

Custom `AbstractArray` object to store an `AbstractArray` `parent` shifted by `shifts` steps (where `shifts` is
a `Tuple` with one `shift` value per dimension of `parent`).
For `s::ShiftedArray`, `s[i...] == s.parent[map(-, i, s.shifts)...]` if `map(-, i, s.shifts)`
is a valid index for `s.parent`, and `s.v[i, ...] == default` otherwise.
Use `copy` to collect the values of a `ShiftedArray` into a normal `Array`.
The recommended constructor is `ShiftedArray(parent, shifts; default = missing)`

# Examples

```jldoctest shiftedarray
julia> v = [1, 3, 5, 4];

julia> s = ShiftedArray(v, (1,))
4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
  missing
 1
 3
 5

julia> v = [1, 3, 5, 4];

julia> s = ShiftedArray(v, (1,))
4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
  missing
 1
 3
 5

julia> copy(s)
4-element Vector{Union{Missing, Int64}}:
  missing
 1
 3
 5

julia> v = reshape(1:16, 4, 4);

julia> s = ShiftedArray(v, (0, 2))
4×4 ShiftedArray{Int64, Missing, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
 missing  missing  1  5
 missing  missing  2  6
 missing  missing  3  7
 missing  missing  4  8

julia> shifts(s)
(0, 2)
```
"""
struct ShiftedArray{T, M, N, S<:AbstractArray} <: AbstractArray{Union{T, M}, N}
    parent::S
    shifts::NTuple{N, Int}
    default::M
end

ShiftedArray(v::AbstractArray{T, N}, n = (); default::M = missing) where {T, N, M} =
     ShiftedArray{T, M, N, typeof(v)}(v, padded_tuple(v, n), default)

"""
    ShiftedVector{T, S<:AbstractArray}

Shorthand for `ShiftedArray{T, 1, S}`.
"""
const ShiftedVector{T, M, S<:AbstractArray} = ShiftedArray{T, M, 1, S}

ShiftedVector(v::AbstractVector, n = (); default = missing) = ShiftedArray(v, n; default = default)

size(s::ShiftedArray) = size(parent(s))
axes(s::ShiftedArray) = axes(parent(s))

# Computing a shifted index (subtracting the offset)
offset(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) where {N} = map(-, inds, offsets)

@inline function getindex(s::ShiftedArray{<:Any, <:Any, N}, x::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(s, x...)
    v, i = parent(s), offset(shifts(s), x)
    return if checkbounds(Bool, v, i...)
        @inbounds v[i...]
    else
        default(s)
    end
end

parent(s::ShiftedArray) = s.parent

"""
    shifts(s::ShiftedArray)

Return amount by which `s` is shifted compared to `parent(s)`.
"""
shifts(s::ShiftedArray) = s.shifts

"""
    default(s::ShiftedArray)

Return default value.
"""
default(s::ShiftedArray) = s.default


text:
import Pkg
Pkg.instantiate()

include("server.jl")

function main()
    @info "starting server .."
    web_server()
end

main()

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
export YaoIR

"""
    YaoIR

The Yao Intermediate Representation. See compilation section for more details.

    YaoIR([m::Module=YaoLang.Compiler], ast::Expr)

Creates a `YaoIR` from Julia AST.
"""
mutable struct YaoIR
    mod::Module
    name::Any
    args::Vector{Any}
    whereparams::Vector{Any}
    body::IR
    pure_quantum::Bool
    qasm_compatible::Bool
end

function YaoIR(m::Module, ast::Expr)
    defs = splitdef(ast; throw = false)
    defs === nothing && throw(ParseError("expect function definition"))

    # potentially we could have code transform pass
    # on frontend AST as well here, but not necessary
    # for now, and all syntax related things should
    # go into to_function transformation
    ex = to_function(m, defs[:body])
    lowered_ast = Meta.lower(m, ex)

    if lowered_ast === nothing
        body = IR()
    else
        body = IR(lowered_ast.args[], 0)
    end

    ir = YaoIR(
        m,
        defs[:name],
        get(defs, :args, Any[]),
        get(defs, :whereparams, Any[]),
        mark_quantum(body),
        false,
        false,
    )
    update_slots!(ir)
    return ir
end

YaoIR(ast::Expr) = YaoIR(@__MODULE__, ast)

"""
    mark_quantum(ir::IR)

swap the statement tag with `:quantum`.
"""
function mark_quantum(ir::IR)
    for (v, st) in ir
        if (st.expr isa Expr) && (st.expr.head === :call) && (st.expr.args[1] isa GlobalRef)
            ref = st.expr.args[1]
            if ref.mod === Compiler && ref.name in RESERVED
                ir[v] = Statement(st; expr = Expr(:quantum, ref.name, st.expr.args[2:end]...))
            end
        end

        # mark quantum meta
        if (st.expr isa Expr) && (st.expr.head === :meta) && (st.expr.args[1] in RESERVED)
            ir[v] = Statement(st; expr = Expr(:quantum, st.expr.args...))
        end
    end
    return ir
end

function update_slots!(ir::YaoIR)
    fn_args = arguements(ir)
    for (v, st) in ir.body
        if st.expr isa Expr
            args = Any[]
            for each in st.expr.args
                if each in fn_args
                    push!(args, IRTools.Slot(each))
                else
                    push!(args, each)
                end
            end
            ir.body[v] = Statement(st; expr = Expr(st.expr.head, args...))
        elseif (st.expr isa Symbol) && (st.expr in fn_args)
            ir.body[v] = Statement(st; expr = IRTools.Slot(st.expr))
        end
    end
    return ir
end


text:
"""
    hubbardMPO(t, U, pchs, N[; T = ComplexF64])
returns a `U1`-symmetric MPO for the Hubbard model
```
    H = -∑_i t (a^†_i a_(i+1) + a_i a^†_(i+1)) + U/2 n_i (n_i-1)
```
on N sites.
"""
function hubbardMPO(t, U, pchs, N; T = ComplexF64)
    pdims = fill(1, length(pchs))
    bulk = DASTensor{T,4}(
        U1(),
        (U1Charges(-1:1), pchs, pchs, U1Charges(-1:1)),
        ([1,2,1], copy(pdims), copy(pdims), [1,2,1]),
        InOut(1,1,-1,-1))
    initwithzero!(bulk)
    preapp  = T[0 0; 0 1] # |0><0|_p ⊗ |1><1|_d
    postapp = T[1 0; 0 0] # |0><0|_p ⊗ |0><0|_d
    app     = T[0 0; 1 0] # |0><0|_p ⊗ |1><0|_d
    mat00to1 = reshape(T[0, 1], 2, 1)
    mat1to01 = reshape(T[1, 0], 1, 2)

    for u1ch in pchs
        ch = u1ch.ch
        peye= Matrix{T}(I, 1, 1)
        #local interaction term
        @tensor deg[1, 2, 3, 4] := preapp[1, 4]  * peye[2, 3] +  # |1><1| * Ι
                                   postapp[1, 4] * peye[2, 3] +  # |0><0| * Ι
                                   U/2 * ch * (ch-1)  * app[1, 4] * peye[2, 3]
                                   # U * ch * (ch-1) * |1><0|
        bulk[DASSector(U1Charge(0), U1Charge(ch), U1Charge(ch), U1Charge(0))] = deg

        #hopping term
        U1Charge(ch-1) in pchs || continue
        #a a^†
        # (0,0) → a_i → (1,0)
        @tensor deg[1, 2, 3, 4] := mat00to1[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(0), U1Charge(ch), U1Charge(ch-1), U1Charge(1))] = -t*deg
        # (1,0) → a^†_{i+1} → (0,1)
        @tensor deg[1, 2, 3, 4] := mat1to01[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(1), U1Charge(ch-1), U1Charge(ch), U1Charge(0))] = deg

        #a^† a
        # (0,0) → a^†_i → (-1,0)
        @tensor deg[1, 2, 3, 4] := mat00to1[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(0), U1Charge(ch-1), U1Charge(ch), U1Charge(-1))] = -t*deg
        # (-1,0) → a_{i+1} → (0,1)
        @tensor deg[1, 2, 3, 4] := mat1to01[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(-1), U1Charge(ch), U1Charge(ch-1), U1Charge(0))] = deg
    end
    # (3)        (4)       (3)
    # [lb](2) (2)[b](3) (2)[rb]
    # (1)        (1)       (1)
    lb = DASTensor{T,1}(U1(), (U1Charges(-1:1),), ([1,2,1],), InOut(-1))
    rb = DASTensor{T,1}(U1(), (U1Charges(-1:1),), ([1,2,1],), InOut( 1))
    lb[DASSector(U1Charge(0))] = [0, 1]
    rb[DASSector(U1Charge(0))] = [1, 0]

    @tensor bulk[1,2,3,4]       := bulk[2,1,4,3]
    @tensor leftbound[1, 2, 3]  := lb[-1] *  bulk[1, -1, 2, 3]
    @tensor rightbound[1, 2, 3] := bulk[1, 2, -1, 3] * rb[-1]

    return MPO{N}(leftbound, bulk,  rightbound)
end

"""
    creationOp(pchs::U1Charges [, T = ComplexF64])
return the creation operator `a†` acting on charges `pchs`
"""
function creationOp(pchs::U1Charges, T = ComplexF64)
    a = DASTensor{T,2}(U1(), (pchs, pchs),
        (fill(1, length(pchs)), fill(1, length(pchs))),
        InOut(1,-1))
    initwithzero!(a)
    for ch in pchs
        nch = ch ⊕ U1Charge(-1)
        nch in pchs || continue
        a[DASSector(nch, ch)] = reshape([sqrt(ch.ch)],1,1)
    end
    return a
end

"""
    annihilationOp(pchs::U1Charges [, T = ComplexF64])
return the annihilation operator `a` acting on charges `pchs`
"""
function annihilationOp(pchs::U1Charges, T = ComplexF64)
    a = DASTensor{T,2}(U1(), (pchs, pchs),
        (fill(1, length(pchs)), fill(1, length(pchs))),
        InOut(1,-1))
    initwithzero!(a)
    for ch in pchs
        nch = ch ⊕ U1Charge(-1)
        nch in pchs || continue
        a[DASSector(ch, nch)] = reshape([sqrt(ch.ch)],1,1)
    end
    return a
end


text:
# Also see related functions in follow.jl and personalization.jl

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-current-users-profile

"""
    user_get_current_profile()

    **Summary**: Get detailed profile information about the current user
                (including the current user's username).

# Example
```julia-repl
julia> Spotify.user_get_current_profile()[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 12 entries:
  :country          => "NL"
  :display_name     => "Itachi"
  :email            => "your_id@gmail.com"
  :explicit_content => {…
  :external_urls    => {…
  :followers        => {…
  :href             => "https://api.spotify.com/v1/users/your_user_id"
  :id               => "your_user_id"
  :images           => Union{}[]
  :product          => "premium"
  :type             => "user"
  :uri              => "spotify:user:your_user_id"
```
"""
function user_get_current_profile()

    return Spotify.spotify_request("me"; scope = "user-read-private")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-users-top-artists-and-tracks

# Implemented in personalization.jl
# top_tracks(offset=0, limit=20, time_range="medium")
# top_artists(offset=0, limit=20, time_range="medium")

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-users-profile

"""
    user_get_profile(user_id::String)

**Summary**: Get public profile information about a Spotify user.

# Arguments
- `user_id::String` : Alphanumeric ID of the user or name (e.g. "smedjan")

# Example
```julia-repl
julia> Spotify.user_get_profile("smedjan")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 8 entries:
  :display_name  => "smedjan"
  :external_urls => {…
  :followers     => {…
  :href          => "https://api.spotify.com/v1/users/smedjan"
  :id            => "smedjan"
  :images        => Union{}[]
  :type          => "user"
  :uri           => "spotify:user:smedjan"
```
"""
function user_get_profile(user_id)

    return Spotify.spotify_request("users/$user_id")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-followed

# Implemented in follow.jl
# follow_artists(type="artist", limit=20)

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-current-user-follows

# Implemented in follow.jl
# follow_check(type, ids)

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-if-user-follows-playlist

# Implemented in follow.jl
# follow_check_playlist(playlist_id, ids)


text:
using AdventOfCode21

AdventOfCode21.run(ARGS[1])


text:
module AppliInvoicing

greet() = print("Invoicing.jl")

export create, process, retrieve_unpaid_invoices, retrieve_paid_invoices, read_bank_statements, report

export id, meta, header, body, stm

# first, link to the model
include("./infrastructure/infrastructure.jl")

# next, submodule Reporting
include("Reporting.jl")

end # module


text:
using Test, bad

pnull = .25
pmcr  = .3
p     = Beta(5, 7)
α, β  = .05, .2
problem = Problem(
    minimise(
        SampleSize(p | (pnull + .2) )
    ),
    Power(p  | pnull) <= α,
    Power(p >= pmcr)  >= 1 - β
)
design = optimise(problem; verbosity = 0)

𝚾 = sample_space(design)

pme1 = PosteriorMean(p)
pme2 = PosteriorMeanPrecalculated(p, design)

@test all( pme1.(𝚾[:,1], 𝚾[:,2], design) ≈ pme2.(𝚾[:,1], 𝚾[:,2], design) )

pme3 = PosteriorMean(JeffreysPrior(design))
pme4 = PosteriorMeanPrecalculated(JeffreysPrior(design), design)
@test all( pme3.(𝚾[:,1], 𝚾[:,2], design) ≈ pme4.(𝚾[:,1], 𝚾[:,2], design) )

mle = MaximumLikelihoodEstimator()
rbe = RaoBlackwellEstimator()

x1_early_stop = early_stop_region(design)
pp = 0:.01:1
@test all( mle.(x1_early_stop, 0, design) .== rbe.(x1_early_stop, 0, design) )
# design is injective on continuation region, rbe = x1/n1?
# allow numerical inaccuracy for rbe!
@test all( rbe.(𝚾[:,1], 𝚾[:,2], design) .≈ 𝚾[:,1]./n1(design) )
@test maximum( abs.( bias.(pp, rbe, design) ) ) ≈ 0 atol = 1e-12
@test all(0 .<= rbe.(𝚾[:,1], 𝚾[:,2], design) .<= 1)

# build shan design
α, β     = .05, .1
pnull       = .6
shan_n1  = 19
shan_n2  = zeros(19 + 1)
shan_n2[(13:17) .+ 1] .= [33, 31, 31, 31, 14]
shan_c   = vcat(
        [Inf for x1 in 0:12],
        [36, 35, 35, 35, 22],
        [-Inf for x1 in 18:shan_n1]
)
design = Design(shan_n2, shan_c .- (0:shan_n1))

@test !mlecompatible(design, pnull, α)["compatible"]

cmle  = CompatibleMLE(design)
@test compatible(EstimatorOrdering(cmle), design, pnull, α)["compatible"]

𝚾     = sample_space(design)
resid = cmle.(𝚾[:,1], 𝚾[:,2], design) .- mle.(𝚾[:,1], 𝚾[:,2], design)
# make sure there is a (minute) difference
@test .001 < maximum(resid)

α, β = .05, .2

get_design(p_sample_size, pnull, palt) = Problem(
        minimise( SampleSize(p | p_sample_size) ),
        Power(p  | pnull) <= α,
        Power(p >= palt)  >= 1 - β;
        maxmultipleonestage = 2.5
    ) |> problem -> optimise(problem; verbosity = 0)

for (pnull, palt) in [(pnull, pnull + .2) for pnull in 0.1:.1:.7]

    for psamplesize in (pnull, palt)
        design = get_design(psamplesize, pnull, palt)
        cmle   = CompatibleMLE(design; lambda = .25)
        if !compatible(EstimatorOrdering(cmle), design, pnull, α)["compatible"]
            global tmp = design, pnull, palt, psamplesize
        end
        @test compatible(EstimatorOrdering(cmle), design, pnull, α)["compatible"]
        𝚾     = sample_space(design)
        resid = cmle.(𝚾[:,1], 𝚾[:,2], design) .- mle.(𝚾[:,1], 𝚾[:,2], design)
        if mlecompatible(design, pnull, α)["compatible"]
            @test maximum(resid) ≈ 0 atol = 1e-2 # numerical inaccuracy
        else
            @test maximum(resid) < 5*1e-1 # shouldnt be tooo different
        end
    end
end

@test true


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

using Markdown
using InteractiveUtils

# ╔═╡ 33bd93ec-77d9-11eb-2da8-87271f8b30c7
begin
  using Pkg
  Pkg.activate("MLJ_env", shared=true)

  # using MLJ
  using CSV
  using DataFrames
  using PlutoUI
  using Test
  using Printf
  using Plots
  using LinearAlgebra
end

# ╔═╡ 6d69b148-77d9-11eb-04f5-034603897c31
begin
  using JSON

  important_words = open("../../ML_UW_Spec/C03/data/important_words.json", "r") do f
      JSON.parse(f)
  end;

  length(important_words)
end

# ╔═╡ 0f6bb770-77df-11eb-0f1f-4dddbd7e2498
begin
	include("./text_proc.jl");
	include("./utils.jl");
end

# ╔═╡ baa8ca1c-77d8-11eb-2e26-3b2445acbfa1
md"""
## Training Logistic Regression via Stochastic Gradient Ascent

The goal of this notebook is to implement a logistic regression classifier using stochastic gradient ascent. We will:
  - Extract features from Amazon product reviews.
  - Convert an SFrame into a NumPy array.
  - Write a function to compute the derivative of log likelihood function with respect to a single coefficient.
  - Implement stochastic gradient ascent.
  - Compare convergence of stochastic gradient ascent with that of batch gradient ascent.

"""

# ╔═╡ 4218be80-7e19-11eb-1706-15ca470f7cf2
begin
	const DF = Union{DataFrame, SubDataFrame}
	const MV = Union{Matrix, Vector}
end

# ╔═╡ 3ad4770e-77d9-11eb-3a93-addf3e427025
md"""
### Load review dataset

For this assignment, we will use a subset of the Amazon product review dataset. The subset was chosen to contain similar numbers of positive and negative reviews, as the original dataset consisted primarily of positive reviews.
"""

# ╔═╡ 6e264178-77d9-11eb-33c1-19462dde66e2
begin
  products = train = CSV.File("../../ML_UW_Spec/C03/data/amazon_baby_subset.csv";
    header=true) |> DataFrame;
  size(products)
end

# ╔═╡ ba7417fe-7e18-11eb-1706-15ca470f7cf2
describe(products, :eltype, :nmissing, :first => first)

# ╔═╡ 36f2c616-77de-11eb-03b4-3b92a532949a
begin
	## how many missing values for the column :review?
	nmiss₀ = sum(ismissing.(products.review))

	## in-place replace(ment) of missing values for review  with empty string
	replace!(products.review, missing =>"");

	nmiss₁ = sum(ismissing.(products.review));
	(missing_before_repl=nmiss₀, missing_after_repl=nmiss₁)
end

# ╔═╡ d7b9c812-77de-11eb-0ebb-6bb75de9725b
## no more missing values for column :review
@test sum(ismissing.(products.review)) == 0

# ╔═╡ 6e09ca34-77d9-11eb-1a99-6578f9141004
md"""
One column of this dataset is 'sentiment', corresponding to the class label with +1 indicating a review with positive sentiment and -1 indicating one with negative sentiment.
"""

# ╔═╡ 6ded8978-77d9-11eb-3d94-e7507a85197e
products.sentiment[1:20]

# ╔═╡ 6dd12c56-77d9-11eb-0350-1dfc08d17dd6
md"""
Let us quickly explore more of this dataset. The 'name' column indicates the name of the product. Here we list the first 10 products in the dataset. We then count the number of positive and negative reviews.
"""

# ╔═╡ 6db87918-77d9-11eb-2082-edecdf8b7576
first(products, 10).name

# ╔═╡ 6d316f86-77d9-11eb-107b-1bf5708655d5
md"""
Now, we will perform 2 simple data transformations:

  1. Remove punctuation using Python's built-in string functionality.
  2. Compute word counts (only for important_words)

"""

# ╔═╡ 734ed4d4-77df-11eb-3dab-e5c10b4cda62
md"""
###### Step 1.
"""

# ╔═╡ 0f528d0e-77df-11eb-1415-2fb0e9d6792d
begin
	products[!, :review_clean] = remove_punctuation.(products[:, :review]);
	select!(products, Not(:review));

	names(products)
end

# ╔═╡ 0f1c6b16-77df-11eb-2197-93c4a82cc135
md"""
###### Step 2.

For each word in `important_words`, we compute a count for the number of times the word occurs in the review. We will store this count in a separate column (one for each word). The result of this feature processing is a single column for each word in `important_words` which keeps a count of the number of times the respective word occurs in the review text.

"""

# ╔═╡ 0057a50e-77e0-11eb-29d6-85c2dac67743
wc(text::String, word::String) = length(findall(word, text))

# ╔═╡ 51a94e8e-77e3-11eb-04c8-d1389b1f5856
function addcols!(df, words)
  for word ∈ words
    df[!, word] = wc.(df.review_clean, word)
  end
end

# ╔═╡ 0f03e474-77df-11eb-2791-f165fd011737
begin
	addcols!(products, important_words);
	length(names(products))
end

# ╔═╡ d39040fe-78d0-11eb-236d-15b506adaffd
md"""
### Train-Validation split

We split the data into a train-validation split with 90% of the data in the training set and 20% of the data in the validation set.
"""

# ╔═╡ 005dc370-78d1-11eb-0fc0-039f1017af8e
begin
	train_data, validation_data = train_test_split(products; split=0.9, seed=42);
	size(train_data), size(validation_data)
end

# ╔═╡ 682a78ba-78d1-11eb-1706-15ca470f7cf2
md"""
### Using Julia Arrays
"""

# ╔═╡ 0e383a04-77df-11eb-2ec9-c5aeffc4b131
function get_data(df::DF, features, output)
  s_features = [Symbol(f) for f ∈ features]
  df[:, :intercept] .= 1.0
  s_features = [:intercept, s_features...]
  X_matrix = convert(Matrix, select(df, s_features))
  y = df[!, output]

  (X_matrix, Vector{eltype(y)}(y))
end

# ╔═╡ 7ad4c752-78d1-11eb-267b-236c60163a33
begin
	# need a copy to add intercept, train_data is actually a view
	feature_matrix_train, sentiment_train = get_data(copy(train_data),
		important_words, :sentiment);
	# same remark as above concerning validation_data
	feature_matrix_valid, sentiment_valid = get_data(copy(validation_data),
		important_words, :sentiment);

	size(feature_matrix_train), size(sentiment_train)
end

# ╔═╡ b093b9ac-78d1-11eb-3a2e-ffb6e1102b1f
size(feature_matrix_valid), size(sentiment_valid)

# ╔═╡ 54c0f8f2-7e1b-11eb-1706-15ca470f7cf2
md"""
**Quiz Question**: In Module 3 assignment, there were 194 features (an intercept + one feature for each of the 193 important words).

In this assignment, we will use stochastic gradient ascent to train the classifier using logistic regression. How does the changing of the solver to stochastic gradient ascent affect the number of features?

  - [ ] Increases
  - [ ] Decreases
  - [x] Stays the same

"""

# ╔═╡ fd23d180-77e6-11eb-2091-af1e84d8fbcf
md"""
#### Building on logistic regression

Recall from lecture that the link function is given by:

$$P(y_i = +1 | \mathbf{x}_i,\mathbf{w}) = \frac{1}{1 + \exp(-\mathbf{w}^T h(\mathbf{x}_i))},$$

where the feature vector $h(\mathbf{x}_i)$ represents the word counts of `important_words` in the review  $\mathbf{x}_i$.

We will use the same code as in Module 3 assignment to make probability predictions, since this part is not affected by using stochastic gradient ascent as a solver. Only the way in which the coefficients are learned is affected by using stochastic gradient ascent as a solver.
"""

# ╔═╡ fd120ad4-77e6-11eb-15ea-bb7d7ea35587
function predict_probability(feature_matrix::MV, coeffs)
	@assert size(coeffs)[2] == 1  # expect a column vector not an array
  	1. ./ (1. .+ exp.(-(feature_matrix * coeffs)))
end

# ╔═╡ fcf1dd2e-77e6-11eb-0619-e7fd2a622877
md"""
#### Derivative of log likelihood with respect to a single coefficient

Let us now work on making minor changes to how the derivative computation is performed for logistic regression.

Recall from the lectures and Module 3 assignment that for logistic regression, **the derivative of log likelihood with respect to a single coefficient** is as follows:

$$\frac{\partial\ell}{\partial w_j} = \sum_{i=1}^N h_j(\mathbf{x}_i)\left(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w})\right)$$

In Module 3 assignment, we wrote a function to compute the derivative of log likelihood with respect to a single coefficient $w_j$. The function accepts the following two parameters:
 * `errors` vector containing $(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w}))$ for all $i$
 * `feature` vector containing $h_j(\mathbf{x}_i)$  for all $i$

Complete the following code block:
"""

# ╔═╡ d6976230-78d2-11eb-1fb7-dd1259c9b9c8
feature_derivative(errors, feature) = dot(errors, feature)

# ╔═╡ 0cea789c-78d4-11eb-36c4-c794c5e1a816
md"""
**Note**. We are not using regularization in this assignment, but, as discussed in the optional video, stochastic gradient can also be used for regularized logistic regression.
"""

# ╔═╡ d6157f9a-7e1b-11eb-267b-236c60163a33
md"""
To verify the correctness of the gradient computation, we provide a function for computing average log likelihood (which we recall from the last assignment was a topic detailed in an advanced optional video, and used here for its numerical stability).

To track the performance of stochastic gradient ascent, we provide a function for computing **average log likelihood**.

$$\ell\ell_A(\mathbf{w}) = {\color{red}{\frac{1}{N}}} \sum_{i=1}^N \Big( (\mathbf{1}[y_i = +1] - 1)\mathbf{w}^T h(\mathbf{x}_i) - \ln\left(1 + \exp(-\mathbf{w}^T h(\mathbf{x}_i))\right) \Big)$$

**Note** that we made one tiny modification to the log likelihood function (called `log_likelihood`) in our earlier assignments. We added a $\color{red}{1/N}$ term which averages the log likelihood accross all data points. The $\color{red}{1/N}$ term makes it easier for us to compare stochastic gradient ascent with batch gradient ascent. We will use this function to generate plots that are similar to those you saw in the lecture.
"""

# ╔═╡ 70a490b2-78d3-11eb-3acd-9bc21f53815e
function avg_log_likelihood(f_matrix::MV, sentiment, coeffs)
	indic = sentiment .== 1
 	scores = f_matrix * coeffs
  	log_exp = log.(1. .+ exp.(-scores))
 	mask = isinf.(log_exp)
  	log_exp[mask] = -scores[mask]
	sum((indic .- 1) .* scores .- log_exp) / length(f_matrix)
end

# ╔═╡ 191a7a84-7e1c-11eb-0389-5978263881b6
md"""
**Quiz Question:** Recall from the lecture and the earlier assignment, the log likelihood (without the averaging term) is given by

$$\ell\ell(\mathbf{w}) = \sum_{i=1}^N \Big( (\mathbf{1}[y_i = +1] - 1)\mathbf{w}^T h(\mathbf{x}_i) - \ln\left(1 + \exp(-\mathbf{w}^T h(\mathbf{x}_i))\right) \Big)$$

How are the functions $\ell\ell(\mathbf{w})$ and $\ell\ell_A(\mathbf{w})$ related?
  - answer: $\ell\ell_A(\mathbf{w}) = \frac{1}{N} \times \ell\ell(\mathbf{w})$
"""

# ╔═╡ f6a97742-7e21-11eb-2a06-cb616b0b43d8
md"""
#### Modifying the derivative for stochastic gradient ascent

Recall from the lecture that the gradient for a single data point $\color{red}{\mathbf{x}_i}$ can be computed using the following formula:

```math
\frac{\partial\ell_{\color{red}{i}}(\mathbf{w})}{\partial w_j} = h_j({\color{red}{\mathbf{x}_i}})

\left(
  \mathbf{1}
  [y_{\color{red}{i}} = +1] - P(y_{\color{red}{i}} = +1|{\color{red}{\mathbf{x}_i}},\mathbf{w})
  \right)
```

**Computing the gradient for a single data point**

Do we really need to re-write all our code to modify $\partial\ell(\mathbf{w})/\partial w_j$ to $\partial\ell_{\color{red}{i}}(\mathbf{w})/{\partial w_j}$?

Thankfully **No!**. we can access $\mathbf{x}_i$ in the training data using `feature_matrix_train[i:i, :]`
and $y_i$ in the training data using `sentiment_train[i:i]`. We can compute $\partial\ell_{\color{red}{i}}(\mathbf{w})/\partial w_j$ by re-using **all the code** written in **feature_derivative** and **predict_probability**.

We compute $\partial\ell_{\color{red}{i}}(\mathbf{w})/\partial w_j$ using the following steps:
* First, compute $P(y_i = +1 | \mathbf{x}_i, \mathbf{w})$ using the **predict_probability** function with `feature_matrix_train[i:i, :]` as the first parameter.
* Next, compute $\mathbf{1}[y_i = +1]$ using `sentiment_train[i:i]`.
* Finally, call the **feature_derivative** function with `feature_matrix_train[i:i,  j]` as one of the parameters.

Let us follow these steps for `j=1` and `i=10`:
"""

# ╔═╡ ba36988a-7e1c-11eb-12f9-1f5e36d0706d
begin
	const N_FEATURES = size(feature_matrix_train)[2]
	ϵ = 1e-8
	j₀ = 1                        # Feature number
	i₀ = 10                       # Data point number
	coeffs₀ = zeros(Float64, (N_FEATURES, 1))

	preds₀ = predict_probability(feature_matrix_train[i₀:i₀, :], coeffs₀)
	indic₀ = sentiment_train[i₀:i₀] .== 1
	errors₀ = indic₀ - preds₀

	∇_single_dp = feature_derivative(errors₀, feature_matrix_train[i₀:i₀, j₀])

	@test abs(∇_single_dp - 0.5) ≤ ϵ
end

# ╔═╡ ba1a42c0-7e1c-11eb-0209-6533fe570748
size(∇_single_dp)  # () means scalar

# ╔═╡ ba03f9ac-7e1c-11eb-1fb7-dd1259c9b9c8
md"""
**Quiz Question:** The code block above computed $\partial\ell_{\color{red}{i}}(\mathbf{w})/{\partial w_j}$ for `j = 1` and `i = 10`.  Is $\partial\ell_{\color{red}{i}}(\mathbf{w})/{\partial w_j}$ a scalar or a 194-dimensional vector?
  - scalar
"""

# ╔═╡ fd15b4ae-7e23-11eb-10cc-0deed0adb371
md"""
#### Modifying the derivative for using a batch of data points

Stochastic gradient estimates the ascent direction using 1 data point, while gradient uses $N$ data points to decide how to update the the parameters.  In an optional video, we discussed the details of a simple change that allows us to use a **mini-batch** of $B \leq N$ data points to estimate the ascent direction. This simple approach is faster than regular gradient but less noisy than stochastic gradient that uses only 1 data point. Although we encorage you to watch the optional video on the topic to better understand why mini-batches help stochastic gradient, in this assignment, we will simply use this technique, since the approach is very simple and will improve your results.

Given a mini-batch (or a set of data points) $\mathbf{x}_{i}, \mathbf{x}_{i+1} \ldots \mathbf{x}_{i+B}$, the gradient function for this mini-batch of data points is given by:

```math
{\color{red}{\sum_{s = i}^{i+B}}} \frac{\partial\ell_{s}}{\partial w_j} = {\color{red}{\sum_{s = i}^{i + B}}} h_j(\mathbf{x}_s)\left(\mathbf{1}[y_s = +1] - P(y_s = +1 | \mathbf{x}_s, \mathbf{w})\right)
```

**Computing the gradient for a "mini-batch" of data points**

Using NumPy, we access the points $\mathbf{x}_i, \mathbf{x}_{i+1} \ldots \mathbf{x}_{i+B}$ in the training data using `feature_matrix_train[i:i+B,:]`
and $y_i$ in the training data using `sentiment_train[i:i+B]`.

We can compute ${\color{red}{\sum_{s = i}^{i+B}}} \partial\ell_{s}/\partial w_j$ easily as follows:
"""

# ╔═╡ fcf0e0b4-7e23-11eb-1499-3558c5ea680c
begin
	B₁ = 10  # Mini-batch size
	coeffs₁ = zeros(Float64, (N_FEATURES, 1))

	preds₁ = predict_probability(feature_matrix_train[i₀:i₀ + B₁, :], coeffs₁)
	indic₁ = sentiment_train[i₀:i₀ + B₁] .== 1
	errors₁ = indic₁ - preds₁

	∇_single_mb = feature_derivative(errors₁, feature_matrix_train[i₀:i₀ + B₁, j₀])

	@test abs(∇_single_mb + 1.5) ≤ ϵ
end

# ╔═╡ b98e0666-7e1c-11eb-3a2e-ffb6e1102b1f
size(∇_single_mb)  # () means scalar

# ╔═╡ f2193cb4-7e24-11eb-27e4-55ec5250ae23
size(feature_matrix_train)[1]

# ╔═╡ e1508b4e-7e24-11eb-36ac-8d6c7b6cdaf9
md"""
**Quiz Question:** The code block above computed
${\color{red}{\sum_{s = i}^{i+B}}}\partial\ell_{s}(\mathbf{w})/{\partial w_j}$
for `j = 10`, `i = 10`, and `B = 10`. Is this a scalar or a 194-dimensional vector?
  - scalar

**Quiz Question:** For what value of `B` is the term
${\color{red}{\sum_{s = 1}^{B}}}\partial\ell_{s}(\mathbf{w})/\partial w_j$
the same as the full gradient $\partial\ell(\mathbf{w})/{\partial w_j}$?

Hint: consider the training set we are using now.
  -  B = size(feature_matrix_train)[1] == 47765
"""

# ╔═╡ 1e90e260-7e25-11eb-37ac-a320a56a628e
md"""
#### Averaging the gradient across a batch

It is a common practice to normalize the gradient update rule by the batch size B:

```math
\frac{\partial\ell_{\color{red}{A}}(\mathbf{w})}{\partial w_j} \approx {\color{red}{\frac{1}{B}}} {\sum_{s = i}^{i + B}} h_j(\mathbf{x}_s)\left(\mathbf{1}[y_s = +1] - P(y_s = +1 | \mathbf{x}_s, \mathbf{w})\right)
```

In other words, we update the coefficients using the **average gradient over data points** (instead of using a summation). By using the average gradient, we ensure that the magnitude of the gradient is approximately the same for all batch sizes. This way, we can more easily compare various batch sizes of stochastic gradient ascent (including a batch size of **all the data points**), and study the effect of batch size on the algorithm as well as the choice of step size.

## Implementing stochastic gradient ascent

Now we are ready to implement our own logistic regression with stochastic gradient ascent. Complete the following function to fit a logistic regression model using gradient ascent:
"""

# ╔═╡ 1e437548-7e25-11eb-3e76-49a4a2f4094f
function shuffle_data(f_matrix::MV, sentiment::Vector; seed=42)
    "Set specific seed to produce consistent results"
    n = size(f_matrix)[1]
	perm = randperm(MersenneTwister(seed), n)
    f_matrix = f_matrix[perm, :]
    sentiment = sentiment[perm]
    (f_matrix, sentiment)
end

# ╔═╡ 970fba24-77ed-11eb-1a3e-e769d6049078
function logistic_regression_SG(f_matrix::MV, sentiment::Vector, init_coeffs;
		η=1e-7, batch_sz=100, max_iter=200)
	#
	coeffs = copy(init_coeffs)
	n = size(f_matrix)[1]
	f_matrix, sentiment = shuffle_data(f_matrix, sentiment)
	ll_all = Float64[]
	ix = 1
  	for it ∈ 1:max_iter
		n_ix = ix + batch_sz - 1

    	preds = predict_probability(f_matrix[ix:n_ix, :], coeffs)
    	indic = sentiment[ix:n_ix] .== 1
    	errors = indic .- preds

    	for jx ∈ 1:length(coeffs)
	  		is_intercept = (jx == 1)
      		∂c = feature_derivative(errors, f_matrix[ix:n_ix, jx])
      		coeffs[jx] += (1. / batch_sz) * η * ∂c
    	end

		ll = avg_log_likelihood(f_matrix[ix:n_ix, :], sentiment[ix:n_ix], coeffs)
        push!(ll_all, ll)

    	if it ≤ 15 || (it ≤ 100 && it % 10 == 0) || (it ≤ 1000 && it % 100 == 0) ||
      	(it ≤ 10000 && it % 1000 == 0) || it % 10000 == 0
      		@printf("iter: %2d: avg log likelihood (of data points in batch [%4d:%4d]): %.8f\n", it, ix, n_ix, ll)
    	end

		ix += batch_sz - 1

		if ix + batch_sz - 1 > n
			f_matrix, sentiment = shuffle_data(f_matrix, sentiment)
			ix = 1
		end
  	end

  	(coeffs, ll_all)
end

# ╔═╡ e6ef9818-7e26-11eb-08c1-13392c483a14
md"""
*Remark*: In practice, the final set of coefficients is rarely used; it is better to use the average of the last K sets of coefficients instead, where K should be adjusted depending on how fast the log likelihood oscillates around the optimum.
"""

# ╔═╡ e6d10c2c-7e26-11eb-05fa-d953b60010ff
md"""
##### Checkpoint

The following cell tests your stochastic gradient ascent function using a toy dataset consisting of two data points. If the test does not pass, make sure you are normalizing the gradient update rule correctly.
"""

# ╔═╡ e6b80b32-7e26-11eb-2b8f-a9b22b5c46fe
begin
	f_matrix₂ = Float64[1. 2. -1.; 1. 0. 1.]
	sentiment₂ = Int[1, -1]
	init_coeffs₂=zeros(Float64, (3, 1))

	coeffs₂, log_ll₂ = logistic_regression_SG(f_matrix₂, sentiment₂, init_coeffs₂;
		η=1., batch_sz=2, max_iter=2)

	ϵ₁ = 1e-3
	exp_coeffs₂ = [-0.0975576 0.682426 -0.779983]
	exp_log_ll₂ = [-0.112582 -0.0781844]

	# coeffs₂, log_ll₂
	@test all(((cₐ, cₑ)=t) -> abs(cₐ - cₑ) ≤ ϵ₁, zip(coeffs₂, exp_coeffs₂))
	@test all(((llₐ, llₑ)=t) -> abs(llₐ - llₑ) ≤ ϵ₁, zip(log_ll₂, exp_log_ll₂))
end

# ╔═╡ e696bce6-7e26-11eb-218d-ff1d2947d397
with_terminal() do
	println("----------------------------------------------------------------")
	println("Coefficients learned: $(coeffs₂)")
	println("Average log likelihood per-iteration: $(log_ll₂)")
	println("----------------------------------------------------------------")
	print("Test passed!")
end

# ╔═╡ 3cbf62b8-7e2d-11eb-3a2e-ffb6e1102b1f
md"""
#### Compare convergence behavior of stochastic gradient ascent

For the remainder of the assignment, we will compare stochastic gradient ascent against batch gradient ascent. For this, we need a reference implementation of batch gradient ascent. But do we need to implement this from scratch?

**Quiz Question:** For what value of batch size `B` above is the stochastic gradient ascent function `logistic_regression_SG` act as a standard gradient ascent algorithm?

Hint: consider the training set we are using now.
"""

# ╔═╡ 54f2b0a6-7e2d-11eb-1fb7-dd1259c9b9c8
# answer B should be:
size(train_data)[1]

# ╔═╡ 7f298886-7e2d-11eb-0209-6533fe570748
md"""
#### Running gradient ascent using the stochastic gradient ascent implementation

Instead of implementing batch gradient ascent separately, we save time by re-using the stochastic gradient ascent function we just wrote — **to perform gradient ascent**, it suffices to set **`batch_size`** to the number of data points in the training data. Yes, we did answer above the quiz question for you, but that is an important point to remember in the future :)

**Small Caveat**. The batch gradient ascent implementation here is slightly different than the one in the earlier assignments, as we now normalize the gradient update rule.

We now **run stochastic gradient ascent** over the `feature_matrix_train` for 10 iterations using:
  - `init_coeffs = zeros(Float64, (194, 1))`
  - `η = 5e-1`
  - `batch_sz = 1`
  - `max_iter = 10`

"""

# ╔═╡ b762f3d6-7e2d-11eb-3906-5ba3566b9ae9
begin
	init_coeffs₃ = zeros(Float64, (N_FEATURES, 1))
	local coeffs₃, log_likelihood₃

	with_terminal() do
		coeffs₃, log_likelihood₃ = logistic_regression_SG(feature_matrix_train,
			sentiment_train, init_coeffs₃,
			η=5e-1, batch_sz=1, max_iter=10);
	end
end

# ╔═╡ b715ba4e-7e2d-11eb-36c4-c794c5e1a816
md"""
**Quiz Question**. When you set `batch_size = 1`, as each iteration passes, how does the average log likelihood in the batch change?
  - [x] Increases
  - [ ] Decreases
  - [ ] Fluctuates
"""

# ╔═╡ b6f8ad50-7e2d-11eb-3acd-9bc21f53815e
md"""
Now run **batch gradient ascent** over the `feature_matrix_train` for 200 iterations using:

  - `init_coeffs = zeros(Float64, (194, 1))`
  - `η = 5e-1`
  - `batch_sz = size(feature_matrix_train)[1]`
  - `max_iter = 200`

"""

# ╔═╡ b6e038ec-7e2d-11eb-13d1-8dafac6a1d6a
begin
	# takes a while ≈ 98.6s
	init_coeffs₄ = zeros(Float64, (N_FEATURES, 1))
	coeffs₄, log_likelihood₄ = nothing, nothing

	with_terminal() do
		global coeffs₄, log_likelihood₄ = logistic_regression_SG(feature_matrix_train,
			sentiment_train, init_coeffs₄,
			η=5e-1, batch_sz=size(feature_matrix_train)[1], max_iter=200);
	end
end

# ╔═╡ b6c44948-7e2d-11eb-12f9-1f5e36d0706d
md"""
**Quiz Question**. When you set `batch_size = len(feature_matrix_train)`, as each iteration passes, how does the average log likelihood in the batch change?
  - [x] Increases
  - [ ] Decreases
  - [ ] Fluctuates
"""

# ╔═╡ 14f788e6-7e30-11eb-27e4-55ec5250ae23
md"""
#### Make "passes" over the dataset

To make a fair comparison betweeen stochastic gradient ascent and batch gradient ascent, we measure the average log likelihood as a function of the number of passes (defined as follows):

$$[\text{\# of passes}] = \frac{[\text{\# of data points touched so far}]}{[\text{size of dataset}]}$$

"""

# ╔═╡ 14dcb3f6-7e30-11eb-36ac-8d6c7b6cdaf9
md"""
**Quiz Question** Suppose that we run stochastic gradient ascent with a batch size of 100. How many gradient updates are performed at the end of two passes over a dataset consisting of 50000 data points?
"""

# ╔═╡ 14c1d310-7e30-11eb-10cc-0deed0adb371
# Answer:
2 * 50_000 ÷ 100

# ╔═╡ 14a70a76-7e30-11eb-1499-3558c5ea680c
md"""
#### Log likelihood plots for stochastic gradient ascent

With the terminology in mind, let us run stochastic gradient ascent for 10 passes. We will use
  - `η=1e-1`
  - `batch_sz=100`
  - `init_coeffs` to all zeros.
"""

# ╔═╡ 148f157e-7e30-11eb-1d8a-1f920b745781
begin
	n_passes₅ = 10
	batch_sz₅ = 100
	n_iter₅ = n_passes₅ * size(feature_matrix_train)[1] ÷ batch_sz₅
	init_coeffs₅ = zeros(Float64, (N_FEATURES, 1))
	coeffs_sgd₅, log_ll_sgd₅ = nothing, nothing

	with_terminal() do
		global coeffs_sgd₅, log_ll_sgd₅ = logistic_regression_SG(feature_matrix_train, 	sentiment_train, init_coeffs₅; η=1e-1, batch_sz=batch_sz₅, max_iter=n_iter₅)
	end
end

# ╔═╡ 09444fda-7e31-11eb-3e76-49a4a2f4094f
# TODO plots
log_ll_sgd₅

# ╔═╡ eae04a48-7e31-11eb-08c1-13392c483a14
md"""
#### Smoothing the stochastic gradient ascent curve

The plotted line oscillates so much that it is hard to see whether the log likelihood is improving. In our plot, we apply a simple smoothing operation using the parameter smoothing_window.
The smoothing is simply a moving average of log likelihood over the last smoothing_window "iterations" of stochastic gradient ascent.

"""

# ╔═╡ a4b3fe7a-7e45-11eb-0209-6533fe570748
function convolve(v₁::AbstractVector, v₂::AbstractVector; mode=:valid)
	if mode != :valid
		# other modes: :same, :full - not required here
		throw(ArgumentError("mode $(mode) not implemented"))
	end
	nv₁, nv₂ = length(v₁), length(v₂)
	@assert nv₁ ≤ nv₂

	vs₁ = sort(v₁, rev=true)
	nv = zeros(eltype(v₁), nv₂ - nv₁ + 1)

	for ix ∈ 1:nv₂ - nv₁ + 1
		nv[ix] = dot(vs₁, v₂[ix:ix+nv₁ - 1])
	end
	nv
end

# ╔═╡ c83d6ece-7e47-11eb-12f9-1f5e36d0706d
@test convolve([1.; 2.; 3.], [0.; 1.; 0.5]) == Float64[2.5]

# ╔═╡ 4f1d3fb4-7e48-11eb-13d1-8dafac6a1d6a
@test convolve([3.; 7.], [1.; 2.; 5.; 7.]) == [13.0, 29.0, 56.0]

# ╔═╡ 498064f0-7e4d-11eb-3a2e-ffb6e1102b1f
function apply_conv(ll_all, swin)
	n = length(ll_all)
	ll_all_ma = convolve(ones(eltype(ll_all), swin) ./ swin, ll_all, mode=:valid)
	(ll_all_ma, n)
end

# ╔═╡ 32e899e0-7e45-11eb-1fb7-dd1259c9b9c8
function make_plot(ll_all, len, batch_sz; swin=1, label="", leg=:topright, over=false)
	ll_all_ma, n = apply_conv(ll_all, swin)
	plot(collect(1:n - swin + 1) * batch_sz / len, ll_all_ma,
		linewidth=3.0, label=label, leg=leg,
		ylabel="Average log likelihood per data point",
		xlabel="# of passes over data")
end

# ╔═╡ 334dece8-7e4d-11eb-0389-5978263881b6
function make_plot!(ll_all, len, batch_sz; swin=1, label="", leg=:topright)
	ll_all_ma, n = apply_conv(ll_all, swin)
	plot!(collect(1:n - swin + 1) * batch_sz / len, ll_all_ma,
		linewidth=3.0, label=label, leg=leg,
		ylabel="Average log likelihood per data point",
		xlabel="# of passes over data")
end

# ╔═╡ d98e2fe0-7e49-11eb-3906-5ba3566b9ae9
make_plot(log_ll_sgd₅, size(feature_matrix_train)[1], batch_sz₅;
	label="stochastic gradient, step_size=1e-1", leg=:topleft)

# ╔═╡ 11f38c2a-7e4c-11eb-1706-15ca470f7cf2
md"""
###### Smoothing the stochastic gradient ascent curve

The plotted line oscillates so much that it is hard to see whether the log likelihood is improving. In our plot, we apply a simple smoothing operation using the parameter smoothing_window.

The smoothing is simply a moving average of log likelihood over the last smoothing_window (`swin`) "iterations" of stochastic gradient ascent.

"""

# ╔═╡ 1a250d88-7e4c-11eb-267b-236c60163a33
make_plot(log_ll_sgd₅, size(feature_matrix_train)[1], batch_sz₅;
          swin=30, label="stochastic gradient, step_size=1e-1", leg=:bottomright)

# ╔═╡ eab39b24-7e31-11eb-2b8f-a9b22b5c46fe
md"""
#### Stochastic gradient ascent vs batch gradient ascent

To compare convergence rates for stochastic gradient ascent with batch gradient ascent, we call `make_plot()` multiple times in the same cell.

We are comparing:
  * **stochastic gradient ascent**: `η = 0.1`, `batch_sz=100`
  * **batch gradient ascent**: `step_size = 0.5`, `batch_size=len(feature_matrix_train)`

Write code to run stochastic gradient ascent for 200 passes using:
  - `η=1e-1`
  - `batch_sz=100`
  - `init_coeffs` to all zeros.
"""

# ╔═╡ 042eab34-7e32-11eb-1f51-2729cf921bc5
begin
	batch_sz₆ = 100
	n_passes₆ = 200
	n_iter₆ = n_passes₆ * size(feature_matrix_train)[1] ÷ batch_sz₆
	coeffs_sgd₆, log_ll_sgd₆ = nothing, nothing

	with_terminal() do
		global coeffs_sgd₆, log_ll_sgd₆ = logistic_regression_SG(feature_matrix_train, sentiment_train, init_coeffs₅; η=1e-1, batch_sz=batch_sz₆, max_iter=n_iter₆)
	end
end

# ╔═╡ ea95b6a4-7e31-11eb-218d-ff1d2947d397
begin
	make_plot(log_ll_sgd₆, size(feature_matrix_train)[1], batch_sz₆;
          swin=30, label="batch, step_size=1e-1", leg=:bottomright)

	make_plot!(log_likelihood₄, size(feature_matrix_train)[1],
		size(feature_matrix_train)[1];
		swin=30, label="stochastic gradient, step_size=5e-1", leg=:bottomright)
end

# ╔═╡ f360444a-7e41-11eb-1706-15ca470f7cf2
md"""
#### Explore the effects of step sizes on stochastic gradient ascent

In previous sections, we chose step sizes for you. In practice, it helps to know how to choose good step sizes yourself.

To start, we explore a wide range of step sizes that are equally spaced in the log space. Run stochastic gradient ascent with `step_size` set to 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, and 1e2. Use the following set of parameters:
  - `init_coeff all zeros`
  - `batch_sz=100`
  - `max_iter` initialized so as to run 10 passes over the data.
"""

# ╔═╡ 22472020-7e42-11eb-267b-236c60163a33
begin
	n_passes₇ = 10
	batch_sz₇ = 100
	n_iter₇ = n_passes₇ * size(feature_matrix_train)[1] ÷ batch_sz₇
	init_coeffs₇ = zeros(Float64, (N_FEATURES, 1))
	c_sgd, log_ll_sgd = Dict{Float64, Any}(), Dict{Float64, Any}()

	const ETA = [1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2]
	with_terminal() do
		for η ∈ ETA
			println("\tη: $(η)")
			global c_sgd[η], log_ll_sgd[η] = logistic_regression_SG( 								feature_matrix_train,
				sentiment_train, init_coeffs₇;
				η=η, batch_sz=batch_sz₇, max_iter=n_iter₇
			)
		end
	end
end

# ╔═╡ b9985a96-7e42-11eb-3a2e-ffb6e1102b1f
md"""
#### Plotting the log likelihood as a function of passes for each step size

Now, we will plot the change in log likelihood using the `make_plot` for each of the following values of η (cf. cell above):

For consistency, we again apply smoothing_window: `swin=30`.
"""

# ╔═╡ 3940bc4a-7e50-11eb-12f9-1f5e36d0706d
begin
	len = size(feature_matrix_train)[1]
	η = 1e-4
	make_plot(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)

	η = 1e-3
	make_plot!(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)

	η = 1e-2
	make_plot!(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)

	η = 1e-1
	make_plot!(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)

	η = 1e0
	make_plot!(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)

	η = 1e1
	make_plot!(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)

	η = 1e2
	make_plot!(log_ll_sgd[η], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η)", leg=:bottomright)
end

# ╔═╡ 848d7b2c-7e53-11eb-2a06-cb616b0b43d8
md"""
Now, let us remove the step size η=1e2 and η=1e1 plot the rest of the curves.
"""

# ╔═╡ 8ce3d06e-7e53-11eb-2873-9fd54b4485ea
begin
	η₁ = 1e-4
	make_plot(log_ll_sgd[η₁], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)

	η₁ = 1e-3
	make_plot!(log_ll_sgd[η₁], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)

	η₁ = 1e-2
	make_plot!(log_ll_sgd[η₁], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)

	η₁ = 1e-1
	make_plot!(log_ll_sgd[η₁], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)

	η₁ = 1e0
	make_plot!(log_ll_sgd[η₁], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)

	# η₁ = 1e1
	# make_plot!(log_ll_sgd[η₁], len, batch_sz₇;
	# 		swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)

	# η₁ = 1e2
	# make_plot!(log_ll_sgd[η₁], len, batch_sz₇;
	# 		swin=30, label="batch, step_size=$(η₁)", leg=:bottomright)
end

# ╔═╡ ebacd9c4-7e53-11eb-1499-3558c5ea680c
begin
	η₂ = 1e1
	make_plot(log_ll_sgd[η₂], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₂)", leg=:bottomright)

	η₂ = 1e2
	make_plot!(log_ll_sgd[η₂], len, batch_sz₇;
			swin=30, label="batch, step_size=$(η₂)", leg=:bottomright)
end

# ╔═╡ dea572a2-7e53-11eb-1d8a-1f920b745781
md"""
**Quiz Question**: Which of the following is the worst step size? Pick the step size that results in the lowest log likelihood in the end.
  1. [ ] 1e-2
  2. [ ] 1e-1
  3. [ ] 1e0
  4. [ ] 1e1
  5. [x] 1e2 (the one we eliminayed...)

**Quiz Question**: Which of the following is the best step size? Pick the step size that results in the highest log likelihood in the end.
  1. [ ] 1e-4
  2. [ ] 1e-2
  3. [x] 1e0
  4. [ ] 1e1
  5. [ ] 1e2
"""

# ╔═╡ Cell order:
# ╟─baa8ca1c-77d8-11eb-2e26-3b2445acbfa1
# ╠═33bd93ec-77d9-11eb-2da8-87271f8b30c7
# ╠═0f6bb770-77df-11eb-0f1f-4dddbd7e2498
# ╠═4218be80-7e19-11eb-1706-15ca470f7cf2
# ╟─3ad4770e-77d9-11eb-3a93-addf3e427025
# ╠═6e264178-77d9-11eb-33c1-19462dde66e2
# ╠═ba7417fe-7e18-11eb-1706-15ca470f7cf2
# ╠═36f2c616-77de-11eb-03b4-3b92a532949a
# ╠═d7b9c812-77de-11eb-0ebb-6bb75de9725b
# ╟─6e09ca34-77d9-11eb-1a99-6578f9141004
# ╠═6ded8978-77d9-11eb-3d94-e7507a85197e
# ╟─6dd12c56-77d9-11eb-0350-1dfc08d17dd6
# ╠═6db87918-77d9-11eb-2082-edecdf8b7576
# ╠═6d69b148-77d9-11eb-04f5-034603897c31
# ╟─6d316f86-77d9-11eb-107b-1bf5708655d5
# ╟─734ed4d4-77df-11eb-3dab-e5c10b4cda62
# ╠═0f528d0e-77df-11eb-1415-2fb0e9d6792d
# ╟─0f1c6b16-77df-11eb-2197-93c4a82cc135
# ╠═0057a50e-77e0-11eb-29d6-85c2dac67743
# ╠═51a94e8e-77e3-11eb-04c8-d1389b1f5856
# ╠═0f03e474-77df-11eb-2791-f165fd011737
# ╟─d39040fe-78d0-11eb-236d-15b506adaffd
# ╠═005dc370-78d1-11eb-0fc0-039f1017af8e
# ╟─682a78ba-78d1-11eb-1706-15ca470f7cf2
# ╠═0e383a04-77df-11eb-2ec9-c5aeffc4b131
# ╠═7ad4c752-78d1-11eb-267b-236c60163a33
# ╠═b093b9ac-78d1-11eb-3a2e-ffb6e1102b1f
# ╟─54c0f8f2-7e1b-11eb-1706-15ca470f7cf2
# ╟─fd23d180-77e6-11eb-2091-af1e84d8fbcf
# ╠═fd120ad4-77e6-11eb-15ea-bb7d7ea35587
# ╟─fcf1dd2e-77e6-11eb-0619-e7fd2a622877
# ╠═d6976230-78d2-11eb-1fb7-dd1259c9b9c8
# ╟─0cea789c-78d4-11eb-36c4-c794c5e1a816
# ╟─d6157f9a-7e1b-11eb-267b-236c60163a33
# ╠═70a490b2-78d3-11eb-3acd-9bc21f53815e
# ╟─191a7a84-7e1c-11eb-0389-5978263881b6
# ╟─f6a97742-7e21-11eb-2a06-cb616b0b43d8
# ╠═ba36988a-7e1c-11eb-12f9-1f5e36d0706d
# ╠═ba1a42c0-7e1c-11eb-0209-6533fe570748
# ╟─ba03f9ac-7e1c-11eb-1fb7-dd1259c9b9c8
# ╟─fd15b4ae-7e23-11eb-10cc-0deed0adb371
# ╠═fcf0e0b4-7e23-11eb-1499-3558c5ea680c
# ╠═b98e0666-7e1c-11eb-3a2e-ffb6e1102b1f
# ╠═f2193cb4-7e24-11eb-27e4-55ec5250ae23
# ╟─e1508b4e-7e24-11eb-36ac-8d6c7b6cdaf9
# ╟─1e90e260-7e25-11eb-37ac-a320a56a628e
# ╠═1e437548-7e25-11eb-3e76-49a4a2f4094f
# ╠═970fba24-77ed-11eb-1a3e-e769d6049078
# ╟─e6ef9818-7e26-11eb-08c1-13392c483a14
# ╟─e6d10c2c-7e26-11eb-05fa-d953b60010ff
# ╠═e6b80b32-7e26-11eb-2b8f-a9b22b5c46fe
# ╠═e696bce6-7e26-11eb-218d-ff1d2947d397
# ╟─3cbf62b8-7e2d-11eb-3a2e-ffb6e1102b1f
# ╠═54f2b0a6-7e2d-11eb-1fb7-dd1259c9b9c8
# ╟─7f298886-7e2d-11eb-0209-6533fe570748
# ╠═b762f3d6-7e2d-11eb-3906-5ba3566b9ae9
# ╟─b715ba4e-7e2d-11eb-36c4-c794c5e1a816
# ╟─b6f8ad50-7e2d-11eb-3acd-9bc21f53815e
# ╠═b6e038ec-7e2d-11eb-13d1-8dafac6a1d6a
# ╟─b6c44948-7e2d-11eb-12f9-1f5e36d0706d
# ╟─14f788e6-7e30-11eb-27e4-55ec5250ae23
# ╟─14dcb3f6-7e30-11eb-36ac-8d6c7b6cdaf9
# ╠═14c1d310-7e30-11eb-10cc-0deed0adb371
# ╟─14a70a76-7e30-11eb-1499-3558c5ea680c
# ╠═148f157e-7e30-11eb-1d8a-1f920b745781
# ╠═09444fda-7e31-11eb-3e76-49a4a2f4094f
# ╟─eae04a48-7e31-11eb-08c1-13392c483a14
# ╠═a4b3fe7a-7e45-11eb-0209-6533fe570748
# ╠═c83d6ece-7e47-11eb-12f9-1f5e36d0706d
# ╠═4f1d3fb4-7e48-11eb-13d1-8dafac6a1d6a
# ╠═498064f0-7e4d-11eb-3a2e-ffb6e1102b1f
# ╠═32e899e0-7e45-11eb-1fb7-dd1259c9b9c8
# ╠═334dece8-7e4d-11eb-0389-5978263881b6
# ╠═d98e2fe0-7e49-11eb-3906-5ba3566b9ae9
# ╟─11f38c2a-7e4c-11eb-1706-15ca470f7cf2
# ╠═1a250d88-7e4c-11eb-267b-236c60163a33
# ╟─eab39b24-7e31-11eb-2b8f-a9b22b5c46fe
# ╠═042eab34-7e32-11eb-1f51-2729cf921bc5
# ╠═ea95b6a4-7e31-11eb-218d-ff1d2947d397
# ╟─f360444a-7e41-11eb-1706-15ca470f7cf2
# ╠═22472020-7e42-11eb-267b-236c60163a33
# ╟─b9985a96-7e42-11eb-3a2e-ffb6e1102b1f
# ╠═3940bc4a-7e50-11eb-12f9-1f5e36d0706d
# ╟─848d7b2c-7e53-11eb-2a06-cb616b0b43d8
# ╠═8ce3d06e-7e53-11eb-2873-9fd54b4485ea
# ╠═ebacd9c4-7e53-11eb-1499-3558c5ea680c
# ╟─dea572a2-7e53-11eb-1d8a-1f920b745781


text:
@reexport module Lists

import FunctionalCollections: append
using Base.Iterators

const ⊚ = Base.map
const car = Base.first

export Cons, List, isnil, ispair, car, cdr, caar, cadr, cddr, nil, lispify,
       append, ++, ⊚, cons, list, islist, Nil

lispify(x) = x
lispify(b::Bool) = b
lispify(i::Integer) = BigInt(i)
lispify(t::Tuple) = isempty(t) ? nil : Cons(lispify(t[1]), lispify(t[2:end]))

struct Nil end
const nil = Nil.instance

struct Cons
    car
    cdr
end
const cons = Cons

car(α::Cons) = α.car
cdr(α::Cons) = α.cdr
caar(α::Cons) = car(car(α))
cadr(α::Cons) = car(cdr(α))
cddr(α::Cons) = cdr(cdr(α))

Base.map(f, ::Nil) = nil
Base.map(f, α::Cons) = Cons(f(car(α)), f ⊚ cdr(α))
Base.filter(p, ::Nil) = nil
Base.filter(p, α::Cons) = let β = filter(p, cdr(α))
    p(car(α)) ? Cons(car(α), β) : β
end

Base.:(==)(α::Cons, β::Cons) = car(α) == car(β) && cdr(α) == cdr(β)

const List = Union{Cons, Nil}

List() = nil
List(xs...) = Cons(lispify(xs[1]), List(xs[2:end]...))
const list = List

isnil(::Nil) = true
isnil(::Cons) = false

ispair(α::List) = !isnil(α)

islist(::Nil) = true
islist(α::Cons) = islist(cdr(α))

Base.iterate(α::List, β=α) = isnil(β) ? nothing : (car(β), cdr(β))

Base.length(::Nil) = 0
Base.length(α::Cons) = 1 + length(cdr(α))

Base.getindex(α::Nil, b) = throw(BoundsError(α, b))
Base.getindex(α::Cons, b) = b == 1 ? car(α) : cdr(α)[b - 1]

unparse(α::List) = "(" * join(unparse ⊚ α, " ") * ")"
unparse(b::Bool) = b ? "#t" : "#f"
unparse(s::Symbol) = string(s)
unparse(s::String) = repr(s)
unparse(i::BigInt) = string(i)
unparse(::Nothing) = "#<void>"

Base.convert(::Type{List}, xs::List) = xs
Base.convert(::Type{List}, xs) = List(xs...)
append(::Nil, β::List) = β
append(α::Cons, β::List) = Cons(car(α), append(cdr(α), β))
append(α::List, β::List, γ::List, γs::List...) = append(append(α, β), γ, γs...)
const (++) = append

include("show.jl")

end  # Lists module


text:
function dq_ri(δ)
    ## Uses the referenceframe of the Kundur page 852 of dq to RI
    return [
        sin(δ) cos(δ)
        -cos(δ) sin(δ)
    ]
end

function ri_dq(δ)
    #Uses the reference frame of the Kundur page 852 of RI to dq
    return [
        sin(δ) -cos(δ)
        cos(δ) sin(δ)
    ]
end


text:
using Plots
rastrigin(x) = 10length(x) + sum(x.^2 .- 10 .* cos.(2π .* x))

dpi = 600
x = y = range(-3, stop=3, length=50)

l = Vector{Float64}()
for i in x
   for j in y
       v = [i j]
       append!(l, 1/rastrigin(v))
   end
end

surface(x,y,l, dpi=dpi)
png("./inverse_rastrigin_surface.png")
contour(x,y,l, dpi=dpi)
png("./inverse_rastrigin_contour.png")


text:
using Statistics

@time begin
    data = [mean(rand(5*10^2)) for _ in 1:10^6]
    println("98% of the means lie in the estimated range: ",
                    (quantile(data,0.01),quantile(data,0.99)) )
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#lc=distinguishable_colors(10);
const _pretty_defaults = Dict(
    :size         =>(700,700),   # size of figure
    :format       =>"png",       # format for output plots
    :plant        =>false,       # bool for simulation of plant as well as MPC
    :X0p          =>false,       # bool for simulation of X0p
    :X0p_marker  =>(:circle,:red,5.0,1.0),
    :plantOnly    =>false,       # bool for plotting only the plant response
    :polyPts      =>true,        # bool for plotting the individual polynomials
    :simulate     =>false,       # bool for simulation
    :L            =>100,         # format for output plots
    :mpc_lines    =>[(2.7,:blue,:solid),(2.7,:red,:solid),(2.7,:darkgreen,:solid),
                     (2.7,:yellow,:solid),(2.7,:black,:solid),(2.7,:firebrick,:solid),
                     (2.7,:gold,:solid),(2.7,:greenyellow,:solid),(2.7,:hotpink,:solid),
                     (2.7,:steelblue,:solid),(2.7,:aqua,:solid),(2.7,:teal,:solid),
                     (2.7,:blue,:dash),(2.7,:red,:dash),(2.7,:darkgreen,:dash),
                     (2.7,:yellow,:dash),(2.7,:black,:dash),(2.7,:firebrick,:dash),
                     (2.7,:gold,:dash),(2.7,:greenyellow,:dash),(2.7,:hotpink,:dash),
                     (2.7,:steelblue,:dash),(2.7,:aqua,:dash),(2.7,:teal,:dash)],
    :mpc_markers  =>(:circle,:blueviolet,0.0,1.0),
    :plant_lines  =>[(4.3,:darkgreen,:solid),(4.3,:darkgreen,:dash),(4.3,:black,:solid),(4.3,:black,:dash)],
    :limit_lines  =>[(2.0,:deeppink1,:solid),(2.0,:deeppink1,:dash),(2.0,:blueviolet,:solid),(2.0,:blueviolet,:dash)],
    :opt_marker  => [(:circle,:green,3.0),(:rect,:red,3.0),(:circle,:black,3.0)],

    # for VehicleModel.jl
    :vehicle_marker  =>(:black,:rect,5.0,1.0),
    :vehicle_fill    =>(0,1,:black),
    :goal_marker     =>[(:circle,:springgreen,10.0,1.0),(:circle,:green,10.0,1.0)],
    :goal_fill       =>[(0,1,:springgreen),(0,1,:green)],
    :goal_line       =>(1.0,0.0,:solid,:green),
    :obstacle_marker =>(:circle,:red,10.0,1.0),
    :obstacle_fill   =>(0,1.0,:red),
    :obstacle_line   =>(1.0,0.0,:solid,:red),
    :tire_force_lims =>(0,12000),
    :ax_lims         =>(-5,2),
    :vehicle_width   =>1.9,
    :vehicle_length  =>3.3,
    :lidar_line      =>(1.0,0.0,:solid,:red),
    :lidar_fill      =>(0,0.2,:yellow),
)

"""
plotSettings(;(:simulate=>true))
plotSettings(;(:mpc_lines   =>[(4.0,:blueviolet,:solid),(4.0,:blueviolet,:dash)]),
              (:plant_lines =>[(3.0,:darkgreen,:solid),(3.0,:darkgreen,:solid)]),
              (:limit_lines =>[(2.0,:turquoise,:solid),(2.0,:violet,:solid),(2.0,:orchid,:solid),(2.0,:darksalmon,:solid)]),
              (:size=>(1000,1000))
             )
plotSettings(;(:obstacle_marker =>(:circle,:red,10.0,1.0)),
                  (:goal_marker=>(:circle,:green,10.0,1.0))
             )
 plotSettings(;(:mpc_lines =>[(4.0,:blueviolet,:solid)]),(:size=>(700,700)))
# this function can be called to modify the default settings in PrettyPlots.jl
# check out Colors.jl for an incredible amount of colors!
"""
function plotSettings(;kwargs...)
  kw = Dict(kwargs);
  for (key,value) in kw
    if haskey(_pretty_defaults,key)
      _pretty_defaults[key]=value
    else
      error("Unknown key: ", kw)
    end
  end
end
"""
currentSettings()
# show the current plot settings
"""
function currentSettings()
  _pretty_defaults
end


text:
with_repo_access("test lookup fails with wrong git obj") do test_repo, path
    # commit
    @test_throws LibGitError{:Invalid,:NotFound} lookup_blob(test_repo,
                                                    Oid("8496071c1b46c854b31185ea97743be6a8774479"))
    # tag
    @test_throws LibGitError{:Invalid,:NotFound} lookup_blob(test_repo,
                                                    Oid("0c37a5391bbff43c37f0d0371823a5509eed5b1d"))
    # tree
    @test_throws LibGitError{:Invalid,:NotFound} lookup_blob(test_repo,
                                                    Oid("c4dc1555e4d4fa0e0c9c3fc46734c7c35b3ce90b"))
end

with_repo_access() do test_repo, path

    context("test read blob data")  do
        id = Oid("fa49b077972391ad58037050f2a75f74e3671e92")
        b = lookup(test_repo, id)
        @test sizeof(b) == 9
        @test bytestring(b) == "new file\n"
        @test isa(b, GitBlob)
        @test Oid(b) == id
        @test text(b) == "new file\n"
    end

    context("test blob sloc") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        @test sloc(b) == 328
    end

    context("test blob content with size") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        c =  rawcontent(b, 10)
        @test c == convert(Vector{Uint8}, "# Rugged\n*")
        @test sizeof(c) == 10
    end

    context("test blob content with size gt file size") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        c =  rawcontent(b, 1000000)
        @test sizeof(b) == sizeof(c)
    end

    context("test blob content with zero size") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        c =  rawcontent(b, 0)
        @test isempty(c)
    end

    context("test blob content with negative size") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        c =  rawcontent(b, -100)
        @test sizeof(b) == sizeof(c)
    end

    context("test blob text with max lines") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        blob = lookup(test_repo, id)
        @test text(blob, 1) == "# Rugged\n"
    end

    context("test blob text with lines gt file lines") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        txt = text(b, 1000000)
        @test length(split(txt, "\n")) - 1 == 464
    end

    context("test blob text with zero lines") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        txt = text(b, 0)
        @test txt == ""
    end

    context("test blob text with negative lines") do
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
        txt = text(b, -100)
        @test length(split(txt, "\n")) - 1 == 464
    end

    #= TODO: blob encoding
    begin :test_blob_text_default_encoding
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
    end

    begin :test_blob_text_set_encoding
        id = Oid("7771329dfa3002caf8c61a0ceb62a31d09023f37")
        b = lookup(test_repo, id)
    end
    =#
end

with_tmp_repo_access("test write blob data") do test_repo, path
    @test blob_from_buffer(test_repo, "a new blob content") ==
          Oid("1d83f106355e4309a293e42ad2a2c4b8bdbe77ae")
end

with_tmp_repo_access("test write blob from workdir") do test_repo, path
    @test Oid("1385f264afb75a56a5bec74243be9b367ba4ca08") ==
      blob_from_workdir(test_repo, "README")
end

with_tmp_repo_access("test write blob from disk") do test_repo, path
    file_path = joinpath(TESTDIR, joinpath("fixtures", "archive.tar.gz"))
    @test isfile(file_path)
    id = blob_from_disk(test_repo, file_path)
    @test isa(id, Oid)
    b = test_repo[id]
    open(file_path, "r") do io
        c,rc  = readbytes(io), rawcontent(b)
        @test length(c) == length(rc)
        for i in length(c)
            @assert c[i] == rc[i]
        end
    end
end

with_tmp_repo_access("test blob is binary") do test_repo, path
    binary_file_path = joinpath(TESTDIR, joinpath("fixtures", "archive.tar.gz"))
    binary_blob = lookup(test_repo, blob_from_disk(test_repo, binary_file_path))
    @test isbinary(binary_blob) == true

    text_file_path = joinpath(TESTDIR, joinpath("fixtures", "text_file.md"))
    text_blob = lookup(test_repo, blob_from_disk(test_repo, text_file_path))
    @test isbinary(text_blob) == false
end

# test blob diff
sandboxed_test("diff") do test_repo, path
    t1 = GitTree(lookup(test_repo, Oid("d70d245ed97ed2aa596dd1af6536e4bfdb047b69")))
    t2 = GitTree(lookup(test_repo, Oid("7a9e0b02e63179929fed24f0a3e0f19168114d10")))

    b = lookup(test_repo, Oid(t1["readme.txt"]))
    o = lookup(test_repo, Oid(t2["readme.txt"]))
    p = diff(test_repo, b, o)

    @test delta(p).status == :modified

    hs = hunks(p)
    @test length(hs) == 3

    @test beginswith(hs[1].header, "@@ -1,4 +1,4 @@")
    @test beginswith(hs[2].header, "@@ -7,10 +7,6 @@")
    @test beginswith(hs[3].header, "@@ -24,12 +20,9 @@")

    ls = lines(hs[1])
    @test length(ls) == 5
    @test :deletion == ls[1].line_origin
    @test "The Git feature that really makes it stand apart from nearly every other SCM\n" == ls[1].content

    @test :addition == ls[2].line_origin
    @test "The Git feature that r3ally mak3s it stand apart from n3arly 3v3ry other SCM\n" == ls[2].content

    @test :context == ls[3].line_origin
    @test "out there is its branching model.\n" == ls[3].content

    @test :context == ls[4].line_origin
    @test "\n" == ls[4].content

    @test :context == ls[5].line_origin
    @test "Git allows and encourages you to have multiple local branches that can be\n" == ls[5].content
end

sandboxed_test("diff", "test diff nothing") do test_repo, path
    t1 = GitTree(lookup(test_repo, Oid("d70d245ed97ed2aa596dd1af6536e4bfdb047b69")))
    b  = lookup(test_repo, Oid(t1["readme.txt"]))
    p  = diff(test_repo, b, nothing)

    @test delta(p).status == :deleted

    hs = hunks(p)
    @test length(hs) == 1

    @test beginswith(first(hs).header, "@@ -1,35 +0,0 @@")
    ls = lines(first(hs))

    @test length(ls) == 35
    for l in ls
      @test l.line_origin == :deletion
    end
end

sandboxed_test("diff", "test diff with paths") do test_repo, path
    t1 = GitTree(lookup(test_repo, Oid("d70d245ed97ed2aa596dd1af6536e4bfdb047b69")))
    t2 = GitTree(lookup(test_repo, Oid("7a9e0b02e63179929fed24f0a3e0f19168114d10")))

    b = lookup(test_repo, Oid(t1["readme.txt"]))
    o = lookup(test_repo, Oid(t2["readme.txt"]))
    p = diff(test_repo, b, o, {:old_path => "old_readme.txt", :new_path => "new_readme.txt"})

    @test "old_readme.txt" == delta(p).old_file.path
    @test "new_readme.txt" == delta(p).new_file.path
end

with_tmp_repo_access("test write blob from io with hintpath") do test_repo, path
    file_path = joinpath(TESTDIR, joinpath("fixtures", "archive.tar.gz"))
    open(file_path, "r") do io
        id = blob_from_stream(test_repo, io, "archive.tar.gz2")
        seekstart(io)
        blob = lookup(test_repo, id)
        c,rc = readbytes(io), rawcontent(blob)
        @test length(c) == length(rc)
        for i in length(rc)
            @assert c[i] == rc[i]
        end
    end
end

with_tmp_repo_access("test write blob from io without hintpath") do test_repo, path
    file_path = joinpath(TESTDIR, joinpath("fixtures", "archive.tar.gz"))
    open(file_path, "r") do io
        id = blob_from_stream(test_repo, io)
        seekstart(io)
        blob = lookup(test_repo, id)
        c,rc = readbytes(io), rawcontent(blob)
        @test length(c) == length(rc)
        for i in length(rc)
            @assert c[i] == rc[i]
        end
    end
end


text:
using LinearAlgebra
using StatsBase

@noinline function compute_cordf(df, columns)
	result = DataFrame(col=columns)
	function compute_cor(name_1, name_2)
		col_1 = df[!, name_1]
		col_2 = df[!, name_2]
		StatsBase.cor(col_1, col_2)
	end
	transformations = map(columns) do name_1
		:col => ByRow(name_2 -> compute_cor(name_1, name_2)) => Symbol(name_1)
	end
	transform(result, transformations...)
end

function find_multicolinearity(cordf, ycol, threshold::Number)
	to_be_removed = String[]
	sizehint!(to_be_removed, size(cordf, 1))
	considers = setdiff(cordf.col, ["SalePrice"])
	name_index = String.(cordf.col)
	df = select(cordf, Not("col"))
	for n1 in considers, n2 in considers
		if n1 in to_be_removed || n2 in to_be_removed || n1 == n2
			continue
		end
		row_ind = findfirst(n1 .== name_index)
		col_ind = findfirst(n2 .== name_index)
		if abs(df[row_ind, col_ind]) > threshold
			append!(to_be_removed, [abs(df[row_ind, ycol]) > abs(df[col_ind, ycol]) ? n2 : n1])
		end
	end
	to_be_removed
end

function find_multicolinearity(df, columns, ycol; threshold_step=0.01, threshold_max=1.0, threshold_min=0.5)
	t_range = reverse(collect(range(start=threshold_min, stop=threshold_max, step=threshold_step)))
	df = select(df, columns)
	cordf = compute_cordf(df, columns)
	for t in t_range
		removed = find_multicolinearity(cordf, ycol, t)
		df2 = select(df, Not(removed))
		if rank(Matrix(df2)) == length(names(df2))
			return removed, t
		end
	end
	String[], threshold_min
end


text:
module ParallelPlanners

using Base: @kwdef
using Parameters: @unpack
using Setfield: @set
using DataStructures: PriorityQueue, OrderedDict, enqueue!, dequeue!
using StatsBase: sample, Weights
using Random, Julog, PDDL, CUDA

include("specifications.jl")
include("solutions.jl")
include("simulators.jl")
include("heuristics.jl")
include("planners.jl")

end # module


text:
## Power iteration-based centrality indices ##

"""
    eigenvector(g::Graph[, eps::Real, maxiter::Int])

Compute eigenvector centrality for all nodes in the graph.

The process does not always converge for directed graphs, in which case
the function returns an error. Oftentimes, it is useful to treat the graph
as undirected, `eigenvector(undirected(g))`.

# Arguments
* `g::Graph`: the graph itself.
* `eps::Real=0.0001`: convergence threshold (when score changes for each node
stop exceeding it, the algorithm terminates).
* `maxiter::Int=100`: maximum number of iterations.
"""
function eigenvector(g::Graph, eps::Real=1e-4, maxiter::Int=100)
    @assert(0 < eps <= 0.1, "convergence threshold (eps) should be positive and small")
    n         = nodecount(g)
    score     = fill(0.,     n)
    prevscore = fill(1. / n, n)
    maxdiff   = 1.
    iter      = 0
    @inbounds while maxdiff > eps
        for i = 1:n             # transfer the scores
            score[i] = 0.
            for j = inneighbors(g, i)
                score[i] += prevscore[j]
            end
        end
        λ = 0                   # normalize the scores
        for i = 1:n
            if abs(score[i]) > λ
                λ = score[i]
            end
        end
        if λ == 0
            error("power iteration cannot converge")
            break
        else
            maxdiff = 0.
            for i = 1:n
                score[i] /= λ
                maxdiff = max(maxdiff, abs(score[i] - prevscore[i]))
            end
        end
        iter += 1
        if iter > maxiter && maxdiff > eps
            warn("power iteration did not converge in $maxiter steps, current convergence value is $maxdiff")
            return score
        end
        prevscore, score = score, prevscore
    end
    return prevscore
end

"""
    pagerank(g::Graph[, d::Real, eps::Real])

Compute PageRank scores for all nodes in the graph.

# Arguments
* `g::Graph`: the graph itself.
* `d::Real=0.85`: damping factor (the lesser it is, the bigger proportion
    of scores is redistributed equally between all nodes at each iteration).
* `eps::Real=0.0001`: convergence threshold (when score changes for each node
    stop exceeding it, the algorithm terminates).
"""
function pagerank(g::Graph, d::Real=0.85, eps::Real=1e-4)
    @assert(0 < d < 1,     "damping factor (d) should be between 0 and 1")
    @assert(0 < eps < 0.1, "convergence threshold (eps) should be positive and small")
    n         = nodecount(g)
    score     = fill(0.,     n)
    prevscore = fill(1. / n, n)
    ratio     = [d / outdegree(g, i) for i = 1:n]
    maxdiff   = 1.
    @inbounds while maxdiff > eps
        for i = 1:n                    # transfer the scores
            score[i] = 0.
            for j = inneighbors(g, i)
                score[i] += prevscore[j] * ratio[j]
            end
        end
        leak = (1. - sum(score)) / n    # put the leaked scores back
        maxdiff = 0.
        for i = 1:n
            score[i] += leak
            maxdiff = max(maxdiff, abs(score[i] - prevscore[i]))
        end
        prevscore, score = score, prevscore
    end
    return prevscore
end


text:
using Documenter, Flux, NNlib

makedocs(modules=[Flux, NNlib],
         doctest = false,
         format = :html,
         analytics = "UA-36890222-9",
         sitename = "Flux",
         assets = ["../flux.css"],
         pages = ["Home" => "index.md",
                  "Building Models" =>
                    ["Basics" => "models/basics.md",
                     "Recurrence" => "models/recurrence.md",
                     "Model Reference" => "models/layers.md"],
                  "Training Models" =>
                    ["Optimisers" => "training/optimisers.md",
                     "Training" => "training/training.md"],
                  "One-Hot Encoding" => "data/onehot.md",
                  "GPU Support" => "gpu.md",
                  "Community" => "community.md"])

deploydocs(
   repo = "github.com/FluxML/Flux.jl.git",
   target = "build",
   osname = "linux",
   julia = "0.6",
   deps = nothing,
   make = nothing)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1SecretReference <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace

    function IoK8sApiCoreV1SecretReference(;name=nothing, namespace=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecretReference, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1SecretReference, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        o
    end
end # type IoK8sApiCoreV1SecretReference

const _property_map_IoK8sApiCoreV1SecretReference = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"))
const _property_types_IoK8sApiCoreV1SecretReference = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("namespace")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1SecretReference }) = collect(keys(_property_map_IoK8sApiCoreV1SecretReference))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecretReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecretReference[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecretReference }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecretReference[property_name]

function check_required(o::IoK8sApiCoreV1SecretReference)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1SecretReference }, name::Symbol, val)
end


text:
using Documenter, DailyMed

makedocs(sitename = "DailyMed Module Documentation", format = Documenter.HTML(prettyurls = false))