text:
#=
HNF algorithm due to Arne Storjohann
"A Fast + Practical + Deterministic Algorithm for Triangularizing Integer Matrices"
=#

# This function only changes CC, zero_cols and N but not T.
function conditioning_with_trafo!(T::MatElem{S}, N::MatElem{S}, CC::Array{S}, zero_cols::BitArray, row::Int, col1::Int) where S <: Nemo.RingElement
  @assert col1 < cols(T)
  @assert T[row, col1] > 0
  R = base_ring(T)
  row1 = row + 1
  # We don't need every field of CC
  for i = (row + 2):rows(T)
    CC[i] = zero!(CC[i])
    zero_cols[i] = true
  end
  N[1, 1] = deepcopy(T[row, col1])
  N[2, 1] = deepcopy(T[row1, col1])
  t = R()
  t1 = R()
  t2 = R()
  col2 = 0
  s = 0
  for c = (col1 + 1):cols(T)
    for i = row1:rows(T)
      t1 = mul!(t1, N[1, 1], T[i, c])
      t2 = mul!(t2, T[row, c], T[i, col1])
      t = t1 - t2
      if !iszero(t)
        col2 = c
        N[1, 2] = deepcopy(T[row, col2])
        N[2, 2] = deepcopy(T[row1, col2])
        # Make sure the rows row and row1 are linearly independent
        if i != row1
          N[2, 1] = addeq!(N[2, 1], T[i, col1])
          N[2, 2] = addeq!(N[2, 2], T[i, col2])
          s = i
        end
        break
      end
    end
    if col2 != 0
      break
    end
  end
  @assert col2 != 0
  if row1 == rows(T)
    return col2
  end
  for l = row + 2:rows(T)
    g = gcd(N[2, 1], T[l, col1])
    if iszero(g)
      continue
    end
    q1, r1 = divrem(divexact(N[2, 1], g), N[1, 1])
    if r1 < 0
      addeq!(r1, N[1, 1])
      q1 = q1 - 1
    end
    q2, r2 = divrem(divexact(T[l, col1], g), N[1, 1])
    if r2 < 0
      addeq!(r2, N[1, 1])
      q2 = q2 - 1
    end
    t1 = mul!(t1, N[1, 1], T[l, col2])
    t2 = mul!(t2, N[1, 2], T[l, col1])
    n = t1 - t2
    if n != 0
      t1 = mul!(t1, N[1, 1], N[2, 2])
      t2 = mul!(t2, N[1, 2], N[2, 1])
      nn = t2 - t1
      @assert nn != 0
      # If nn/n in Z and positive:
      if sign(n) == sign(nn) && rem(nn, n) == 0
        r2 = -r2
      end
    end
    t1 = zero!(t1)
    t2 = deepcopy(r1)
    while gcd(t2, N[1, 1]) != 1
      t2 = addeq!(t2, r2)
      t1 += 1
    end
    if !iszero(t1)
      if r2 < 0
        t1 = -t1
      end
      CC[l] = deepcopy(t1)
      zero_cols[l] = false
      t = mul!(t, T[l, col1], t1)
      N[2, 1] = addeq!(N[2, 1], t)
      t = mul!(t, T[l, col2], t1)
      N[2, 2] = addeq!(N[2, 2], t)
    end
  end
  if s > row1
    CC[s] = CC[s] + R(1)
    zero_cols[s] = false
  end
  return col2
end

# This function only changes N and QC, but not T.
function column_reduction_with_trafo!(T::MatElem{S}, N::MatElem{S}, QC::MatElem{S}, row::Int, col1::Int, col2::Int) where S <: Nemo.RingElement
  @assert T[ row, col1] > 0
  @assert rows(T) == rows(QC)
  R = base_ring(T)
  row1 = row + 1
  t = R()
  t1 = R()
  t2 = R()
  g, a, b = gcdx(N[1, 1], N[2, 1])
  QC[row, row] = a
  QC[row, row1] = b
  QC[row1, row] = -divexact(N[2, 1], g)
  QC[row1, row1] = divexact(N[1, 1], g)
  t = deepcopy(N[1, 2])
  t1 = mul!(t1, N[1, 2], a)
  t2 = mul!(t2, N[2, 2], b)
  N[1, 2] = add!(N[1, 2], t1, t2)
  t1 = mul!(t1, QC[row1, row], t)
  t2 = mul!(t2, QC[row1, row1], N[2, 2])
  N[2, 2] = add!(N[2, 2], t1, t2)
  if N[2, 2] < 0
    QC[row1, row] = -QC[row1, row]
    QC[row1, row1] = -QC[row1, row1]
    N[2, 2] = -N[2, 2]
  end
  N[1, 1] = g
  # Make 0 below T[row, col1] and reduce below T[row1, col2]
  for i = (row + 2):rows(T)
    if iszero(T[i, col1])
      QC[i, row] = zero!(QC[i, row])
      QC[i, row1] = zero!(QC[i, row1])
    else
      g, a, b = gcdx(N[1, 1], T[i, col1])
      g1 = divexact(N[1, 1], g)
      g2 = -divexact(T[i, col1], g)
      N[1, 1] = g
      t = deepcopy(N[1, 2])
      t1 = mul!(t1, a, N[1, 2])
      t2 = mul!(t2, b, T[i, col2])
      N[1, 2] = add!(N[1, 2], t1, t2)
      t1 = mul!(t1, g1, T[i, col2])
      t2 = mul!(t2, g2, t)
      t = add!(t, t1, t2) # this would be T[i, col2], but we don't change T here
      for c = row:row1
        QC[i, c] = mul!(QC[i, c], g2, QC[row, c])
        QC[row, c] = mul!(QC[row, c], a, QC[row, c])
      end
    end
    if t < 0 || t >= N[2, 2]
      q = -fdiv(t, N[2, 2])
      for c = row:row1
        t1 = mul!(t1, q, QC[row1, c])
        QC[i, c] = addeq!(QC[i, c], t1)
      end
    end
  end
  # Reduce above T[row, col1] and above T[row1, col2]
  for i = 1:row
    t = deepcopy(T[i, col2])
    if i != row && (T[i, col1] < 0 || T[i, col1] >= N[1, 1])
      q = -fdiv(T[i, col1], N[1, 1])
      t1 = mul!(t1, q, N[1, 2])
      t = addeq!(t, t1) # this would be T[i, col2]
      for c = row:row1
        QC[i, c] = mul!(QC[i, c], q, QC[row, c])
      end
    else
      if i != row
        QC[i, row] = zero!(QC[i, row])
        QC[i, row1] = zero!(QC[i, row1])
      end
    end
    if t < 0 || t >= N[2, 2]
      q = -fdiv(t, N[2, 2])
      for c = row:row1
        t1 = mul!(t1, q, QC[row1, c])
        QC[i, c] = addeq!(QC[i, c], t1)
      end
    end
  end
  return nothing
end

