text:
export ambient_space, rank, gram_matrix, inner_product, involution,
       islocal_square, isequivalent, isrationally_equivalent, quadratic_space,
       hermitian_space, diagonal, invariants, hasse_invariant, witt_invariant, orthogonal_basis, fixed_field

################################################################################
#
#  Basic invariants
#
################################################################################

@doc Markdown.doc"""
    isquadratic(V::AbsSpace) -> Bool

Returns whether $V$ is quadratic.
"""
isquadratic(V::AbsSpace)

@doc Markdown.doc"""
    ishermitian(V::AbsSpace) -> Bool

Returns whether $V$ is hermitian.
"""
ishermitian(V::AbsSpace)

@doc Markdown.doc"""
    rank(V::AbsSpace) -> Int

Return the rank of the space `V`.
"""
rank(L::AbsSpace) = nrows(L.gram)

@doc Markdown.doc"""
    dim(V::AbsSpace) -> Int

Return the dimension of the space `V`.
"""
dim(V::AbsSpace) = rank(V)

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace) -> MatElem

Return the Gram matrix of `V`.
"""
gram_matrix(V::AbsSpace) = V.gram

# Once we have quaternion spaces the following makes more sense

@doc Markdown.doc"""
    base_ring(V::AbsSpace) -> NumField

Return the base field of `V`.
"""
base_ring(V::AbsSpace) = _base_algebra(V)

@doc Markdown.doc"""
    fixed_field(V::AbsSpace) -> NumField

Return the fixed field of `V`.
"""
fixed_field(::AbsSpace)

@doc Markdown.doc"""
    involution(V::AbsSpace) -> NumField

Return the involution of `V`.
"""
involution(V::AbsSpace)

# Maybe cache this

@doc Markdown.doc"""
    isregular(V::AbsSpace) -> Bool

Return whether `V` is regular, that is, if the Gram matrix
has full rank.
"""
function isregular(V::AbsSpace)
  G = gram_matrix(V)
  return rank(G) == nrows(G)
end

################################################################################
#
#  Determinant and discriminant
#
################################################################################

function det(V::AbsSpace)
  d = det(gram_matrix(V))
  return fixed_field(V)(d)
end

@doc Markdown.doc"""
    det(V::AbsSpace) -> FieldElem

Returns the determinant of the space `V` as an element of the fixed field.
"""
det(::AbsSpace)

@doc Markdown.doc"""
    discriminant(V::AbsSpace) -> FieldElem

Returns the discriminant of the space `V` as an element of the fixed field.
"""
function discriminant(V::AbsSpace)
  d = det(V)
  n = mod(rank(V), 4)
  if n == 0 || n == 1
    return d
  else
    return -d
  end
end

function _discriminant(G)
  d = det(G)
  n = nrows(G)
  if mod(n, 4) == 0 || mod(n, 4) == 1
    return d
  else
    return -d
  end
end

################################################################################
#
#  Computing inner products
#
################################################################################

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace, M::MatElem) -> MatElem

Returns the gram matrix of the rows of `M`.
"""
function gram_matrix(V::AbsSpace{T}, M::MatElem{S}) where {S, T}
  if S === elem_type(T)
    return M * gram_matrix(V) * transpose(_map(M, involution(V)))
  else
    Mc = change_base_ring(base_ring(V), M)
    return Mc * gram_matrix(V) * transpose(_map(Mc, involution(V)))
  end
end

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace, S::Vector{Vector}) -> MatElem

Returns the gram matrix of the sequence `S`.
"""
function gram_matrix(V::AbsSpace{T}, S::Vector{Vector{U}}) where {T, U}
  m = zero_matrix(base_ring(V), length(S), rank(V))
  for i in 1:length(S)
    if length(S[i]) != rank(V)
      throw(error("Vectors must be of length $(rank(V))"))
    end
    for j in 1:rank(V)
      m[i, j] = S[i][j]
    end
  end
  return gram_matrix(V, m)
end

@doc Markdown.doc"""
    inner_product(V::AbsSpace, v::Vector, w::Vector) -> FieldElem

Returns the inner product of `v` and `w` with respect to `V`.
"""
inner_product(V::AbsSpace, v::Vector, w::Vector)

################################################################################
#
#  Diagonalization
#
################################################################################

@doc Markdown.doc"""
    orthogonal_basis(V::AbsSpace) -> MatElem

Returns a matrix `M`, such that the rows of `M` form an orthgonal basis of `V`.
"""
function orthogonal_basis(V::AbsSpace)
  _, B = _gram_schmidt(gram_matrix(V), involution(V))
  return B
end

@doc Markdown.doc"""
    diagonal(V::AbsSpace) -> Vector{FieldElem}

Returns a vector of elements $a_1,\dotsc,a_n$ such that `V` is isometric to
the diagonal space $\langle a_1,\dotsc,a_n \rangle$.

The elements will be contained in the fixed field of `V`.
"""
diagonal(V::AbsSpace)

################################################################################
#
#  Gram-Schmidt
#
################################################################################

# Clean this up
function _gram_schmidt(M::MatElem, a, nondeg = true)
  F = deepcopy(M)
  K = base_ring(F)
  n = nrows(F)
  S = identity_matrix(K, n)
  okk = isdiagonal(F)
  if !okk
    for i in 1:n
      if iszero(F[i,i])
        T = identity_matrix(K, n)
        ok = 0
        for j in (i + 1):n
          if !iszero(F[j, j])
            ok = j
            break
          end
        end
        #ok = findfirst(j -> !iszero(F[j, j]), (i + 1):n)
        if ok != 0 # ok !== nothing
          #j = ok + i # findfirst gives the index
          j = ok
          T[i,i] = 0
          T[j,j] = 0
          T[i,j] = 1
          T[j,i] = 1
        else
          ok = 0
          for j in (i + 1):n
            if !iszero(F[i, j])
              ok = j
              break
            end
          end
          if ok == 0
            if nondeg
              error("Matrix is not of full rank")
            end
          else
            j = ok
            T[i, j] = 1 // (2 * F[j, i])
          end
        end
        S = T * S
        F = T * F * transpose(_map(T, a))
      end
      T = identity_matrix(K, n)
      for j in (i + 1):n
        T[j, i] = divexact(-F[j, i], F[i, i])
      end
      F = T * F * transpose(_map(T, a))
      S = T * S
    end
    @assert isdiagonal(F)
  end
  return F, S
end

################################################################################
#
#  Local equivalence
#
################################################################################

@doc Markdown.doc"""
    isequivalent(L::AbsSpace, M::AbsSpace, p::Union{InfPlc, NfOrdIdl}) -> Bool