function hnf_storjohann_with_trafo(A::MatElem{S}) where S <: Nemo.RingElement
  #timeRed = 0
  #timeProdC = 0
  #timeBuildQC = 0
  #timeProdQC = 0
  #timeU = 0
  R = base_ring(A)
  T = similar(A, rows(A) + 2, cols(A) + 2)
  n = rows(T)
  m = cols(T)
  T[1, 1] = R(1)
  T[n, m] = R(1)
  n1 = n - 1
  m1 = m - 1
  for i = 2:n1
    for j = 2:m1
      T[i, j] = deepcopy(A[i - 1, j - 1])
    end
  end
  #=
      1 0 0
  T = 0 A 0
      0 0 1
  =#
  Q = eye(T, n)
  C = eye(T, n)
  QC = eye(T, n)
  #=
  We use QC at the same time for Q_k and the product Q_k * CC.
    Q_k shoud be
  1 0 ... * * ... 0
  0 1 ... * * ... 0
  . . . . * * . . .
  0 0 ... * * ... 1
  Only the columns row and row + 1 (with the *'s) are saved in QC.
  The columns before row will not be reset since they are never needed again.
  =#
  N = similar(T, 2, 2)
  col1 = 1
  col2 = 0
  row = 1
  CC = Array{S, 1}(n)
  for i = 1:n
    CC[i] = R()
  end
  #=
  CC should be
  1 0 ..... 0
  0 1 ..... 0
  . . . . . .
  0 0 1 * * *
  . . . . . .
  0 0 0 0 0 1
  Only the "row + 1 row" (the one with the *'s) is saved here.
  =#
  zero_cols = BitArray(n) # zero_cols[i] == true iff CC[i] == 0
  non_zero_entries = Array{Array{Int, 1}, 1}(n)
  # If i is in non_zero_entries[j] then C[i, j] != 0 (but we don't
  # save j in non_zero_entries[j] although there is always a 1 on
  # the diagonal).
  for i = 1:n
    non_zero_entries[i] = Array{Int, 1}()
  end
  temp = Array{S, 1}(n)
  for i = 1:n
    temp[i] = R()
  end
  t = R()
  while col2 < m
    row1 = row + 1
    row2 = row + 2
    #tic()
    col2 = conditioning_with_trafo!(T, N, CC, zero_cols, row, col1)
    column_reduction_with_trafo!(T, N, QC, row, col1, col2)
    col1 = col2
    #timeRed += toq();
    # C = CC*C
    #tic()
    for i = row2:n
      if zero_cols[i]
        continue
      end
      for j = row2:i # C is upper triangular
        if iszero(C[i, j])
          continue
        end
        t = mul!(t, CC[i], C[i, j])
        C[row1, j] = addeq!(C[row1, j], t)
        push!(non_zero_entries[j], row1)
      end
    end
    #timeProdC += toq();
    # "QC = Q_k*CC"
    #tic()
    for j = row2:n
      if zero_cols[j]
        continue
      end
      for i = 1:n
        QC[i, j] = mul!(QC[i, j], QC[i, row1], CC[j])
        if i == j
          QC[i, j] += 1
        end
      end
    end
    #timeBuildQC += toq();
    zero_cols[row] = false
    zero_cols[row1] = false
    # Q = QC*Q
    # We only need to change the columns 1 to row1 since we already know,
    # that the columns row2 to n will eventually be an "indentity matrix".
    # In particular, we don't need to compute Q*inv(CC).
    #tic()
    for j = 1:row1
      for i = row:n
        temp[i] = zero!(temp[i])
      end
      for l = row:n
        if zero_cols[l]
          temp[l] = addeq!(temp[l], Q[l, j])
          continue
        end
        for i = 1:row - 1
          t = mul!(t, QC[i, l], Q[l, j])
          Q[i, j] = addeq!(Q[i, j], t)
        end
        for i = row:n
          t = mul!(t, QC[i, l], Q[l, j])
          temp[i] = addeq!(temp[i], t)
        end
      end
      for i = row:n
        Q[i, j] = deepcopy(temp[i])
      end
    end
    # T = QC*T
    for j = row:m
      for i = row:n
        temp[i] = zero!(temp[i])
      end
      for l = row:n
        if zero_cols[l]
          temp[l] = addeq!(temp[l], T[l, j])
          continue
        end
        for i = 1:row - 1
          t = mul!(t, QC[i, l], T[l, j])
          T[i, j] = addeq!(T[i, j], t)
        end
        for i = row:n
          t = mul!(t, QC[i, l], T[l, j])
          temp[i] = addeq!(temp[i], t)
        end
      end
      for i = row:n
        T[i, j] = deepcopy(temp[i])
      end
    end
    #timeProdQC += toq();
    row += 1
  end
  #tic()
  U = deepcopy(Q)
  for j = 1:n
    if length(non_zero_entries[j]) == 0
      continue
    end
    for i = 1:n
      for l in non_zero_entries[j]
        t = mul!(t, Q[i, l], C[l, j])
        U[i, j] = addeq!(U[i, j], t)
      end
    end
  end
  #timeU += toq();
  #println("Time reduction: $timeRed")
  #println("Time products with CC: $timeProdC")
  #println("Time to compute QC: $timeBuildQC")
  #println("Time products with QC: $timeProdQC")
  #println("Time to compute U: $timeU")
  return sub(T, 2:n1, 2:m1), sub(Q, 2:n1, 2:n1), sub(C, 2:n1, 2:n1), sub(U, 2:n1, 2:n1)
end

function conditioning!(T::MatElem{S}, row::Int, col1::Int) where S <: Nemo.RingElement
  @assert col1 < cols(T)
  @assert T[row, col1] > 0
  R = base_ring(T)
  row1 = row + 1
  t = R()
  t1 = R()
  t2 = R()
  col2 = 0
  for c = (col1 + 1):cols(T)
    for i = row1:rows(T)
      t1 = mul!(t1, T[row, col1], T[i, c])
      t2 = mul!(t2, T[row, c], T[i, col1])
      t = t1 - t2
      if !iszero(t)
        col2 = c
        # Make sure the rows row and row1 are linearly independent
        if i != row1
          for j = col1:cols(T)
            T[row1, j] = addeq!(T[row1, j], T[i, j])
          end
        end
        break
      end
    end
    if col2 != 0
      break
    end
  end
  @assert col2 != 0
  if row1 == rows(T)
    return col2
  end
  for l = row + 2:rows(T)
    g = gcd(T[row1, col1], T[l, col1])
    if iszero(g)
      continue
    end
    q1, r1 = divrem(divexact(T[row1, col1], g), T[row, col1])
    if r1 < 0
      addeq!(r1, T[row, col1])
      q1 = q1 - 1
    end
    q2, r2 = divrem(divexact(T[l, col1], g), T[row, col1])
    if r2 < 0
      addeq!(r2, T[row, col1])
      q2 = q2 - 1
    end
    t1 = mul!(t1, T[row, col1], T[l, col2])
    t2 = mul!(t2, T[row, col2], T[l, col1])
    n = t1 - t2
    if n != 0
      t1 = mul!(t1, T[row, col1], T[row1, col2])
      t2 = mul!(t2, T[row, col2], T[row1, col1])
      nn = t2 - t1
      @assert nn != 0
      # If nn/n in Z and positive:
      if sign(n) == sign(nn) && rem(nn, n) == 0
        r2 = -r2
      end
    end
    t1 = zero!(t1)
    t2 = deepcopy(r1)
    while gcd(t2, T[row, col1]) != 1
      t2 = addeq!(t2, r2)
      t1 += 1
    end
    if !iszero(t1)
      if r2 < 0
        t1 = -t1
      end
      for j = col1:(cols(T) - 1)
        t = mul!(t, T[l, j], t1)
        T[row1, j] = addeq!(T[row1, j], t)
      end
    end
  end
  return col2
end

function column_reduction!(T::MatElem{S}, row::Int, col1::Int, col2::Int) where S <: Nemo.RingElement
  @assert T[ row, col1] > 0
  R = base_ring(T)
  row1 = row + 1
  col11 = col1 + 1
  n = cols(T) - 1
  t = R()
  t1 = R()
  t2 = R()
  # Make 0 below T[row, col1]
  for i = row1:rows(T)
    if i == row1 || !iszero(T[i, col1])
      g, a, b = gcdx(T[row, col1], T[i, col1])
      c = divexact(T[row, col1], g)
      d = -divexact(T[i, col1], g)
      if i == row1
        t1 = mul!(t1, d, T[row, col2])
        t2 = mul!(t2, c, T[row1, col2])
        t = add!(t, t1, t2)
        if t < 0
          c = -c
          d = -d
        end
      end
      T[row, col1] = deepcopy(g)
      T[i, col1] = zero!(T[i, col1])
      for j = col11:n
        t = deepcopy(T[row, j])
        t1 = mul!(t1, T[row, j], a)
        t2 = mul!(t2, T[i, j], b)
        T[row, j] = add!(T[row, j], t1, t2)
        t1 = mul!(t1, T[i, j], c)
        t2 = mul!(t2, t, d)
        T[i, j] = add!(T[i, j], t1, t2)
      end
    end
  end
  # Reduce below T[row1, col2]
  for i = (row + 2):rows(T)
    if T[i, col2] < 0 || T[i, col2] >= T[row1, col2]
      q = -fdiv(T[i, col2], T[row1, col2])
      for j = col11:n
        t = mul!(t, q, T[row1, j])
        T[i, j] = addeq!(T[i, j], t)
      end
    end
  end
  # Reduce above T[row, col1]
  for i = 1:(row - 1)
    if T[i, col1] < 0 || T[i, col1] >= T[row, col1]
      q = -fdiv(T[i, col1], T[row, col1])
      for j = col1:n
        t = mul!(t, q, T[row, j])
        T[i, j] = addeq!(T[i, j], t)
      end
    end
  end
  # Reduce above T[row1, col2]
  for i = 1:row
    if T[i, col2] < 0 || T[i, col2] >= T[row1, col2]
      q = -fdiv(T[i, col2], T[row1, col2])
      for j = col11:n
        t = mul!(t, q, T[row1, j])
        T[i, j] = addeq!(T[i, j], t)
      end
    end
  end
  return nothing
end

function hnf_storjohann(A::MatElem{S}) where S <: Nemo.RingElement
  R = base_ring(A)
  T = similar(A, rows(A) + 2, cols(A) + 2)
  n = rows(T)
  m = cols(T)
  T[1, 1] = R(1)
  T[n, m] = R(1)
  n1 = n - 1
  m1 = m - 1
  for i = 2:n1
    for j = 2:m1
      T[i, j] = deepcopy(A[i - 1, j - 1])
    end
  end
  #=
      1 0 0
  T = 0 A 0
      0 0 1
  =#
  col1 = 1
  col2 = 0
  row = 1
  while col2 < m
    col2 = conditioning!(T, row, col1)
    column_reduction!(T, row, col1, col2)
    col1 = col2
    row += 1
  end
  return sub(T, 2:n1, 2:m1)
end


text:
@testset "datenum.jl" begin
    using Dates
    # see https://gist.github.com/RobBlackwell/913b129d8d91e72e9879b20cb38fd529
    MATLAB_EPOCH = Dates.DateTime(-0001,12,31)

    """
         datenum(d::Dates.DateTime)
    Converts a Julia DateTime to a MATLAB style DateNumber.
    MATLAB represents time as DateNumber, a double precision floating
    point number being the the number of days since January 0, 0000
    Example
        datenum(now())
    """
    function datenum_RobBlackwell(d::Dates.DateTime)
        Dates.value(d - MATLAB_EPOCH) /(1000 * 60 * 60 * 24)
    end
    dt = DateTime(2006,11,14,9,15,37);
    @test isequal(datenum(dt), datenum_RobBlackwell(dt))
    dt = DateTime(2033,5,1,23,59,59);
    @test isequal(datenum(dt), datenum_RobBlackwell(dt))
    dt = DateTime(1,2,28,0,13,59);
    @test isequal(datenum(dt), datenum_RobBlackwell(dt))
    dt = DateTime(1999,12,31,23,59,59);
    @test isequal(datenum(dt), datenum_RobBlackwell(dt))
    @test isequal(toordinal(1970,1,1), 719163);
end

ext: jl


text:
"""
    simulate!(rng::AbstractRNG, m::MixedModel{T}; β=m.β, σ=m.σ, θ=T[])
    simulate!(m::MixedModel; β=m.β, σ=m.σ, θ=m.θ)

Overwrite the response (i.e. `m.trms[end]`) with a simulated response vector from model `m`.
"""
function simulate!(
    rng::AbstractRNG,
    m::LinearMixedModel{T};
    β = coef(m),
    σ = m.σ,
    θ = T[],
) where {T}
    length(β) == length(fixef(m)) ||
        length(β) == length(coef(m)) ||
            throw(ArgumentError("You must specify all (non-singular) βs"))

    β = convert(Vector{T},β)
    σ = T(σ)
    θ = convert(Vector{T},θ)
    isempty(θ) || setθ!(m, θ)

    if length(β) ≠ length(coef(m))
        padding = length(coef(m)) - length(β)
        for ii in 1:padding
            push!(β, -0.0)
        end
    end

    y = randn!(rng, response(m))      # initialize y to standard normal

    for trm in m.reterms              # add the unscaled random effects
        unscaledre!(rng, y, trm)
    end
                    # scale by σ and add fixed-effects contribution
    mul!(y, m.X, β, one(T), σ)

    m
end

function simulate!(
    rng::AbstractRNG,
    m::GeneralizedLinearMixedModel{T};
    β = coef(m),
    σ = m.σ,
    θ = T[],
) where {T}
    length(β) == length(fixef(m)) ||
        length(β) == length(coef(m)) ||
            throw(ArgumentError("You must specify all (non-singular) βs"))

    dispersion_parameter(m) || ismissing(σ) ||
        throw(ArgumentError("You must not specify a dispersion parameter for model families without a dispersion parameter"))

    β = convert(Vector{T},β)
    if σ !== missing
        σ = T(σ)
    end
    θ = convert(Vector{T},θ)

    d = m.resp.d

    if length(β) ≠ length(coef(m))
        padding = length(coef(m)) - length(β)
        for ii in 1:padding
            push!(β, -0.0)
        end
    end

    fast = (length(m.θ) == length(m.optsum.final))
    setpar! = fast ? setθ! : setβθ!
    params = fast ? θ : vcat(β, θ)
    setpar!(m, params)

    lm = m.LMM

    # note that these m.resp.y and m.LMM.y will later be sychronized in (re)fit!()
    # but for now we use them as distinct scratch buffers to avoid allocations

    # the noise term is actually in the GLM and not the LMM part so no noise
    # at the LMM level
    η = fill!(m.LMM.y, zero(T))
    y = m.resp.y

    # assemble the linear predictor

    # add the unscaled random effects
    # note that unit scaling may not be correct for
    # families with a dispersion parameter
    @inbounds for trm in m.reterms
        unscaledre!(rng, η, trm)
    end

    # add fixed-effects contribution
    # note that unit scaling may not be correct for
    # families with a dispersion parameter
    mul!(η, lm.X, β, one(T), one(T))

    # from η to μ
    GLM.updateμ!(m.resp, η)

    # convert to the distribution / add in noise
    @inbounds for (idx, val) in enumerate(m.resp.mu)
        n = isempty(m.wt) ? 1 : m.wt[idx]
        y[idx] = _rand(rng, d, val, σ, n)
    end

    m
end

"""
    _rand(rng::AbstractRNG, d::Distribution, location, scale=missing, n=1)

A convenience function taking a draw from a distrbution.

Note that `d` is specified as an existing distribution, such as
from the `GlmResp.d` field. This isn't vectorized nicely because
for distributions where the scale/dispersion is dependent on the
location (e.g. Bernoulli, Binomial, Poisson), it's not really
possible to avoid creating multiple `Distribution` objects.

Note that `n` is the `n` parameter for the Binomial distribution,
*not* the number of draws from the RNG. It is then used to change the
random draw (an integer in [0, n]) into a probability (a float in [0,1]).
"""
function _rand(rng::AbstractRNG, d::Distribution, location, scale=NaN, n=1)
    if !ismissing(scale)
        throw(ArgumentError("Families with a dispersion parameter not yet supported"))
    end

    if d isa Binomial
        dist = Binomial(Int(n), location)
    else
        dist = typeof(d)(location)
    end

    rand(rng, dist) / n
end

function simulate!(m::MixedModel{T}; β = coef(m), σ = m.σ, θ = T[]) where {T}
    simulate!(Random.GLOBAL_RNG, m, β = β, σ = σ, θ = θ)
end

"""
    unscaledre!(y::AbstractVector{T}, M::ReMat{T}, b) where {T}
    unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, M::ReMat{T}) where {T}

Add unscaled random effects defined by `M` and `b` to `y`.  When `rng` is present the `b`
vector is generated as `randn(rng, size(M, 2))`
"""
function unscaledre! end

function unscaledre!(y::AbstractVector{T}, A::ReMat{T,1}, b::AbstractVector{T}) where {T}
    m, n = size(A)
    length(y) == m && length(b) == n || throw(DimensionMismatch(""))
    z = A.z
    @inbounds for (i, r) in enumerate(A.refs)
        y[i] += b[r] * z[i]
    end
    y
end

unscaledre!(y::AbstractVector{T}, A::ReMat{T,1}, B::AbstractMatrix{T}) where {T} =
    unscaledre!(y, A, vec(B))

function unscaledre!(y::AbstractVector{T}, A::ReMat{T,S}, b::AbstractMatrix{T}) where {T,S}
    Z = A.z
    k, n = size(Z)
    l = nlevs(A)
    length(y) == n && size(b) == (k, l) || throw(DimensionMismatch(""))
    @inbounds for (i, ii) in enumerate(A.refs)
        for j = 1:k
            y[i] += Z[j, i] * b[j, ii]
        end
    end
    y
end

function unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,S}) where {T,S}
    rng_nums = randn(rng, S, nlevs(A))

    unscaledre!(y, A, lmul!(A.λ, rng_nums))
end

function unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,1}) where {T}
    rng_nums = randn(rng, 1, nlevs(A))

    unscaledre!(y, A, lmul!(first(A.λ), rng_nums))
end

unscaledre!(y::AbstractVector, A::ReMat) = unscaledre!(Random.GLOBAL_RNG, y, A)


text:
module TestAsInterval

using Base.Test
using ArbFloats

afIvl1 = midpoint_radius(ArbFloat(1.0), ArbFloat(1.0e-8))
afIvl2 = midpoint_radius(ArbFloat(1.0), ArbFloat(1.0e-10))
afIvl3 = nextfloat(afIvl2)

@test contains(afIvl1, afIvl2)
@test !contains(afIvl2, afIvl1)
@test overlap(afIvl2, afIvl3)

@test midpoint(afIvl3) > midpoint(afIvl2)
@test radius(afIvl3) == radius(afIvl2)

end # module


text:
using PkgTemplates

t = Template(;
            user="Shoichiro-Tsutsui", # GitHubのユーザー名
            dir=".",
            plugins=[
                License(; name="MIT"), # MITライセンスであることを明示する。
                TravisCI(), # CIサービスとしてTravisCIを使う
                Codecov(),
                Documenter{TravisCI}(), # ドキュメントのテンプレを作成する。
            ],
            )

t("MyPkg")

max_stars_repo_name: Shoichiro-Tsutsui/LearnTravisCI

max_stars_repo_licenses: FTL RSA-MD

max_issues_repo_name: Shoichiro-Tsutsui/LearnTravisCI

max_issues_repo_licenses: FTL RSA-MD

max_forks_repo_name: Shoichiro-Tsutsui/LearnTravisCI

max_forks_repo_licenses: FTL RSA-MD