Returns whether `L` and `M` are equivalent over the completion at `p`.
"""
isequivalent(L::AbsSpace, M::AbsSpace, p)

################################################################################
#
#  Definitness
#
################################################################################

# Returns 0 if V is not definite
# Returns an element a != 0 such that a * canonical_basis of V has
# positive Gram matrix
function _isdefinite(V::AbsSpace)
  K = fixed_field(V)
  R = maximal_order(K)
  if !istotally_real(K) || (ishermitian(V) && !istotally_complex(K))
    return zero(K)
  end
  D = diagonal(V)
  signs_to_consider = Tuple{InfPlc, Int}[]
  for v in real_places(K)
    S = Int[sign(d, v) for d in D]
    if length(unique(S)) != 1
      return zero(K)
    else
      push!(signs_to_consider, (v, S[1]))
    end
  end
  if length(signs_to_consider) == 1
    return K(signs_to_consider[1][2])
  else
    return element_with_signs(K, signs_to_consider)
  end
end

function ispositive_definite(V::AbsSpace)
  E = base_ring(V)
  if isquadratic(V)
    K = E
  else
    K = base_field(E)
  end
  if !istotally_real(K)
    return false
  end
  D = diagonal(V)
  for d in D
    if !istotally_positive(d)
      return false
    end
  end
  return true
end

function isnegative_definite(V::AbsSpace)
  E = base_ring(V)
  if isquadratic(V)
    K = E
  else
    K = base_field(E)
  end

  if !istotally_real(K)
    return false
  end
  D = diagonal(V)
  for d in D
    if !istotally_positive(-d)
      return false
    end
  end
  return true
end

function isdefinite(V::AbsSpace)
  return ispositive_definite(V) || isnegative_definite(V)
end

################################################################################
#
#  Isotropic
#
################################################################################

isisotropic(V::AbsSpace, p::InfPlc) = _isisotropic(V, p)

function _isisotropic(V::AbsSpace, p::InfPlc)
  n = rank(V)
  d = det(V)
  E = base_ring(V)
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif iscomplex(p)
    return true
  else
    D = diagonal(V)
    return length(unique!(Int[sign(d, p) for d in D])) == 2
  end
end

################################################################################
#
#  Helper functions (sort them later)
#
################################################################################

function image(f::NfRelToNfRelMor{T, T}, I::NfRelOrdIdl{T, S}) where {T, S}
  #f has to be an automorphism!!!!
  O = order(I)
  @assert ismaximal(O) # Otherwise the order might change
  K = nf(O)

  B = absolute_basis(I)

  if I.is_prime == 1
    lp = prime_decomposition(O, minimum(I))
    for (Q, e) in lp
      if I.splitting_type[2] == e
        if all(b -> f(b) in Q, B)
          return Q
        end
      end
    end
  end

  pb = pseudo_basis(I)
  pm = basis_pmatrix(I)

  m = zero(matrix(pm))

  c = coefficient_ideals(pm)

  for i in 1:length(pb)
    cc = coordinates(O(f(pb[i][1])))
    for j in 1:length(cc)
      m[i, j] = cc[j]
    end
  end

  J = ideal(O, pseudo_matrix(m, c))

  if isdefined(I, :minimum)
    J.minimum = I.minimum
  end

  J.has_norm = I.has_norm

  if isdefined(I, :norm)
    J.norm = I.norm
  end

  if isdefined(I, :is_prime)
    J.is_prime = I.is_prime
  end

  if isdefined(I, :splitting_type)
    J.splitting_type = I.splitting_type
  end

  return J
end

function image(f::NfRelToNfRelMor{T, T}, I::NfRelOrdFracIdl{T, S}) where {T, S}
  #S has to be an automorphism!!!!
  O = order(I)
  @assert ismaximal(O) # Otherwise the order might change
  K = nf(O)

  pb = pseudo_basis(I)

  z = sum(b * (f(a) * O) for (a, b) in pb)
  return z
end

# An element is locally a square if and only if the quadratic defect is 0, that is
# the valuation is inf.
# (see O'Meara, Introduction to quadratic forms, 3rd edition, p. 160)
function islocal_square(a, p)
  return quadratic_defect(a, p) isa PosInf
end

function _map(a::AbstractAlgebra.MatrixElem, f)
  z = similar(a)
  for i in 1:nrows(a)
    for j in 1:ncols(a)
      z[i, j] = f(a[i, j])
    end
  end
  return z
end

# I think I need a can_change_base_ring version

function element_with_signs(K, D::Dict{InfPlc, Int})
  return _element_with_signs(K, D)
end

function _element_with_signs(K, D)
  OK = maximal_order(K)
  G, mG = infinite_primes_map(OK, real_places(K), 1*OK)
  r = real_places(K)
  z = id(G)
  for (v, s) in D
    for i in 1:length(r)
      if s == 1
        ss = 0
      else
        ss = 1
      end

      if v == r[i]
        z = z + ss * G[i]
      end
    end
  end
  zz = elem_in_nf(mG(z))::elem_type(K)
  @assert all(u -> sign(zz, u[1]) == u[2], D)
  return zz
end

function element_with_signs(K, P::Vector{InfPlc}, S::Vector{Int})
  return _element_with_signs(K, zip(P, S))
end

function element_with_signs(K, A::Vector{Tuple{InfPlc, Int}})
  return _element_with_signs(K, A)
end

function prime_ideals_up_to(O::NfRelOrd, n::Union{Int, fmpz})
  p = 2
  z = ideal_type(O)[]
  while p <= n
    lp = prime_decomposition(base_ring(O), p)
    for q in lp
      if norm(q[1]) > n
        continue
      else
        lq = prime_decomposition(O, q[1])
        for Q in lq
          if absolute_norm(Q[1]) <= n
            push!(z, Q[1])
          end
        end
      end
    end
    p = next_prime(p)
  end
  return sort!(z, by = a -> absolute_norm(a))
end

################################################################################
#
#  Restriction of scalars
#
################################################################################

function restrict_scalars(V::AbsSpace, K::Field = FlintQQ, alpha = one(base_ring(V)))
  E = base_ring(V)
  n = rank(V)
  d = absolute_degree(E)
  B = absolute_basis(E)
  v = elem_type(E)[zero(E) for i in 1:n]
  w = elem_type(E)[zero(E) for i in 1:n]
  G = zero_matrix(FlintQQ, d * n, d * n)
  r = 1
  c = 1
  indices = Vector{Tuple{Int, Int}}(undef, n * d)
  for i in 1:n
    for bi in 1:length(B)
      v[i] = B[bi]
      c = 1
      for j in 1:n
        for bj in 1:length(B)
          w[j] = B[bj]
          G[r, c] = trace(alpha * inner_product(V, v, w), FlintQQ)
          w[j] = zero(E)
          c = c + 1
        end
      end
      v[i] = zero(E)
      indices[r] = (i, bi)
      r = r + 1
    end
  end

  VabstoV = function(v)
    @assert length(v) == d * n
    z = Vector{elem_type(E)}(undef, n)
    for i in 1:n
      z[i] = zero(E)
    end
    for k in 1:(d * n)
      (i, j) = indices[k]
      z[i] = z[i] + v[k] * B[j]
    end
    return z
  end

  VtoVabs = function(w)
    z = Vector{fmpq}(undef, d * n)
    k = 1
    for i in 1:n
      y = w[i]
      for j in 1:d
        z[k] = absolute_coeff(y, j - 1)
        k = k + 1
      end
    end
    return z
  end

  return quadratic_space(FlintQQ, G), VabstoV, VtoVabs
end

# Careful, starts at 0!
function absolute_coeff(z::nf_elem, i)
  return coeff(z, i)
end

function absolute_coeff(z::NumFieldElem, i)
  d = absolute_degree(base_field(parent(z)))
  rowindex = fld(i, d)
  colindex = (i % d)
  return absolute_coeff(coeff(z, rowindex), colindex)
end

function absolute_basis(K::NumField)
  k = base_field(K)
  kabs = absolute_basis(k)
  res = elem_type(K)[]
  for b in basis(K)
    for bb in kabs
      push!(res, bb * b)
    end
  end
  return res
end

function absolute_basis(K::NumField{fmpq})
  return basis(K)
end

#

istotally_real(::FlintRationalField) = true

istotally_positive(x::fmpq) = x > 0

# This function is really slow...
function denominator(M::fmpq_mat)
  d = one(FlintZZ)
  for i in 1:nrows(M)
    for j in 1:ncols(M)
      d = lcm!(d, d, denominator(M[i, j]))
    end
  end
  return d
end

function _weak_approximation_coprime(IP, S, M)
  R = order(M)
  A, _exp, _log = infinite_primes_map(R, IP, M)

  t = (1 + _exp(A([ S[j] == 1 ? 0 : -1 for j in 1:length(IP)])))
  @assert all(i -> sign(t, IP[i]) == S[i], 1:length(IP))
  return t
end


text:
function init_ped_grid(env::UrbanEnv, pos_res::Float64, vel_res::Float64)
    d = Dict{LaneTag, RectangleGrid{3}}()
    phi_space = SVector{2, Float64}(0., float(pi))
    v_space = get_ped_vspace(env, vel_res)
    for lane in env.crosswalks
        l_space = get_discretized_lane(lane.tag, env.roadway, pos_res)
        grid = RectangleGrid(l_space, v_space,  phi_space)
        d[lane.tag] = grid
    end
    return d
end

function init_car_grid(env::UrbanEnv, pos_res::Float64, vel_res::Float64)
    d = Dict{LaneTag, RectangleGrid{2}}()
    v_space = get_car_vspace(env, vel_res)
    for i=1:length(env.roadway.segments)
        for lane in env.roadway.segments[i].lanes
            l_space = get_discretized_lane(lane.tag, env.roadway, pos_res)
            grid = RectangleGrid(l_space, v_space)
            d[lane.tag] = grid
        end
    end
    return d
end
function init_l_grid(env::UrbanEnv, pos_res::Float64)
    d = Dict{LaneTag, RectangleGrid{1}}()
    for i=1:length(env.roadway.segments)
        for lane in env.roadway.segments[i].lanes
            l_space = get_discretized_lane(lane.tag, env.roadway, pos_res)
            grid = RectangleGrid(l_space)
            d[lane.tag] = grid
        end
    end
    return d
end

function init_v_grid(env::UrbanEnv, vel_res::Float64)
    v_space = get_car_vspace(env, vel_res)
    return RectangleGrid(v_space)
end

# a bunch of interpolation helpers
function AutomotivePOMDPs.interpolate_state(mdp::PedCarMDP, state::VehicleState)
    # interpolate longitudinal position and velocity
    if state == mdp.off_grid
        return VehicleState[state], Float64[1.0]
    end
    lane = get_lane(mdp.env.roadway, state)
    l_space = get_discretized_lane(lane.tag, mdp.env.roadway, mdp.pos_res)
    grid = mdp._car_grid[lane.tag]
    real_state = SVector{2, Float64}(state.posF.s, state.v)
    idx, weights = interpolants(grid, real_state)
    n_pts = length(idx)
    states = Vector{VehicleState}(undef, n_pts)
    probs = zeros(n_pts)
    for i=1:n_pts
        sg, vg = ind2x(grid, idx[i])
        states[i] = VehicleState(Frenet(lane, sg), mdp.env.roadway, vg)
        probs[i] = weights[i]
    end
    return states, probs
end

# take into account heading as well
function AutomotivePOMDPs.interpolate_pedestrian(mdp::PedCarMDP, state::VehicleState)
    # interpolate longitudinal position and velocity
    if state == mdp.off_grid
        return (state,), (1.0,)
    end
    lane = get_lane(mdp.env.ped_roadway, state)
    grid = mdp._ped_grid[lane.tag]
    real_state = SVector{3, Float64}(state.posF.s, state.v, abs(state.posF.ϕ))
    idx, weights = interpolants(grid, real_state)
    n_pts = length(idx)
    states = Vector{VehicleState}(undef, n_pts)
    probs = zeros(n_pts)
    for i=1:n_pts
        sg, vg, phig = ind2x(grid, idx[i])
        states[i] = VehicleState(Frenet(lane, sg, 0., phig), mdp.env.ped_roadway, vg)
        probs[i] = weights[i]
    end
    return states, probs
end

function AutomotivePOMDPs.interpolate_state(mdp::PedCarMDP, s::PedCarMDPState)
    # interpolate s in the grid
    vspace = get_car_vspace(mdp.env, mdp.vel_res)
    v_ped_space = get_ped_vspace(mdp.env, mdp.vel_ped_res)
    itp_ped, itp_ped_w = interpolate_pedestrian(mdp, s.ped)
    itp_car, itp_car_w = interpolate_state(mdp, s.car)
    itp_ego, itp_ego_w = interpolate_state(mdp, s.ego)
    itp_states = Vector{PedCarMDPState}(undef, length(itp_ego)*length(itp_car)*length(itp_ped))
    itp_w = Vector{Float64}(undef, length(itp_states))
    l = 1
    for (i, p) in enumerate(itp_ped)
        for (j, c) in enumerate(itp_car)
            for (k, e) in enumerate(itp_ego)
                itp_states[l] = PedCarMDPState(crash(mdp, s), e, p, c, s.route)
                itp_w[l] = itp_ped_w[i]*itp_car_w[j]*itp_ego_w[k]
                l += 1
            end
        end
    end
    @assert sum(itp_w) ≈ 1.
    return itp_states, itp_w
end

function AutomotivePOMDPs.get_mdp_state(mdp::PedCarMDP, pomdp::UrbanPOMDP, s::Scene, ped_id, car_id)
    return AutomotivePOMDPs.get_mdp_state(mdp, pomdp.models, s, ped_id, car_id)
end

function AutomotivePOMDPs.get_mdp_state(mdp::PedCarMDP, models::Dict{Int64, DriverModel}, s::Scene, ped_id, car_id)
    @assert length(s) <= 3
    car_i = findfirst(CAR_ID, s)
    car = Vehicle(mdp.off_grid, mdp.car_type, car_id)
    if car_i != nothing
        car = s[car_i]
    end
    ped_i = findfirst(PED_ID, s)
    ped = Vehicle(mdp.off_grid, mdp.ped_type, ped_id)
    if ped_i != nothing
        ped = s[ped_i]
    end
    ego = get_ego(s)
    # find route
    sroute = nothing
    if haskey(models, car_id) && car_i != nothing
        # find the exact route from the list of routes
        curr_route = [l.tag for l in models[car_id].navigator.route]
        for route in get_car_routes(mdp.env)
            tags = intersect(Set(curr_route), Set(route))
            if length(tags) >= 2
                sroute = SVector{2, LaneTag}(route[1], route[end])
            elseif length(curr_route) == 1 && curr_route[1] ∈ route
                sroute = SVector{2, LaneTag}(route[1], route[end])
            end
        end
        if sroute == nothing
            println(curr_route)
        end
    else
        sroute = OFF_ROUTE
    end
    if isapprox(ped.state.posG.θ, float(pi))
        posG = ped.state.posG
        lane = mdp.env.ped_roadway[LaneTag(19,1)]
        posF = Frenet(proj(posG, lane, mdp.env.roadway, move_along_curves=false), mdp.env.roadway)
        p_state = VehicleState(posG, posF, ped.state.v)
    else
        p_state = VehicleState(ped.state.posG, mdp.env.ped_roadway, ped.state.v)
    end
    e_state = VehicleState(ego.state.posG, car_roadway(mdp.env), ego.state.v)
    c_state = VehicleState(car.state.posG, car_roadway(mdp.env), car.state.v)
    return PedCarMDPState(is_colliding(ego, car), e_state, p_state, c_state, sroute)
end

ext: jl


text:
using Distributed
@testset "drop_repeats" begin
    Signals.async_mode(false)
    A = Signal(1)
    B = Signal(droprepeats(A); state = 0) do a, state
        state.x += 1
        a + 1
    end
    C = Signal(B; state = 0) do b, state
        state.x += 1
        b + 1
    end

    @test B.state.x == 1
    @test C.state.x == 1
    A(1); A(1);
    @test B.state.x == 1
    @test C.state.x == 1
    A(3.14)
    @test B.state.x == 2
    @test C.state.x == 2
end

@testset "when do" begin
    Signals.async_mode(false)
    A = Signal(1)
    cond = Signal(A) do a
        a < 10
    end
    B = when(cond, A) do a
        a*2
    end
    A(5)
    @test B[] == 10
    A(15)
    @test B[] == 10

    A[] = 2
    @test B() == 4
    A[] = 15
    @test B() == 4
end

@testset "sampleon" begin
    Signals.async_mode(false)
    A = Signal(1)
    B = Signal(10)
    C = sampleon(A, B)

    @test C[] == 10
    B(0)
    @test C[] == 10
    A(2)
    @test C[] == 0

    B[] = 200
    @test C() == 0
    A[] = 9
    @test C() == 200
end

@testset "filter" begin
    Signals.async_mode(false)
    A = Signal(100000)
    B = filter(x -> x<10, 1, A)

    @test B[] == 1
    A(100)
    @test B[] == 1
    A(9)
    @test B[] == 9

    A[] = 1
    @test B() == 1
    A[] = 100
    @test B() == 1
    A[] = 9
    @test B() == 9
end

@testset "foldp" begin
    Signals.async_mode(false)
    A = Signal(1)
    B = foldp(+, 0, A)

    A(2)
    A(3)
    @test B[] == 6

    A[] = 6
    @test B() == 12
end

@testset "merge" begin
    Signals.async_mode(false)
    A = Signal(1)
    B = Signal(1)
    C = merge(A, B)

    A(10)
    @test C[] == 10
    B[] = "aa"
    @test C() == "aa"

    A[] = 0; B[] = 100
    @test C() == 100

    A = Signal(1)
    B = Signal(1)
    C = droprepeats(A)
    D = merge(B, C)

    A(100)
    @test D[] == 100
    B[] = "hh"; A[] = 100
    @test D() == "hh"
end

@testset "async_signal" begin
    addprocs(1)
    @everywhere using Signals
    Signals.async_mode(true)

    A = Signal(1)
    B = async_signal(x -> x+1, A)
    C = remote_signal(x -> x+1, A)

    A(10)
    sleep(2)
    @test B() == 11
    sleep(2)
    @test C() == 11
end

@testset "recursion_free" begin
    A = Signal(1)
    C = count(A)

    B = recursion_free(A) do a
        A(a+1)
    end

    Signals.async_mode(true)
    A(10)
    sleep(0.1)
    @test C() == 3
    @test A() == 11

    Signals.async_mode(false)
    A(100)
    @test C() == 4
    @test A() == 101
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

"""
    onthread(f::F, tid::Int; wait::Bool=true) where {F<:Function}

Execute a function `f` on thread `tid`.

To execute `f` on a thread other than 1 can speed it up
significantly if it depends on asynchronous tasks.

# Arguments
- `f::Function`:     function to execute
- `tid::Int`:        thread id
- `wait::Bool=true`: if true, it waits for function to finish

# Example
```jldoctest
julia> using DiscreteEvents, .Threads

julia> onthread(2) do; threadid(); end
2
```
"""
function onthread(f::F, tid::Int; wait::Bool=true) where {F<:Function}
    t = Task(nothing)
    @assert tid in 1:nthreads() "thread $tid not available!"
    @threads for i in 1:nthreads()
        if i == tid
            t = @async f()
        end
    end
    wait && fetch(t)
end

"""
```
pseed!(s::Int)
```
Seed each of the thread local RNGs with `s*threadid()` to get
reproducible, but different random number sequences on each thread.
"""
function pseed!(s::Int)
    @threads for i in 1:nthreads()
        Random.seed!(s*i)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Openresty"
version = v"1.15.8"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://openresty.org/download/openresty-1.15.8.3.tar.gz", "b68cf3aa7878db16771c96d9af9887ce11f3e96a1e5e68755637ecaff75134a8"),
    ArchiveSource("https://ftp.pcre.org/pub/pcre/pcre-8.43.tar.gz", "0b8e7465dc5e98c757cc3650a20a7843ee4c3edf50aaf60bb33fd879690d2c73"),
    ArchiveSource("https://www.openssl.org/source/openssl-1.0.2t.tar.gz", "14cb464efe7ac6b54799b34456bd69558a749a4931ecfd9cf9f71d7881cac7bc"),
    ArchiveSource("https://www.zlib.net/zlib-1.2.11.tar.gz", "c3e5e9fdd5004dcb542feda5ee4f0ff0744628baf8ed2dd5d66f8ca1197cb1a1")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/openresty-*/
export SUPER_VERBOSE=1
./configure --prefix=${prefix} \
    --with-cc=$CC \
    --with-zlib=$WORKSPACE/srcdir/zlib-1.2.11 \
    --with-openssl=$WORKSPACE/srcdir/openssl-1.0.2t \
    --with-pcre=$WORKSPACE/srcdir/pcre-8.43 \
    --with-pcre-jit
make
make install
rm ${bindir}/openresty
ln -s ../nginx/sbin/nginx ${bindir}/openresty
install_license $prefix/COPYRIGHT
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("openresty", :openresty),
    ExecutableProduct("opm", :opm),
    ExecutableProduct("resty", :resty),
    FileProduct("luajit", :luajit_dir),
    FileProduct("lualib", :lualib_dir),
    FileProduct("nginx", :nginx_dir),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
struct KnapsackGenerator <: AbstractModelGenerator
    nb_items::Int
    max_weight::Int
    correlation::Real
end

KnapsackGenerator(nb_items::Int, max_weight::Int) = KnapsackGenerator(nb_items, max_weight, 1)

"""
    fill_with_generator!(cpmodel::CPModel, gen::KnapsackGenerator)::CPModel

Fill the cpmodel with variables and constraints for the knapsack problem.
The weights are uniformly distributed between 1 and `gen.max_weight`, the values are uniformly distributed
between (for each `weight`) `weight - max_weight/(10*correlation)` and `weight + max_weight/(10*correlation)`.
It is possible to give `Inf` as the `gen.correlation` to have a strict equality between the weights and their values.
`gen.correlation` must be strictly positive.
This method is from the following paper:
https://www.researchgate.net/publication/2548374_Core_Problems_in_Knapsack_Algorithms
"""
function fill_with_generator!(cpmodel::CPModel, gen::KnapsackGenerator; seed=nothing)
    if !isnothing(seed)
        Random.seed!(seed)
    end

    correlation = gen.correlation
    max_weight = gen.max_weight
    nb_items = gen.nb_items

    @assert correlation > 0

    # create values, weights and capacity
    weights_distr = DiscreteUniform(1, max_weight)
    weights = rand(weights_distr, nb_items)

    deviation = floor(max_weight/(10*correlation))
    value_distr = truncated.(DiscreteUniform.(weights .- deviation, weights .+ deviation), 1, Inf)
    values = rand.(value_distr)

    c = floor(nb_items * max_weight/2 / 4)
    capacity = rand(DiscreteUniform(c, c*4))

    ### Variables
    x = SeaPearl.IntVar[]
    for i in 1:nb_items
        push!(x, SeaPearl.IntVar(0, 1, "x[" * string(i) * "]", cpmodel.trailer))
        addVariable!(cpmodel, last(x))
    end

    ### Constraints

    # Creating the totalWeight variable
    varsWeight = SeaPearl.AbstractIntVar[]
    maxWeight = 0
    for i in 1:nb_items
        wx_i = SeaPearl.IntVarViewMul(x[i], weights[i], "w["*string(i)*"]*x["*string(i)*"]")
        push!(varsWeight, wx_i)
        maxWeight += weights[i]
    end
    totalWeight = SeaPearl.IntVar(0, maxWeight, "totalWeight", cpmodel.trailer)
    minusTotalWeight = SeaPearl.IntVarViewOpposite(totalWeight, "-totalWeight")
    SeaPearl.addVariable!(cpmodel, totalWeight;branchable=false)
    SeaPearl.addVariable!(cpmodel, minusTotalWeight;branchable=false)
    push!(varsWeight, minusTotalWeight)
    weightEquality = SeaPearl.SumToZero(varsWeight, cpmodel.trailer)
    SeaPearl.addConstraint!(cpmodel, weightEquality)

    # Making sure it is below the capacity
    weightConstraint = SeaPearl.LessOrEqualConstant(totalWeight, capacity, cpmodel.trailer)
    SeaPearl.addConstraint!(cpmodel, weightConstraint)

    ### Objective ### minimize: -sum(v[i]*x_a[i])

    # Creating the sum
    varsValue = SeaPearl.AbstractIntVar[]
    maxValue = 0
    for i in 1:nb_items
        vx_i = SeaPearl.IntVarViewMul(x[i], values[i], "v["*string(i)*"]*x["*string(i)*"]")
        push!(varsValue, vx_i)
        maxValue += values[i]
    end
    totalValue = SeaPearl.IntVar(-maxValue, 0, "totalValue", cpmodel.trailer)
    SeaPearl.addVariable!(cpmodel, totalValue;branchable=false)
    push!(varsValue, totalValue)
    valueEquality = SeaPearl.SumToZero(varsValue, cpmodel.trailer)
    SeaPearl.addConstraint!(cpmodel, valueEquality)

    # Setting it as the objective
    SeaPearl.addObjective!(cpmodel,totalValue)

    nothing
end

feature_length(gen::KnapsackGenerator,  t::Type{DefaultStateRepresentation{F}}) where {F} = 10


text:
using BNPRegressionGGA2021
using CSV
using DataFrames
using Random
const BNP = BNPRegressionGGA2021

# Simulate the data
begin
    include("dgp-normal-normal1.jl");
    Random.seed!(1);
    N0, N1 = 500, 50;
    dy, y0, x0, X0, y1, x1, X1 = simulate_sample(N0, N1);
    m0, o0 = mean(y0), std(y0)
    y0 = (y0 .- m0) ./ o0
end;
# did not work with 1000
# did not work with 2000
# did not work with 3000
# did not work with 4000
# worked better with 5000
# worked better with 6000
# worked with 10000

# Fit the model
begin
    update_γ = false;
    mapping = [[1], collect(2:size(X0, 2))];
    m = BNP.DGSBPNormal(; y0, X0, y1, X1, mapping, update_γ);
    chainf, chainβ, chaing = BNP.sample!(m; mcmcsize = 20000, burnin = 10000);
end;

# Save the results as a CSV file
begin
    df = DataFrame(
        y = y1,
        x = x1,
        fh = mean(chainf),
        f0 = pdf.(dy.(x1), m0 .+ o0 .* y1) .* o0
    )
    filename = "data/dgp-normal-normal-example1.csv"
    CSV.write(filename, df)
end


text:
function Mn_series_coeff!(t::Transit_Struct{T}) where {T <: Real}
# Use series expansion to compute M_n. (don't need its derivative, dM_n/dk.)
# Computing leading coefficient (n=0).  Need to compute for
# n_max-3 to n_max.  Also, need to keep track of k^2 < 1
# and k^2 >= 1.  So, Mn_coeff is 3-dimensional array:
# Mn_coeff[2,4,jmax], where first dimension is for k^2 < 1 (1)
# and k^2 >= 1 (2); second dimension is for n_max-3 (1) to n_max (4),
# while third dimension is for the series coefficients.
# Computing leading coefficient (n=0):
for k2 = 0:1
coeff = zero(T)
# Loop over n_max to n_max-3:
for j=1:4
  # m is n_max to n_max-3 (we need both since downward recursion requires two):
  m = t.n_max+j-4
  mhalf = 0.5*m
  if k2 < 1
    coeff = sqrt(pi)*exp(lgamma(mhalf+1.0)-lgamma(mhalf+1.5))
    # Add leading term to M_n:
    t.Mn_coeff[1,j,1] = coeff
    # Now, compute higher order terms until desired precision is reached:
    @inbounds for i=1:t.jmax-1
      coeff *= (2*i-1)^2/(2*i*(1+m+2*i))
      t.Mn_coeff[1,j,i+1] = coeff
    end
  else # k^2 >= 1
    coeff = convert(T,pi)
    # Store leading terms:
    t.Mn_coeff[2,j,1] = coeff
    # Loop over higher order terms:
    if iseven(m)
      # If m is even, then series truncates:
      jmax = div(m,2)
    else
      jmax = t.jmax-1
    end
    for i = 1:jmax
      coeff *= (2+m-2*i)*(1-2*i)/(4*i^2)
      t.Mn_coeff[2,j,i+1] = coeff
    end
  end
end
end
end

max_stars_repo_name: dfm/limbdark

max_issues_repo_name: dfm/limbdark

max_forks_repo_name: dfm/limbdark


text:
# This file was generated, do not modify it. # hide
r = report(spca).reports[1]
cumsum(r.principalvars ./ r.tvar)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
module Subsystems

using LinearAlgebra
using BlockArrays

abstract type AbstractSubsystem{T} end
abstract type LinearSubsystem{T} <: AbstractSubsystem{T} end

export LdtiSubsystem,
    LinearSubsystem,
    OutputMatrix,
    AbstractSubsystem,
    UIO,
    DLUE,
    UMVE,
    addNeighbour,
    setNeighbour,
    setNeighbourCoupling,
    removeNeighbour,
    getNeighbourIndices,
    getNeighbourStates,
    getOutboundWeights,
    outputMap,
    setF,
    updateState,
    setInitialState

OutputMatrix{T} = Union{typeof(I), AbstractMatrix{T}}

function _string_mat_with_headers(X::AbstractVecOrMat)
    p = (io, m) -> Base.print_matrix(io, m)
    return replace(sprint(p, X), "\"" => " ")
end

include("LdtiSubsystem.jl")
include("Observers/UIO.jl")
include("Observers/DLuenberger.jl")
include("Observers/UMVE.jl")
include("LSS.jl")

export _getdimarrays,
    tomonolith

end


text:
using Statistics
using StatsBase

# Absolute energy of the time series, the sum over the squared values
abs_energy(x::AbstractArray{T} where T<:Real) = mapreduce(x -> x * x, +, x) #using reduce(map()) is faster but with more mem allocation

# Highest absolute value of the time series x
absolute_maximum(x::AbstractArray{T} where T<:Real) = maximum(abs, x)

# Sum over the abs value of consecutive changes in the series x
function absolute_sum_of_changes(x::AbstractArray{T} where T<:Real)
    diff = Real[]

    for i in 1:length(x)-1
        push!(diff,x[i+1] - x[i])
    end

    return sum(map(abs,diff))
end

"""
Descriptive statistics on the autocorrelation of the time series x.
Calculates the value of an aggregation function fₐ over
    the autocorrelation for different lags.
Input:  x -> time series to calculate the feature of
        param -> dictionaries {"f_agg":x, "maxlag":n},
            with x, a string, the name of the aggregation funcion(e.g. mean, variance, ...),
            and n, an integer, the maximal number of lags to consider
"""
#TODO functions that calls function must have parameter of type "Function" not "String"
function agg_autocorrelation(x::AbstractArray{T} where T<:Real, f::Function, maxlag::Int)

    R::AbstractArray = zeros(maxlag)
    μ = mean(x)
    σ²= var(x)

    for i in 1:maxlag
        aux = 0.

        for t in 1:(length(x) - i)
            aux += (x[t] - μ) * (x[t+i] - μ)
        end

        R[i] = ( 1 / ( (length(x) - i) * σ² ) ) * aux
    end

    return f(R)

end

"""
Implements a vectorized approximate entropy algorithm.
Input: x-> time series to calculate the feature of
       m-> length of compared run of data (int)
       r-> filtering level, must be positive (float)
"""

"""
function approximate_entropy(x::AbstractArray{T} where T<:Real, m::Int, r::Real)
    dist(x_i, x_j) =

    n = length(x)

    x_re = Array([[ x[j] for j in i:(i+m-1) ] for i in 1:(n-m+1)])

    # GO AHEAD
end
"""

# augmented_dickey_fuller TODO

# Calculates the autocorrelation of the specified lag.
function autocorrelation(x::AbstractArray{T} where T<:Real, lag::Int)
    n = length(x)
    μ = mean(x)
    σ²= var(x)

    sum = 0.
    for t in 1:(n-lag)
        sum += (x[t] - μ) * (x[t+lag] - μ)
    end

    return (1 / ((n-lag) * σ²)) * sum
end

"""
Returns the correlation from first digit distribution when compared
 to the Newcomb-Benford's Law fistribution