text:
module JMPS
    import Base.getindex,Base.setindex!,Base.*,Base.transpose
    using LinearAlgebra,TensorOperations,OMEinsum,CUDA

    include("./abstract_types.jl")
    include("./Bdim.jl")
    include("./MPS.jl")
    include("./normalization.jl")
    include("./MPO.jl")
    include("./measure.jl")
    include("./subroutine.jl")
    include("./mvcuda.jl")
    include("./randMPS.jl")
    # debug_tools
    include("./utils.jl")
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text: include("brentmin.jl")

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
function AdaptiveStepsizeFieldIntegrator(arg0::Field, arg1::JString, arg2::jdouble, arg3::jdouble, arg4::Vector{jdouble}, arg5::Vector{jdouble})
    return AdaptiveStepsizeFieldIntegrator((Field, JString, jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5)
end

function AdaptiveStepsizeFieldIntegrator(arg0::Field, arg1::JString, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return AdaptiveStepsizeFieldIntegrator((Field, JString, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function integrate(obj::FieldODEIntegrator, arg0::FieldExpandableODE, arg1::FieldODEState, arg2::RealFieldElement)
    return jcall(obj, "integrate", FieldODEStateAndDerivative, (FieldExpandableODE, FieldODEState, RealFieldElement), arg0, arg1, arg2)
end


text:
@testset "server + client" begin
    script = joinpath(@__DIR__, "..", "src", "servers", "standalone.jl")
    cmd = `$(joinpath(JULIA_HOME, "julia")) -- $script --zmq-url="tcp://127.0.0.1:5568"`
    if !haskey(ENV, "CI")
        cmd = `$cmd --open`
    end
    stream, proc = open(cmd)
    try
        line = readline(stream)
        @show line
        @test match(r"^Listening on 127.0.0.1\:[0-9]*", line) !== nothing
        line = readline(stream)
        @show line
        m = match(r"^zmq_url=(.*)$", line)
        @test m !== nothing
        zmq_url = m[1]
        @test zmq_url == "tcp://127.0.0.1:5568"
        line = readline(stream)
        m = match(r"^web_url=(.*)$", line)
        @show line
        @test m !== nothing
        web_url = m[1]

        vis = Visualizer(zmq_url)
        @test url(vis) == web_url

        if !haskey(ENV, "CI")
            wait(vis)
        end
        setobject!(vis, HyperSphere(Point(0., 0, 0), 0.5))
        settransform!(vis, Translation(0, 0, 1))
    finally
        kill(proc)
    end
end


text:
import Base.rand

import StatsBase.quantile
import StatsBase.median
import Distributions.support

import Distributions.ecdf
import Base:
    minimum, maximum,
    max, min

"""
    UncertainScalarKDE(d::KernelDensity.UnivariateKDE, values::AbstractVector{T}, range, pdf)

An empirical value represented by a distribution estimated from actual data.

## Fields

- **`distribution`**: The `UnivariateKDE` estimate for the distribution of `values`.
- **`values`**: The values from which `distribution` is estimated.
- **`range`**: The values for which the pdf is estimated.
- **`pdf`**: The values of the pdf at each point in `range`.
"""
struct UncertainScalarKDE{T} <: AbstractUncertainScalarKDE{T}
    distribution::KernelDensity.UnivariateKDE
    values::AbstractVector{T}
    range
    pdf::StatsBase.Weights
end

"""
    TruncatedUncertainScalarKDE

A truncated [`UncertainScalarKDE`](@ref).
"""
struct TruncatedUncertainScalarKDE{T} <: AbstractUncertainScalarKDE{T}
    distribution::KernelDensity.UnivariateKDE
    values::AbstractVector{T}
    range
    pdf::StatsBase.Weights
end

function summarise(uv::AbstractUncertainScalarKDE{T}) where {T}
    dist = uv.distribution
    range = uv.range
    dist = typeof(uv.distribution)
    _type = typeof(uv)
    "$_type($dist, range = $range)"
end

Base.show(io::IO, uv::AbstractUncertainScalarKDE{T}) where {T} = print(io, summarise(uv))

"""
    rand(uv::AbstractUncertainScalarKDE)

Sample a random number from an uncertain value represented by a kernel
density estimate.
"""
function rand(uv::AbstractUncertainScalarKDE)
    # Box width
    δ = step(uv.range)

    # Sample a box
    sampled_val = sample(uv.range, uv.pdf)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    rand(uv::AbstractUncertainScalarKDE, n::Int)

Sample a random number from an uncertain value represented by a kernel
density estimate.
"""
function rand(uv::AbstractUncertainScalarKDE, n::Int)
    # Box width
    δ = step(uv.range)

    # Sample n boxes according to estimated pdf
    sampled_vals = Vector{Float64}(undef, n)
    sample!(uv.range, uv.pdf, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

"""
    ecdf(uv::AbstractUncertainScalarKDE)

Empirical cumulative distribution function for an uncertain value approximated
by kernel density estimation.
"""
ecdf(uv::AbstractUncertainScalarKDE) = cumsum(uv.pdf)

"""
The mode (most likely value) of an uncertain value represented by a
kernel density estimate.
"""
mode(uv::AbstractUncertainScalarKDE) = uv.range(findmax(uv.distribution.density)[2])

"""
    quantile(uv::UncertainScalarKDE, q)

Return the `q`-th quantile of the distribution furnishing the uncertain value.
"""
function quantile(uv::AbstractUncertainScalarKDE{T}, q) where T
    uv.range[findfirst(ecdf(uv) .> q)]
end

median(uv::AbstractUncertainScalarKDE{T}) where T = quantile(uv, 0.5)

"""
    support(uv::UncertainScalarKDE)

Return the support of an uncertain value furnished by a kernel density
estimate.
"""
support(uv::AbstractUncertainScalarKDE{T}) where T =
    interval(minimum(uv.range), maximum(uv.range))

"""
    getquantileindex(uv::AbstractUncertainScalarKDE, q::Float64)

Return the index of the range/density value corresponding to the `q`-th quantile
of an uncertain value furnished by a kernel density estimate.
"""
function getquantileindex(uv::AbstractUncertainScalarKDE{T}, q::Float64) where T
    findfirst(ecdf(uv) .> q)
end

minimum(uv::AbstractUncertainScalarKDE) = minimum(uv.range)
maximum(uv::AbstractUncertainScalarKDE) = maximum(uv.range)

min(uv::AbstractUncertainScalarKDE) = minimum(uv.range)
max(uv::AbstractUncertainScalarKDE) = maximum(uv.range)

export
AbstractUncertainScalarKDE,
UncertainScalarKDE
ecdf,
support,
getquantileindex,
UnivariateKDE,
minimum,
maximum


text:
"""
    iterate(cplx::AbstractComplex) -> Union{Nothing, AbstractCell}

Advance the iterator to obtain the next cell in the complex `cplx`
"""
function iterate(cplx::AbstractComplex, (tid, d, did)=(0, -1, 1))
    tid >= length(cplx) && return nothing
    if d < 0
        d = 0
    end
    if did > size(cplx, d)
        d += 1
        did = 1
    end
    return cells(cplx, d)[did], (tid+1, d, did+1)
end

"""
    iterate(flt::Filtration{C,FI}) -> Union{Nothing, Tuple{FI,Vector{Tuple{Int,UInt}}}}

Advance the filtration iterator `flt` to obtain the next filtration value of type `FI` with
a corresponding collection of cell parameters: dimension & identifier.
"""
function iterate(flt::Filtration{C,FI}, state=nothing) where {C<:AbstractComplex, FI<:AbstractFloat}
    ord = order(flt)
    if state === nothing # calculate initial state
        idx = 1
        fval = ord[idx][3]
        incr = (maximum(flt)-minimum(flt)) / flt.divisions
    else
        idx, fval, incr = state
    end
    idx > length(ord) && return nothing # done
    splxs = Tuple{Int,UInt}[] #simplex dim & index
    while idx <= length(ord) && (fval+incr) >= ord[idx][3]
        push!(splxs, ord[idx][1:2])
        idx += 1
    end
    nextfval = fval+incr
    if idx <= length(ord) && isinf(flt.divisions)
        nextfval = ord[idx][3]
    end
    return (fval, splxs), (idx, nextfval, incr)
end

"""
Persistent homology group iterator for a filtration
"""
struct PersistentHomology{R<:AbstractPersistenceReduction, T} <: AbstractHomology
    diagram::Dict{Int, PersistenceDiagram{T}}
end
function persistenthomology(::Type{R}, flt::Filtration{C,FI};
                            kwargs...) where {R<:AbstractPersistenceReduction, C<:AbstractComplex, FI}
    dgm = diagram(R, flt; kwargs...)
    return PersistentHomology{R,FI}(dgm)
end
persistenthomology(flt::Filtration; kwargs...) = persistenthomology(TwistReduction, flt; kwargs...)

diagram(ph::PersistentHomology) = ph.diagram
show(io::IO, h::PersistentHomology{R}) where {R <: AbstractPersistenceReduction} =
    print(io, "Persistent Homology with $R")

#
# Iterator methods
#
values(ph::PersistentHomology) = Iterators.flatten( [birth.(dgm); death.(dgm)] for (d,dgm) in ph.diagram) |> unique |> sort!
length(ph::PersistentHomology) = length(values(ph))

"""Return persistent homology group type: filtraton value & Betti numbers"""
eltype(ph::PersistentHomology{R,T}) where {R, T} = Tuple{T, NTuple}

function iterate(ph::PersistentHomology,
                 p=(1, sort!(collect(keys(ph.diagram))), values(ph)) )
    fidx, dims, fvals = p
    fidx > length(fvals) && return nothing
    maxd = length(dims)
    fv = fvals[fidx]
    βₚ = ntuple(maxd) do i
        if isinf(fv)
            count(v->isinf(death(v)), ph.diagram[i-1])
        else
            count(v->fv ∈ v, ph.diagram[i-1])
        end
    end
    return (fv, βₚ), (fidx+1, dims, fvals)
end

betti(g::PersistentHomology) = map(last, g)

generators(ph::PersistentHomology{PersistentCocycleReduction{R}, T}) where {R, T} =
    Dict( d=>getfield.(cyc, :g) for (d, cyc) in ph.diagram)


text:
@testset "Galileo E1B" begin

    @test @inferred(get_center_frequency(GalileoE1B)) == 1.57542e9Hz
    @test @inferred(get_code_length(GalileoE1B)) == 4092
    @test @inferred(get_secondary_code_length(GalileoE1B)) == 1
    @test @inferred(get_code(GalileoE1B, 0, 1)) == 1
    @test @inferred(get_code(GalileoE1B, 0.0, 1)) == 1
    @test @inferred(get_code(GalileoE1B, 0.5, 1)) == -1
    @test @inferred(get_code(GalileoE1B, 1.0, 1)) == 1
    @test @inferred(get_code(GalileoE1B, 1.5, 1)) == -1
    @test @inferred(get_code_unsafe(GalileoE1B, 0.0, 1)) == 1
    @test @inferred(get_data_frequency(GalileoE1B)) == 250Hz
    @test @inferred(get_code_frequency(GalileoE1B)) == 1023e3Hz

end


text:
module FESpacesTestsSeq
using PartitionedArrays
include("../FESpacesTests.jl")
prun(FESpacesTests.main,sequential,(2,2))
end # module


text:
using CSTParser: KEYWORD

"Strip a trailing line of purely whitespace (if present)."
function maybe_strip_trailing_trivia_line(body)
    ws = trailing_trivia(body)
    idx = findlast(c->c=='\n', ws)
    idx == 0 && return body
    any(c->!isspace(c), ws[idx:end]) && return body
    TriviaReplacementNode(nothing, body, leading_trivia(body), ws[1:prevind(ws, idx)])
end

function last_line(ws)
    idx = findlast(c->c=='\n', ws)
    idx == 0 && return ""
    ws[idx+1:end]
end

function except_first_line(ws)
    idx = findfirst(ws, '\n')
    idx == 0 && return ""
    ws[idx+1:end]
end

function except_leading_ws(ws)
    idx = findfirst(!isspace, ws)
    idx == 0 && return ""
    ws[idx:end]
end

function compute_resolve_inline_body(expr, inline_body, replace_expr)
    # Inline else body
    _last_line = last_line(trailing_trivia(children(expr)[2]))
    indent = (isempty(_last_line) ? 0 : sum(charwidth, _last_line)) - line_pos(replace_expr, first(replace_expr.span))
    # Keep any trivia before the first statement of the body
    ws = string(prev_node_ws(inline_body), leading_trivia(inline_body))
    body = format_addindent_body(TriviaReplacementNode(nothing, inline_body,
        ws, trailing_trivia(inline_body)), -indent)
    body = maybe_strip_trailing_trivia_line(body)
    # Strip the first line (the content after the if/else)
    body = TriviaReplacementNode(nothing, body, except_leading_ws(except_first_line(leading_trivia(body))), trailing_trivia(body))
    buf = IOBuffer()
    print_replacement(buf, body, true, true)
    String(take!(buf))
end

function resolve_inline_body(dep, resolutions, expr, replace_expr)
    push!(resolutions, TextReplacement(dep, replace_expr.span, compute_resolve_inline_body(expr, children(expr)[3], replace_expr)))
end

function resolve_delete_expr(dep, resolutions, expr, replace_expr, prefix="")
    if length(children(expr)) <= 4
        # Find our current indentation
        prev = findprev_str(i -> !isspace(i), replace_expr.buffer, prevind(replace_expr.buffer, first(replace_expr.fullspan)))
        if prev == 0
            prev = 1
        else
            # Essentially, we drop everything except for the content up until
            # the first newline from the leading whitespace and do the opposite
            # for the trailing whitespace
            while true
                prev = nextind(replace_expr.buffer, prev)
                c = replace_expr.buffer[prev]
                !isspace(c) && break
                if c == '\n'
                    prev = nextind(replace_expr.buffer, prev)
                    break
                end
            end
        end
        push!(resolutions, TextReplacement(dep, prev:last(replace_expr.fullspan),
                                           string(prefix, except_first_line(trailing_trivia(replace_expr)))))
    elseif isexpr(children(expr)[4], KEYWORD, Tokens.ELSE)
        # Inline else body
        push!(resolutions, TextReplacement(dep, replace_expr.span, string(prefix,
            compute_resolve_inline_body(expr, children(expr)[5], replace_expr))))
    else
        # elseif -> if
        indent = sum(charwidth, trailing_trivia(children(expr)[2]))
        repl = ChildReplacementNode(nothing, Any[], expr)
        eif = children(expr)[4]
        @assert isexpr(eif, KEYWORD, Tokens.ELSEIF)
        push!(repl.children, ReplacementNode("if", leading_trivia(eif), trailing_trivia(eif)))
        append!(repl.children, children(expr)[5:end])
        buf = IOBuffer()
        print_replacement(buf, repl, false, true)
        # Here we actually replace expr, rather than replace_expr, since we're not removing the
        # expr, entirely, just removing a branch.
        push!(resolutions, TextReplacement(dep, expr.fullspan, string(prefix, String(take!(buf)))))
    end
end

    # Could be implemented in the future
    function is_statically_effect_free(expr)
        if isexpr(expr, CSTParser.Call)
            is_identifier(children(expr)[1], "isdefined") || return false
            all(is_statically_effect_free, children(expr)[2:end]) || return false
            return true
        end
        # TODO: Only if guaranteed to be defined
        isexpr(expr, CSTParser.IDENTIFIER) && return true
        isexpr(expr, CSTParser.PUNCTUATION) && return true
        isexpr(expr, CSTParser.Quotenode) && return true
        return false
    end

    function resolve_boolean(dep, resolutions, context, expr, alwaystruefalse, effects = Any[])
        p = parent(expr)
        if isexpr(p, CSTParser.BinarySyntaxOpCall)
            context.in_macrocall && return
            is_and = isexpr(children(p)[2], OPERATOR, Tokens.LAZY_AND)
            is_or = isexpr(children(p)[2], OPERATOR, Tokens.LAZY_OR)
            (is_and || is_or) || @goto out
            if expr == children(p)[1]
                if is_and ? alwaystruefalse : !alwaystruefalse
                    !isempty(effects) && return
                    push!(resolutions,
                        TextReplacement(dep,
                        first(p.span):(first(children(p)[3].span)-1), ""))
                    return
                elseif is_and ? !alwaystruefalse : alwaystruefalse
                    resolve_boolean(dep, resolutions, context, p, is_and ? false : true, effects)
                end
            else
                @assert expr == children(p)[3]
                if is_and ? alwaystruefalse : !alwaystruefalse
                    !isempty(effects) && return
                    push!(resolutions,
                        TextReplacement(dep),
                        first(children(p)[2].span):last(children(p)[3].span),
                        "")
                elseif is_and ? !alwaystruefalse : alwaystruefalse
                    fcond = children(p)[1]
                    if !is_statically_effect_free(fcond)
                        length(effects) >= 2 && return
                        if length(effects) == 1
                            effects[1] = ChildReplacementNode(nothing, Any[children(p)[1:2]..., effects[1]], p)
                        else
                            pushfirst!(effects, TriviaReplacementNode(
                                nothing,
                                fcond,
                                leading_trivia(fcond), ""))
                        end
                    end
                    resolve_boolean(dep, resolutions, context, p, is_and ? false : true, effects)
                end
            end
            return
        elseif isexpr(p, CSTParser.ConditionalOpCall)
                        is_ternary = isexpr(children(p)[2], OPERATOR, Tokens.CONDITIONAL)
            if alwaystruefalse
                push!(resolutions,
                    TextReplacement(dep,
                    first(p.span):(first(children(p)[3].span)-1), ""))
                push!(resolutions,
                    TextReplacement(dep,
                    (last(children(p)[3].span)+1):last(p.span), ""))
            else
                push!(resolutions,
                    TextReplacement(dep,
                    first(p.span):(first(children(p)[end].span)-1), ""))
            end
            return
        elseif isexpr(p, CSTParser.If)
            replace_expr = p
            if context.in_macrocall
                context.top_macrocall == parent(p) || return nothing
                is_macroname(context.top_macrocall, "static") || return nothing
                replace_expr = context.top_macrocall
            end
            indent = 0
            buf = IOBuffer()
            anythingleft = alwaystruefalse || length(children(p)) > 4
            ind = indentation(replace_expr)
            for (i, node) in enumerate(effects)
                print_replacement(buf, node, true, true)
                if i !== length(effects) || anythingleft
                    print(buf, "\n", ind)
                end
            end
            if alwaystruefalse
                repl = compute_resolve_inline_body(p, children(p)[3], replace_expr)
                push!(resolutions, TextReplacement(dep, replace_expr.span, string(String(take!(buf)), repl)))
            else
                resolve_delete_expr(dep, resolutions, p, replace_expr, String(take!(buf)))
            end
            return
        end

        @label out
        if in_statement_position(expr)
            resolve_delete_expr(dep, resolutions, expr, expr)
        else
            push!(resolutions,
                TextReplacement(dep,
                expr.span, alwaystruefalse ? "true" : "false"))
        end
    end

function replace_node(tree, node, replacement, current=tree, parent=nothing)
    (current == node) && return replacement
    isempty(children(current)) && return current
    newtree = ChildReplacementNode(parent, Any[], current)
    for c in children(current)
        if c == node
            push!(newtree.children, replacement)
        else
            push!(newtree.children, replace_node(tree, node, replacement, c, newtree))
        end
    end
    newtree
end

function indentation(tree)
    all_leading_trivia = string(prev_node_ws(tree), leading_trivia(tree))
    # Find the amount the number characters from the last newline to the start
    # expr's span.
    lastn = rsearch(all_leading_trivia, '\n')
    all_leading_trivia[lastn+1:end]
end

"""
    Aligns new lines to the opening parenthesis of a function call's argument list
"""
function format_align_arguments(tree)

    call = function_def_call(tree)
    lparen = children(call)[2]
    nindent = mapreduce(charwidth, +, 0, tree.buffer[1 + (first(tree.span):first(lparen.span))]) +
              mapreduce(charwidth, +, 0, indentation(tree))
    ftree = format_setindent_body(call, nindent)
    rtree = replace_node(tree, call, ftree)
end

function apply_formatter(f, tree)
    rtree = f(tree)
    buf = IOBuffer()
    print_replacement(buf, rtree)
    TextReplacement(nothing, tree.span, String(take!(buf)))
end


text:
simple_loss(W::AbstractMatrix, b::AbstractVector, x::AbstractArray) = sum(W * x .+ b)

@testset "simplegrad" begin
    W, b, x = rand(128, 784), rand(128), rand(784, 100)

    val1, g = grad(simple_loss, W, b, x)

    simple_loss_grad = simplegrad(simple_loss, W, b, x)
    val2, dW, db, dx = simple_loss_grad(W, b, x)

    @test val1 == val2
    @test g[1] == dW
    @test g[2] == db
    @test g[3] == dx
end


text:
using Random; Random.seed!(1234);

include("src/data.jl")
using .DataUtil
include("src/stepwise/SelectiveInference.jl")
using .SelectiveInference
include("src/selective_p.jl")
using .MultiTruncatedDistributions
include("src/intervals.jl")
using .Intersection
using LinearAlgebra:I
using Distributions
using IntervalSets
using Parameters
using ProgressMeter

using Plots
default(:bglegend, plot_color(default(:bg), 0.8))
default(:fglegend, plot_color(ifelse(isdark(plot_color(default(:bg))), :white, :black), 0.6))

function selective_p1(test_prob::SelectiveHypothesisTest; isOC=false)
    @unpack alg, alg_result, test_stats, 𝐗, 𝐲_obs, 𝚺 = test_prob
    M_obs = alg_result.model.M
    1 ∈ M_obs || return NaN
    for (i, j) in enumerate(M_obs)
        # β₁ ≠ 0
        j == 1 || continue
        test_stat = test_stats[i]
        z_obs, σ²_z = test_stat.z_obs, test_stat.σ²_z

        # compute Z = {z | selection_event}
        z_min = -abs(z_obs) - 10 * √(σ²_z); z_max = abs(z_obs) + 10 * √(σ²_z)
        Z = isOC ?
        SelectiveInference.region(𝐗, 𝐲_obs, 𝚺, test_stat, alg)[2] ∩ [z_min..z_max] :
        SelectiveInference.truncated_interval(alg_result.model, 𝐗, 𝚺, test_stat, alg;
            z_min=z_min, z_max=z_max)

        𝜋 = cdf(TruncatedDistribution(Normal(0, √(σ²_z)), Z), z_obs)
        return 2 * min(𝜋, 1 - 𝜋)
    end
end

# TPR experiment
    function TPR_experiment(param::DataParam, alg::StepwiseFeatureSelection;
    isOC=false,
    α=0.05,
    denominator_upper=100,
    iter_size=5)
    @assert param.𝛃[1] ≠ 0
    TPR = fill(NaN, iter_size)
    @showprogress "[isOC:$(isOC),N:$(param.N),p:$(length(param.𝛃)),β:$(param.𝛃[1])]" for l = 1:iter_size
        numerator = 0
        denominator = 0
        for _ = 1:denominator_upper
            try
                𝐗, 𝐲_obs = make_dataset(param)
                alg_result = stepwise(𝐗, 𝐲_obs, I, alg)
                p_value = selective_p1(SelectiveHypothesisTest(alg, alg_result, 𝐗, 𝐲_obs, I); isOC=isOC)
                isnan(p_value) && continue
                numerator += p_value < α ? 1 : 0
                denominator += 1
            catch e
                println(e)
                continue
            end
        end
        TPR[l] = numerator / denominator
    end
    return TPR
end

function main(alg::StepwiseFeatureSelection)
    # parameters
    β₀ = 0.0
    Ns = (50, 100, 150)
    ps = (10, 20, 50)
    βs = (0.01, 0.25, 0.5, 1.)

    iter_size = 10
    TPR = Array{Float64}(undef, iter_size, length(Ns), length(ps), length(βs))
    TPR_OC = similar(TPR)
    α = 0.05

    # run experiment
    for (i, N) in enumerate(Ns)
        for (j, p) in enumerate(ps)
            for (k, β) in enumerate(βs)
                𝛃 = [fill(β, floor(Int, p / 2)); zeros(ceil(Int, p / 2))]
                param = DataParam(N=N, β₀=β₀, 𝛃=𝛃)
                TPR[:,i,j,k] = TPR_experiment(param, alg; isOC=false, α=α, iter_size=iter_size)
                TPR_OC[:,i,j,k] = TPR_experiment(param, alg; isOC=true, α=α, iter_size=iter_size)
            end
        end
    end

    # plot results
    TPR_mean = reshape(mean(TPR, dims=1), size(TPR, 2), size(TPR, 3), size(TPR, 4))
    TPR_SE = reshape(std(TPR, dims=1) ./ √size(TPR, 1), size(TPR, 2), size(TPR, 3), size(TPR, 4))
    TPR_OC_mean = reshape(mean(TPR_OC, dims=1), size(TPR_OC, 2), size(TPR_OC, 3), size(TPR_OC, 4))
    TPR_OC_SE = reshape(std(TPR_OC, dims=1) ./ √size(TPR_OC, 1), size(TPR_OC, 2), size(TPR_OC, 3), size(TPR_OC, 4))
    p = Array{typeof(plot())}(undef, length(βs))
    h_cs = (:red2, :violetred2, :darkred)
    o_cs = (:blue2, :skyblue, :darkblue)
    for k = eachindex(βs)
        p[k] = plot()
        for j = eachindex(ps) # p
            # Homotopy
            plot!(p[k], 1:length(Ns), (TPR_mean - TPR_SE)[:,j,k]; fillrange=(TPR_mean + TPR_SE)[:,j,k], fillalpha=0.2, α=0, label=nothing, color=h_cs[j])
            plot!(p[k], 1:length(Ns), TPR_mean[:,j,k], label="Homotopy, p=$(ps[j])", color=h_cs[j], markershape=:circle)
            # Quadratic
            plot!(p[k], 1:length(Ns), (TPR_OC_mean - TPR_OC_SE)[:,j,k]; fillrange=(TPR_OC_mean + TPR_OC_SE)[:,j,k], fillalpha=0.2, α=0, label=nothing, color=o_cs[j])
            plot!(p[k], 1:length(Ns), TPR_OC_mean[:,j,k], label="Quadratic, p=$(ps[j])", color=o_cs[j], markershape=:rect)
        end
        # plot setting
        plot!(p[k]; ylims=(0, 1))
        plot!(p[k]; legend=:best)
        plot!(p[k]; xlabel="Sample size", ylabel="True Positive Rate (TPR)")
        plot!(p[k]; guidefontsize=14, legendfontsize=9, tickfontsize=9)
        plot!(p[k]; xticks=(1:length(Ns), Ns))
        plot!(p[k]; title="signal:$(βs[k])")
    end
    fig = plot(p..., layout=(1, 4), titlefontsize=16, size=(900, 900))
    savefig(fig, "img/TPR_forward_backward.pdf")
end

main(StepwiseFeatureSelection())


text:
using IRTools: @dynamo, argument!, IR

# Implementation of `invoke(f, argtypes::Type, args...`, where `args` don't have
# to conform to `argtypes`.  i.e., it allows you to call `f` on `args`, but using the
# method for `argtypes`.
# See https://julialang.slack.com/archives/CJ357CY2Y/p1579598240011500.

@dynamo function sneakyinvoke(f, ::Type{T}, args...) where T<:Tuple
  ir = IR(f, T.parameters...)
  argument!(ir, at = 2)
  return ir
end

f(x::Float64) = x^2

sneakyinvoke(f, Tuple{Float64}, 2)


text:
using EconPDEs, Distributions

Base.@kwdef struct HaddadModel
  # consumption process parameters
  μbar::Float64 = 0.018
  vbar::Float64 = 0.00052
  κμ::Float64 = 0.3
  νμ::Float64 = 0.456
  κv::Float64 = 0.012
  νv::Float64 = 0.00472

  # active capital
  αbar::Float64 = 1.2
  λ::Float64 = 0.018

  # utility parameters
  ρ::Float64 = 0.0132
  γ::Float64 = 10.0
  ψ::Float64 = 1.5
end

function initialize_stategrid(m::HaddadModel; μn = 30, vn = 30)
  μbar = m.μbar; vbar = m.vbar ; κμ = m.κμ ; νμ = m.νμ ; κv = m.κv ; νv = m.νv ; αbar = m.αbar ; λ = m.λ ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ

  σ = sqrt(νμ^2 * vbar / (2 * κμ))
  μmin = quantile(Normal(μbar, σ), 0.001)
  μmax = quantile(Normal(μbar, σ), 0.999)
  μs = range(μmin,  μmax, length = μn)

  α = 2 * κv * vbar / νv^2
  β = νv^2 / (2 * κv)
  vmin = quantile(Gamma(α, β), 0.001)
  vmax = quantile(Gamma(α, β), 0.999)
  vs = range(vmin,  vmax, length = vn)

  OrderedDict(:μ => μs, :v => vs)
end

function (m::HaddadModel)(state::NamedTuple, y::NamedTuple)
  (; μbar, vbar, κμ, νμ, κv, νv, αbar, λ, ρ, γ, ψ) = m
  (; μ, v) = state
  (; p, pμ_up, pμ_down, pv_up, pv_down, pμμ, pμv, pvv) = y

  # drift and volatility of μ, ν, p
  μc = μ
  σc = sqrt(v)
  μμ = κμ * (μbar - μ)
  σμ = νμ * sqrt(v)
  μv = κv * (vbar - v)

  pμ = (μμ >= 0) ? pμ_up : pμ_down
  pv = (νμ >= 0) ? pv_up : pv_down
  σv = νv * sqrt(v)
  σp_Zμ = pμ / p * σμ
  σp_Zv = pv / p * σv
  σp2 = σp_Zμ^2 + σp_Zv^2
  μp = pμ / p * μμ + pv / p * μv + 0.5 * pμμ / p * σμ^2 + 0.5 * pvv / p * σv^2

  # Market Price of Risk κ
  αstar = αbar - sqrt(λ * αbar / (0.5 * γ * (σc^2 + σp_Zμ^2 + σp_Zv^2)))
  αstar = clamp(αstar, 0.0, 1.0)
  κ_Zc = αstar * γ * σc
  κ_Zμ = (αstar - (1 / γ - 1) / (ψ - 1)) * γ * σp_Zμ
  κ_Zv = (αstar - (1 / γ - 1) / (ψ - 1)) * γ * σp_Zv
  κ2 = κ_Zc^2 + κ_Zμ^2 + κ_Zv^2

  # Interest rate r
  #r = ρ + μc / ψ - (1 - 1 / ψ) / (2 * γ) * κ2 - (1 / ψ - γ) / (2 * γ * (ψ - 1)) * σp2 - (1 - γ) / (γ * ψ) * (κμ * σμ + κv * σv) + 1 / ψ  * (κ_Zc * σc + κμ * σμ + κv * σv)
  #out = p * (1 / p + μc + μp - r - κ_Zc * σc - κ_Zμ * σp_Zμ - κ_Zv * σp_Zv)

  pt = - p * (1 / p - ρ + (1 - 1 / ψ) * (μc - 0.5 * γ * σc^2 * (1 - (αstar - 1)^2)) + μp + (0.5 * (1 / ψ - γ) / (1 - 1 / ψ) + 0.5 * γ * (1 - 1 / ψ) * (αstar - 1)^2) * σp2)
  return (; pt)
end

m = HaddadModel()
μn = 30
vn = 30
stategrid = initialize_stategrid(m; μn = μn, vn = vn)
yend = OrderedDict(:p => ones(μn, vn))
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5


text:
using CUDA

# Each block is assigned a row.
# Each thread checks  the main row
# with one of the element of the columns pointed by the row.
function countTri_elem_per_th(csr_rows, col_indx, out_sum, max_nnz)
    bid::Int32 = blockIdx().x + (blockIdx().y-1)*gridDim().x + (blockIdx().z-1)*gridDim().x*gridDim().y
    rows::Int32 = size(csr_rows)[1]
    nnz::Int32 = size(col_indx)[1]
    tid::Int32 = threadIdx().x + (threadIdx().y-1)*blockDim().x
    stride::Int32 = blockDim().x*blockDim().y
    if bid > rows
        return nothing
    end
    sh_row = @cuDynamicSharedMem(Int32, stride)
    sh_sum = @cuDynamicSharedMem(Int32, stride, stride*sizeof(Int32))
    sh_cols = @cuDynamicSharedMem(Int32, max_nnz*max_nnz, sizeof(Int32)*2*stride)

    @inbounds sh_sum[tid] = 0
    @inbounds row_start::Int32 =  csr_rows[bid]
    @inbounds row_end::Int32 = (bid == rows) ?  nnz+1 : csr_rows[bid + 1]
    len::Int32 = row_end-row_start
    if tid <= len
        @inbounds sh_row[tid] = col_indx[row_start + tid - 1]
    end

    sync_threads()

    # #Read all the other rows
    start_save = 0
    for i::Int32 = 1:len
        @inbounds col::Int32 = sh_row[i]
        @inbounds row_start = csr_rows[col]
        @inbounds row_end = (col == rows) ?  nnz+1 : csr_rows[col + 1]
        templen::Int32 =  row_end-row_start
        if tid  <= templen
            @inbounds sh_cols[start_save+tid] = col_indx[row_start+tid-1]
        end
        start_save += templen
    end
    sync_threads()
    sum::Int32 = 0
    for i::Int32 = 0:stride:(start_save-1)
        if tid + i <= start_save
            for j = 1:len
                @inbounds sum += sh_row[((tid+j-2)%len)+1] == sh_cols[i + tid]
            end
        end
    end
    sync_threads()
    w::Int32 = 16
    while w > 0
        sum += shfl_down_sync(0xffffffff, sum, w)
        w = floor(w/2)
    end
    if threadIdx().x == 1
        @inbounds sh_sum[threadIdx().y] = sum
    end
    sync_threads()
    if tid == 1
        for i::Int32 = 2:blockDim().y
            @inbounds sum += sh_sum[i]
        end
        @inbounds out_sum[bid] = sum
    end
    return nothing
end

function CountTriangles_ElemPerTh(csr_rows, col_indx, max_per_row)
    dev_csr_rows = CuArray(csr_rows)
    dev_col_indx = CuArray(col_indx)
    dev_out_sum = CUDA.zeros(Int32, size(csr_rows)[1])
    th_groups = 1
    threads = 0
    while 32*th_groups < max_per_row
        th_groups += 1
    end
    threads = 32*th_groups
    blocks = Int64(floor((size(csr_rows)[1])/(512*512)+ 1))
    println("2 Total blocks: ", 512*512*blocks, " vs rows: ", size(csr_rows)[1])
    println("\t Threads: ", threads)
    @cuda(blocks=(512,512,blocks), threads=(32,th_groups,1),
        shmem=sizeof(Int32)*(threads*2 + max_per_row*max_per_row),
        countTri_elem_per_th(dev_csr_rows, dev_col_indx, dev_out_sum, max_per_row))
    out_sum = Array(dev_out_sum)
    # println(dev_out_sum)
    return reduce(+, out_sum, init=0)
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(@__DIR__), "src", "CodingTheory.jl"))
using .CodingTheory

using DataFrames, CSV
using DataFramesMeta: @where
using Formatting
# using Lazy: @>

const ruled_out = Tuple[
        # (3, 11, 5), # n words ≈ 186 ≠ 729
        # (21, 4, 3), # n words ≈ 359 ≠ 2401
        # (6, 7, 3), # n words ≈ 2995 ≠ 7776
        # (2, 23, 7), # this is one of the Golay codes, silly me
    ]

function load_df(pathname::AbstractString)
    # use bounds.jl to get the hamming bounds
    df = DataFrame(CSV.read(pathname))
    println("The complete list has $(format(nrow(df), commas = true)) observations.")
    insert!(df, ncol(df) + 1, CodingTheory.sizeof_perfect_code.(df.q, df.n, df.d), :size_of_perfect_code)
    df = round.(BigInt, df)

    return df
end

function filter_df(
    df::DataFrame;
    q_not_one::Bool=true,
    d_not_one::Bool=true,
    d_not_two::Bool=true,
    d_not_n::Bool=true,
    d_leq_n::Bool=true,
    not_hamming_perfect::Bool=true,
    d_not_even::Bool=true,
    q_prime_power::Bool=false,
    q_not_prime_power::Bool=true,
    hamming_bound_isinteger::Bool=true,
    hamming_bound_is_lower_bound::Bool=true
    )

    if q_not_one
        df = @where(df, :q .≠ 1)
    end

    if d_not_one
        df = @where(df, :d .≠ 1)
    end

    if d_not_two
        df = @where(df, :d .≠ 2)
    end

    if d_not_n
        df = @where(df, :d .≠ :n)
    end

    if d_leq_n
        df = @where(df, :d .≤ :n)
    end

    if not_hamming_perfect
        df = @where(df, .! ishammingperfect.(:q, :n, :d))
    end

    if d_not_even
        df = @where(df, .! iseven.(:d))
    end

    if q_prime_power
        df = @where(df, CodingTheory.isprimepower.(:q))
    end

    if q_not_prime_power
        df = @where(df, .! CodingTheory.isprimepower.(:q))
    end

    if hamming_bound_isinteger
        df = @where(df, isinteger.(:hamming_bound))
    end

    if hamming_bound_is_lower_bound
		df = @where(df, isequal.(:hamming_bound, :smallest_bound))
        # df = @where(df, isequal.(:hamming_bound, min.(:hamming_bound, :singleton_bound)))
    end

    # Filter out those (q, n, d) configurations already ruled out.
    # df = @where(df, tuple.(:q, :n, :d) .∉ Ref(ruled_out))

    return df
end

function view_df(
    df::DataFrame;
    nrows::Integer=10,
    q_not_one::Bool=true,
    d_not_one::Bool=true,
    d_not_two::Bool=true,
    d_not_n::Bool=true,
    d_leq_n::Bool=true,
    not_hamming_perfect::Bool=true,
    d_not_even::Bool=true,
    q_prime_power::Bool=false,
    q_not_prime_power::Bool=true,
    hamming_bound_isinteger::Bool=true,
    hamming_bound_is_lower_bound::Bool=true
    )
    df = filter_df(
        df;
        q_not_one=q_not_one,
        d_not_one=d_not_one,
        d_not_two=d_not_two,
        d_not_n=d_not_n,
        d_leq_n=d_leq_n,
        not_hamming_perfect=not_hamming_perfect,
        d_not_even=d_not_even,
        q_prime_power=q_prime_power,
        q_not_prime_power=q_not_prime_power,
        hamming_bound_isinteger=hamming_bound_isinteger,
        hamming_bound_is_lower_bound=hamming_bound_is_lower_bound
        )

    # sort last!
    df = sort(df, :hamming_bound)
    # df = sort(df, :size_of_perfect_code)

    println("The filtered list has $(format(nrow(df), commas = true)) observations matching the filtered criteria.  Now showing the 10 smallest values.")
	# return select(df, Not([:smallest_bound]))[1:nrows, :]
	return df[1:nrows, :]
    # return select(df, Not([:singleton_bound]))[1:nrows, :]
end

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

datafile = "/Users/jakeireland/projects/CodingTheory.jl/other/research/bounds/bound_integers_100000_very_loose_q_neq_1.csv" # <=====
# datafile = "/Users/jakeireland/projects/CodingTheory.jl/other/bound_integers_1000000.csv"

# df = load_df("/Users/jakeireland/projects/CodingTheory.jl/other/research/bounds/bound_integers_100000_very_loose_q_neq_1.csv")
# println(df[end, :])

# df = load_df("/Users/jakeireland/projects/CodingTheory.jl/other/bound_integers_1000000.csv")

# println(df[end, :])

df = view_df(
    load_df(datafile);
    nrows=10,
    q_not_one=true,
    d_not_one=true, # distance equal to one is trivially the list of all words with q and n
    d_not_two=true, # distance equal to two is trivial as the words' coordinates sum to zero modulo q
    d_not_n=true, # distance equal to two is trivially the repetition code
    d_leq_n=true,
    not_hamming_perfect=true,
    d_not_even=false, # cannot get a perfect code if d is even, because of how hamming bound is calculated
    q_prime_power=false, # inverse of the below
    q_not_prime_power=true, # we are interested in alphabet size not prime power
    hamming_bound_isinteger=true,
    hamming_bound_is_lower_bound=true # cannnot obtain perfect code if singleton bound is less than hamming bound
    )

println(df)

# df2 = filter_df(
#     load_df(datafile);
#     nrows=10,
#     q_not_one=true,
#     d_not_one=true, # distance equal to one is trivially the list of all words with q and n
#     d_not_two=true, # distance equal to two is trivial as the words' coordinates sum to zero modulo q
#     d_not_n=true, # distance equal to two is trivially the repetition code
#     d_leq_n=true,
#     not_hamming_perfect=true,
#     d_not_even=false, # cannot get a perfect code if d is even, because of how hamming bound is calculated
#     q_prime_power=false, # inverse of the below
#     q_not_prime_power=false, # we are interested in alphabet size not prime power
#     hamming_bound_isinteger=true,
#     hamming_bound_is_lower_bound=false # cannnot obtain perfect code if singleton bound is less than hamming bound
#     )
#
# df2 = select(df2, Not([:smallest_bound]))
#
# D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[], size_of_perfect_code = Number[])
#
# for row in eachrow(df2)
#     if isinteger(log(row.q, hamming_bound(row.q, row.n, row.d)))
#         push!(D, row)
#     end
# end
#
# println(size(D))
# println(sort(D, :size_of_perfect_code))

## original code:
# df = round.(BigInt, DataFrame(CSV.read(datafile)))
# df = @where(df, .! isprimepower.(:q))
# df = sort(df, :hamming_bound)[1:10, :]
# println(df)


text: using Test, JuliaReachTemplatePkg

max_stars_repo_name: JuliaReach/JuliaReachTemplatePkg

max_issues_repo_name: JuliaReach/JuliaReachTemplatePkg

max_forks_repo_name: JuliaReach/JuliaReachTemplatePkg


text:
# NOTE: create symlinks called "results_old.json" and "results_new.json" in this directory
# to the results files you want to compare

using BenchmarkTools, PyPlot, Unitful

const tref = 1000

function xycmp(results_old, results_new, filterstrs...)
    x, y = Float64[], Float64[]
    for (k, v) in results_new
        passes = true
        for str in filterstrs
            passes &= occursin(str, k)
        end
        passes || continue
        if haskey(results_old, k)
            push!(x, minimum(results_old[k]).time/tref)
            push!(y, minimum(v).time/tref)
        end
    end
    x, y
end

function plotcmp(results_old, results_new, keys, filterstrs...)
    for key in keys
        results_old = results_old[key]
        results_new = results_new[key]
    end
    x, y = xycmp(results_old, results_new, filterstrs...)
    maxxy = max(maximum(x), maximum(y))
    scatter(x, y)
    plot([0, maxxy], [0, maxxy], "--")
    # title(string(keys, filterstrs))
    title(string(keys))
end

function plotcmppanels(results_old, results_new, keys, filterstrs...)
    nrows = ceil(Int, sqrt(length(keys)))
    ncols = ceil(Int, length(keys)/nrows)
    figure()
    for k = 1:length(keys)
        subplot(nrows, ncols, k)
        plotcmp(results_old, results_new, keys[k], filterstrs...)
    end
    suptitle(string(filterstrs))
end

results_old = BenchmarkTools.load("results_old.json")[1]
results_new = BenchmarkTools.load("results_new.json")[1]

results_old = results_old["bsplines"]
results_new = results_new["bsplines"]

for filt in (("use",), ("construct",))
    # for keys in (("constant",), ("linear",), ("quadratic",), ("cubic",))
    #     plotcmp(results_old, results_new, keys, filt...)
    # end
    keys = (("constant",), ("linear",), ("quadratic",), ("cubic",))
    plotcmppanels(results_old, results_new, keys, filt...)
end


text:
module NoTitle
function ui()
end
module Server
end
end

module NoUI
title = "a"
module Server
end
end

module NoServer
title = "a"
function ui()
end
end

@test_throws ErrorException("You have not defined the `ui` function in your app module") run_app(NoUI)
@test_throws ErrorException("You have not set a `title` in your app module") run_app(NoTitle)
@test_throws ErrorException("You have not defined the `Server` module in your app module") run_app(NoServer)


text:
module AudioDeclipping

using BenchmarkTools
using StructuredOptimization, AbstractOperators
using DSP, WAV
using PyPlot

function set_up()

	Nl = 2^10                           # frames length
	Fs = 16e3                           # sampling frequency
	x0, Fs0 = wavread("demo.wav")       # load input signal
	x0 = x0[:,1]
	x0 = resample(x0,Int(Fs)//Int(Fs0))
	normalize!(x0,Inf)
	Nt = length(x0)                     # signal length
	x0 .*= 0.8

	C = 0.02                                       # clip threshold
	xc = copy(x0)
	xc[abs.(xc).>=C] = C.*sign.(xc[abs.(xc).>=C])  # clip original signal

	#time win options
	K = 2 # K = 2 hanning
	win = sqrt.(tukey(Nl+1,2/K)[1:Nl])
	overlap = div(Nl,K)

	xd = zeros(x0)                    # allocate output
	x, y = Variable(Nl), Variable(Nl) # optimization variables
	cf = ls( idct(x) - y )            # cost function
	yw = zeros(Nl)

	return Fs, x, x0, xd, xc, y, yw, Nl, Nt, C, win, overlap, cf
end

function run_demo()
	slv = PANOC(tol = 1e-5, verbose = 0)
	setup = set_up()
	@time solve_problem!(slv, setup..., true)

	Fs, x, x0, xd, xc, y, yw, Nl, Nt, C, overlap, cf = setup
	Irp = find(xc.>=C)  #positive clipping indices
	Irn = find(xc.<=-C) #negative clipping indices

	SDRy  = 20*log10(norm(x0[[Irp;Irn]])/norm(x0[[Irp;Irn]]-xc[[Irp;Irn]]))
	SDRyd = 20*log10(norm(x0[[Irp;Irn]])/norm(x0[[Irp;Irn]]-xd[[Irp;Irn]]))

	println("SDRy: $(SDRy) SDRyd: $(SDRyd) diff: $(SDRyd-SDRy)")
	return setup
end

function solve_problem!(slv, Fs, x, x0, xd, xc, y, yw, Nl, Nt, C, win, overlap, cf, verb)
	its = Int64[]
	fit_tol = 1e-5
	z = 0                             # frame initial intex
	counter = 0
    winC = C.*win
	#weighted Overlap-Add
	while z+Nl < Nt
		counter += 1
		push!(its,0)

		fill!(x.x,0.) # initialize variables
		fill!(y.x,0.)

		Irp = sort(find(     xc[z+1:z+Nl]  .>=  C) ) #positive clipping indices
		Irn = sort(find(     xc[z+1:z+Nl]  .<= -C) ) #negative clipping indices
		Im  = sort(find(abs.(xc[z+1:z+Nl]) .<   C)) #non clipped indices
		yw .= xc[z+1:z+Nl].*win

		Mp = GetIndex((Nl,),Irp)
		Mn = GetIndex((Nl,),Irn)
		M  = GetIndex((Nl,),Im )

		if isempty(Irp) && isempty(Irn)
			xd[z+1:z+Nl] .+= yw.*win
			if verb @printf("%7d / %7d |           no_clip          | \n", z, Nt) end
		else
			N = 0 # number of active components in DCT
			for N = 30:30:30*div(Nl,30)
				cstr = (norm(x,0) <= N,
					norm(M*y-M*yw) <= sqrt(fit_tol),
					Mp*y >=  Mp*winC,
                    Mn*y <= -Mn*winC)
                it, = @minimize cf st cstr with slv
				its[counter] += it
                if norm(idct(~x)- ~y) <= sqrt(fit_tol) break end
			end
			@views xd[z+1:z+Nl] .+= (~y).*win
			if verb @printf("%7d / %7d | N: %7d /  %7d | \n", z, Nt, N, Nl) end
		end

		z += Nl-overlap                      # update index
	end
	return mean(its)
end

function benchmark(;verb = 0, samples = 1, seconds = 100)

	suite = BenchmarkGroup()

	tol = 1e-4
	solvers = [
               "PANOC",
               "ZeroFPR",
               "PG"
          ]
	slv_opt = ["(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)"]

	its = Dict([(sol,0.) for sol in solvers])
	for i in eachindex(solvers)

		setup = set_up()
		solver = eval(parse(solvers[i]*slv_opt[i]))

		suite[solvers[i]] =
		@benchmarkable(it = solve_problem!(solver, setup..., false),
			       setup = (
					it = 0;
					setup = deepcopy($setup);
					solver = deepcopy($solver) ),
			       teardown = (
					  $its[$solvers[$i]] = it;
					  ),
			       evals = 1, samples = samples, seconds = seconds)
	end

	results = run(suite, verbose = (verb != 0))
	println("AudioDecliping (mean of per-frame) its")
	println(its)
	return results
end

function show_results(Fs, x, x0, xd, xc, y, yw, Nl, Nt, C, win, overlap, cf)
	idxs = 2*Nl+1:4*Nl
	t = linspace(0,length(idxs)/Fs,length(idxs))
	figure()
	plot(t,x0[idxs])
	plot(t,xd[idxs])
	plot(t,xc[idxs])
end

function save_wav(Fs, x, x0, xd, xc, y, yw, Nl, Nt, C, win, overlap, cf)
    wavwrite(x0 ,"unclipped.wav"; Fs = Fs, nbits = 16, compression=WAVE_FORMAT_PCM)
    wavwrite(xc ,"clipped.wav";   Fs = Fs, nbits = 16, compression=WAVE_FORMAT_PCM)
    wavwrite(0.8*normalize(xd,Inf) ,"declipped.wav"; Fs = Fs, nbits = 16, compression=WAVE_FORMAT_PCM)
end

end


text:
"""
    learn_compose(;
        nvars, dom_size, param=nothing, icn=ICN(nvars, dom_size, param),
        X, X_sols, global_iter=100, local_iter=100, metric=hamming, popSize=200
    )
Create an ICN, optimize it, and return its composition.
"""
function learn_compose(
    solutions,
    non_sltns,
    dom_size,
    param=nothing;
    global_iter=Threads.nthreads(),
    local_iter=64,
    metric=:hamming,
    pop_size=64,
    sampler=nothing,
    memoize=false,
)
    icn = ICN(; param=!isnothing(param))
    _, weigths = optimize!(
        icn,
        solutions,
        non_sltns,
        global_iter,
        local_iter,
        dom_size,
        param,
        metric,
        pop_size;
        sampler,
        memoize
    )
    compositions = Dictionary{Composition,Int}()
    for (bv, occurences) in pairs(weigths)
        set!(compositions, compose(deepcopy(icn), bv), occurences)
    end

    return compose(icn), icn, compositions
end

"""
    explore_learn_compose(concept; domains, param = nothing, search = :complete, global_iter = 10, local_iter = 100, metric = hamming, popSize = 200, action = :composition)

Explore a search space, learn a composition from an ICN, and compose an error function.

# Arguments:
- `concept`: the concept of the targeted constraint
- `domains`: domains of the variables that define the training space
- `param`: an optional parameter of the constraint
- `search`: either `flexible`,`:partial` or `:complete` search. Flexible search will use `search_limit` and `solutions_limit` to determine if the search space needs to be partially or completely explored
- `global_iter`: number of learning iteration
- `local_iter`: number of generation in the genetic algorithm
- `metric`: the metric to measure the distance between a configuration and known solutions
- `popSize`: size of the population in the genetic algorithm
- `action`: either `:symbols` to have a description of the composition or `:composition` to have the composed function itself
"""
function explore_learn_compose(
    domains,
    concept,
    param=nothing;
    global_iter=Threads.nthreads(),
    local_iter=64,
    metric=:hamming,
    pop_size=64,
    search=:flexible,
    complete_search_limit=1000,
    solutions_limit=100,
    sampler=nothing,
    configurations=explore(
        domains, concept, param; search, complete_search_limit, solutions_limit
    ),
    memoize=false,
)
    dom_size = maximum(length, domains)
    solutions, non_sltns = configurations
    return learn_compose(
        solutions,
        non_sltns,
        dom_size,
        param;
        local_iter,
        global_iter,
        metric,
        pop_size,
        sampler,
        memoize,
    )
end

"""
    compose_to_file!(concept, name, path; domains, param = nothing, language = :Julia, search = :complete, global_iter = 10, local_iter = 100, metric = hamming, popSize = 200)

Explore, learn and compose a function and write it to a file.

# Arguments:
- `concept`: the concept to learn
- `name`: the name to give to the constraint
- `path`: path of the output file
# Keywords arguments:
- `domains`: domains that defines the search space
- `param`: an optional paramater of the constraint
- `language`: the language to export to, default to `:julia`
- `search`: either `:partial` or `:complete` search
- `global_iter`: number of learning iteration
- `local_iter`: number of generation in the genetic algorithm
- `metric`: the metric to measure the distance between a configuration and known solutions
- `popSize`: size of the population in the genetic algorithm
"""
function compose_to_file!(
    concept,
    name,
    path;
    configurations=explore(domains, concept, param; search, search_limit, solutions_limit),
    domains,
    param=nothing,
    global_iter=Threads.nthreads(),
    language=:Julia,
    local_iter=64,
    metric=hamming,
    pop_size=64,
    search=:flexible,
    search_limit=1000,
    solutions_limit=100,
    sampler=nothing,
    memoize=false,
)
    compo, icn, _ = explore_learn_compose(
        domains,
        concept,
        param;
        configurations,
        global_iter,
        local_iter,
        metric,
        pop_size,
        search,
        search_limit,
        solutions_limit,
        sampler,
        memoize,
    )

    composition_to_file!(compo, path, name, language)
    return icn
end


text:
function updateWd(pMis::Array{RemoteChannel},Dc::Array{RemoteChannel})
@sync begin
	@async begin
		for k=1:length(pMis)
			pMis[k] = remotecall_fetch(updateWd,pMis[k].where,pMis[k],Dc[k]);
		end
	end
end
return pMis;
end

function updateWd(pMisRF::RemoteChannel,Dc::RemoteChannel)
pMis  = take!(pMisRF)
Dc = fetch(Dc);
pMis.Wd = 1./(real(Dc - pMis.dobs) + 1e-3*mean(abs(pMis.dobs[:]))) + 1im./(imag(Dc - pMis.dobs) + 1e-3*mean(abs(pMis.dobs[:])));
put!(pMisRF,pMis)
return pMisRF;
end

function multWd(pMis::Array{RemoteChannel},beta::Float64)
@sync begin
	@async begin
		for k=1:length(pMis)
			pMis[k] = remotecall_fetch(multWd,pMis[k].where,pMis[k],beta);
		end
	end
end
return pMis;
end

function multWd(pMisRF::RemoteChannel,beta::Float64)
pMis  = take!(pMisRF)
pMis.Wd *= beta;
put!(pMisRF,pMis)
return pMisRF;
end


text:
module Tests
using AverageShiftedHistograms, StatsBase, Base.Test, Distributions

info("Messy Output")
show(ash(randn(1000)))
println("\n\n")
info("Begin Tests")

@testset "Kernels" begin
    for kernel in [Kernels.biweight,
                   Kernels.cosine,
                   Kernels.epanechnikov,
                   Kernels.triangular,
                   Kernels.tricube,
                   Kernels.triweight,
                   Kernels.uniform,
                   Kernels.gaussian,
                   Kernels.logistic]
        for u in [0.0, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0]
            @test kernel(-u) ≈ kernel(u)
        end
    end
end

@testset "Ash" begin
    y = randn(10_000)
    o = ash(y, -4:.1:4; warnout = false)
    @test pdf(o, 0.0) > 0
    @test pdf(o, -10) ≈ 0.0
    @test cdf(o, -10) ≈ 0.0
    @test cdf(o, 10) ≈ 1.0

    for f in [mean, var, std, x -> quantile(x, .4), x -> quantile(x, .4:.1:.6)]
        @test_approx_eq_eps f(o) f(y) .1
    end
    AverageShiftedHistograms.histdensity(o)
    @test quantile(o, 1e-20) ≈ o.rng[1]

    # check that histogram is correct
    h = fit(Histogram, y, (-4:.1:4.1)-.05)
    @test h.weights == o.v

    fit!(o, y; warnout = false)
    @test nobs(o) == 20_000

    o = ash([.1, .1], -1:.1:1)
    @test nout(o) == 0
end

@testset "MVAsh" begin
    x = randn(1000)
    y = x + randn(1000)
    o = ash(x, y)
    show(o)

    mean(o)
    var(o)
    std(o)
    fit!(o, x, y; mx = 3)
    xyz(o)
end

end #module


text:
#=----------------------------------------------------------------------------

Copyright (c) 2015 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

----------------------------------------------------------------------------=#

using PyPlot
export viewlabspace

"""
viewlabspace:  Visualisation of L*a*b* space

```
Usage:    viewlabspace(L = 50, figNo = 1)

Arguments:     L - Lightness level in which to display slice of L*a*b* space
           figNo - PyPlot figure to use
```
Function allows interactive viewing of a sequence of images corresponding to
different slices of lightness in L*a*b* space.  Lightness varies from 0 to
100.  Initially a slice at a lightness of 50 is displayed.
Pressing 'l' or arrow up/right will increase the lightness by dL.
Pressing 'd' or arrow down/left will darken by dL.
Press 'x' to exit.

To Do:
The CIELAB colour coordinates of the cursor position within the slice images
should be updated continuously.  This is useful for determining suitable control
points for the definition of colourmap paths through CIELAB space in cmap().

See also: colourmappath, cmap
"""
function viewlabspace(L=50, figNo=1)
    # March 2013
    # November 2013  Interactive CIELAB coordinate feedback from mouse position
    # December 2015  Ported from MATLAB to Julia but loses interactivity...

    # Define some reference colours in rgb
    rgb = [1 0 0
           0 1 0
           0 0 1
           1 1 0
           0 1 1
           1 0 1]

    colours = ["red    "
               "green  "
               "blue   "
               "yellow "
               "cyan   "
               "magenta"]

    # ... and convert them to lab
    labv = srgb2lab(rgb)

    # Obtain cylindrical coordinates in lab space
    labradius = sqrt(labv[:,2].^2+labv[:,3].^2)
    labtheta = atan(labv[:,3], labv[:,2])

    # Define a*b* grid for image
    scale = 2;
    (a, b) = meshgrid(-127:1/scale:127)
    (rows,cols) = size(a)

    # Scale and offset lab coords to fit image coords
    labc = zeros(size(labv))
    labc[:,1] = round(labv[:,1])
    labc[:,2] = round(scale*labv[:,2] + cols/2)
    labc[:,3] = round(scale*labv[:,3] + rows/2)

    # Print out lab values
    labv = round(labv)
    @printf("\nCoordinates of standard colours in L*a*b* space\n\n");
    for n = 1:length(colours)
        @printf("%s  L%3d   a %4d  b %4d    angle %4.1f  radius %4d\n",
                colours[n], labv[n,1], labv[n,2], labv[n,3],
                labtheta[n], round(labradius[n]))
    end

    @printf("\n\n")

    # Generate axis tick values
    tickval = [-100, -50, 0, 50, 100]
    tickcoords = scale*tickval + cols/2
    ticklabels = ("-100", "-50", "0", "50", "100")

    renderlabslice(L, a, b, figNo, labv, labc, colours)

    @printf("Type in a desired Lightness value or type 'x' to exit\n\n")
    resp = ""
    while lowercase(resp) != 'x'
        resp = lowercase(readline())

        if resp[1] == 'x'
            break
        else
            try
                L = parse(Float64, resp)
                renderlabslice(L, a, b, figNo, labv, labc, colours)
            catch
                @printf("Type in a desired Lightness value or type 'x' to exit\n")
            end
        end
    end

end

#--------------------------------------------------------

function renderlabslice(L, a, b, figNo, labv, labc, colours)

    # Build image in lab space
    (rows,cols) = size(a)

    lab = zeros(rows,cols,3)
    lab[:,:,1] = L
    lab[:,:,2] = a
    lab[:,:,3] = b

    # Generate rgb values from lab
    rgb = lab2srgb(lab)

    # Invert to reconstruct the lab values
    lab2 = srgb2lab(rgb)

    # Where the reconstructed lab values differ from the specified values is
    # an indication that we have gone outside of the rgb gamut.  Apply a
    # mask to the rgb values accordingly
    mask = squeeze(maximum(abs(lab-lab2),3),3)

    for n = 1:3
        rgb[:,:,n] = rgb[:,:,n].*(mask.<2)  # tolerance of 2
    end

    figure(figNo)
    imshow(rgb, origin="lower", aspect="equal")
    title(@sprintf("Lab space:  Lightness %d", L))

    # Generate axis tick values
    tickval = [-100, -50, 0, 50, 100]
    scale = 2
    tickcoords = scale*tickval + cols/2
    ticklabels = ("-100", "-50", "0", "50", "100")
    xticks(tickcoords, ticklabels)
    yticks(tickcoords, ticklabels)

    xlabel("a*")
    ylabel("b*")

    hold(true),
    plot(cols/2, rows/2, "r+")   # Centre point for reference

    # Plot reference colour positions
    for n = 1:length(colours)
        plot(labc[n,2], labc[n,3], "w+")
        text(labc[n,2], labc[n,3],
             @sprintf("   %s\n  %d %d %d  ",colours[n],
                     labv[n,1], labv[n,2], labv[n,3]),
             color = [1, 1, 1])

    end

    hold(false)

end


text:
export PDDoc,
       pdDocOpen,
       pdDocClose,
       pdDocGetCatalog,
       pdDocGetNamesDict,
       pdDocGetInfo,
       pdDocGetCosDoc,
       pdDocGetPage,
       pdDocGetPageCount,
       pdDocGetPageRange,
       pdDocHasPageLabels

using ..Common

"""
```
    PDDoc
```
An in memory representation of a PDF document. Once created this type has to be
used to access a PDF document.
"""
abstract type PDDoc end

"""
```
    pdDocOpen(filepath::AbstractString) -> PDDoc
```
Opens a PDF document and provides the PDDoc document object for subsequent query
into the PDF file. `filepath` is the path to the PDF file in the relative or
absolute path format.

Remember to release the document with `pdDocClose`, once the object is used.
"""
function pdDocOpen(filepath::AbstractString)
    doc = PDDocImpl(filepath)
    update_page_tree(doc)
    update_structure_tree!(doc)
    return doc
end

"""
```
    pdDocClose(doc::PDDoc, num::Int) -> PDDoc
```
Reclaim the resources associated with a `PDDoc` object. Once called the `PDDoc`
object cannot be further used.
"""
function pdDocClose(doc::PDDoc)
  cosDocClose(doc.cosDoc)
end

"""
```
    pdDocGetPageCount(doc::PDDoc) -> Int
```
Returns the number of pages associated with the document.
"""
function pdDocGetPageCount(doc::PDDoc)
  return Cos.get_internal_pagecount(doc.pages)
end

"""
```
    pdDocGetCatalog(doc::PDDoc) -> CosObject
```
`Catalog` is considered the topmost level object in  PDF document that is
subsequently used to traverse and extract information from a PDF document. To be
used for accessing PDF internal objects from document structure when no direct
API is available.
"""
pdDocGetCatalog(doc::PDDoc) = doc.catalog

"""
```
    pdDocGetCosDoc(doc::PDDoc) -> CosDoc
```
PDF document format is developed in two layers. A logical PDF document
information is represented over a physical file structure called COS. `CosDoc` is
an access object to the physical file structure of the PDF document. To be used
for accessing PDF internal objects from document structure when no direct API is
available.

One can access any aspect of PDF using the COS level APIs alone. However, they
may require you to know the PDF specification in details and it is not the most
intuititive.
"""
pdDocGetCosDoc(doc::PDDoc)= doc.cosDoc

"""
```
    pdDocGetPage(doc::PDDoc, num::Int) -> PDPage
```
Given a document absolute page number, provides the associated page object.
"""
pdDocGetPage

"""
```
    pdDocGetPageRange(doc::PDDoc, nums::AbstractRange{Int}) -> Vector{PDPage}
    pdDocGetPageRange(doc::PDDoc, label::AbstractString) -> Vector{PDPage}
```
Given a range of page numbers or a label returns an array of pages associated
with it.
For a detailed explanation on page labels, refer to the method `pdDocHasPageLabels`.
"""
function pdDocGetPageRange(doc::PDDoc, nums::AbstractRange{Int})
    pages = []
    for i in nums
        push!(pages, pdDocGetPage(doc, i))
    end
    return pages
end

function pdDocGetPageRange(doc::PDDoc, label::AbstractString)
    catalog = pdDocGetCatalog(doc)
    pr = cosDocGetPageNumbers(doc.cosDoc, catalog, label)
    return pdDocGetPageRange(doc, pr)
end

"""
```
    pdDocHasPageLabels(doc::PDDoc) -> Bool
```
Returns `true` if the document has page labels defined.

As per PDF Specification 1.7 Section 12.4.2, a document may optionally define page
labels (PDF 1.3) to identifyeach page visually on the screen or in print. Page labels
and page indices need not coincide: the indices shallbe fixed, running consecutively
through the document starting from 0 for the first page, but the labels may be
specified in any way that is appropriate for the particular document.
"""
function pdDocHasPageLabels(doc::PDDoc)
    catalog = pdDocGetCatalog(doc)
    return get(catalog, cn"PageLabels") !== CosNull
end

"""
```
    pdDocGetInfo(doc::PDDoc) -> Dict
```
Given a PDF document provides the document information available in the `Document
Info` dictionary. The information typically includes *creation date, modification
date, author, creator* used etc. However, all information content are not
mandatory. Hence, all information needed may not be available in a document.
If document does not have Info dictionary at all this method returns `nothing`.

Please refer to the PDF specification for further details.
"""
function pdDocGetInfo(doc::PDDoc)
    obj = cosDocGetInfo(doc.cosDoc)
    obj === CosNull && return nothing
    dInfo = Dict{CDTextString, Union{CDTextString, CDDate, CosObject}}()
    for (key, val) in get(obj)
        skey = CDTextString(key)
        try
            dInfo[skey] = (skey == "CreationDate") ||
                          (skey == "ModDate") ? CDDate(val) :
                          (skey == "Trapped") ? val : CDTextString(val)
        catch
            # no op: we skipp the key that cannot be properly decoded
        end
    end
    return dInfo
end

"""
```
    pdDocGetNamesDict(doc::PDDoc) -> CosObject
```
Some information in PDF is stored as name and value pairs not essentially a
dictionary. They are all aggregated and can be accessed from one `names`
dictionary object in the document catalog. This method provides access to such
values in a PDF file. Not all PDF document may have a names dictionary. In such
cases, a `CosNull` object may be returned.

Please refer to the PDF specification for further details.
"""
function pdDocGetNamesDict(doc::PDDoc)
    catalog = pdDocGetCatalog(doc)
    ref = get(catalog, CosName("Names"))
    obj = cosDocGetObject(doc.cosDoc, ref)
end


text:
# Based on PWSCF documentation (version 6.2)

function gen_lattice_cubic( a::Float64 )
    v1 = a*[1,0,0]
    v2 = a*[0,1,0]
    v3 = a*[0,0,1]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

gen_lattice_sc(a::Float64) = gen_lattice_cubic(a)

function gen_lattice_fcc( a::Float64 )
    v1 = 0.5*a*[-1,0,1]
    v2 = 0.5*a*[0,1,1]
    v3 = 0.5*a*[-1,1,0]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

function gen_lattice_bcc( a::Float64 )
    v1 = 0.5*a*[1,1,1]
    v2 = 0.5*a*[-1,1,1]
    v3 = 0.5*a*[-1,-1,1]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

# more symmetric axis:
function gen_lattice_bcc_v2( a::Float64 )
    v1 = 0.5*a*[-1,1,1]
    v2 = 0.5*a*[1,-1,1]
    v3 = 0.5*a*[1,1,-1]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

# also for trigonal P
function gen_lattice_hexagonal( a::Float64, c::Float64 )
    v1 = a*[1,0,0]
    v2 = a*[-1/2,sqrt(3)/2,0]
    v3 = [0,0,c]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

# 5          Trigonal R, 3fold axis c        celldm(4)=cos(gamma)
# The crystallographic vectors form a three-fold star around
# the z-axis, the primitive cell is a simple rhombohedron:
# v1 = a(tx,-ty,tz),   v2 = a(0,2ty,tz),   v3 = a(-tx,-ty,tz)
# where c=cos(gamma) is the cosine of the angle gamma between
# any pair of crystallographic vectors, tx, ty, tz are:
#   tx=sqrt((1-c)/2), ty=sqrt((1-c)/6), tz=sqrt((1+2c)/3)
function gen_lattice_trigonal( a::Float64, gamma_degree::Float64 )
    c = cos( gamma_degree*pi/180 )
    tx = sqrt((1-c)/2)
    ty = sqrt((1-c)/6)
    tz = sqrt((1+2c)/3)
    #
    v1 = a*[tx,-ty,tz]
    v2 = a*[0,2*ty,tz]
    v3 = a*[-tx,-ty,tz]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

#   -5          Trigonal R, 3fold axis <111>    celldm(4)=cos(gamma)
# The crystallographic vectors form a three-fold star around
# <111>. Defining a' = a/sqrt(3) :
# v1 = a' (u,v,v),   v2 = a' (v,u,v),   v3 = a' (v,v,u)
# where u and v are defined as
#   u = tz - 2*sqrt(2)*ty,  v = tz + sqrt(2)*ty
# and tx, ty, tz as for case ibrav=5
# Note: if you prefer x,y,z as axis in the cubic limit,
#       set  u = tz + 2*sqrt(2)*ty,  v = tz - sqrt(2)*ty
#       See also the note in Modules/latgen.f90
function gen_lattice_trigonal_v2( a::Float64, gamma_degree::Float64 )
    c = cos( gamma_degree*pi/180 )
    tx = sqrt((1-c)/2)
    ty = sqrt((1-c)/6)
    tz = sqrt((1+2c)/3)
    u = tz - 2*sqrt(2)*ty
    v = tz + sqrt(2)*ty
    ap = a/sqrt(3)
    #
    v1 = ap*[u,v,v]
    v2 = ap*[v,u,v]
    v3 = ap*[v,v,u]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

#Tetragonal P (st)               celldm(3)=c/a
#v1 = a(1,0,0),  v2 = a(0,1,0),  v3 = a(0,0,c/a)
function gen_lattice_tetragonal_P( a::Float64, c::Float64)
    v1 = a*[1,0,0]
    v2 = a*[0,1,0]
    v3 = [0,0,c]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

#7          Tetragonal I (bct)              celldm(3)=c/a
#v1=(a/2)(1,-1,c/a),  v2=(a/2)(1,1,c/a),  v3=(a/2)(-1,-1,c/a)
function gen_lattice_tetragonal_I( a::Float64, c::Float64 )
    v1 = 0.5*[a,-a,c]
    v2 = 0.5*[a,a,c]
    v3 = 0.5*[-a,-a,c]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

# 8          Orthorhombic P                  celldm(2)=b/a
#                                              celldm(3)=c/a
#       v1 = (a,0,0),  v2 = (0,b,0), v3 = (0,0,c)
function gen_lattice_orthorhombic_P( a::Float64, b::Float64, c::Float64 )
    v1 = [a,0,0]
    v2 = [0,b,0]
    v3 = [0,0,c]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

# 12          Monoclinic P, unique axis c     celldm(2)=b/a
#                                              celldm(3)=c/a,
#                                              celldm(4)=cos(ab)
#       v1=(a,0,0), v2=(b*cos(gamma),b*sin(gamma),0),  v3 = (0,0,c)
#       where gamma is the angle between axis a and b.
function gen_lattice_monoclinic_P( a::Float64, b::Float64, c::Float64, gamma_degree::Float64 )
    gamma = gamma_degree*pi/180
    v1 = [a,0,0]
    v2 = [b*cos(gamma), b*sin(gamma), 0]
    v3 = [0,0,c]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end

# 14          Triclinic                       celldm(2)= b/a,
#                                             celldm(3)= c/a,
#                                             celldm(4)= cos(bc),
#                                             celldm(5)= cos(ac),
#                                             celldm(6)= cos(ab)
#       v1 = (a, 0, 0),
#       v2 = (b*cos(gamma), b*sin(gamma), 0)
#       v3 = (c*cos(beta),  c*(cos(alpha)-cos(beta)cos(gamma))/sin(gamma),
#            c*sqrt( 1 + 2*cos(alpha)cos(beta)cos(gamma)
#                      - cos(alpha)^2-cos(beta)^2-cos(gamma)^2 )/sin(gamma) )
#       where alpha is the angle between axis b and c
#              beta is the angle between axis a and c
#             gamma is the angle between axis a and b

function gen_lattice_triclinic(a::Float64, b::Float64, c::Float64,
                               alpha_degree::Float64, beta_degree::Float64, gamma_degree::Float64)

    if alpha_degree + beta_degree + gamma_degree <= 180.0
        println("ERROR in gen_lattice_triclinic")
        println("sum of angles must be larger than 180°")
        exit()
    end

    alpha = alpha_degree*pi/180
    beta = beta_degree*pi/180
    gamma = gamma_degree*pi/180

    #
    v1 = [a, 0, 0]
    v2 = [b*cos(gamma), b*sin(gamma), 0]
    t1 = c*(cos(alpha)-cos(beta)cos(gamma))/sin(gamma)
    t2 = c*sqrt( 1 + 2*cos(alpha)cos(beta)cos(gamma) -
                 cos(alpha)^2-cos(beta)^2-cos(gamma)^2 )/sin(gamma)
    v3 = [c*cos(beta), t1, t2]
    #
    LL = zeros(3,3)
    LL[:,1] = v1
    LL[:,2] = v2
    LL[:,3] = v3
    return LL
end


text:
####################################################################################################
function normalize(x)
	out = copyto!(similar(x), x)
	rmul!(out, norm(x))
	return out
end
####################################################################################################
"""
This function checks whether the solution with eigenvalues `eigvalues` is stable and also compute the number of unstable eigenvalues with nonzero imaginary part
"""
function is_stable(contparams::ContinuationPar, eigvalues)::Tuple{Bool, Int64, Int64}
	# the return type definition above is to remove type instability in continuation
	# numerical precision for deciding if an eigenvalue is above a threshold
	precision = contparams.precisionStability

	# update number of unstable eigenvalues
	n_unstable = mapreduce(x -> real(x) > precision, +, eigvalues)

	# update number of unstable eigenvalues with nonzero imaginary part
	n_imag = mapreduce(x -> (abs(round(imag(x), digits = 15)) > precision) * (round(real(x), digits = 15) > precision), +, eigvalues)

	return n_unstable == 0, n_unstable, n_imag
end
####################################################################################################
interval(a, b) = (min(a, b), max(a, b))

# Test function for Fold bifurcation
@inline detectFold(p1, p2, p3) = (p3 - p2) * (p2 - p1) < 0

function locateFold!(contparams::ContinuationPar, contres::ContResult, z, tau, normC, printsolution, verbosity)
	branch = contres.branch
	# Fold point detection based on continuation parameter monotony
	if contparams.detectFold && size(branch)[2] > 2 && detectFold(branch[1, end-2], branch[1, end-1], branch[1, end])
		(verbosity > 0) && printstyled(color=:red, "!! Fold bifurcation point in $(interval(branch[1, end-1],branch[1, end])) \n")
		push!(contres.foldpoint, (
				type = :fold,
				idx = length(branch)-1,
				param = branch[1, end-1],
				norm = normC(z.u),
				printsol = branch[2, end-1], #printsolution(z.u, z.p),
				x = _copy(z.u), tau = normalize(tau.u),
				ind_bif = 0,
				step = length(branch)-1,
				status = :guess,
				δ = (0, 0)))
		detected = true
	end
end

locateFold!(contres::ContResult, iter::PALCIterable, state::PALCStateVariables) = locateFold!(iter.contParams, contres, solution(state), state.tau, iter.normC, iter.printSolution, iter.verbosity)

####################################################################################################
"""
Function for coarse detection of bifurcation points.
"""
function getBifurcationType(contparams::ContinuationPar{T,S,E}, state::PALCStateVariables, normC, printsolution, verbosity, status::Symbol) where {T,S,E}

	detected = false

	# get current number of unstable eigenvalues and
	# unstable eigenvalues with nonzero imaginary part
	n_unstable = state.n_unstable[1]
	n_imag = state.n_imag[1]

	#previous values
	n_unstable_prev = state.n_unstable[2]
	n_imag_prev = state.n_imag[2]

	# computation of the index of the bifurcating eigenvalue
	ind_bif = n_unstable
	if n_unstable < n_unstable_prev; ind_bif += 1; end

	tp = :none

	δn_unstable = abs(n_unstable - n_unstable_prev)
	δn_imag		= abs(n_imag - n_imag_prev)

	# codim 1 bifurcation point detection based on eigenvalue distribution

	if δn_unstable == 1
	# In this case, only a single eigenvalue crossed the imaginary axis
	# Either it is a Branch Point
		if δn_imag == 0
			tp = :bp
		elseif δn_imag == 1
	# Either it is a Hopf bifurcation for a Complex valued system
			tp = contparams.newtonOptions.eigsolver isa AbstractFloquetSolver ? :pd : :hopf
		else # I dont know what bifurcation this is
			tp = :bp
		end
		detected = true
	elseif δn_unstable == 2
		if δn_imag == 2
			tp = contparams.newtonOptions.eigsolver isa AbstractFloquetSolver ? :ns : :hopf
		else
			tp = :bp
		end
		detected = true
	elseif δn_unstable > 2
		tp = :nd
		detected = true
	end

	if δn_unstable < δn_imag
		@error "Error in eigenvalues computation. It seems an eigenvalue is missing, probably conj(λ) for some already computed eigenvalue λ. This makes the identification (but not the detection) of bifurcation points erroneous. You should increase the number of requested eigenvalues."
		tp = :nd
		detected = true
	end

	# rule out initial condition where we populate n_unstable = (-1,-1) and n_imag = (-1,-1)
	if prod(state.n_unstable) < 0 || prod(state.n_imag) < 0
		tp = :nd
		detected = true
	end

	if detected
		# record information about the bifurcation point
		param_bif = (
			type = tp,
			idx = state.step+1,			# this is the index in br.eig
			param = getp(state),
			norm = normC(getx(state)),
			printsol = printsolution(getx(state), getp(state)),
			x = _copy(getx(state)),
			tau = normalize(state.tau.u),
			ind_bif = ind_bif,
			step = state.step,
			status = status,
			δ = (n_unstable - n_unstable_prev, n_imag - n_imag_prev))
		(verbosity>0) && printstyled(color=:red, "!! $(tp) Bifurcation point around p ≈ $(getp(state)), δn_unstable = $δn_unstable, δn_imag = $δn_imag \n")
	end

	return detected, param_bif
end

function closesttozero(ev)
	I = sortperm(abs.(real.(ev)))
	return ev[I]
end

"""
Function to locate precisely bifurcation points using a bisection algorithm.
"""
function locateBifurcation!(iter::PALCIterable, _state::PALCStateVariables, verbose::Bool = true)
	@assert detectBifucation(_state) "No bifucation detected for the state"
	verbose && println("----> Entering [Loc-Bif], state.n_unstable = ", _state.n_unstable)

	# number of unstable eigenvalues
	(n2, n1) = _state.n_unstable
	if n1 == -1 || n2 == -1 return :none end

	# we create a new state for stepping through the continuation routine
	state = copy(_state)

	# iter = @set iter.contParams.newtonOptions.verbose = false

	verbose && println("----> [Loc-Bif] initial ds = ", _state.ds)

	# the bifurcation point is before the current state
	# so we want to first iterate backward with half step size
	# we turn off stepsizecontrol because it would not make a
	# bisection otherwise
	state.ds /= -1
	state.step = 0
	state.stepsizecontrol = false

	next = (state, state)

	if abs(state.ds) < iter.contParams.dsmin; return :none; end

	# record sequence of unstable eigenvalue number
	nunstbls = [n2]
	nimags   = [state.n_imag[1]]

	verbose && println("----> [Loc-Bif] state.ds = ", state.ds)

	# we put this to be able to reference it at the end of this function
	# we don't know its type yet
	eiginfo = nothing

	# we compute the number of changes in n_unstable
	n_inversion = 0
	status = :guess

	biflocated = false

	# emulate a do-while
	while true

		if state.isconverged == false
			@error "----> Newton failed when locating bifurcation!"
			break
		 end

		# we get the current state
		(i, state) = next

		eiginfo, _, n_unstable, n_imag = computeEigenvalues(iter, state)
		updatestability!(state, n_unstable, n_imag)
		push!(nunstbls, n_unstable)
		push!(nimags, n_imag)

		if nunstbls[end] == nunstbls[end-1]
			# bifurcation point still after current state, keep going
			state.ds /= 2
		else
			# we passed the bifurcation point, reverse continuation
			state.ds /= -2
			n_inversion += 1
		end

		verbose &&	printstyled(color=:blue, "----> $(state.step) - [Loc-Bif] (n1, nc, n2) = ",(n1, nunstbls[end], n2), ", ds = $(state.ds), p = ", getp(state), ", #reverse = ", n_inversion,"\n Eigenvalues:\n")
		verbose && Base.display(closesttozero(eiginfo[1])[1:min(5, length(getx(state)))])

		biflocated = abs(real.(closesttozero(eiginfo[1]))[1]) < iter.contParams.tolBisectionEigenvalue

		!(next !== nothing &&
				abs(state.ds) >= iter.contParams.dsminBisection &&
				state.step < iter.contParams.maxBisectionSteps &&
				n_inversion <= iter.contParams.nInversion &&
				biflocated == false) && break

		next = iterate(iter, state; _verbosity = 0)
	end
	verbose && printstyled(color=:red, "----> Found at p = ", getp(state), ", δn = ", abs(2nunstbls[end]-n1-n2),", δim = ",abs(2nimags[end]-sum(state.n_imag))," from p = ",getp(_state),"\n")

	# update current state
	if isodd(n_inversion) || biflocated
		status = :converged
		copyto!(_state.z_pred, state.z_pred)
		copyto!(_state.z_old,  state.z_old)
		copyto!(_state.tau, state.tau)

		_state.eigvals = eiginfo[1]
		if iter.contParams.saveEigenvectors
			_state.eigvecs = eiginfo[2]
		end

		# to prevent bifurcation detection, update the following numbers carefully
		_state.n_unstable = (state.n_unstable[2], _state.n_unstable[2])
		_state.n_imag = (state.n_imag[2], _state.n_imag[2])
	else
		@warn "Bisection failed to locate bifurcation point precisely around p = $(getp(_state)). Fall back to original guess for the bifurcation point. Number of Bisections = $n_inversion"
	end
	verbose && println("----> Leaving [Loc-Bif]")
	return status
end