"""
function benford_correlation(x::AbstractArray{T} where T<:Real)
    # Take first digit from data
    x = Array([last(digits(trunc(Int, abs(x[i])))) for i in 1:length(x)])

    # Benford distribution
    P = Array([log10(1 + 1/n) for n in 1:9])

    # Data distribution
    D = Array([count(i -> i==n, x)/length(x) for n in 1:9])

    # Return the correlation between benford distribution and data distribution
    return cor(P,D)
end

"""
First bins the values of x into max_bins equidistant bins.
Then calculates the value of binned entropy
"""
function binned_entropy(x::AbstractArray{T} where T<:Real, max_bins::Int)
    # Create equidistant max_bins
    edges = LinRange(minimum(x), maximum(x), max_bins)
    # Fit function requires column vector
    edges = collect(permutedims(edges))

    # Create histogram with the calculated bins
    h = fit(Histogram, vec(x), vec(edges))
    hist = h.weights

    # Vector pₖ with percentages of samples in bin k
    p = hist / length(x)

    # Where the percentage is 0, set to 1 to avoid Nan values?
    p[p.==0] .= 1

    # Return Entropy = - ∑ᵐᵢ₌₁ pᵢ * log(pᵢ)
    return - reduce(+, map(x->x*log(x), p))
end

# Uses c3 statistics to measure non linearity in the time series
function c3(x::AbstractArray{T} where T<:Real, lag::Int)
    n = length(x)
    if 2*lag ≥ n
        return 0
    end

    sum=0.
    for i in 1:(n - 2*lag)
        sum += x[i + 2*lag] * x[i + lag] * x[i]
    end

    return (1 / (n - 2*lag)) * sum
end

"""
First takes a corridor given by the quatiles ql and qh of
the distribution of x. Then calculates the average, absolute
value of consecutive changes of the series x inside this corridor.
"""
function change_quantiles(x::AbstractArray{T} where T<:Real, ql::Float64, qh::Float64, f_agg::Function)
        # Compute the quantiles
        q = quantile(vec(x), [ql, qh], sorted=true)

        # Corridor between the 2 quantiles
        corr = Real[]
        for i in 1:length(x)
            if x[i] >= q[1] && x[i] <= q[2]
                push!(corr, x[i])
            end
        end

        # Return the aggregation function of corridor
        return f_agg(corr)
end

#This function calculator is an estimate for a time series complexity
function cid_ce(x::AbstractArray{T} where T<:Real, normalize::Bool)
    n = length(x)

    if normalize
        s = std(x)
        if s != 0
            x = (x .- mean(x)) / s
        else
            return 0
        end
    end

    sum = 0.
    for i in 2:n
        sum += (x[i] - x[i-1])^2
    end

    return sqrt(sum)
end

# Returns the percentage of values in x that are higher than t
count_above(x::AbstractArray{T} where T<:Real, t::Float64) = count(i -> i > t, x) / length(x)

# Returns the number of values in x that are higher than the mean of x
count_above_mean(x::AbstractArray{T} where T<:Real) = count(i -> i > mean(x), x)

# Returns the percentage of values in x that are lower than t
count_below(x::AbstractArray{T} where T<:Real, t::Float64) = count(i -> i < t, x) / length(x)

# Returns the number of values in x that are lower than the mean of x
count_below_mean(x::AbstractArray{T} where T<:Real) = count(i -> i < mean(x), x)

"""
Return the first location of the maximum value of x.
The position is calculated relatively to the length of x
"""
first_location_of_maximum(x::AbstractArray{T} where T<:Real) = findmax(x)[2] / length(x)

"""
Return the first location of the minimum value of x.
The position is calculated relatively to the length of x
"""
first_location_of_minimum(x::AbstractArray{T} where T<:Real) = findmin(x)[2] / length(x)

# Checks if any value in x occurs more than once
has_duplicate(x::AbstractArray{T} where T<:Real) = return length(x) != length(unique(x))

# Checks if the maximum value of x is observed more than once
has_duplicate_max(x::AbstractArray{T} where T<:Real) = if count(i -> i == maximum(x), x) > 1 return true else return false end

# Checks if the minimum value of x is observed more than once
has_duplicate_min(x::AbstractArray{T} where T<:Real) = if count(i -> i == minimum(x), x) > 1 return true else return false end

"""
Does time series have large standard deviation?
Return boolean denoting if the standard deviation of x
is higher than 'r' times the range.
"""
large_standard_deviation(x::AbstractArray{T} where T<:Real, r::Float64) = return std(x) > r * (maximum(x) - minimum(x))

# Returns the relative last location of tha maximum value of x
last_location_of_maximum(x::AbstractArray{T} where T<:Real) = findmax(reverse(x))[2] / length(x)

# Returns the relative last location of tha minimum value of x
last_location_of_minimum(x::AbstractArray{T} where T<:Real) = findmin(reverse(x))[2] / length(x)

# Returns the length of the longest consecutive subsequence in x that is bigger than the mean of x
function longest_strike_above_mean(x::AbstractArray{T} where T<:Real)
    aux = x.>mean(x)
    max_strike = 0
    i=1

    while i <= length(aux)
        if aux[i] == 1
            c = 1
            i+=1
            while i <= length(aux)
                if aux[i] == 1
                    c+=1
                else
                    break
                end
                i+=1
            end
            if c > max_strike
                max_strike = c
            end
        else
            i+=1
        end
    end

    return max_strike
end

# Returns the length of the longest consecutive subsequence in x that is lower than the mean of x
function longest_strike_below_mean(x::AbstractArray{T} where T<:Real)
    aux = x.<mean(x)
    max_strike = 0
    i=1

    while i <= length(aux)
        if aux[i] == 1
            c = 1
            i+=1
            while i <= length(aux)
                if aux[i] == 1
                    c+=1
                else
                    break
                end
                i+=1
            end
            if c > max_strike
                max_strike = c
            end
        else
            i+=1
        end
    end

    return max_strike
end

# Average over first differences
function mean_abs_change(x::AbstractArray{T} where T<:Real)
    diff = Real[]

    for i in 1:length(x)-1
        push!(diff,x[i+1] - x[i])
    end

    return mean(map(abs, diff))
end

# Average over time sereis differences
mean_change(x::AbstractArray{T} where T<:Real) = (x[end] - x[1]) / (length(x) - 1)

# Calculates the arithmetic mean of the n absolute maximum values of the time series
function mean_n_absolute_max(x::AbstractArray{T} where T<:Real, n::Int)
    aux = sort(map(abs, x))
    return mean(aux[length(aux) - (n - 1) : end])
end

# Returns the mean value of a central approximation of the secondo derivative
mean_second_derivative_central(x::AbstractArray{T} where T<:Real) = (x[end] - x[end - 1] - x[2] - x[1]) / (2 * (length(x) - 2))

# Calculates the number of crossing of x on m.
function number_crossing_m(x::AbstractArray{T} where T<:Real, m::Float64)
    positive = x .> m
    diff = Int[]

    for i in 1:length(x) - 1
        push!(diff, positive[i+1] - positive[i])
    end

    return count(i -> i != 0, diff)
end

# TODO Gabriele
# # Returns the percentage of non-unique data points.
# function percentage_of_reoccurring_datapoints_to_all_datapoints(x::AbstractArray{T} where T<:Real)
#     value_counts = countmap(x)
#     reoccurring_values = 0
#
#     for item in value_counts
#         if item[2] > 1
#             reoccurring_values += 1
#         end
#     end
#
#     return reoccurring_values / length(x)
# end

# TODO Gabriele
# # Returns the percentage of values that are present in the time series more than once
# function percentage_of_reoccurring_values_to_all_values(x::AbstractArray{T} where T<:Real)
#     value_counts = countmap(x)
#     reoccurring_values = 0
#     unique_values = 0
#
#     for item in value_counts
#         if item[2] > 1
#             reoccurring_values += 1
#         else
#             unique_values += 1
#         end
#     end
#
#     return reoccurring_values / unique_values
# end

# Count observed values within the interval [min, max)
range_count(x::AbstractArray{T} where T<:Real, min::Float64, max::Float64) = count(i -> i>=min && i<max, x)

# Ratio of values that are mroe than r * std(x) away from the mean of x
function ratio_beyond_r_sigma(x::AbstractArray{T} where T<:Real, r::Float64)
    μ = mean(x)
    r_sigma = r * std(x)

    aux = map(abs, map(x -> x - μ, x))

    aux = Array([aux[i] > r_sigma ? aux[i] : 0 for i in 1:length(aux)])

    return sum(aux) / length(x)
end

# TODO Gabriele
# """
# Returns a factow which is 1 if all values in the time series occur only once,
# and below one if this is not the case
# """
# function ratio_value_number_to_time_series_length(x::AbstractArray{T} where T<:Real)
#     value_counts = countmap(x)
#     unique = Real[]
#
#     for item in value_counts
#         if item[2] == 1
#             push!(unique, item[2])
#         end
#     end
#
#     return length(unique) / length(x)
# end

# Returns the root mean square of the time series
root_mean_square(x::AbstractArray{T} where T<:Real) = sqrt(mean(map(x -> x^2, x)))

# Calculate and return sample entropy of x
function sample_entropy(x::AbstractArray{T} where T<:Real)
    m = 2
    tolerance = 0.2 * std(x)

    xmi = Array([x[i : i+m-1] for i in 1:length(x)-m])
    xmj = Array([x[i : i+m-1] for i in 1:length(x)-m+1])

    xmi = maximum(collect(hcat(xmi...)'), dims=2)
    xmj = maximum(collect(hcat(xmj...)'), dims=2)

    B = sum( Array([sum(x -> x <= tolerance, map(abs, xmi_i .- xmj)) - 1 for xmi_i in xmi]))

    m += 1
    xm = Array([x[i:i+m-1] for i in 1:length(x)-m+1])
    xm = maximum(collect(hcat(xm...)'), dims=2)

    A = sum( Array([sum(x -> x <= tolerance, map(abs, xm_i .- xm)) - 1 for xm_i in xm]))

    return -log(A/B)
end

# # Returns the sum of all data points, that are present in the time series more than once
# function sum_of_reoccurring_data_points(x::AbstractArray{T} where T<:Real)
#     value_counts = countmap(x)
#     reoc = Real[]
#
#     for item in value_counts
#         if item[2] > 1
#             for i in 1:item[2]
#                 push!(reoc, item[1])
#             end
#         end
#     end
#
#     return sum(reoc)
# end

function window_slice(x::AbstractArray{T} where T<:Real;
        from::Union{Missing,Int}=missing,
        to::Union{Missing,Int}=missing)
        ismissing(from) && ismissing(to) && return x
        ismissing(to) && return x[from:end]
        ismissing(from) && return x[1:to]
        return x[from:to]
end

function paa(x::AbstractArray{T} where T <: Real;
        n_chunks::Union{Missing,Int}=missing,
        f::Function=mean,
        kwargs...)
    if ismissing(n_chunks) return x
    else
        N = length(x)
        @assert 0 ≤ n_chunks && n_chunks ≤ N "The number of chunks must be in [0,$(N)]"

        y = Array{Float64}(undef, n_chunks) # TODO Float64?
        for i in 1:n_chunks
            l = Int(ceil((N*(i-1)/n_chunks) + 1))
            h = Int(ceil(N*i/n_chunks))
            y[i] = f(x[l:h]; kwargs...)
        end
        return y
    end
end


text:
function neighborhood(A::Array, I::CartesianIndex)
    xmax, ymax = size(A)
    view(
        A,
        max(1, I[1] - 1):min(ymax, I[1] + 1),
        max(1, I[2] - 1):min(ymax, I[2] + 1),
    )
end

function liferule(A, I::CartesianIndex)
    num_ones = sum(neighborhood(A, I)) - A[I]
    if num_ones == 3
        return 1
    elseif num_ones == 2
        return A[I]
    else
        return 0
    end
end

function next(A::Array)
    B = zeros(size(A))
    for I in CartesianIndices(A)
        B[I] = liferule(A, I)
    end
    return B
end

max_stars_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_issues_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_forks_repo_name: VasanthManiVasi/Conways-Game-of-Life


text:
# Themelis, Stella, Patrinos, "Forward-backward envelope for the sum of two
# nonconvex functions: Further properties and nonmonotone line-search
# algorithms", SIAM Journal on Optimization, vol. 28, no. 3, pp. 2274–2303
# (2018).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    ZeroFPRIteration(; <keyword-arguments>)

Instantiate the ZeroFPR algorithm (see [1]) for solving optimization problems
of the form

    minimize f(Ax) + g(x),

where `f` is smooth and `A` is a linear mapping (for example, a matrix).

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `A=I`: linear operator (e.g. a matrix).
- `g=Zero()`: proximable objective term.
- `Lf=nothing`: Lipschitz constant of the gradient of x ↦ f(Ax).
- `gamma=nothing`: stepsize to use, defaults to `1/Lf` if not set (but `Lf` is).
- `adaptive=false`: forces the method stepsize to be adaptively adjusted.
- `minimum_gamma=1e-7`: lower bound to `gamma` in case `adaptive == true`.
- `max_backtracks=20`: maximum number of line-search backtracks.
- `directions=LBFGS(5)`: strategy to use to compute line-search directions.

# References
1. Themelis, Stella, Patrinos, "Forward-backward envelope for the sum of two nonconvex functions: Further properties and nonmonotone line-search algorithms", SIAM Journal on Optimization, vol. 28, no. 3, pp. 2274–2303 (2018).
"""
Base.@kwdef struct ZeroFPRIteration{R,Tx,Tf,TA,Tg,TLf,Tgamma,D}
    f::Tf = Zero()
    A::TA = I
    g::Tg = Zero()
    x0::Tx
    alpha::R = real(eltype(x0))(0.95)
    beta::R = real(eltype(x0))(0.5)
    Lf::TLf = nothing
    gamma::Tgamma = Lf === nothing ? nothing : (alpha / Lf)
    adaptive::Bool = gamma === nothing
    minimum_gamma::R = real(eltype(x0))(1e-7)
    max_backtracks::Int = 20
    directions::D = LBFGS(5)
end

Base.IteratorSize(::Type{<:ZeroFPRIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct ZeroFPRState{R,Tx,TAx,TH}
    x::Tx             # iterate
    Ax::TAx           # A times x
    f_Ax::R           # value of smooth term
    grad_f_Ax::TAx    # gradient of f at Ax
    At_grad_f_Ax::Tx  # gradient of smooth term
    gamma::R          # stepsize parameter of forward and backward steps
    y::Tx             # forward point
    xbar::Tx          # forward-backward point
    g_xbar::R         # value of nonsmooth term (at xbar)
    res::Tx           # fixed-point residual at iterate (= x - xbar)
    H::TH             # variable metric
    tau::R = zero(gamma)
    Axbar::TAx = similar(Ax)
    grad_f_Axbar::TAx = similar(Ax)
    At_grad_f_Axbar::Tx = similar(x)
    xbarbar::Tx = similar(x)
    res_xbar::Tx = similar(x)
    xbar_prev::Tx = similar(x)
    res_xbar_prev::Tx = similar(x)
    is_prev_set::Bool = false
    d::Tx = similar(x)
    Ad::TAx = similar(Ax)
end

f_model(iter::ZeroFPRIteration, state::ZeroFPRState) = f_model(state.f_Ax, state.At_grad_f_Ax, state.res, iter.alpha / state.gamma)

function Base.iterate(iter::ZeroFPRIteration{R}) where R
    x = copy(iter.x0)
    Ax = iter.A * x
    grad_f_Ax, f_Ax = gradient(iter.f, Ax)
    gamma = iter.gamma === nothing ? iter.alpha / lower_bound_smoothness_constant(iter.f, iter.A, x, grad_f_Ax) : iter.gamma
    At_grad_f_Ax = iter.A' * grad_f_Ax
    y = x - gamma .* At_grad_f_Ax
    xbar, g_xbar = prox(iter.g, y, gamma)
    state = ZeroFPRState(
        x=x, Ax=Ax, f_Ax=f_Ax, grad_f_Ax=grad_f_Ax, At_grad_f_Ax=At_grad_f_Ax,
        gamma=gamma, y=y, xbar=xbar, g_xbar=g_xbar, res=x - xbar, H=initialize(iter.directions, x),
    )
    return state, state
end

set_next_direction!(::QuasiNewtonStyle, ::ZeroFPRIteration, state::ZeroFPRState) = mul!(state.d, state.H, -state.res_xbar)
set_next_direction!(::NoAccelerationStyle, ::ZeroFPRIteration, state::ZeroFPRState) = state.d .= .-state.res
set_next_direction!(iter::ZeroFPRIteration, state::ZeroFPRState) = set_next_direction!(acceleration_style(typeof(iter.directions)), iter, state)

update_direction_state!(::QuasiNewtonStyle, ::ZeroFPRIteration, state::ZeroFPRState) = update!(state.H, state.xbar - state.xbar_prev, state.res_xbar - state.res_xbar_prev)
update_direction_state!(::NoAccelerationStyle, ::ZeroFPRIteration, state::ZeroFPRState) = return
update_direction_state!(iter::ZeroFPRIteration, state::ZeroFPRState) = update_direction_state!(acceleration_style(typeof(iter.directions)), iter, state)

reset_direction_state!(::QuasiNewtonStyle, ::ZeroFPRIteration, state::ZeroFPRState) = reset!(state.H)
reset_direction_state!(::NoAccelerationStyle, ::ZeroFPRIteration, state::ZeroFPRState) = return
reset_direction_state!(iter::ZeroFPRIteration, state::ZeroFPRState) = reset_direction_state!(acceleration_style(typeof(iter.directions)), iter, state)

function Base.iterate(iter::ZeroFPRIteration{R}, state::ZeroFPRState) where R
    f_Axbar_upp, f_Axbar = if iter.adaptive == true
        gamma_prev = state.gamma
        state.gamma, state.g_xbar, f_Axbar, f_Axbar_upp = backtrack_stepsize!(
            state.gamma, iter.f, iter.A, iter.g,
            state.x, state.f_Ax, state.At_grad_f_Ax, state.y, state.xbar, state.g_xbar, state.res,
            state.Axbar, state.grad_f_Axbar,
            alpha = iter.alpha, minimum_gamma = iter.minimum_gamma,
        )
        if state.gamma != gamma_prev
            reset_direction_state!(iter, state)
        end
        f_Axbar_upp, f_Axbar
    else
        mul!(state.Axbar, iter.A, state.xbar)
        f_model(iter, state), gradient!(state.grad_f_Axbar, iter.f, state.Axbar)
    end

    # compute FBE
    FBE_x = f_Axbar_upp + state.g_xbar

    # compute residual at xbar
    mul!(state.At_grad_f_Axbar, iter.A', state.grad_f_Axbar)
    state.y .= state.xbar .- state.gamma .* state.At_grad_f_Axbar
    g_xbarbar = prox!(state.xbarbar, iter.g, state.y, state.gamma)
    state.res_xbar .= state.xbar .- state.xbarbar

    if state.is_prev_set == true
        update_direction_state!(iter, state)
    end

    copyto!(state.xbar_prev, state.xbar)
    copyto!(state.res_xbar_prev, state.res_xbar)
    state.is_prev_set = true

    set_next_direction!(iter, state)

    # Perform line-search over the FBE
    tau = R(1)
    mul!(state.Ad, iter.A, state.d)

    sigma = iter.beta * (0.5 / state.gamma) * (1 - iter.alpha)
    tol = 10 * eps(R) * (1 + abs(FBE_x))
    threshold = FBE_x - sigma * norm(state.res)^2 + tol

    for _ in 1:iter.max_backtracks
        state.x .= state.xbar_prev .+ tau .* state.d
        state.Ax .= state.Axbar .+ tau .* state.Ad
        # TODO: can precompute most of next line in case f is quadratic
        state.f_Ax = gradient!(state.grad_f_Ax, iter.f, state.Ax)
        mul!(state.At_grad_f_Ax, iter.A', state.grad_f_Ax)
        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_xbar = prox!(state.xbar, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.xbar
        FBE_x = f_model(iter, state) + state.g_xbar

        if FBE_x <= threshold
            state.tau = tau
            return state, state
        end

        tau *= 0.5
    end

    @warn "stepsize `tau` became too small ($(tau)), stopping the iterations"
    return nothing
end

# Solver

struct ZeroFPR{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::ZeroFPR)(x0; kwargs...)
    stop(state::ZeroFPRState) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it,
        state.gamma,
        norm(state.res, Inf) / state.gamma,
        state.tau,
    )
    iter = ZeroFPRIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.xbar, num_iters
end

ZeroFPR(; maxit=1_000, tol=1e-8, verbose=false, freq=10, kwargs...) =
    ZeroFPR(maxit, tol, verbose, freq, kwargs)


text:
using SIMD
using Base.Test

macro showtest(expr)
    if length(expr.args) == 3
        lhs = expr.args[1]
        rhs = expr.args[3]
    else
        lhs = expr.args[1]
        rhs = nothing
    end
    esc(quote
        # @show $lhs
        # @show $rhs
        @test $expr
        # println()
    end)
end

info("Basic definitions")

# The vector we are testing. Ideally, we should be able to use any vector size
# anywhere, but LLVM codegen bugs prevent us from doing so -- thus we make this
# a parameter.
const nbytes = 32

const L8 = nbytes÷4
const L4 = nbytes÷8

typealias V8I32 Vec{L8,Int32}
typealias V4F64 Vec{L4,Float64}

info("Type properties")

@showtest eltype(V8I32) === Int32
@showtest eltype(V4F64) === Float64
@showtest length(V8I32) == L8
@showtest length(V4F64) == L4
@showtest ndims(V8I32) == 1
@showtest ndims(V4F64) == 1
@showtest size(V8I32,1) == L8
@showtest size(V4F64,1) == L4
@showtest size(V8I32) == (L8,)
@showtest size(V4F64) == (L4,)

info("Type conversion")

const v8i32 = ntuple(i->Int32(ifelse(isodd(i), i, -i)), L8)
const v4f64 = ntuple(i->Float64(ifelse(isodd(i), i, -i)), L4)

@showtest string(V8I32(v8i32)) == "Int32<" * string(v8i32)[2:end-1] * ">"
@showtest string(V4F64(v4f64)) == "Float64<" * string(v4f64)[2:end-1] * ">"

@showtest V8I32(v8i32).elts === v8i32
@showtest V4F64(v4f64).elts === v4f64

@showtest V8I32(9).elts === ntuple(i->Int32(9), L8)
@showtest V4F64(9).elts === ntuple(i->9.0, L4)
@showtest V8I32(ntuple(i->Float32(v8i32[i]), L8)).elts === v8i32
@showtest V4F64(ntuple(i->Int64(v4f64[i]), L4)).elts === v4f64

@showtest NTuple{L8,Int32}(V8I32(v8i32)) === v8i32
@showtest NTuple{L4,Float64}(V4F64(v4f64)) === v4f64

info("Element-wise access")

for i in 1:L8
    @showtest setindex(V8I32(v8i32), Val{i}, 9.0).elts ===
        ntuple(j->Int32(ifelse(j==i, 9, v8i32[j])), L8)
    @showtest setindex(V8I32(v8i32), i, 9.0).elts ===
        ntuple(j->Int32(ifelse(j==i, 9, v8i32[j])), L8)

    @showtest V8I32(v8i32)[Val{i}] === v8i32[i]
    @showtest V8I32(v8i32)[i] === v8i32[i]
end
@test_throws BoundsError setindex(V8I32(v8i32), Val{0}, 0)
@test_throws BoundsError setindex(V8I32(v8i32), Val{L8+1}, 0)
@test_throws BoundsError setindex(V8I32(v8i32), 0, 0)
@test_throws BoundsError setindex(V8I32(v8i32), L8+1, 0)
@test_throws BoundsError V8I32(v8i32)[Val{0}]
@test_throws BoundsError V8I32(v8i32)[Val{L8+1}]
@test_throws BoundsError V8I32(v8i32)[0]
@test_throws BoundsError V8I32(v8i32)[L8+1]

for i in 1:L4
    @showtest setindex(V4F64(v4f64), Val{i}, 9).elts ===
        ntuple(j->Float64(ifelse(j==i, 9.0, v4f64[j])), L4)
    @showtest setindex(V4F64(v4f64), i, 9).elts ===
        ntuple(j->Float64(ifelse(j==i, 9.0, v4f64[j])), L4)

    @showtest V4F64(v4f64)[Val{i}] === v4f64[i]
    @showtest V4F64(v4f64)[i] === v4f64[i]
end

info("Integer arithmetic functions")

const v8i32b = map(x->Int32(x+1), v8i32)
const v8i32c = map(x->Int32(x*2), v8i32)

notbool(x) = !(x>=typeof(x)(0))
for op in (~, +, -, abs, notbool, sign, signbit)
    # @show op
    @showtest op(V8I32(v8i32)).elts === map(op, v8i32)
end

for op in (
        +, -, *, ÷, %, ==, !=, <, <=, >, >=,
        copysign, div, flipsign, max, min, rem)
    # @show op
    @showtest op(V8I32(v8i32), V8I32(v8i32b)).elts === map(op, v8i32, v8i32b)
end

ifelsebool(x,y,z) = ifelse(x>=typeof(x)(0),y,z)
for op in (ifelsebool, muladd)
    # @show op
    @showtest op(V8I32(v8i32), V8I32(v8i32b), V8I32(v8i32c)).elts ===
        map(op, v8i32, v8i32b, v8i32c)
end

info("Floating point arithmetic functions")

const v4f64b = map(x->Float64(x+1), v4f64)
const v4f64c = map(x->Float64(x*2), v4f64)

logabs(x) = log(abs(x))
log10abs(x) = log10(abs(x))
log2abs(x) = log2(abs(x))
powi4(x) = x^4
sqrtabs(x) = sqrt(abs(x))
for op in (
        +, -,
        abs, ceil, inv, isinf, isinf, isnan, issubnormal, floor, powi4,
        round, sign, signbit, sqrtabs, trunc)
    # @show op
    @showtest op(V4F64(v4f64)).elts === map(op, v4f64)
end
function Base.isapprox(t1::Tuple,t2::Tuple)
    length(t1)==length(t2) &&
        all(Bool[isapprox(t1[i], t2[i]) for i in 1:length(t1)])
end
for op in (cos, exp, exp10, exp2, logabs, log10abs, log2abs, sin)
    # @show op
    rvec = op(V4F64(v4f64)).elts
    rsca = map(op, v4f64)
    @showtest typeof(rvec) === typeof(rsca)
    @showtest isapprox(rvec, rsca)
end

@showtest isfinite(V4F64(0.0))[1]
@showtest isfinite(V4F64(-0.0))[1]
@showtest isfinite(V4F64(nextfloat(0.0)))[1]
@showtest isfinite(V4F64(-nextfloat(0.0)))[1]
@showtest isfinite(V4F64(1.0))[1]
@showtest isfinite(V4F64(-1.0))[1]
@showtest !isfinite(V4F64(Inf))[1]
@showtest !isfinite(V4F64(-Inf))[1]
@showtest !isfinite(V4F64(NaN))[1]

@showtest !isinf(V4F64(0.0))[1]
@showtest !isinf(V4F64(-0.0))[1]
@showtest !isinf(V4F64(nextfloat(0.0)))[1]
@showtest !isinf(V4F64(-nextfloat(0.0)))[1]
@showtest !isinf(V4F64(1.0))[1]
@showtest !isinf(V4F64(-1.0))[1]
@showtest isinf(V4F64(Inf))[1]
@showtest isinf(V4F64(-Inf))[1]
@showtest !isinf(V4F64(NaN))[1]

@showtest !isnan(V4F64(0.0))[1]
@showtest !isnan(V4F64(-0.0))[1]
@showtest !isnan(V4F64(nextfloat(0.0)))[1]
@showtest !isnan(V4F64(-nextfloat(0.0)))[1]
@showtest !isnan(V4F64(1.0))[1]
@showtest !isnan(V4F64(-1.0))[1]
@showtest !isnan(V4F64(Inf))[1]
@showtest !isnan(V4F64(-Inf))[1]
@showtest isnan(V4F64(NaN))[1]

@showtest !issubnormal(V4F64(0.0))[1]
@showtest !issubnormal(V4F64(-0.0))[1]
@showtest issubnormal(V4F64(nextfloat(0.0)))[1]
@showtest issubnormal(V4F64(-nextfloat(0.0)))[1]
@showtest !issubnormal(V4F64(1.0))[1]
@showtest !issubnormal(V4F64(-1.0))[1]
@showtest !issubnormal(V4F64(Inf))[1]
@showtest !issubnormal(V4F64(-Inf))[1]
@showtest !issubnormal(V4F64(NaN))[1]

@showtest !signbit(V4F64(0.0))[1]
@showtest signbit(V4F64(-0.0))[1]
@showtest !signbit(V4F64(nextfloat(0.0)))[1]
@showtest signbit(V4F64(-nextfloat(0.0)))[1]
@showtest !signbit(V4F64(1.0))[1]
@showtest signbit(V4F64(-1.0))[1]
@showtest !signbit(V4F64(Inf))[1]
@showtest signbit(V4F64(-Inf))[1]
@showtest !signbit(V4F64(NaN))[1]

for op in (
        +, -, *, /, %, ^, ==, !=, <, <=, >, >=,
        copysign, flipsign, max, min, rem)
    # @show op
    @showtest op(V4F64(v4f64), V4F64(v4f64b)).elts === map(op, v4f64, v4f64b)
end

for op in (fma, ifelsebool, muladd)
    # @show op
    @showtest op(V4F64(v4f64), V4F64(v4f64b), V4F64(v4f64c)).elts ===
        map(op, v4f64, v4f64b, v4f64c)
end

info("Type promotion")

for op in (
        ==, !=, <, <=, >, >=,
        &, |, $, +, -, *, copysign, div, flipsign, max, min, rem)
    @showtest op(42, V8I32(v8i32)) === op(V8I32(42), V8I32(v8i32))
    @showtest op(V8I32(v8i32), 42) === op(V8I32(v8i32), V8I32(42))
end
@showtest ifelse(signbit(V8I32(v8i32)), 42, V8I32(v8i32)) ===
    ifelse(signbit(V8I32(v8i32)), V8I32(42), V8I32(v8i32))
@showtest ifelse(signbit(V8I32(v8i32)), V8I32(v8i32), 42) ===
    ifelse(signbit(V8I32(v8i32)), V8I32(v8i32), V8I32(42))
for op in (muladd,)
    @showtest op(42, 42, V8I32(v8i32)) ===
        op(V8I32(42), V8I32(42), V8I32(v8i32))
    @showtest op(42, V8I32(v8i32), V8I32(v8i32)) ===
        op(V8I32(42), V8I32(v8i32), V8I32(v8i32))
    @showtest op(V8I32(v8i32), 42, V8I32(v8i32)) ===
        op(V8I32(v8i32), V8I32(42), V8I32(v8i32))
    @showtest op(V8I32(v8i32), V8I32(v8i32), 42) ===
        op(V8I32(v8i32), V8I32(v8i32), V8I32(42))
    @showtest op(42, V8I32(v8i32), 42) ===
        op(V8I32(42), V8I32(v8i32), V8I32(42))
    @showtest op(V8I32(v8i32), 42, 42) ===
        op(V8I32(v8i32), V8I32(42), V8I32(42))
end

for op in (
        ==, !=, <, <=, >, >=,
        +, -, *, /, ^, copysign, flipsign, max, min, rem)
    @showtest op(42, V4F64(v4f64)) === op(V4F64(42), V4F64(v4f64))
    @showtest op(V4F64(v4f64), 42) === op(V4F64(v4f64), V4F64(42))
end
@showtest ifelse(signbit(V4F64(v4f64)), 42, V4F64(v4f64)) ===
    ifelse(signbit(V4F64(v4f64)), V4F64(42), V4F64(v4f64))
@showtest ifelse(signbit(V4F64(v4f64)), V4F64(v4f64), 42) ===
    ifelse(signbit(V4F64(v4f64)), V4F64(v4f64), V4F64(42))
for op in (fma, muladd)
    @showtest op(42, 42, V4F64(v4f64)) ===
        op(V4F64(42), V4F64(42), V4F64(v4f64))
    @showtest op(42, V4F64(v4f64), V4F64(v4f64)) ===
        op(V4F64(42), V4F64(v4f64), V4F64(v4f64))
    @showtest op(V4F64(v4f64), 42, V4F64(v4f64)) ===
        op(V4F64(v4f64), V4F64(42), V4F64(v4f64))
    @showtest op(V4F64(v4f64), V4F64(v4f64), 42) ===
        op(V4F64(v4f64), V4F64(v4f64), V4F64(42))
    @showtest op(42, V4F64(v4f64), 42) ===
        op(V4F64(42), V4F64(v4f64), V4F64(42))
    @showtest op(V4F64(v4f64), 42, 42) ===
        op(V4F64(v4f64), V4F64(42), V4F64(42))
end

info("Reduction operations")

for op in (maximum, minimum, sum, prod)
    # @show op
    @showtest op(V8I32(v8i32)) === op(v8i32)
end
@showtest all(V8I32(v8i32)) == reduce(&, v8i32)
@showtest any(V8I32(v8i32)) == reduce(|, v8i32)

for op in (maximum, minimum, sum, prod)
    # @show op
    @showtest op(V4F64(v4f64)) === op(v4f64)
end

info("Load and store functions")

const arri32 = Int32[i for i in 1:(2*L8)]
for i in 1:length(arri32)-(L8-1)
    @showtest vload(V8I32, arri32, i) === V8I32(ntuple(j->i+j-1, L8))
end
for i in 1:L8:length(arri32)-(L8-1)
    @showtest vloada(V8I32, arri32, i) === V8I32(ntuple(j->i+j-1, L8))
end
vstorea(V8I32(0), arri32, 1)
vstore(V8I32(1), arri32, 2)
for i in 1:length(arri32)
    @showtest arri32[i] == if i==1 0 elseif i<=(L8+1) 1 else i end
end

const arrf64 = Float64[i for i in 1:(4*L4)]
for i in 1:length(arrf64)-(L4-1)
    @showtest vload(V4F64, arrf64, i) === V4F64(ntuple(j->i+j-1, L4))
end
for i in 1:4:length(arrf64)-(L4-1)
    @showtest vloada(V4F64, arrf64, i) === V4F64(ntuple(j->i+j-1, L4))
end
vstorea(V4F64(0), arrf64, 1)
vstore(V4F64(1), arrf64, 2)
for i in 1:length(arrf64)
    @showtest arrf64[i] == if i==1 0 elseif i<=(L4+1) 1 else i end
end

info("Real-world examples")

function vadd!{N,T}(xs::Vector{T}, ys::Vector{T}, ::Type{Vec{N,T}})
    @assert length(ys) == length(xs)
    @assert length(xs) % N == 0
    @inbounds for i in 1:N:length(xs)
        xv = vload(Vec{N,T}, xs, i)
        yv = vload(Vec{N,T}, ys, i)
        xv += yv
        vstore(xv, xs, i)
    end
end

let xs = Float64[i for i in 1:(4*L4)];
    ys = Float64[1 for i in 1:(4*L4)]
    vadd!(xs, ys, V4F64)
    @showtest xs == Float64[i+1 for i in 1:(4*L4)]
    # @code_native vadd!(xs, ys, V4F64)
end

function vsum{N,T}(xs::Vector{T}, ::Type{Vec{N,T}})
    @assert length(xs) % N == 0
    sv = Vec{N,T}(0)
    @inbounds for i in 1:N:length(xs)
        xv = vload(Vec{N,T}, xs, i)
        sv += xv
    end
    sum(sv)
end

let xs = Float64[i for i in 1:(4*L4)]
    s = vsum(xs, V4F64)
    @showtest s === (x->(x^2+x)/2)(Float64(4*L4))
    # @code_native vsum(xs, V4F64)
end

function vadd_masked!{N,T}(xs::Vector{T}, ys::Vector{T}, ::Type{Vec{N,T}})
    @assert length(ys) == length(xs)
    limit = length(xs) - (N-1)
    vlimit = Vec{N,Int}(let l=length(xs); (l:l+N-1...) end)
    @inbounds for i in 1:N:length(xs)
        xv = vload(Vec{N,T}, xs, i)
        yv = vload(Vec{N,T}, ys, i)
        xv += yv
        if i <= limit
            vstore(xv, xs, i)
        else
            mask = Vec{N,Int}(i) <= vlimit
            vstore(xv, xs, i, mask)
        end
    end
end

let xs = Float64[i for i in 1:13],
    ys = Float64[1 for i in 1:13]
    vadd_masked!(xs, ys, V4F64)
    @showtest xs == Float64[i+1 for i in 1:13]
    # @code_native vadd!(xs, ys, V4F64)
end

function vsum_masked{N,T}(xs::Vector{T}, ::Type{Vec{N,T}})
    vlimit = Vec{N,Int}(let l=length(xs); (l:l+N-1...) end)
    sv = Vec{N,T}(0)
    @inbounds for i in 1:N:length(xs)
        mask = Vec{N,Int}(i) <= vlimit
        xv = vload(Vec{N,T}, xs, i, mask)
        sv += xv
    end
    sum(sv)
end

let xs = Float64[i for i in 1:13]
    s = vsum_masked(xs, V4F64)
    @showtest s === sum(xs)
    # @code_native vsum(xs, V4F64)
end


text:
push!(LOAD_PATH, "../src/")
include("../src/sail_route.jl")

using Documenter, sail_route
ENV["DOCUMENTER_DEBUG"]=true

makedocs(
    modules = [sail_route]
)

deploydocs(
    repo   = "github.com/TAJD/sail_route.jl.git",
    julia = "nightly",
    osname = "osx"
)


text:
using LinearAlgebra
using CUDA
using Distributions
using Random

struct OptimizerCmaEsCfg
    population_size::Int
    sigma::Float64

    function OptimizerCmaEsCfg(configuration::OrderedDict)
        new(configuration["population_size"], configuration["sigma"])
    end
end

mutable struct OptimizerCmaEs
    lambda_::Any
    dim::Any
    chiN::Any
    mu::Any
    weights::Any
    mueff::Any
    cc::Any
    cs::Any
    ps::Any
    pc::Any
    centroid::Any
    update_count::Any
    ccov1::Any
    ccovmu::Any
    C::Any
    sigma::Any
    damps::Any
    diagD::Any
    B::Any
    BD::Any
    genomes::Any

    function OptimizerCmaEs(individual_size::Int, optimizer_configuration::OrderedDict; test = false, eigenvectors1 = Nothing, indx1 = Nothing)

        config = OptimizerCmaEsCfg(optimizer_configuration)

        centroid = zeros(individual_size)

        dim = individual_size
        pc = zeros(dim)
        ps = zeros(dim)
        chiN = sqrt(dim) * (1 - 1 / (4 * dim) + 1 / (21 * dim^2))

        C = Matrix(1.0I, dim, dim)

        C_GPU = CuArray(C)
        val_GPU, vec_GPU = CUDA.CUSOLVER.syevd!('V', 'U', C_GPU)
        diagD = Array(val_GPU)
        B = Array(vec_GPU)

        indx = sortperm(diagD)

        if test == true
            eigenvectors_test = copy(B)
            B = copy(eigenvectors1)
            indx_test = copy(indx)
            indx = copy(indx1)
            diagD_test = copy(diagD)
        else
            eigenvectors_test = Nothing
            indx_test = Nothing
            diagD_test = Nothing
        end

        diagD = diagD[indx] .^ 0.5
        B = B[:, indx]
        BD = B .* diagD'

        lambda_ = config.population_size
        update_count = 0

        # Compute params
        mu = Int(lambda_ / 2)
        weights = log(mu + 0.5) .- log.(collect(1:mu))
        weights /= sum(weights)
        mueff = 1 / sum(weights .^ 2)
        cc = 4 / (dim + 4)
        cs = (mueff + 2) / (dim + mueff + 3)
        ccov1 = 2 / ((dim + 1.3)^2 + mueff)
        ccovmu = 2 * (mueff - 2 + 1 / mueff) / ((dim + 2)^2 + mueff)
        ccovmu = min(1 - ccov1, ccovmu)
        damps = 1 + 2 * max(0, sqrt((mueff - 1) / (dim + 1)) - 1) + cs

        genomes = zeros(lambda_, individual_size)

        optimizer = new(lambda_, dim, chiN, mu, weights, mueff, cc, cs, ps, pc, centroid, update_count, ccov1, ccovmu, C, config.sigma, damps, diagD, B, BD, genomes)

        if test == true
            return optimizer, eigenvectors_test, indx_test, diagD_test
        else
            return optimizer
        end
    end
end

function ask(optimizer::OptimizerCmaEs; test = false, randoms1 = Nothing)

    randoms = rand(Normal(), size(optimizer.genomes))

    if test == true
        randoms_test = copy(randoms)
        randoms = copy(randoms1)
    else
        randoms_test = Nothing
    end

    optimizer.genomes = optimizer.centroid' .+ (optimizer.sigma .* (randoms * optimizer.BD'))

    if test == true
        return optimizer.genomes, randoms_test
    else
        return optimizer.genomes
    end

end

function tell(optimizer::OptimizerCmaEs, rewards_training; test = false, eigenvectors1 = Nothing, indx1 = Nothing)

    genomes_sorted = optimizer.genomes[sortperm(rewards_training, rev = true), :]

    old_centroid = copy(optimizer.centroid)
    optimizer.centroid = genomes_sorted[1:optimizer.mu, :]' * optimizer.weights

    c_diff = optimizer.centroid - old_centroid

    # Cumulation : update evolution path
    optimizer.ps =
        (1 - optimizer.cs) .* optimizer.ps +
        sqrt(optimizer.cs * (2 - optimizer.cs) * optimizer.mueff) ./ optimizer.sigma *
        optimizer.B *
        ((1 ./ optimizer.diagD) .* optimizer.B' * c_diff)

    hsig = float(
        norm(optimizer.ps) /
        sqrt(1.0 - (1 - optimizer.cs)^(2 * (optimizer.update_count + 1))) /
        optimizer.chiN < (1.4 + 2 / (optimizer.dim + 1)),
    )

    optimizer.update_count += 1

    optimizer.pc =
        (1 - optimizer.cc) * optimizer.pc +
        hsig * sqrt(optimizer.cc * (2 - optimizer.cc) * optimizer.mueff) / optimizer.sigma *
        c_diff

    # Update covariance matrix
    artmp = genomes_sorted[1:optimizer.mu, :]' .- old_centroid

    optimizer.C =
        (
            1 - optimizer.ccov1 - optimizer.ccovmu +
            (1 - hsig) * optimizer.ccov1 * optimizer.cc * (2 - optimizer.cc)
        ) * optimizer.C +
        optimizer.ccov1 * optimizer.pc * optimizer.pc' +
        optimizer.ccovmu * (optimizer.weights' .* artmp) * artmp' / optimizer.sigma^2

    optimizer.sigma *=
        exp((norm(optimizer.ps) / optimizer.chiN - 1) * optimizer.cs / optimizer.damps)

    C_GPU = CuArray(optimizer.C)
    val_GPU, vec_GPU = CUDA.CUSOLVER.syevd!('V', 'U', C_GPU)
    optimizer.diagD = Array(val_GPU)
    optimizer.B = Array(vec_GPU)

    indx = sortperm(optimizer.diagD)

    if test == true
        eigenvectors_test = copy(optimizer.B)
        optimizer.B = copy(eigenvectors1)
        indx_test = copy(indx)
        indx = copy(indx1)
        diagD_test = copy(optimizer.diagD)
    else
        eigenvectors_test = Nothing
        indx_test = Nothing
        diagD_test = Nothing
    end

    optimizer.diagD = optimizer.diagD[indx] .^ 0.5
    optimizer.B = optimizer.B[:, indx]
    optimizer.BD = optimizer.B .* optimizer.diagD'

    if test == true
        return eigenvectors_test, indx_test, diagD_test
    end

end


text:
module OpenDiagrams
export OpenDiagram, oapply, draw

using Catlab, Catlab.CategoricalAlgebra, Catlab.WiringDiagrams, Catlab.Programs
import Catlab.CategoricalAlgebra: apex, legs, feet
import Catlab.WiringDiagrams: oapply
using Catlab.Programs.DiagrammaticPrograms: NamedGraph # FIXME: Should export?
using Catlab.Graphics

""" Open diagram as a structured multicospan in R-form.

An open diagram is a diagram, represented as a `FinDomFunctor`, together with a
cospan of finite sets whose apex is the object set of the diagram's indexing
category.
"""
struct OpenDiagram{F<:FinDomFunctor,Cosp<:Multicospan}
  functor::F
  cospan::Cosp
end

apex(diagram::OpenDiagram) = diagram.functor
legs(diagram::OpenDiagram) = legs(diagram.cospan)
feet(diagram::OpenDiagram) = feet(diagram.cospan)

function OpenDiagram(F::FinDomFunctor, names::AbstractVector{Symbol})
  g = graph(dom(F))
  legs = map(name -> FinFunction([incident(g, name, :vname)], nv(g)), names)
  OpenDiagram(F, Multicospan(FinSet(1), legs))
end

const OpenFreeDiagramOb, OpenFreeDiagram = OpenACSetTypes(FreeDiagram, :V)

OpenFreeDiagram(d::FreeDiagram{Ob,Hom}, args...) where {Ob,Hom} =
  OpenFreeDiagram{Ob,Hom}(d, args...)

function oapply(uwd::RelationDiagram, diagrams::AbstractVector{<:OpenDiagram})
  open_free_diagram = oapply(uwd, map(diagrams) do d
    OpenFreeDiagram(FreeDiagram(apex(d)), d.cospan)
  end)
  free_diagram = apex(open_free_diagram)
  leg_functions = map(leg -> leg[:V], legs(open_free_diagram))

  g = NamedGraph{Symbol,Union{Symbol,Nothing}}()
  copy_parts!(g, free_diagram)
  g[:vname] = collect(apex(colimit_vertex_names(uwd, diagrams)))
  g[:ename] = nothing # FIXME: Coproduct edge names, allowing for nulls.

  diagram = FinFunctor(ob(free_diagram), hom(free_diagram),
                       FinCat(g), codom(apex(first(diagrams))))
  OpenDiagram(diagram, Multicospan(leg_functions))
end

function colimit_vertex_names(uwd::RelationDiagram,
                              diagrams::AbstractVector{<:OpenDiagram})
  @assert nboxes(uwd) == length(diagrams)
  bfd = BipartiteFreeDiagram{FinSet{<:Any,Symbol},FinFunction}()
  add_vertices₁!(bfd, njunctions(uwd), ob₁=map(uwd[:variable]) do var
    FinSet([var])
  end)
  add_vertices₂!(bfd, nboxes(uwd), ob₂=map(diagrams) do diagram
    FinSet(graph(dom(diagram.functor))[:vname])
  end)
  for (b, diagram) in zip(boxes(uwd), diagrams)
    g = graph(dom(diagram.functor))
    for (p, leg) in zip(ports(uwd, b), legs(diagram.cospan))
      j, v = junction(uwd, p), only(collect(leg))
      f = FinFunction(Dict(uwd[j, :variable] => g[v, :vname]), ob₂(bfd, b))
      add_edge!(bfd, j, b, hom=f)
    end
  end
  colimit(bfd)
end

####################
# Graphics Tooling #
####################

#TODO: Rampant type piracy in here, fix
draw(d) = to_graphviz(d, box_labels=:name, junction_labels=:variable,
graph_attrs=Dict(:start => "2", :overlap=>"scale"), port_labels=true)

using Catlab.Graphs
using Catlab.Graphs.BasicGraphs
function to_graph(J::FinCat)
    g = BasicGraphs.Graph()
    copy_parts!(g, graph(J))
    return g
end

Graphics.to_graphviz(F::FinFunctor; kw...) =
to_graphviz(GraphvizGraphs.to_graphviz_property_graph(F; kw...))

function GraphvizGraphs.to_graphviz_property_graph(F::FinFunctor; kw...)
    simplify_vname(s) = begin
      table = Dict(
    "Form0(X)" => "Ω₀",
    "Form1(X)" => "Ω₁",
    "Form2(X)" => "Ω₂",
    "DualForm0(X)" => "Ω̃₀",
    "DualForm1(X)" => "Ω̃₁",
    "DualForm2(X)" => "Ω̃₂",
    "otimes(Form0(X),Form0(X))" => "Ω₀²",
    "otimes(Form1(X),Form1(X))" => "Ω₁²",
    "otimes(Form1(X),DualForm2(X))" => "Ω₁×Ω̃₂",
    "otimes(Form1(X),Form1(X),Form1(X))" => "Ω₁³",
    "otimes(Form1(X),Form1(X),Form1(X),Form1(X))" => "Ω₁⁴"
    )
    if string(s) in keys(table)
        return table[string(s)]
    else
        println(string(s))
        return string(s)
    end
    end

    simplify_ename(s) = begin
        b = IOBuffer()
        show_unicode(b, s)
        return replace(String(take!(b)), r"{X}"=>"")
    end

    J = dom(F)
    G = graph(J)
    pg = GraphvizGraphs.to_graphviz_property_graph(to_graph(J); kw...)
    for v in vertices(G)
        lᵥ = G[v, :vname]
        tᵥ = simplify_vname(F.ob_map[v])
        set_vprops!(pg, v, Dict(:label => "$(lᵥ):$tᵥ"))
    end
    for e in edges(G)
        tₑ = F.hom_map[e]
        set_eprops!(pg, e, Dict(:label => "$(simplify_ename(tₑ))"))
    end
    pg
end

Graphics.to_graphviz(cosp::OpenDiagram; kw...) =
to_graphviz(GraphvizGraphs.to_graphviz_property_graph(cosp; kw...))

function GraphvizGraphs.to_graphviz_property_graph(cosp::OpenDiagram; kw...)
    pg = GraphvizGraphs.to_graphviz_property_graph(cosp.functor; kw...)
    label(I, l, i) = length(dom(l)) > 1 ? "$(i):$I" : "$I"
    for (I,l) in enumerate(legs(cosp.cospan))
        for i in dom(l)
            v = add_vertex!(pg)
            set_vprops!(pg, v, Dict(:label=>label(I, l,i), :shape=>"circle", :color=>"blue"))
            e = add_edge!(pg, v, l(i))
            set_eprops!(pg, e, Dict(:style=>"dashed"))
        end
    end
    return pg
end

end


text:
function ReferenceEllipsoid(arg0::jdouble, arg1::jdouble, arg2::Frame, arg3::jdouble, arg4::jdouble)
    return ReferenceEllipsoid((jdouble, jdouble, Frame, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function get_a(obj::Ellipsoid)
    return jcall(obj, "getA", jdouble, ())
end

function get_b(obj::Ellipsoid)
    return jcall(obj, "getB", jdouble, ())
end

function get_body_frame(obj::OneAxisEllipsoid)
    return jcall(obj, "getBodyFrame", Frame, ())
end

function get_c(obj::Ellipsoid)
    return jcall(obj, "getC", jdouble, ())
end

function get_c2n0(obj::ReferenceEllipsoid, arg0::jint)
    return jcall(obj, "getC2n0", jdouble, (jint,), arg0)
end

function get_cartesian_intersection_point(obj::OneAxisEllipsoid, arg0::FieldLine, arg1::FieldVector3D, arg2::Frame, arg3::FieldAbsoluteDate)
    return jcall(obj, "getCartesianIntersectionPoint", FieldVector3D, (FieldLine, FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_cartesian_intersection_point(obj::OneAxisEllipsoid, arg0::Line, arg1::Vector3D, arg2::Frame, arg3::AbsoluteDate)
    return jcall(obj, "getCartesianIntersectionPoint", Vector3D, (Line, Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_ellipsoid(obj::ReferenceEllipsoid)
    return jcall(obj, "getEllipsoid", ReferenceEllipsoid, ())
end

function get_equatorial_radius(obj::OneAxisEllipsoid)
    return jcall(obj, "getEquatorialRadius", jdouble, ())
end

function get_flattening(obj::OneAxisEllipsoid)
    return jcall(obj, "getFlattening", jdouble, ())
end

function get_frame(obj::Ellipsoid)
    return jcall(obj, "getFrame", Frame, ())
end

function get_gm(obj::ReferenceEllipsoid)
    return jcall(obj, "getGM", jdouble, ())
end

function get_grs80(::Type{ReferenceEllipsoid}, arg0::Frame)
    return jcall(ReferenceEllipsoid, "getGrs80", ReferenceEllipsoid, (Frame,), arg0)
end

function get_iers2003(::Type{ReferenceEllipsoid}, arg0::Frame)
    return jcall(ReferenceEllipsoid, "getIers2003", ReferenceEllipsoid, (Frame,), arg0)
end

function get_iers2010(::Type{ReferenceEllipsoid}, arg0::Frame)
    return jcall(ReferenceEllipsoid, "getIers2010", ReferenceEllipsoid, (Frame,), arg0)
end

function get_iers96(::Type{ReferenceEllipsoid}, arg0::Frame)
    return jcall(ReferenceEllipsoid, "getIers96", ReferenceEllipsoid, (Frame,), arg0)
end

function get_intersection_point(obj::OneAxisEllipsoid, arg0::FieldLine, arg1::FieldVector3D, arg2::Frame, arg3::FieldAbsoluteDate)
    return jcall(obj, "getIntersectionPoint", FieldGeodeticPoint, (FieldLine, FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_intersection_point(obj::OneAxisEllipsoid, arg0::Line, arg1::Vector3D, arg2::Frame, arg3::AbsoluteDate)
    return jcall(obj, "getIntersectionPoint", GeodeticPoint, (Line, Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_normal_gravity(obj::ReferenceEllipsoid, arg0::jdouble)
    return jcall(obj, "getNormalGravity", jdouble, (jdouble,), arg0)
end

function get_plane_section(obj::Ellipsoid, arg0::Vector3D, arg1::Vector3D)
    return jcall(obj, "getPlaneSection", Ellipse, (Vector3D, Vector3D), arg0, arg1)
end

function get_polar_radius(obj::ReferenceEllipsoid)
    return jcall(obj, "getPolarRadius", jdouble, ())
end

function get_spin(obj::ReferenceEllipsoid)
    return jcall(obj, "getSpin", jdouble, ())
end

function get_wgs84(::Type{ReferenceEllipsoid}, arg0::Frame)
    return jcall(ReferenceEllipsoid, "getWgs84", ReferenceEllipsoid, (Frame,), arg0)
end

function is_inside(obj::Ellipsoid, arg0::Vector3D)
    return jcall(obj, "isInside", jboolean, (Vector3D,), arg0)
end

function point_on_limb(obj::Ellipsoid, arg0::Vector3D, arg1::Vector3D)
    return jcall(obj, "pointOnLimb", Vector3D, (Vector3D, Vector3D), arg0, arg1)
end

function project_to_ground(obj::OneAxisEllipsoid, arg0::TimeStampedPVCoordinates, arg1::Frame)
    return jcall(obj, "projectToGround", TimeStampedPVCoordinates, (TimeStampedPVCoordinates, Frame), arg0, arg1)
end

function project_to_ground(obj::OneAxisEllipsoid, arg0::Vector3D, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "projectToGround", Vector3D, (Vector3D, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function set_angular_threshold(obj::OneAxisEllipsoid, arg0::jdouble)
    return jcall(obj, "setAngularThreshold", void, (jdouble,), arg0)
end

function transform(obj::OneAxisEllipsoid, arg0::FieldGeodeticPoint)
    return jcall(obj, "transform", FieldVector3D, (FieldGeodeticPoint,), arg0)
end

function transform(obj::OneAxisEllipsoid, arg0::FieldVector3D, arg1::Frame, arg2::FieldAbsoluteDate)
    return jcall(obj, "transform", FieldGeodeticPoint, (FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2)
end

function transform(obj::OneAxisEllipsoid, arg0::GeodeticPoint)
    return jcall(obj, "transform", Vector3D, (GeodeticPoint,), arg0)
end

function transform(obj::OneAxisEllipsoid, arg0::PVCoordinates, arg1::Frame, arg2::AbsoluteDate)
    return jcall(obj, "transform", FieldGeodeticPoint, (PVCoordinates, Frame, AbsoluteDate), arg0, arg1, arg2)
end

function transform(obj::OneAxisEllipsoid, arg0::Vector3D, arg1::Frame, arg2::AbsoluteDate)
    return jcall(obj, "transform", GeodeticPoint, (Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2)
end


text:
using Turing, TuringBenchmarks.TuringTools
using LinearAlgebra

include("lr_helper.jl")

X, Y = readlrdata()

@model lr_nuts(x, y, σ) = begin

    N,D = size(x)

    α ~ Normal(0, σ)
    β ~ MvNormal(zeros(D), ones(D)*σ)

    for n = 1:N
        y[n] ~ BinomialLogit(1, dot(x[n,:], β) + α)
    end
end

# Sampling parameter settings
n_adapts = 1_000

# Sampling
@tbenchmark chain = sample(lr_nuts(x, y, 100), NUTS(0.65), n_samples);


text:
"""
---
Structure describing window parameters
# --- Syntax
- winFunc	: Name of window
- winLength : Size of window
- window	: Window coefficients
# ---
# v 1.0 - Robin Gerzaguet.
"""
struct Window
	winFunc::String;
	winLength::Int;
	window::Array{Float64};
end

"""
---
Returns the desired window used in transmitter or receiver weigthed overlap and add methods (WOLA). See [1] for window reference design.
# --- Syntax
window = getWolaWindow(winFunc,nFFT,nCp,winlength)
# --- Input parameters
- winFunc	: Window type (see supported format below) [String]
- nFFT		: FFT size [Int]
- nCP		: CP size [Int]
- winlength : length of window [Int]
# --- Output parameters
window	: Window [Window]
# ---
# Supported window
- "Triangle"		: Triangle window
- "srrc"			: Square Root Raised Cosine
- "Meyer"			: Meyer window (See [1])
# ---
# References
- [1] R. Zayani, Y. Medjahdi, H. Shaiek and D. Roviras, "WOLA-OFDM: A Potential Candidate for Asynchronous 5G," 2016.
# ---
# v 1.0 - Robin Gerzaguet.
"""
function getWolaWindow(winFunc,nFFT,nCP,winlength)
	# --- Overall parameters
	sizeSymb		= nFFT + nCP;
	# --- Switch to lowercase to avoid case mistakes.
	winName = lowercase(winFunc);
	if winFunc == "triangle"
		# ----------------------------------------------------
		# --- Triangle window
		# ----------------------------------------------------
		# Classic ramp up // ramp down window
		filterTimeUp  = 1/(winlength-1) * collect(0:1:winlength-1);
		filterTimeDown= 1/(winlength-1) * collect(winlength-1:-1:0);
		window      = [filterTimeUp ; ones(Float64,sizeSymb-winlength) ; filterTimeDown];
	elseif 	winFunc == "meyer"
		# ----------------------------------------------------
		# --- Meyer window
		# ----------------------------------------------------
		ww      = Int.(floor(winlength)/2);
		# --- Up
		n       = collect(0:ww-1) ./(ww-1);
		V       = n.^4 .*(35 .- 84n .+ (70 .*n.^2) .- (20 .*n.^3));
		wUp     = 1/2 .- 1/2 .* cos.( pi .* V);
		# --- Down
		wDown     = collect(wUp[ end:-1:1 ]);
		# --- Complete filter
		window                = [wUp;ones(Float64,sizeSymb) ;wDown];
	else
		error("Unsupported apoodisation window");
	end
	return Window(winFunc,winlength,window)
end


text:
test_list = ["CalendarTest.jl";"YearFractionTest.jl"]

println("Running tests:\n")
i=1;
for current_test in test_list
	println("------------------------------------------------------------")
    println("  * $(current_test) *")
    include(current_test)
	println("------------------------------------------------------------")
	if (i<length(test_list))
		println("")
	end
	i+=1;
end

println("DateLib Test Passed")


text:
module NerdHelperBouncyJellyfish

using ..Ahorn, Maple

@mapdef Entity "NerdHelper/BouncyJellyfish" BouncyJellyfish(x::Integer, y::Integer, bubble::Bool=false, tutorial::Bool=false)

const placements = Ahorn.PlacementDict(
    "Bouncy Jellyfish (Nerd Helper)" => Ahorn.EntityPlacement(
        BouncyJellyfish
    ),
    "Bouncy Jellyfish (Floating) (Nerd Helper)" => Ahorn.EntityPlacement(
        BouncyJellyfish,
        "point",
        Dict{String, Any}(
            "bubble" => true
        )
    )
)

sprite = "objects/glider/idle0"

function Ahorn.selection(entity::BouncyJellyfish)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BouncyJellyfish, room::Maple.Room)
    Ahorn.drawSprite(ctx, sprite, 0, 0)

    if get(entity, "bubble", false)
        curve = Ahorn.SimpleCurve((-7, -1), (7, -1), (0, -6))
        Ahorn.drawSimpleCurve(ctx, curve, (1.0, 1.0, 1.0, 1.0), thickness=1)
    end
end

end

max_stars_repo_name: Kalobi/NerdHelper

max_issues_repo_name: Kalobi/NerdHelper

max_forks_repo_name: Kalobi/NerdHelper


text:
using Cytof5
using Distributions
using Flux, Flux.Tracker
include("util.jl")
include("logprob_lam.jl")

function sample_lam(state::Cytof5.VB.StateMP, y::Vector{Matrix{Float64}}, c::Cytof5.VB.Constants)
  _, s, ys, _ = Cytof5.VB.rsample(state, y, c)
  lp_lam= logprob_lam(data(s), Tracker.data.(ys), c)
  lam = [[Cytof5.MCMC.wsample_logprob(lp_lam[i][n, :])
          for n in 1:c.N[i]] for i in 1:c.I]

  # Relabel noisy class as 0
  for i in 1:c.I
    for n in 1:c.N[i]
      if lam[i][n] > c.K
        lam[i][n] = 0
      end
    end
  end

  return Vector{Int8}.(lam)
end

VV = Vector{Vector{I}} where I
VVV = Vector{VV{I}}  where I

function lam_f(lam_samps::VVV{T}, i::Integer, n::Integer, f::Function) where {T <: Integer}
  B = length(lam_samps)
  return f(lam_samps[b][i][n] for b in 1:B)
end

function lam_f(lam_samps::VVV{T}, f::Function) where {T <: Integer}
  I = length(lam_samps[1])
  N = length.(lam_samps[1])
  B = length(lam_samps)

  return [[lam_f(lam_samps, i, n, f) for n in 1:N[i]] for i in 1:I]
end

#= TEST -------------------------------------------------------
# read data
using BSON
output_path = "../results/vb-sim-paper/test/0/output.bson"
out = BSON.load(output_path)
simdat = out[:simdat]
simdat[:y] = Matrix{Float64}.(simdat[:y])

c = out[:c]
elbo = out[:metrics][:elbo] / sum(c.N)
state = out[:state]
state_hist = out[:state_hist]
metrics = out[:metrics]
m = [isnan.(yi) for yi in simdat[:y]]

# Test stuff
@time lam_samps = [sample_lam(state, simdat[:y], c) for b in 1:10]
@time lam_mode = lam_f(lam_samps, mode)
=#


text:
for (nm, mmf, nnl, nl) in ((:AsympReg, :AsympRegmmf, 1, 2),
                           (:AsympOff, :AsympOffmmf, 2, 1),
                           (:AsympOrig, :AsympOrigmmf, 1, 1),
                           (:Biexp, :Biexpmmf, 2, 2),
                           (:BolusSD1, :BolusSD1mmf, 1, 1),
                           (:Chwirut, :Chwirutmmf, 2, 1),
                           (:Gompertz, :Gompertzmmf, 2, 1),
                           (:Logis3P, :Logis3Pmmf, 2, 1),
                           (:Logis4P, :Logis4Pmmf, 2, 2),
                           (:MicMen, :MicMenmmf, 1, 1))
    @eval begin
        immutable $nm{T<:FP} <: PLregModF{T}
            x::Matrix{T}
            y::Vector{T}
            mu::Vector{T}
            resid::Vector{T}
            tgrad::Matrix{T}
            MMD::Array{T,3}
            mmf::Function
        end
        function $nm{T<:FP}(x::Vector{T},y::Vector{T})
            (n = length(x)) == length(y)|| throw(DimensionMismatch(""))
            $nm(reshape(x,(1,n)),y,similar(y),similar(y),ones(T,($(nl + nnl),n)),
                zeros(T,($nnl,$nl,n)),$mmf)
        end
        $nm(x::DataVector,y::DataVector) = $nm(float(x),float(y))
        function $nm(f::Formula,dat::AbstractDataFrame)
            mf = ModelFrame(f,dat)
            mat = ModelMatrix(mf).m
            rr = model_response(mf)
            T = promote_type(eltype(mat),eltype(rr))
            $nm(convert(Vector{T},mat[:,end]),convert(Vector{T},rr))
        end
        $nm(ex::Expr,dat::AbstractDataFrame) = $nm(Formula(ex),dat)
    end
end

### Michaelis-Menten model for enzyme kinetics
formula(::MicMen) = :(Vm*x/(K+x))
pnames(::MicMen) = [:Vm, :K, :x]
dim(::MicMen) = (1,1,1)

function MicMenmmf(nlp,x,tg,MMD)
    @inbounds begin
        x1 = x[1]
        invde = inv(nlp[1] + x1)
        MMD[1,1] = -(tg[1] =  x1*invde)*invde
    end
end

function initpars(m::MicMen)
    x = vec(m.x); y = m.y
    length(y) < 2 && return [one(eltype(y))]
    cc = linreg(map(inv,x), map(inv,y))
    [cc[2] / cc[1]]
end

### Asymptotic Regression model
formula(::AsympReg) = :(Asym + (Asym - R0)*exp(-rc*x))
pnames(::AsympReg) = [:Asym,:R0,:rc,:x]
dim(::AsympReg) = (2,1,1) # conditionally linear, nonlinear, covariates

function AsympRegmmf(nlp,x,tg,MMD)
    x1 = x[1]
    tg[1] = 1 - (tg[2] = eKx = exp(-nlp[1]*x1))
    MMD[2,1] = -(MMD[1,1] = x1 * eKx)
end

function initpars(m::AsympReg)
    y = m.y
    mn,mx = extrema(y)
    [abs(linreg(vec(m.x), log(y - mn + (mx-mn)/4.))[2])]
end

### Asymptotic Regression model constrained to pass through the origin
formula(::AsympOrig) = :(Asym*(1 - exp(-rc*x)))
pnames(::AsympOrig) = [:Asym,:rc,:x]
dim(::AsympOrig) = (1,1,1)

function AsympOrigmmf(nlp,x,tg,MMD)
    x1 = x[1]
    tg[1] =  1 - (ex = exp(-nlp[1]*x1))
    MMD[1,1] = x1*ex
end

function initpars(m::AsympOrig)
    y = m.y
    mn,mx = extrema(y)
    A0 = mx + (mx-mn)/4.
    [abs(mean(log(1 - y/A0) ./ vec(m.x)))]
end

### Asymptotic Regression model expressed with an offset in x
formula(::AsympOff) = :(Asym*(1-exp(-rc*(x-c0))))
pnames(::AsympOff) = [:Asym,:c0,:rc,:x]
dim(::AsympOff) = (1,2,1)

function AsympOffmmf(nlp,x,tg,MMD)
    tg[1] = 1 - (eKx = exp(-(K = nlp[2])*(x1 = x[1] - nlp[1])))
    MMD[1,1] = -K * eKx
    MMD[2,1] = x1 * eKx
end

function initpars(m::AsympOff)
    cc = coef(fit(AsympReg(vec(m.x),m.y))) # fit the model with AsympReg pars
    rc = cc[3]; A = cc[1]
    [-log(A/(A-cc[2]))/rc, rc]
end

### Bolus single dose in measured compartment
formula(::BolusSD1) = :(V*exp(-K*x))
pnames(::BolusSD1) = [:V,:K,:x]
dim(::BolusSD1) = (1,1,1)

function BolusSD1mmf(nlp,x,tg,MMD)
    x1 = x[1]
    nKx1 = -nlp[1] * x1            # negative K*x[1] where nlp[1] is K
    MMD[1,1] = -x1 * (tg[1] = exp(nKx1))
end

function initpars{T<:FP}(m::BolusSD1{T})
    (n = length(m.y)) < 2 && return [-one(T)]
    cc = linreg(vec(m.x[1,:]),log(m.y))
    [cc[2] < 0. ? -cc[2] : -one(T)]
end

### 3-parameter Logistic
formula(::Logis3P) = :(Asym/(1 + exp(-(x-xmid)/scal)))
pnames(::Logis3P) = [:Asym,:xmid,:scal,:x]
dim(::Logis3P) = (1,2,1)

function Logis3Pmmf(nlp,x,tg,MMD)
    scal = nlp[2]
    nd = nlp[1] - x[1]                # negative difference from xmid
    ed = exp(nd/scal)                 # exp of standardized difference
    oped = 1 + ed
    tg[1] = 1/oped
    MMD[2,1] = -nd*(MMD[1,1] = -(ed/scal)/abs2(oped))/scal
end

function initpars(m::Logis3P)
    z = copy(m.y)
    if (minz = minimum(z)) < 0.       # ensure minimum(z) is positive
        z -= 1.05 * minz
    end
    z /= (1.05 * maximum(z))          # all z values in (0,1)
    linreg(log(z ./ (one(eltype(z)) - z)),vec(m.x))
end

### 4-parameter Logistic
formula(::Logis4P) = :(A + (B-A)/(1 + exp(-(x-xmid)/scal)))
pnames(::Logis4P) = [:A,:B,:xmid,:scal,:x]
dim(::Logis3P) = (2,2,1)

function Logis4Pmmf(nlp,x,tg,MMD)
    scal = nlp[2]
    nd = nlp[1] - x[1]                # negative difference from xmid
    ed = exp(nd/scal)                 # exp of standardized difference
    oped = 1 + ed
    tg[1] = ed*(tg[2] = 1/oped)
    MMD[1,2] = -(bb = MMD[1,1] = ed/(scal * abs2(oped)))
    MMD[2,1] = -(MMD[2,2] = nd * bb/scal)
end

function initpars(m::Logis4P)
    y = m.y
    x = vec(m.x)
    rg = range(y)
    z = (y - minimum(y) + rg/20)/(1.1rg)
    linreg(log(z ./ (1 - z)),x)
end

### Biexponential
formula(::Biexp) = :(A1 * exp(-K1*x) + A2 * exp(-K2*x))
pnames(::Biexp) = [:A1,:A2,:K1,:K2,:x]
dim(::Biexp) = (2,2,1)
function Biexpmmf(nlp,x,tg,MMD)
    x1 = x[1]
    tg[1] = eK1 = exp(-nlp[1]*x1)
    tg[2] = eK2 = exp(-nlp[2]*x1)
    MMD[1,1] = -x1 * eK1
    MMD[2,2] = -x1 * eK2
end

function initpars(m::Biexp)
    x = copy(vec(m.x))
    perm = sortperm(x)
    x = x[perm]
    y = copy(m.y)[perm]
    n = length(y)
    nlast = max(3, n>>1)
    ll = (n - nlast + 1):n      # "last half" of data but min of 3 els
    cc = linreg(x[ll],log(y[ll]))
    rc2 = abs(cc[2])
    [abs(linreg(x,log(abs(y - exp(cc[1] - rc2 * x))))[2]),rc2]
end

### Gompertz

formula(::Gompertz) = :(Asym*exp(-b2*b3^x))
pnames(::Gompertz) = [:Asym,:a,:b,:x]
dim(::Gompertz) = (1,2,1)
function Gompertzmmf(nlp,x,tg,MMD)
    x1 = x[1]
    a = nlp[1]
    b = nlp[2]
    p1 = b^x1
    tg[1] = e1 = exp(-a*p1)
    MMD[1,1] = -p1 * e1
    MMD[2,1] = -e1 * a * b^(x1 - 1) * x1
end

formula(::Chwirut) = :(R0*exp(-rc*x)/1+m*x)
pnames(::Chwirut) = [:R0,:rc,:m,:x]
dim(::Chwirut) = (1,2,1)
function Chwirutmmf(nlp,x,tg,MMD)
    x1 = x[1]
    et = exp(-nlp[1]*x1)
    denom = 1 + nlp[2] * x1
    tg[1] = et/denom
    etx = et * x1
    MMD[2,1] = (MMD[1,1] = -etx/denom)/denom
end

initpars(m::Chwirut) = [-(linreg(vec(m.x),log(m.y))[2]), mean(m.x)]

immutable Logsd1{T<:FP} <: NLregModF{T}
    x::Matrix{T}
    y::Vector{T}
    mu::Vector{T}
    resid::Vector{T}
    tgrad::Matrix{T}
    f::Function
end
function Logsd1f(p::StridedVector,x::StridedVector,tg::StridedVector)
    V = exp(p[1]); nKx1 = -exp(p[2])*x[1] # negative K * x[1]
    tg[2] = nKx1 * (mm = V*exp(nKx1))
    tg[1] = V*mm
    mm
end
function Logsd1{T<:FP}(t::Vector{T},y::Vector{T})
    n = length(y); length(t) == n || error("Dimension mismatch")
    Logsd1(reshape(t,(1,n)),y,similar(y),similar(y),Array(T,(2,n)),Logsd1f)
end
Logsd1{T<:FP}(t::DataArray{T,1},y::DataArray{T,1}) = Logsd1(vector(t),vector(y))
function Logsd1(f::Formula,dat::AbstractDataFrame)
    mf = ModelFrame(f,dat)
    mm = ModelMatrix(mf)
    Logsd1(mm.m[:,end],model_response(mf))
end
Logsd1(ex::Expr,dat::AbstractDataFrame) = Logsd1(Formula(ex),dat)

function initpars{T<:FP}(m::Logsd1{T})
    (n = length(m.y)) < 2 && return [zero(T),-one(T)]
    cc = hcat(ones(n),vec(m.x[1,:]))\log(m.y)
    [cc[2] < 0. ? [cc[1],log(-cc[2])] : [cc[1],-one(T)]]
end

pnames(m::Logsd1) = [:logV,:logK,:x]