text:
function rrule(::typeof(MPS), x::Vector{<:ITensor}; kwargs...)
  y = MPS(x; kwargs...)
  function MPS_pullback(ȳ)
    ȳtensors = ȳ.data
    n = length(ȳtensors)
    envL = [ȳtensors[1] * dag(x[1])]
    envR = [ȳtensors[n] * dag(x[n])]
    for j in 2:(n - 1)
      push!(envL, envL[j - 1] * ȳtensors[j] * dag(x[j]))
      push!(envR, envR[j - 1] * ȳtensors[n + 1 - j] * dag(x[n + 1 - j]))
    end

    x̄ = ITensor[]
    push!(x̄, ȳtensors[1] * envR[n - 1])
    for j in 2:(n - 1)
      push!(x̄, envL[j - 1] * ȳtensors[j] * envR[n - j])
    end
    push!(x̄, envL[n - 1] * ȳtensors[n])
    return (NoTangent(), x̄)
  end
  return y, MPS_pullback
end

function rrule(::typeof(inner), x1::MPS, x2::MPS; kwargs...)
  if !hassameinds(siteinds, x1, x2)
    error(
      "Taking gradients of `inner(::MPS, ::MPS)` is not supported if the site indices of the input MPS don't match. If you input `inner(x, Ay)` where `Ay` is the result of something like `contract(A::MPO, y::MPS)`, try `inner(x', Ay)` or `inner(x, replaceprime(Ay, 1 => 0))`instead.",
    )
  end
  y = inner(x1, x2)
  function inner_pullback(ȳ)
    x̄1 = dag(ȳ) * x2
    # `dag` of `x1` gets reversed by `inner`
    x̄2 = x1 * ȳ
    return (NoTangent(), x̄1, x̄2)
  end
  return y, inner_pullback
end

# TODO: Define a more general version in ITensors.jl
function _contract(::Type{ITensor}, ψ::Union{MPS,MPO}, ϕ::Union{MPS,MPO}; kwargs...)
  T = ITensor(1)
  for n in 1:length(ψ)
    T = T * ψ[n] * ϕ[n]
  end
  return T
end

function _contract(::Type{MPO}, ψ::MPS, ϕ::MPS; kwargs...)
  ψmat = convert(MPO, ψ)
  ϕmat = convert(MPO, ϕ)
  return contract(ψmat, ϕmat; kwargs...)
end

function rrule(::typeof(apply), x1::Vector{ITensor}, x2::Union{MPS,MPO}; kwargs...)
  N = length(x1) + 1
  apply_dag = x2 isa MPO ? kwargs[:apply_dag] : nothing

  # Apply circuit and store intermediates in the forward direction
  x1x2 = Vector{typeof(x2)}(undef, N)
  x1x2[1] = x2
  for n in 2:N
    x1x2[n] = apply(x1[n - 1], x1x2[n - 1]; move_sites_back=true, kwargs...)
  end
  y = x1x2[end]

  function apply_pullback(ȳ)
    x1x2dag = dag.(x1x2)
    x1dag = [swapprime(dag(x), 0 => 1) for x in x1]

    # Apply circuit and store intermediates in the reverse direction
    x1dag_ȳ = Vector{typeof(x2)}(undef, N)
    x1dag_ȳ[end] = ȳ
    for n in (N - 1):-1:1
      x1dag_ȳ[n] = apply(x1dag[n], x1dag_ȳ[n + 1]; move_sites_back=true, kwargs...)
    end

    x̄1 = similar(x1)
    for n in 1:length(x1)
      # check if it's not a noisy gate (rank-3 tensor)
      if iseven(length(inds(x1[n])))
        gateinds = inds(x1[n]; plev=0)
        if x2 isa MPS
          ξ̃ = prime(x1dag_ȳ[n + 1], gateinds)
          ϕ̃ = x1x2dag[n]
        else
          # apply U on one side of the MPO
          if apply_dag
            ϕ̃ = swapprime(x1x2dag[n], 0 => 1)
            ϕ̃ = apply(x1[n], ϕ̃; move_sites_back=true, apply_dag=false)
            ϕ̃ = mapprime(ϕ̃, 1 => 2, 0 => 1)
            ϕ̃ = replaceprime(ϕ̃, 1 => 0; inds=gateinds')
            ξ̃ = 2 * dag(x1dag_ȳ[n + 1])'
          else
            ϕ̃ = mapprime(x1x2dag[n], 0 => 2)
            ϕ̃ = replaceprime(ϕ̃, 1 => 0; inds=gateinds')
            ξ̃ = mapprime(x1dag_ȳ[n + 1], 0 => 2)
          end
        end
        x̄1[n] = _contract(ITensor, ξ̃, ϕ̃)
      else
        s = inds(x1[n])
        x̄1[n] = itensor(zeros(dim.(s)), s...)
      end
    end
    x̄2 = x1dag_ȳ[end]
    return (NoTangent(), x̄1, x̄2)
  end
  return y, apply_pullback
end


text:
function log_timeseries!(lf::Logfile, step::Int; args...)
    log_timeseries!(lf, step, args)
end

function log_timeseries!(lf::Logfile, step::Int, args)
    h5open(lf.path, "r+") do root
        # create "/timeseries" super group if it does not yet exist
        supergroup = g_create_or_get(root, GROUP_TIMESERIES)
        # compute time stamp for current entry
        time = div(Int(Dates.value(now())), 1000)
        for (name_symb, value) in args
            @assert value isa Number
            name = string(name_symb)
            # check if the "/timeseries/$name" group for current
            # series already exists or not
            subgroup = g_create_or_get(supergroup, name)
            # push the values to the data arrays
            d_create_and_push(subgroup, DATA_TIMESERIES_TIME, time)
            d_create_and_push(subgroup, DATA_TIMESERIES_STEPS, step)
            d_create_and_push(subgroup, DATA_TIMESERIES_VALUE, value)
        end
    end
end

ext: jl


text:
module day15

function play_game(input, turns)
    seen = Dict(x=>i for (i,x) in enumerate(input[1:end-1]))
    lastnum = input[end]
    for i in length(input):(turns-1)
        next = lastnum in keys(seen) ? i - seen[lastnum] : 0
        seen[lastnum] = i
        lastnum = next
    end
    lastnum
end

part1(input) = play_game(input, 2020)
part2(input) = play_game(input, 30_000_000)

@assert part1([0,3,6]) == 436
@assert part1([3,1,2]) == 1836
@assert part2([3,1,2]) == 362

input = [6,4,12,1,20,0,16]
@info "Day 15, Part 1 answer: $(part1(input))"
@info "Day 15, Part 2 answer: $(part2(input))"

end


text:
export encode_String, decode_String,
       encode_string, decode_string,
       encode_string_array,
       Guid, Guids,
       encode_Guid, decode_Guid,
       encode_Guid_array, decode_Guid_array,
       is_empty_guid,
       encode_bool, decode_bool,
       encode_byte, decode_byte,
       encode_int, decode_int, decode_int_or_error, decode_int_or_nothing,
       encode_long, decode_long, decode_long_or_error,
       encode_int_array, decode_int_array, decode_int_or_error_array,
       encode_float, decode_float,
       encode_float_array, decode_float_array,
       encode_double, decode_double,
       encode_double_array, decode_double_array,
       encode_Point3d, decode_Point3d,
       encode_Point3d_array, decode_Point3d_array,
       encode_Point2d, decode_Point2d,
       encode_Point2d_array, decode_Point2d_array,
       encode_Vector3d, decode_Vector3d,
       encode_Frame3d, decode_Frame3d,
       decode_void,
       encode_object, decode_object,
       encode_object_array, decode_object_array

encode_String(c::IO, v::String) = begin
  str = string(v)
  size = length(str)
  array = UInt8[]
  while true
    byte = size & 0x7f
    size >>= 7
    if size > 0
      push!(array, byte | 0x80)
    else
      push!(array, byte)
      break
    end
  end
  write(c, array)
  write(c, str)
end

decode_String(c::IO) = begin
  loop(size::Int, shift::Int) = begin
    b = convert(Int, read(c, UInt8))
    size = size | ((b & 0x7f) << shift)
    if (b & 0x80) == 0
      String(read(c, size))
    else
      loop(size, shift + 7)
    end
  end
  loop(0, 0)
end

# C# uses two different names for strings
encode_string = encode_String
decode_string = decode_String

Guid = Vector{UInt8}
Guids = Vector{Guid}

is_empty_guid(g::Guid) = all(v -> v == 0, g)

encode_Guid(c::IO, v::Guid) = write(c, v)
decode_Guid(c::IO) =
  let guid = read(c, 16)
      if is_empty_guid(guid)
        error("Backend Error: $(decode_String(c))")
      else
        guid
      end
  end

encode_Guid_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_Guid(c, e) end
end
decode_Guid_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{Guid}(undef, len)
  for i in 1:len
    r[i] = decode_Guid(c)
  end
  r
end

encode_string = encode_String
decode_string = decode_String

encode_bool(c::IO, v::Bool) = encode_byte(c, v ? UInt8(1) : UInt8(0))
decode_bool(c::IO) =
  let i = decode_byte(c)
    if i == 127
      error("Backend Error: $(decode_String(c))")
    else
      i == 1
    end
  end

encode_byte(c::IO, v::UInt8) = write(c, convert(UInt8, v))
decode_byte(c::IO) = convert(UInt8, read(c, UInt8))

encode_int(c::IO, v::Int) = write(c, convert(Int32, v))
decode_int(c::IO) = convert(Int, read(c, Int32))

encode_long(c::IO, v::Int) = write(c, convert(Int64, v))
decode_long(c::IO) = convert(Int, read(c, Int64))

encode_double(c::IO, v::Real) = write(c, convert(Float64, v))
decode_double(c::IO) =
    let d = read(c, Float64)
        if isnan(d)
            error("Backend Error: $(decode_String(c))")
        else
            d
        end
    end
encode_double3(c::IO, v0::Real, v1::Real, v2::Real) = begin
    encode_double(c, v0)
    encode_double(c, v1)
    encode_double(c, v2)
end

encode_float(c::IO, v::Real) = write(c, convert(Float32, v))
decode_float(c::IO) =
    let d = read(c, Float32)
        if isnan(d)
            error("Backend Error: $(decode_String(c))")
        else
            convert(Float64, d)
        end
    end
encode_float3(c::IO, v0::Real, v1::Real, v2::Real) = begin
    encode_float(c, v0)
    encode_float(c, v1)
    encode_float(c, v2)
end

decode_int_or_error_numbered(err_num) = (c::IO) ->
  let i = decode_int(c)
    if i == err_num
      error("Backend Error: $(decode_String(c))")
    else
      i
    end
  end

decode_int_or_error = decode_int_or_error_numbered(-1)

decode_int_or_nothing(c::IO) =
  let i = decode_int_or_error(c)
    if i == -2
      nothing
    else
      i
    end
  end

decode_long_or_error_numbered(err_num) = (c::IO) ->
  let i = decode_long(c)
    if i == err_num
      error("Backend Error: $(decode_String(c))")
    else
      i
    end
  end

decode_long_or_error = decode_long_or_error_numbered(-1)

encode_string_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_string(c, e) end
end

encode_double_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_double(c, e) end
end
decode_double_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{Float64}(undef, len)
  for i in 1:len
    r[i] = decode_double(c)
  end
  r
end

encode_int_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_int(c, e) end
end
decode_int_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{Int}(undef, len)
  for i in 1:len
    r[i] = decode_int(c)
  end
  r
end

decode_int_or_error_numbered_array(err_num) =
    let decode_int_or_err_num = decode_int_or_error_numbered(err_num)
        (c::IO) -> begin
            len = decode_int_or_error(c)
            r = Vector{Int}(undef, len)
            for i in 1:len
                r[i] = decode_int_or_err_num(c)
            end
            r
        end
    end

decode_int_or_error_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{Int}(undef, len)
  for i in 1:len
    r[i] = decode_int_or_error(c)
  end
  r
end

# Must convert from local to world coordinates
encode_Point3d(c::IO, v::XYZ) = begin
  v = in_world(v)
  encode_double3(c, v.x, v.y, v.z)
end
decode_Point3d(c::IO) = xyz(decode_double(c), decode_double(c), decode_double(c), world_cs)
encode_Point3d_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_Point3d(c, e) end
end
decode_Point3d_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{XYZ}(undef, len)
  for i in 1:len
    r[i] = decode_Point3d(c)
  end
  r
end
encode_Point2d(c::IO, v::XYZ) = begin
  v = in_world(v)
  @assert v.z == zero(v.z) "Non zero Z coordinate";
  encode_double(c, v.x)
  encode_double(c, v.y)
end
decode_Point2d(c::IO) = xyz(decode_double(c), decode_double(c), 0.0, world_cs)
encode_Point2d_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v
    encode_Point2d(c, e)
  end
end
decode_Point2d_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{XYZ}(undef, len)
  for i in 1:len
    r[i] = decode_Point2d(c)
  end
  r
end

encode_Vector3d(c::IO, v::VXYZ) = begin
  v = in_world(v)
  encode_double3(c, v.x, v.y, v.z)
end
decode_Vector3d(c::IO) = vxyz(decode_double(c), decode_double(c), decode_double(c), world_cs)

encode_Vector3d_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_Vector3d(c, e) end
end
decode_Vector3d_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{VXYZ}(undef, len)
  for i in 1:len
    r[i] = decode_Vector3d(c)
  end
  r
end

encode_Frame3d(c::IO, v::XYZ) = begin
  encode_Point3d(c, v)
  t = v.cs.transform
  encode_double3(c, t[1,1], t[2,1], t[3,1])
  encode_double3(c, t[1,2], t[2,2], t[3,2])
  encode_double3(c, t[1,3], t[2,3], t[3,3])
end

decode_Frame3d(c::IO) =
    u0(cs_from_o_vx_vy_vz(decode_Point3d(c), decode_Vector3d(c), decode_Vector3d(c), decode_Vector3d(c)))

decode_void(c::IO) = begin
  v = decode_byte(c)
  if v == 127
    error("Backend Error: $(decode_String(c))")
  end
end

encode_object(c::IO, v::Bool) = (encode_byte(c, 0x0); encode_bool(c, v))
encode_object(c::IO, v::UInt8) = (encode_byte(c, 0x1); encode_byte(c, v))
encode_object(c::IO, v::Int32) = (encode_byte(c, 0x2); encode_int(c, v))
encode_object(c::IO, v::Int64) = (encode_byte(c, 0x3); encode_long(c, v))
encode_object(c::IO, v::Float32) = (encode_byte(c, 0x4); encode_float(c, v))
encode_object(c::IO, v::Float64) = (encode_byte(c, 0x5); encode_double(c, v))

encode_object_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_object(c, e) end
end

create_backend_connection(backend::String, port::Integer) =
  for i in 1:10
    try
      return connect(port)
    catch e
      @info("Please, start/restart $(backend).")
      sleep(5)
      if i == 9
        throw(e)
      end
    end
  end


text:
using NDTensors, Test

@static if VERSION >= v"1.5"
  using Pkg
  Pkg.add("Octavian")
  using Octavian
end

@testset "Dense Tensors" begin
  @testset "DenseTensor basic functionality" begin
    A = Tensor(3, 4)

    for I in eachindex(A)
      @test A[I] == 0
    end

    randn!(A)

    for I in eachindex(A)
      @test A[I] != 0
    end

    for I in eachindex(A)
      @test A[I] != 0
    end

    @test ndims(A) == 2
    @test dims(A) == (3, 4)
    @test inds(A) == (3, 4)

    A[1, 1] = 11

    @test A[1, 1] == 11

    Aview = A[2:3, 2:3]

    @test dims(Aview) == (2, 2)
    @test A[2, 2] == Aview[1, 1]

    @test A * 2.0 == 2.0 * A

    Asim = similar(data(A), 10)
    @test eltype(Asim) == Float64
    @test length(Asim) == 10

    B = Tensor(undef, 3, 4)
    randn!(B)

    C = A + B

    for I in eachindex(C)
      @test C[I] == A[I] + B[I]
    end

    Ap = permutedims(A, (2, 1))

    for I in eachindex(A)
      @test A[I] == Ap[NDTensors.permute(I, (2, 1))]
    end

    t = Tensor(ComplexF64, 100, 100)
    randn!(t)
    @test conj(data(store(t))) == data(store(conj(t)))
    @test typeof(conj(t)) <: DenseTensor

    @test Dense(ComplexF64) == Dense{ComplexF64}()
    @test Dense(ComplexF64) == complex(Dense(Float64))

    D = Tensor(ComplexF64, (100, 100))
    @test eltype(D) == ComplexF64
    @test ndims(D) == 2
    @test dim(D) == 100^2

    E = Tensor(ComplexF64, undef, (100, 100))
    @test eltype(E) == ComplexF64
    @test ndims(E) == 2
    @test dim(E) == 100^2

    F = Tensor((100, 100))
    @test eltype(F) == Float64
    @test ndims(F) == 2
    @test dim(F) == 100^2

    G = Tensor(undef, (100, 100))
    @test eltype(G) == Float64
    @test ndims(G) == 2
    @test dim(G) == 100^2

    H = Tensor(ComplexF64, undef, 100, 100)
    @test eltype(H) == ComplexF64
    @test ndims(H) == 2
    @test dim(H) == 100^2

    I_arr = rand(10, 10, 10)
    I = Tensor(I_arr, (10, 10, 10))
    @test eltype(I) == Float64
    @test dim(I) == 1000
    @test Array(I) == I_arr

    J = Tensor(2, 2)
    K = Tensor(2, 2)
    @test Array(J * K) ≈ Array(J) * Array(K)
  end

  @testset "Random constructor" begin
    T = randomTensor(2, 2)
    @test dims(T) == (2, 2)
    @test eltype(T) == Float64
    @test T[1, 1] ≉ 0
    @test norm(T) ≉ 0

    Tc = randomTensor(ComplexF64, 2, 2)
    @test dims(Tc) == (2, 2)
    @test eltype(Tc) == ComplexF64
    @test Tc[1, 1] ≉ 0
    @test norm(Tc) ≉ 0
  end

  @testset "Complex Valued Tensors" begin
    d1, d2, d3 = 2, 3, 4
    T = randomTensor(ComplexF64, d1, d2, d3)

    rT = real(T)
    iT = imag(T)
    cT = conj(T)

    for n1 in 1:d1, n2 in 1:d2, n3 in 1:d3
      @test rT[n1, n2, n3] ≈ real(T[n1, n2, n3])
      @test iT[n1, n2, n3] ≈ imag(T[n1, n2, n3])
      @test cT[n1, n2, n3] ≈ conj(T[n1, n2, n3])
    end
  end

  @testset "Custom inds types" begin
    struct MyInd
      dim::Int
    end
    NDTensors.dim(i::MyInd) = i.dim

    T = Tensor((MyInd(2), MyInd(3), MyInd(4)))
    @test store(T) isa Dense
    @test eltype(T) == Float64
    @test norm(T) == 0
    @test dims(T) == (2, 3, 4)
    @test ndims(T) == 3
    @test inds(T) == (MyInd(2), MyInd(3), MyInd(4))
    T[2, 1, 2] = 1.21
    @test T[2, 1, 2] == 1.21
    @test norm(T) == 1.21

    T = randomTensor(ComplexF64, (MyInd(4), MyInd(3)))
    @test store(T) isa Dense
    @test eltype(T) == ComplexF64
    @test norm(T) > 0
    @test dims(T) == (4, 3)
    @test ndims(T) == 2
    @test inds(T) == (MyInd(4), MyInd(3))

    T2 = 2 * T
    @test eltype(T2) == ComplexF64
    @test store(T2) isa Dense
    @test norm(T2) > 0
    @test norm(T2) / norm(T) ≈ 2
    @test dims(T2) == (4, 3)
    @test ndims(T2) == 2
    @test inds(T2) == (MyInd(4), MyInd(3))
  end

  @testset "generic contraction" begin
    # correctness of _gemm!
    for alpha in [0.0, 1.0, 2.0]
      for beta in [0.0, 1.0, 2.0]
        for tA in ['N', 'T']
          for tB in ['N', 'T']
            A = randn(4, 4)
            B = randn(4, 4)
            C = randn(4, 4)
            A = BigFloat.(A)
            B = BigFloat.(B)
            C2 = BigFloat.(C)
            NDTensors._gemm!(tA, tB, alpha, A, B, beta, C)
            NDTensors._gemm!(tA, tB, alpha, A, B, beta, C2)
            @test C ≈ C2
          end
        end
      end
    end
  end

  @testset "Contraction with size 1 block and NaN" begin
    @testset "No permutation" begin
      R = Tensor(ComplexF64, 2, 2, 1)
      fill!(R, NaN)
      @test any(isnan, R)
      T1 = randomTensor(2, 2, 1)
      T2 = randomTensor(ComplexF64, 1, 1)
      NDTensors.contract!(R, (1, 2, 3), T1, (1, 2, -1), T2, (-1, 1))
      @test !any(isnan, R)
      @test convert(Array, R) ≈ convert(Array, T1) * T2[1]
    end

    @testset "Permutation" begin
      R = Tensor(ComplexF64, 2, 2, 1)
      fill!(R, NaN)
      @test any(isnan, R)
      T1 = randomTensor(2, 2, 1)
      T2 = randomTensor(ComplexF64, 1, 1)
      NDTensors.contract!(R, (2, 1, 3), T1, (1, 2, -1), T2, (-1, 1))
      @test !any(isnan, R)
      @test convert(Array, R) ≈ permutedims(convert(Array, T1), (2, 1, 3)) * T2[1]
    end
  end

  @testset "change backends" begin
    a, b, c = [randn(5, 5) for i in 1:3]
    backend_auto()
    @test NDTensors.gemm_backend[] == :Auto
    @test NDTensors.auto_select_backend(typeof.((a, b, c))...) ==
          NDTensors.GemmBackend(:BLAS)
    res1 = NDTensors._gemm!('N', 'N', 2.0, a, b, 0.2, copy(c))
    backend_blas()
    @test NDTensors.gemm_backend[] == :BLAS
    res2 = NDTensors._gemm!('N', 'N', 2.0, a, b, 0.2, copy(c))
    backend_generic()
    @test NDTensors.gemm_backend[] == :Generic
    res3 = NDTensors._gemm!('N', 'N', 2.0, a, b, 0.2, copy(c))
    @test res1 == res2
    @test res1 ≈ res3
    backend_auto()
  end

  @static if VERSION >= v"1.5"
    @testset "change backends" begin
      a, b, c = [randn(5, 5) for i in 1:3]
      backend_auto()
      @test NDTensors.gemm_backend[] == :Auto
      @test NDTensors.auto_select_backend(typeof.((a, b, c))...) ==
            NDTensors.GemmBackend(:BLAS)
      res1 = NDTensors._gemm!('N', 'N', 2.0, a, b, 0.2, copy(c))
      backend_octavian()
      @test NDTensors.gemm_backend[] == :Octavian
      res4 = NDTensors._gemm!('N', 'N', 2.0, a, b, 0.2, copy(c))
      @test res1 ≈ res4
      backend_auto()
    end
  end
end

nothing


text:
"""
this will be invoked after first pass kernel so we already have a working queue  and we will work on it
we will finish when mainActiveCounterNext will evaluate 0 at the end of iteration
as we will looop here intensively and after each iteration we will need to synchronize all blocks we will use cooperative goups and sync grid function

    reducedArrays =  [reducedSegmA,reducedGoldB,reducedGoldA,reducedSegmB]
        - list of 3 dimensional boolean arrays with main data that we will work on

    metaData - basic data about data blocks - stored in 3 dimensional array
    metadataDims - dimensions of metaData array
    resArrays = [resArrayA, resArrayB]- 3 dimensional array of size of reduced source array where entry will be in place where we had covered the voxel, and value will mark in which itwration
        - we have two result arrays separate for each pass
        - when we will cover some point  we need to report it in the resArray
    resArraysCounters= [resArrayAcounter, resArrayBcounter] - stores the numbers of results related to dilatation of each array - those will be needed to establish that
                            - we can do early finish of the algorithm

    datablockDim - the edge length  describing the block of data we are working on  - defoult is 32 - which would produce blocks of size 32x32xzDim
        - we will have  the plane of threads of size 32x32 and in order to preserve memory coalescence we should have a warp have the  set of next x positions

    IMPORTANT we will alternate  between A nd B entries  and one will be used as working queue and other as a queue for next iteration
    mainQuesCounterArr=[mainQuesCounterA,mainQuesCounterB] -array of two counters that we will update atomically and will be usefull to populate the work queue
    mainWorkQueueArr = [mainWorkQueueA,mainWorkQueueB]-2 element array of lists of the indicies of  data blocks in metadata with additional information is it referencing the goldpass or second one

    isSecondCounter - Bool when true will mean that if true we should currently take data aout blocks wchich needs processed from second entry
                        of mainQuesCounterArr and  mainWorkQueueArr and write activated blocks to the first ones - in case it is false - we do it in opposite mainQuesCounter
    tailCounter - variable accessed atomically by each block (if needed) and on each access increased - we will finish as this tail counter will equal mainQuesCounter
    tailParts - number controlling how much of active blocks will be left to tail (part that is available to all blocks) - so the bigger this number the more of active blosk will go to
                tail queue
    numberOfThreadBlocks - how many thread blocks we can squeeeze on our GPU  at once
    isOddPass - boolean that will mark that we have an odd pass will be initialized to false as we had already first pass in other kernel
                    will be important to point out which        mainQuesCounter and mainWorkQueue is currently a source of data and where we should put active and activated blocks
    iterationNumber - initialized to 1 and increased atomically after every dilation cycle (this variable is global)
    fp,fn - amount of false positives and false negatives - used in order to be able to get early termination
                    """
module MainLoopIterations
using CUDA, ..CUDAGpuUtils, Logging,StaticArrays
using ..HFUtils

function mainLoopIterationsKernel(reducedArrays
                                    ,metaData
                                    ,metadataDims::Tuple{UInt8,UInt8,UInt8}
                                    ,resArrays
                                    ,resArraysCounters
                                    ,datablockDim
                                    ,mainQuesCounterArr
                                    ,mainWorkQueueArr
                                    ,tailCounter
                                    ,tailParts::UInt8
                                    ,numberOfThreadBlocks::UInt16
                                    ,isOddPass
                                    ,iterationNumber
                                    ,fp,fn)

    #needed to manage cooperative groups functions
    grid_handle = this_grid()
    #storing intermidiate results +2 in order to get the one padding
    resShmem =  @cuStaticSharedMem(Bool,(34,34,34))
    #storing values loaded from analyzed array ...
    sourceShmem =  @cuStaticSharedMem(Bool,(34,34,34))
    #for storing sums for reductions
    shmemSum =  @cuStaticSharedMem(Float32,35,14) # we need this additional 33th an 34th spots

    #coordinates of data in main array
    #we will use this to establish weather we should mark  the data block as empty or full ...
    isMaskFull= false
    #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
    locArr= Int32(0)
    offsetIter = UInt16(0)

    #we will store here the data about the blocks that we want to process
    toIterWorkQueueShmem =  @cuStaticSharedMem(UInt8,32,4)
    #indicates where we are in general work queue in given moment if we are iterating over part of work queue owned by this thread block
    positionInLocalWorkQueaue =  @cuStaticSharedMem(UInt16,1)
    #boolean usefull in the iterating over private part of work queue
    isAnyBiggerThanZero =  @cuStaticSharedMem(Bool,1)
    #used for iterating over a tail
    tailCounterInShmem= @cuStaticSharedMem(UInt16,1)
    workCounterInshmem= @cuStaticSharedMem(UInt16,1)

    #loading data to shared memory from global about is it even or odd pass
    isOddPassShmem =  @cuStaticSharedMem(Bool,1)
    iterationNumberShmem =  @cuStaticSharedMem(UInt16,1)
    #current spot in tail - usefull to save where we need to access the tail to get proper block
    currentTailPosition =  @cuStaticSharedMem(UInt16,1)
    #shared memory variables needed to marks wheather we are already finished with  any dilatation step
    goldToBeDilatated =  @cuStaticSharedMem(Bool,1)
    segmToBeDilatated =  @cuStaticSharedMem(Bool,1)
    #true when we have more than 0 blocks to analyze in next iteration
    workCounterBiggerThan0 =  @cuStaticSharedMem(Bool,1)

    #resetting
    @ifXY 1 1 tailCounterInShmem[1]=0
    @ifXY 2 1 workCounterInshmem[1]=0

    """
    just making easier to invoke it without passing arguments all the time
    """
    macro  innersingleDataBlockPass(ispassGoldd::Bool,currMatadataBlockX::UInt8,currMatadataBlockY::UInt8 ,currMatadataBlockZ::UInt8 )
        singleDataBlockPass(reducedArrays[ispassGoldd*2+1]
        ,reducedArrays[ispassGoldd*2+2]
        ,iterationNumberShmem[]
        ,((currMatadataBlockX-1)*32)+1
        ,((currMatadataBlockY-1)*32)+1
        ,((currMatadataBlockZ-1)*32)+1
        ,isMaskFull ,resShmem ,locArr ,metaData ,metadataDims ,isPassGold
        ,currMatadataBlockX
        ,currMatadataBlockY
        ,currMatadataBlockZ
        ,mainQuesCounterArr[!isOddPassShmem[]+1]# we negate those as we want to pass the data about new active blocks to queue of next pass
        ,mainWorkQueueArr[!isOddPassShmem[]+1]
        ,resArrays[ispassGoldd+1]
        ,resArraysCounters[ispassGoldd+1]  )

    loadDataAtTheBegOfDilatationStep(isOddPassShmem,iterationNumberShmem,iterationNumber,positionInMainWorkQueaue,workCounterInshmem,mainQuesCounterArr,isAnyBiggerThanZero,goldToBeDilatated,segmToBeDilatated, resArraysCounters  )
    sync_threads()
    #we check first wheather next dilatation step should be done or not we also establish some shared memory variables to know wheather both passes should continue or just one
    # checking weather we already finished so we need to check
    # - is amount of results related to gold mask dilatations is equal to false positives
    # - is amount of results related to other  mask dilatations is equal to false negatives
    # - is amount of workQueue that we will want to analyze now is bigger than 0
    while(goldToBeDilatated[1] && goldToBeDilatated[1] && workCounterBiggerThan0[1])
        while(isAnyBiggerThanZero[])
            loadOwnedWorkQueueIntoShmem(mainWorkQueue,mainQuesCounter,toIterWorkQueueShmem,positionInMainWorkQueaue ,numberOfThreadBlocks,tailParts)
            #at this point if we have anything in the private part of the  work queue we will have it in toIterWorkQueueShmem, in case  we will find some 0 inside it means that queue is exhousted and we need to loo into tail
            @unroll for i in UInt16(1):32# most outer loop is responsible for z dimension
                if(shmemIter[i,1]>0)
                    #we need to check also wheather given dilatation step is already finished - for example it is possible that it do not make sense to dilatate gold mask but still we need dilate other
                    if( (goldToBeDilatated[1]&&shmemIter[i,4]==1) || (segmToBeDilatated[1]&&shmemIter[i,4]==0)  )
                        innersingleDataBlockPass(shmemIter[i,4]==1,shmemIter[i,1],shmemIter[i,2] ,shmemIter[i,3] )
                    end
                    sync_threads()
                else    #at this point we got some 0 in the shmemIter
                    if(threadIdxX()==1)# sadly only threads with this id are managing the work queue
                        isAnyBiggerThanZero[]=false
                    end#inner if
                end#if
            end#for
        end#while

        # in this moment we have nothing left in private part of work queue and we need to check is there sth in tail to process
        #first load data of tail counter once

        loadDataNeededForTailAnalysisToShmem(currentTailPosition,tailCounter )
        sync_threads()
        while (currentTailPosition[1]< workCounterInshmem[1])
            #below we access tail of working queue in a way that will be atomic
            if(threadIdxY()==1 && threadIdxX()==1 )
                shmemIter[1]= mainWorkQueueArr[isOddPassShmem[]+1][currentTailPosition[1]]
            end
            sync_threads()
            #we need to check also wheather given dilatation step is already finished - for example it is possible that it do not make sense to dilatate gold mask but still we need dilate other
            if( (goldToBeDilatated[1]&&shmemIter[1,4]==1) || (segmToBeDilatated[1]&&shmemIter[1,4]==0)  )
                innersingleDataBlockPass(shmemIter[1,4]==1,shmemIter[1,1],shmemIter[1,2] ,shmemIter[1,3] )
            end
            loadDataNeededForTailAnalysisToShmem(currentTailPosition,tailCounter )
            sync_threads()
        end
        #we set it before sync grid - as we do not need it for this one
        if(threadIdxX()==10 && threadIdxY()==10)
            #we are just negiting it so it should be the same in all blocks
            isOddPassShmem[1]= !isOddPassShmem[1]
        end

        #if we are here it means we had covered all blocks that were marked as active and we need to prepare to next dilatation step
        sync_grid(grid_handle)
            #we clear  and add negation to !isOddPassShmem becouse we want to set the previously updated counter to 0
            clearBeforeNextDilatation(locArr,resShmem,mainQuesCounterArr[!isOddPassShmem[1]+1])
            prepareForNextDilatationStep(iterationNumber,tailCounter,numberOfThreadBlocks,tailParts,mainQuesCounterArr,isOddPassShmem)
            loadDataAtTheBegOfDilatationStep(isOddPassShmem,iterationNumberShmem,iterationNumber,positionInMainWorkQueaue,workCounterInshmem,mainQuesCounterArr,isAnyBiggerThanZero,goldToBeDilatated,segmToBeDilatated, resArraysCounters)

        sync_grid(grid_handle)
    end#while isNotFinished

end#mainLoopIterationsKernel

"""
checking weather we already finished so we need to check
    - is amount of results related to gold mask dilatations is equal to false positives
    - is amount of results related to other  mask dilatations is equal to false negatives
    - is amount of workQueue that we will want to analyze now is bigger than 0
"""
function isNotFinished(resArraysCounters,fp,fn,workCounterInshmem,goldToBeDilatated,segmToBeDilatated)::Bool

    resArraysCounters[2] == fp[1]
    resArraysCounters[1] == fn[1]

end#isNotFinished

"""
clearing data in order to enable their reuse in next iteration
"""
function clearBeforeNextDilatation(locArr,resShmem,oldBlockCounter)
    #resetting the counter so when we will add to this new items it will count from 0
    if(threadIdxX()==7 && threadIdxY()==7)
        oldBlockCounter[1]=0
    end

    #probably we do not need to clear work queue as we will just overwrite it
    #workQueue[((blockDimX()-1)*1024)+ (((threadIdxX()-1)*32)+1)+threadIdxY()]
    clearLocArr(locArr)
    clearMainShmem(resShmem)
    clearPadding(resShmem)
end#clearBeforeNextDilatation

"""
before next loop after grid sync we need to
    - increase iteration Number
    - set tail counter to correct value (1 bigger than last value from private queue of last block)
those should be invoked only once per grid hence it will  be invoked only in first thread block (first chosen arbitrarly)
"""
function prepareForNextDilatationStep(iterationNumber,tailCounter,numberOfThreadBlocks,tailParts,mainQuesCounterArr,isOddPassShmem)
    if(blockIdx().x==1)
        @IfYX 8 8 iterationNumber[1]= CUDA.atomic_inc!(pointer(iterationNumber), UInt16(1))+1
        @IfYX 9 9 tailCounter[1]=cld(mainQuesCounterArr[isOddPassShmem[1]+1],(numberOfThreadBlocks+tailParts))*numberOfThreadBlocks +1
    end#if
end#prepareForNextDilatationStep

"""
load data needed for tail analysis into shared memory
"""
function loadDataNeededForTailAnalysisToShmem(currentTailPosition,tailCounter )
    @IfYX 4 4 currentTailPosition[1]= CUDA.atomic_inc!(pointer(tailCounter), UInt16(1))+1
end#loadDataNeededForTailAnalysisToShmem

"""
loads data at the begining of each dilatation step
    we need to set some variables in shared memory ro initial values
"""
function loadDataAtTheBegOfDilatationStep(isOddPassShmem,iterationNumberShmem,iterationNumber,positionInMainWorkQueaue,workCounterInshmem,mainQuesCounterArr,isAnyBiggerThanZero,goldToBeDilatated,segmToBeDilatated, resArraysCounters  )

    @IfYX 2 2 iterationNumberShmem[1]= iterationNumber[1]
    @IfYX 3 3 positionInMainWorkQueaue[1]=0
    @IfYX 4 4 begin
        workCounterInshmem[1]= mainQuesCounterArr[isOddPassShmem[1]+1]
        workCounterBiggerThan0[1]= (workCounterInshmem[1]>0)
                    end
    @IfYX 5 5 isAnyBiggerThanZero[]=true
    @IfYX 6 6 goldToBeDilatated[1]=(resArraysCounters[2] < fp[1])
    @IfYX 7 7  segmToBeDilatated[1]=(resArraysCounters[1] < fn[1])

end

"""
load data from current queue (part related to this block) to shared memory  so we can iterate over effectively
    we can load no more than 32 data points at a time

    mainWorkQueue - data we are currently iterating over - list of the indicies of  data blocks in metadata with additional information is it referencing the goldpass or second one
    mainQuesCounter - wil give us information how many active blocks we have in our work queue
    positionInLocalWorkQueaue- indicates where we are in locla shared memory work queue in given moment if we are iterating over part of work queue owned by this thread block
    toIterWorkQueueShmem - we will store here the data about the blocks that we want to process
    numberOfThreadBlocks - how many thread blocks we can squeeeze on our GPU  at once
    tailParts - number controlling how much of active blocks will be left to tail (part that is available to all blocks) - so the bigger this number the more of active blosk will go to
                tail queue

"""
function loadOwnedWorkQueueIntoShmem(mainWorkQueue
                                    ,mainQuesCounter
                                    ,toIterWorkQueueShmem
                                    ,positionInMainWorkQueaue
                                    ,numberOfThreadBlocks
                                    ,tailParts)

 privatePartLength=cld(mainQuesCounter[1],(numberOfThreadBlocks+tailParts))
 localPosition= blockIdx().x *privatePartLength + positionInMainWorkQueaue[1]+threadIdxY()
    #TODO() check for coalescence
    #as we want to load only 32 at a time we use 1 warp
    if(threadIdxX()==1)
        toIterWorkQueueShmem[threadIdxY(),:]= [UInt8(0),UInt8(0),UInt8(0),UInt8(0)]# so it is 0 and if we will not go through loop below it will remain 0
    end
    if(threadIdxX()==1 &&  localPosition < ((blockIdx().x+1) *privatePartLength) )
    toIterWorkQueueShmem[threadIdxY(),:]=mainWorkQueue[localPosition]#we add short 4 element arrays here
    end
    #now we need to mark that we loaded 32  elemnts

    positionInLocalWorkQueaue[1]+=32

end

"""
collects multiple functions that will be invoked over a single data blocks
analyzedArr - array that we dilatate
refAray - array we are comparing with
iterationNumber - in first pass it is 1
blockBeginingX,blockBeginingY,blockBeginingZ - coordinates where our block is begining - will be used as offset by our threads
isMaskFull,isMaskEmpty - needed for establishing state  of a current block after processing
resShmem- shared memory
locArr - thread local array with values loaded from global memory
metaData - 3 dim array with metadata of data blocks
metadataDims - dimensions of metadataDims
currMatadataBlockX,currMatadataBlockY, currMatadataBlockZ - cartesian coordinates of current block in metadaa!!!

IMPORTANT mainQuesCounter and mainWorkQueue will alternate as we will have 2 of them and every iteration we will swith between one to another
mainQuesCounter - counter that we will update atomically and will be usefull to populate the work queue
mainWorkQueue - the list of the indicies of  data blocks in metadata with additional information is it referencing the goldpass or second one
resArraysCounter - counter needed to add keep track on how many results we have in our result array
                    - we have separate counter for gold pass and other pass
"""
function singleDataBlockPass(analyzedArr
                ,refAray
                ,iterationNumber
                ,blockBeginingX
                ,blockBeginingY
                ,blockBeginingZ
                ,isMaskFull
                ,resShmem
                ,locArr
                ,metaData
                ,metadataDims::Tuple{UInt8,UInt8,UInt8}
                ,isPassGold::Bool
                ,currMatadataBlockX::UInt8
                ,currMatadataBlockY::UInt8
                ,currMatadataBlockZ::UInt8
                ,mainQuesCounter
                ,mainWorkQueue
                ,resArray
                ,resArraysCounter )

            #clear shared memory
            clearMainShmem(resShmem,32)
            clearPadding(resShmem,32)# we separately clear padding
            ############### execution
            executeDataIterOtherPasses(analyzedArr, refAray,iterationNumber,blockBeginingX,blockBeginingY,blockBeginingZ,isMaskFull,resShmem,locArr,resArraysCounter)
            #for futher processing we need to have space in main shmem
            clearMainShmem(shmem)
            #now we need to deal with padding in shmem res
            processAllPaddingPlanes(x,y,z,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArray,metaData,metadataDims,isPassGold,locArr,resArraysCounter)
            # now let's check weather block is eligible for futher processing - for this we need sums ...
            isActiveForNormalPass(isMaskFull, isMaskEmpty,resShmem,currMatadataBlockX,currMatadataBlockY,currMatadataBlockZ,isPassGold,metaData,mainQuesCounter,mainWorkQueue,resArraysCounter)
end#singleDataBlockPass

end


text:
##############################
#
# obj-Files
#
##############################

function load{MT <: AbstractMesh}(io::Stream{format"OBJ"}, MeshType::Type{MT}=GLNormalMesh)
    io           = stream(io)
    lineNumber   = 1
    Tv,Tn,Tuv,Tf = vertextype(MT), normaltype(MT), texturecoordinatetype(MT), facetype(MT)
    v,n,uv,f     = Tv[], Tn[], Tuv[], Tf[]
    last_command = ""
    attrib_type  = nothing
    for line in eachline(io)
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(line))
        !isvalid(line) && error("non valid ascii in obj")

        if !startswith(line, "#") && !isempty(line) && !iscntrl(line) #ignore comments
            lines   = split(line)
            command = shift!(lines) #first is the command, rest the data

            if "v" == command # mesh always has vertices
                push!(v, Point{3, Float32}(lines)) # should automatically convert to the right type in vertices(mesh)
            elseif "vn" == command && hasnormals(MT)
                push!(n, Normal{3, Float32}(lines))
            elseif "vt" == command && hastexturecoordinates(MT)
                length(lines) == 2 && push!(lines, "0.0") # length can be two, but obj normaly does three coordinates with the third equals 0.
                push!(uv, UVW{Float32}(lines))
            elseif "f" == command #mesh always has faces
                if any(x->contains(x, "/"), lines)
                    fs = process_face_uv_or_normal(lines)
                elseif any(x->contains(x, "//"), lines)
                    fs = process_face_normal(lines)
                else
                    push!(f, triangulated_faces(Tf, lines)...)
                    continue
                end
                push!(f, triangulated_faces(Tf, map(first, fs))...)
            else
                #TODO
            end
        end
        # read next line
        lineNumber += 1
    end
    attributes = Dict{Symbol, Any}()
    !isempty(f) && (attributes[:faces] = f)
    !isempty(v) && (attributes[:vertices] = v)
    if !isempty(n)
        if length(v) != length(n) # these are not per vertex normals, which we
            # can't deal with at the moment
            n = normals(v, f, Tn)
        end
        attributes[:normals] = n
    end
    !isempty(uv) && (attributes[:texturecoordinates] = uv)

    return MT(HomogenousMesh(attributes))
end

immutable SplitFunctor
    seperator::Compat.UTF8String
end
@compat (s::SplitFunctor)(array) = split(array, s.seperator)

# of form "f v1 v2 v3 ....""
process_face{S <: AbstractString}(lines::Vector{S}) = (lines,) # just put it in the same format as the others
# of form "f v1//vn1 v2//vn2 v3//vn3 ..."
process_face_normal{S <: AbstractString}(lines::Vector{S}) = map(SplitFunctor("//"), lines)
# of form "f v1/vt1 v2/vt2 v3/vt3 ..." or of form "f v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...."
process_face_uv_or_normal{S <: AbstractString}(lines::Vector{S}) = map(SplitFunctor("/"), lines)

immutable ParseFunctor{T}
    T::Type{T}
end
@compat (::ParseFunctor{T}){T}(x) = parse(T, x)

function triangulated_faces{Tf}(::Type{Tf}, vertex_indices::Vector{SubString{Compat.ASCIIString}})
    poly_face = Face{length(vertex_indices), UInt32, 0}(map(ParseFunctor(UInt32), vertex_indices))
    decompose(Tf, poly_face)
end


text:
using HelloWorld
using Base.Test

# write your own tests here
@test 1 == 1
@test 1 + 1 ==2
@test 2 == 2


text:
function discriminants(LDA::LinearDiscriminantModel{T}, X::Matrix{T}) where T
    dims = LDA.Θ.dims
    n = size(X, dims)
    m = LDA.Θ.m
    Δ = dims == 1 ? Matrix{T}(undef, n, m) : Matrix{T}(undef, m, n)
    return discriminants!(Δ, LDA, X)
end

function _posteriors!(Δ::Matrix{T}, LDA::LinearDiscriminantModel{T}) where T
    Π = broadcast!(exp, Δ, Δ)
    alt_dims = LDA.Θ.dims == 1 ? 2 : 1
    return broadcast!(/, Π, Π, sum(Π, dims=alt_dims))
end

function posteriors!(Π::Matrix{T}, LDA::LinearDiscriminantModel{T}, X::Matrix{T}) where T
    Δ = discriminants!(Π, LDA, X)
    return _posteriors!(Δ, LDA)
end

function posteriors(LDA::LinearDiscriminantModel{T}, X::Matrix{T}) where T
    Δ = discriminants(LDA, X)
    return _posteriors!(Δ, LDA)
end

function classify(LDA::LinearDiscriminantModel{T}, X::Matrix{T}) where T
    Π = posteriors(LDA, X)
    alt_dims = LDA.Θ.dims == 1 ? 2 : 1
    return [y[alt_dims] for y in argmax(Π, dims=alt_dims)]
end


text:
"""
Fast tokenization function.
"""
tokenize_for_conceptnet(doc::Vector{S}, splitter::Regex=DEFAULT_SPLITTER
                       ) where S<:AbstractString = begin doc end

tokenize_for_conceptnet(doc::S, splitter::Regex=DEFAULT_SPLITTER
                       ) where S<:AbstractString = begin
    # First, split
    tokens = strip.(split(doc, splitter))
    # Filter out empty strings
    filter!(!isempty, tokens)
end

"""
Retrieves the embedding matrix for a given `document`.
"""
function embed_document(conceptnet::ConceptNet{L,K,E},
                        document::AbstractString;
                        language=Languages.English(),
                        keep_size::Bool=true,
                        compound_word_separator::String="_",
                        max_compound_word_length::Int=1,
                        wildcard_matching::Bool=false,
                        print_matched_words::Bool=false
                       ) where {L<:Language, K, E<:Real}
    # Split document into tokens and embed
    return embed_document(conceptnet,
                          tokenize_for_conceptnet(document),
                          language=language,
                          keep_size=keep_size,
                          compound_word_separator=compound_word_separator,
                          max_compound_word_length=max_compound_word_length,
                          wildcard_matching=wildcard_matching,
                          print_matched_words=print_matched_words)
end

function embed_document(conceptnet::ConceptNet{L,K,E},
                        document_tokens::Vector{S};
                        language=Languages.English(),
                        keep_size::Bool=true,
                        compound_word_separator::String="_",
                        max_compound_word_length::Int=1,
                        wildcard_matching::Bool=false,
                        print_matched_words::Bool=false
                       ) where {L<:Language, K, E<:Real, S<:AbstractString}
    # Initializations
    embeddings = conceptnet.embeddings[language]
    # Get positions of words that can be used for indexing (found)
    # and those of words that can be searched (not_found)
    found_positions = token_search(conceptnet,
                                   document_tokens;
                                   language=language,
                                   separator=compound_word_separator,
                                   max_length=max_compound_word_length,
                                   wildcard_matching=wildcard_matching)
    # Get found words
    found_words = Vector{String}()
    for pos in found_positions
        word = make_word_from_tokens(document_tokens,
                                     pos,
                                     separator=compound_word_separator,
                                     separator_last=compound_word_separator)
        push!(found_words, word)
    end
    # Get best matches for not found words
    not_found_positions = setdiff(1:length(document_tokens),
                                  collect.(found_positions)...)
    # Return
    if print_matched_words
        words_not_found = document_tokens[not_found_positions]
        println("Embedded words: $found_words")
        println("Mismatched words: $words_not_found")
    end
    default = zeros(E, conceptnet.width)
    _embdoc = get(conceptnet.embeddings[language],
                  found_words,
                  default,
                  conceptnet.fuzzy_words[language],
                  n=conceptnet.width,
                  wildcard_matching=wildcard_matching)
    if keep_size
        embedded_document = hcat(_embdoc, zeros(E, conceptnet.width,
                                                length(not_found_positions)))
    else
        embedded_document = _embdoc
    end
    return embedded_document, not_found_positions
end

# Small function that builds a compound word
function make_word_from_tokens(tokens::Vector{S},
                               pos;
                               separator::String="_",
                               separator_last::String="_") where
        S<:AbstractString
    if length(pos) == 1
        return join(tokens[pos])
    else
        return join(tokens[pos], separator, separator_last)
    end
end

# Function that searches subphrases (continuous token combinations)
# from a document in a the embedded words and returns the positions of matched
# subphrases/words
# Example:
#   - for a vector: String[a, simpler, world, would, be, more, complicated],
#     max_length=7 and separator='_', it would generate:
#       String[a_simple_world_..._complicated,
#              a_simple_world_..._more,
#              ...
#              a_simple,
#              a,
#              simple_world_..._complicated,
#              simple_world_..._more,
#              ...
#              simple_world,
#              simple,
#              ...
#              ...
#              more_complicated,
#              complicated]
function token_search(conceptnet::ConceptNet{L,K,E},
                      tokens::Vector{S};
                      language::L=Languages.English(),
                      separator::String="_",
                      max_length::Int=3,
                      wildcard_matching::Bool=false) where
        {L<:Language, K, E<:Real, S<:AbstractString}
    # Initializations
    found = Vector{UnitRange{Int}}()
    n = length(tokens)
    i = 1
    j = n
    while i <= n
        if j-i+1 <= max_length
            token = join(tokens[i:j], separator, separator)
            is_match = !isempty(get(conceptnet[language], token,
                                    Vector{E}(),
                                    conceptnet.fuzzy_words[language],
                                    wildcard_matching=wildcard_matching))
            if is_match
                push!(found, i:j)
                i = j + 1
                j = n
                continue
            end
        end
        if i == j
            j = n
            i+= 1
        else
            j-= 1
        end
    end
    return found
end


text:
import GLAbstraction: bind, draw, color_attachment, depth_attachment

@render_program PeelingCompositingProgram
@render_program PeelingProgram
@render_program BlendProgram
@render_program InstancedPeelingProgram

@render_target ColorBlendTarget
@render_target PeelTarget

@vao PeelingVao
@instanced_vao InstancedPeelingVao

# Using the shared uploader system inside uploading.jl
PeelingUploader() = Uploader{PeelingProgram}()
InstancedPeelingUploader() = Uploader{InstancedPeelingProgram}()

@with_kw struct DepthPeelingRenderer <: AbstractRenderSystem
    num_passes::Int = 4
end

function Overseer.requested_components(::DepthPeelingRenderer)
    return (PeelingVao, PeelingProgram, InstancedPeelingVao, InstancedPeelingProgram,
            BlendProgram, PeelingCompositingProgram, CompositingProgram, ModelMat, Material,
            PointLight, UniformColor, BufferColor, Spatial, Camera3D, PeelTarget,
            ColorBlendTarget, IOTarget)
end

function Overseer.prepare(::DepthPeelingRenderer, dio::Diorama)
    if isempty(dio[BlendProgram])
        dio[Entity(1)] = BlendProgram(Program(blending_shaders()))
    end
    if isempty(dio[PeelingCompositingProgram])
        dio[Entity(1)] = PeelingCompositingProgram(Program(peeling_compositing_shaders()))
    end
    if isempty(dio[CompositingProgram])
        dio[Entity(1)] = CompositingProgram(Program(compositing_shaders()))
    end
    c = singleton(dio, Canvas)
    wh = size(c)
    while length(dio[PeelTarget]) < 2
        Entity(dio.ledger, DioEntity(),
               PeelTarget(GLA.FrameBuffer(wh,
                                          GLA.Texture(RGBAf0, wh; internalformat = GL_RGBA),
                                          GLA.Texture(RGBAf0, wh; internalformat = GL_RGBA),
                                          GLA.Texture(GLA.Depth{Float32}, wh)),
                          c.background))
    end
    if isempty(dio[ColorBlendTarget])
        dio[Entity(1)] = ColorBlendTarget(GLA.FrameBuffer(wh,
                                                          map(x -> GLA.Texture(x, wh),
                                                              (RGBAf0, RGBAf0,
                                                               GLA.Depth{Float32}))...),
                                          c.background)
    end
end

function Overseer.update(renderer::DepthPeelingRenderer, m::AbstractLedger)
    glDisableCullFace()
    vao = m[PeelingVao]
    ivao = m[InstancedPeelingVao]
    if isempty(vao) && isempty(ivao)
        return
    end
    rem1(x, y) = (x - 1) % y + 1

    spatial  = m[Spatial]
    material = m[Material]
    modelmat = m[ModelMat]
    ucolor   = m[UniformColor]
    bcolor   = m[BufferColor]
    light    = m[PointLight]
    camera   = m[Camera3D]
    alpha    = m[Alpha]

    peeling_program = m[PeelingProgram][1]
    ipeeling_program = m[InstancedPeelingProgram][1]

    peel_comp_program   = m[PeelingCompositingProgram][1]
    blending_program    = m[BlendProgram][1]
    compositing_program = m[CompositingProgram][1]

    colorblender    = m[ColorBlendTarget][1]
    peeling_targets = m[PeelTarget].data[1:2]
    iofbo           = m[IOTarget][1]
    fullscreenvao   = m[FullscreenVao][1]
    vis             = m[Visible]

    set_light_camera_uniforms = (prog) -> begin
        for e in @entities_in(light && ucolor && spatial)
            gluniform(prog, light[e], ucolor[e], spatial[e])
        end
        for e in @entities_in(spatial && camera)
            gluniform(prog, spatial[e], camera[e])
        end
    end

    set_model_material = (e_modelmat, e_material) -> begin
        gluniform(peeling_program, :material, Vec2(e_material.specpow, e_material.specint))
        gluniform(peeling_program, :modelmat, e_modelmat.modelmat)
    end

    bind(colorblender)
    draw(colorblender)
    clear!(colorblender)
    #     #TODO change this nonsense
    canvas_width, canvas_height = Float32.(size(iofbo))

    resize!(colorblender, (Int(canvas_width), Int(canvas_height)))
    resize!(peeling_targets[1], (Int(canvas_width), Int(canvas_height)))
    resize!(peeling_targets[2], (Int(canvas_width), Int(canvas_height)))

    glEnable(GL_DEPTH_TEST)
    glDepthFunc(GL_LEQUAL)
    glDisable(GL_BLEND)

    # # first pass: Render the previous opaque stuff first
    bind(peel_comp_program)
    gluniform(peel_comp_program, :first_pass, true)
    gluniform(peel_comp_program, :color_texture, 0, color_attachment(iofbo, 1))
    gluniform(peel_comp_program, :depth_texture, 1, depth_attachment(iofbo))
    gluniform(peel_comp_program, :color_id_texture, 3, color_attachment(iofbo, 2))
    bind(fullscreenvao)
    draw(fullscreenvao)
    gluniform(peel_comp_program, :first_pass, false)

    it1 = @entities_in(vao && modelmat && material && bcolor)

    renderall_separate = () -> begin
        set_light_camera_uniforms(peeling_program)
        for e in it1
            evao = vao[e]
            if vis[e].visible
                set_model_material(modelmat[e], material[e])
                gluniform(peeling_program, :alpha, alpha[e].α)
                GLA.bind(evao)
                GLA.draw(evao)
            end
        end
    end

    renderall_instanced = () -> begin
        set_light_camera_uniforms(ipeeling_program)
        for evao in ivao
            GLA.bind(evao)
            GLA.draw(evao)
        end
    end

    function render_start(prog, renderfunc)
        bind(prog)
        set_light_camera_uniforms(prog)
        gluniform(prog, :first_pass, true)
        gluniform(prog, :canvas_width, canvas_width)
        gluniform(prog, :canvas_height, canvas_height)
        renderfunc()
        return gluniform(prog, :first_pass, false)
    end

    # first pass: Render all the transparent stuff
    # separate
    render_start(peeling_program, renderall_separate)
    render_start(ipeeling_program, renderall_instanced)

    #start peeling passes
    for layer in 1:renderer.num_passes
        currid  = rem1(layer, 2)
        currfbo = peeling_targets[currid]
        previd  = 3 - currid
        prevfbo = layer == 1 ? colorblender : peeling_targets[previd]
        bind(currfbo)
        draw(currfbo)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
        glEnable(GL_DEPTH_TEST)
        glDisable(GL_BLEND)

        #         # peel: Render all opaque stuff
        bind(peel_comp_program)
        gluniform(peel_comp_program, :color_texture, 0, color_attachment(iofbo, 1))
        gluniform(peel_comp_program, :depth_texture, 1, depth_attachment(iofbo))
        gluniform(peel_comp_program, :prev_depth, 2, depth_attachment(prevfbo))
        bind(fullscreenvao)
        draw(fullscreenvao)

        bind(peeling_program)
        gluniform(peeling_program, :depth_texture, 0, depth_attachment(prevfbo))
        renderall_separate()

        bind(ipeeling_program)
        gluniform(ipeeling_program, :depth_texture, 0, depth_attachment(prevfbo))
        renderall_instanced()

        bind(colorblender)
        draw(colorblender, 1)

        glDisable(GL_DEPTH_TEST)
        glEnable(GL_BLEND)
        glBlendEquation(GL_FUNC_ADD)
        glBlendFuncSeparate(GL_DST_ALPHA, GL_ONE, GL_ZERO, GL_ONE_MINUS_SRC_ALPHA)

        bind(blending_program)
        gluniform(blending_program, :color_texture, 0, color_attachment(currfbo, 1))

        bind(fullscreenvao)
        draw(fullscreenvao)
    end
    bind(iofbo)
    draw(iofbo)
    glDisable(GL_BLEND)

    bind(compositing_program)
    gluniform(compositing_program, :color_texture, 0, color_attachment(colorblender, 1))
    gluniform(compositing_program, :color_id_texture, 1, color_attachment(colorblender, 2))
    bind(fullscreenvao)
    return draw(fullscreenvao)
end


text:
include("2.jl")
include("5.jl")
function obolochkaper(r::Robot)
    angleofsquare(r)
    kub(r)
end

function angleofsquare(r::Robot)
    while ((!isborder(r,Sud)) || (!isborder(r,Ost)))
        if (!isborder(r,Sud))
            move!(r,Sud)
        end
        if (!isborder(r,Ost))
            move!(r,Ost)
        end
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
module DynamicProgramming

# Dependencies
using LinearAlgebra
using SparseArrays
using MarkovChains
using ProgressMeter

import Base:
	show

# Include statements
include("hjb.jl")

# Export types
export
	HJBEquation

# Export methods
export
	solve

"""
A Julia package for solving Dynamic Programming problems.

"""

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Emoslib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Emoslib")
JLLWrappers.@generate_main_file("Emoslib", UUID("497af919-c005-5197-85cc-b1f3ea720854"))
end  # module Emoslib_jll


text:
using MPIReco

@testset "system matrix recovery" begin
  # patterns of system matrix
  Random.seed!(1234)
  bSF = MPIFile(joinpath(datadir, "calibrations", "12.mdf"))
  shape = calibSize(bSF)
  numFreqs = 3
  f = filterFrequencies(bSF, sortBySNR=true)
  f = f[shuffle(1:100)[1:numFreqs]]
  S = getSystemMatrix(bSF,bgCorrection=true)
  S = S[:,f]

  # undersampled measurement
  numSamples = Int64(div(prod(shape),2))
  samplingIdx = shuffle(1:prod(shape))[1:numSamples]
  y = S[samplingIdx,:]

  # solver parameters
  params = Dict{Symbol,Any}()
  params[:shape] = Tuple(shape)
  params[:ρ_l1] = 2.0
  params[:λ_l1] = 2.e-2
  params[:iterationsInner] = 50
  params[:iterations] = 10
  params[:relTol] = 1.e-2
  params[:absTol] = 1.e-4

  # DCT based recovery
  @info "recovery using DCT"
  S_dct = smRecovery(y,samplingIdx,params)
  for k=1:numFreqs
    @test norm(S[:,k] - S_dct[:,k]) / norm(S[:,k]) ≈ 0 atol=0.15
  end

  # DCT and low rank
  @info "recovery using DCT and LR"
  params[:λ_lr] = 2.e-2
  params[:ρ_lr] = 2.e-1
  S_dctlr = smRecovery(y,samplingIdx,params)
  for k=1:numFreqs
    @test norm(S[:,k] - S_dct[:,k]) / norm(S[:,k]) ≈ 0 atol=0.15
  end
end


text:
using Tricks
using InteractiveUtils
using Test

struct Iterable end; struct NonIterable end;
iterableness_static(::Type{T}) where T = static_hasmethod(iterate, Tuple{T}) ? Iterable() : NonIterable()

struct Foo end

@testset "static_hasmethod" begin
    @testset "positive: $(typeof(data))" for data in (
        "abc", [1,2,3], (2,3), ones(4,10,2), 'a',  1:100
    )
        T = typeof(data)
        @test iterableness_static(T) === Iterable()
        code_typed_ir = (@code_typed iterableness_static(T))[1].code
        @test code_typed_ir == [:(return $(QuoteNode(Iterable())))]
    end

    @testset "negative: $(typeof(data))" for data in (
        :a, rand, Int
    )
        T = typeof(data)
        @test iterableness_static(T) === NonIterable()
        code_typed = (@code_typed iterableness_static(T))
        @test code_typed[2] === NonIterable  # return type
        @test length(code_typed[1].code) <= 3  # currently has dead conditions on statements 1 and 2, but LLVM removed them.
    end

    @testset "add method" begin
        @test iterableness_static(Foo) === NonIterable()

        Base.iterate(::Foo) = ("Foo", nothing);
        Base.iterate(::Foo, ::Nothing) = nothing;
        Base.length(::Foo) = 1;
        @test collect(Foo()) == ["Foo"]

        @test iterableness_static(Foo) === Iterable()
    end

    @testset "delete method" begin
        @test iterableness_static(Foo) === Iterable()
        meth = first(methods(iterate, Tuple{Foo}))
        Base.delete_method(meth)
        @test_throws MethodError collect(Foo())

        @test iterableness_static(Foo) === NonIterable()
    end
end


text:
"""
    VanillaNetStats()

Network statistics for `VanillaConfig`.
"""
mutable struct VanillaNetStats <: Statistics
    # network traffic
    uplink::Vector{Int}
    downlink::Vector{Int}

    delta_earned::Int

    VanillaNetStats() = new(Vector{Int}(undef, 0), Vector{Int}(undef, 0), 0)
end

"""
    update_stats!(stats, round_num)

Updates the network traffic stats of a `VanillaConfig`.
"""
function update_stats!(stats::VanillaNetStats, round_num::Int)
    num_comm = round_num * STATS.base.num_clients_per_round
    updown = num_comm * STATS.base.num_weights * sizeof(STATS.T)
    push!(stats.uplink, updown)
    push!(stats.downlink, updown - stats.delta_earned)
end


text:
module GraphSignal

using Distances                 # This is necessary in Ggrid
using Distributions             # This is necessary in AddNoise
using SparseArrays
using Statistics
using LinearAlgebra

export GraphSig, replace_data!, snr, gpath, ggrid, Adj2InvEuc, AddNoise

# GraphSig data structure and constructors
"""
    G = GraphSig(W, xy, f, name, plotspecs, length, dim, size)

is a data structure for a GraphSig object containing the following fields:
* `W::SparseMatrixCSC{Float64,Int}`: the edge weight matrix in a sparse format
* `xy::Matrix{Float64}`: the coordinates of the vertices (could be >= 3D)
* `f::Matrix{Float64}`: the values of the function/data at the vertices; each column is an input signal vector, i.e., number of rows == number of nodes.
* `name::String`: a title/name for the data
* `plotspecs::Vector{Symbol}`: specifications for plotting:
    (most likely may be dropped later thanks to julia/Plots.jl)
                   - symm: use a symmetric colorscale
                   - gray: use grayscale
                   - gray255: plot a grayscale image with color bounds [0,255]
                   - copper: use a copper color scale
                   - notitle: don't display the title
                   - nocolorbar: don't display a colorbar
                   - stem: use a stem plot
                   - CLim[cmin,cmax]: set the dynamic display range to [cmin,cmax]
                   - size25: set the size of the nodes to 25 (or whatever value is specified)
                   - LineWidth2: set the width of the lines to 2 (or whatever value is specified)
                   - LineColorb: set the color of the lines/graph edges to color 'b' (or whatever color is specified)
                   - red/black/blue: make the nodes red, black, or blue
                   - verbatim{{...}}: execute the command in the brackets
"""
mutable struct GraphSig
    W::SparseMatrixCSC{Float64,Int} # edge weight matrix
    xy::Matrix{Float64}               # spatial coordinates (could be >= 3D)
    f::Matrix{Float64}                # the data sequences (signals) on the nodes
    name::String                      # name of the GraphSig object
    plotspecs::String        # special instructions for plotting
    # The following three are deriable from W, xy, and f:
    length::Int                     # length of the data sequence
    dim::Int                        # the dimension of the node coordinates
    size::Tuple{Int,Int}          # size of all the data specified

    # An inner constructor here.
    function GraphSig(W::SparseMatrixCSC{Float64,Int};
                      xy::Matrix{Float64} = Matrix{Float64}(undef, 0, 0),
                      f::Matrix{Float64} = zeros(Base.size(W, 1), 1),
                      name::String = "",
                      plotspecs :: String = "",
                      #plotspecs::Vector{Symbol} = Vector{Symbol}(0),
                      length::Int = Base.size(W, 1),
                      dim::Int = Base.size(xy, 2),
                      size::Tuple{Int,Int} = (Base.size(W, 1), Base.size(W, 2) + Base.size(xy, 2) + Base.size(f, 2)))
#
        if isempty(W)
            error("Cannot setup a graph without an edge weight matrix W!")
        elseif Base.size(W, 1) != Base.size(f, 1)
            error("Signal length != # Nodes!!")
        elseif Base.size(W, 1) != Base.size(xy, 1) && Base.size(xy, 1) != 0
            @warn("# Node coordinates != # Nodes. Hence, use null coordinates.")
            new( W, Matrix{Float64}(0, 0), f, name, plotspecs, Base.size(W, 1),
                 0, (Base.size(W, 1), Base.size(W, 2) + Base.size(f, 2)) )
        else # Finally, the full constructor!
            new( W, xy, f, name, plotspecs, Base.size(W, 1), Base.size(xy, 2),
                 (Base.size(W, 1), Base.size(W, 2) + Base.size(xy, 2) + Base.size(f, 2)) )
        end
    end # of the inner constructor
end # of type GraphSig

#
# Other GraphSig related functions
#
"""
    replace_data!(G::GraphSig, newdata::Matrix{Float64})

Replace the data in a GraphSig object

### Input Arguments
* `G::GraphSig`: the input GraphSig object whose data is to be replaced; various fields are updated after execution of this function
* `newdata::Matrix{Float64}`: the new data (must be a matrix)
"""
function replace_data!(G::GraphSig, newdata::Matrix{Float64})

    (rows, cols) = Base.size(newdata)
    if G.length == rows         # update the relevant fields accordingly
        G.f = newdata
        G.size = (rows, Base.size(G.W, 2) + Base.size(G.xy, 2) + cols)
    else
        error("The length of each signal (", rows, ") does not match the number of nodes (", G.length, ")!")
    # elseif G.length == cols # This is a dangerous assumption; so prohibit it.
    #    G.f = newdata'
    # elseif isscalar(newdata) # This does not work anymore unlike MATLAB
    #    G.f = newdata*ones(G.length, 1)
    end
end # of function replace_data!

"""
    (value, sigma) = snr(G1::GraphSig, G2::GraphSig, stdp::Bool = false)

    Compute the SNR between G1 (original signal) and G2 (noisy signal)

### Input Arguments
* `G1::GraphSig`: Original reference graph signal
* `G2::GraphSig`: Restored or noisy graph signal
* `stdp::Bool`: A flag to specify to compute the standard deviation of the difference

### Output Arguments
* `value`: Signal/Noise Ratio in the dB unit
* `sigma`: the standard deviation of the noise (if stdp == true)
"""
function snr(G1::GraphSig, G2::GraphSig, stdp=false)
  tmp = norm(G2.f - G1.f)
  if tmp < 10 * eps()
    @warn("Two signals are the same, hence SNR = Inf")
    value = Inf
  else
    value = 20 * log10(norm(G1.f) / tmp)
  end

  if stdp
    return value, std(G2.f - G1.f)
  else
    return value
  end

end # of snr

"""
    G = gpath(N, f, name)

Generate a GraphSig object for a 1-D path of length `N`

### Input Arguments
* `N::Int`: the length of the path
* `f::Matrix{Float64}`: the signal to be used
* `name::String`: the name for the GraphSig object

### Output Argument
* `G::GraphSig`: the GraphSig ojbect representing the 1-D path of length `N`
"""
function gpath(N::Int,
               f::Matrix{Float64} = ones(N, 1),
               name::String = "Path of length $(N)")

    # the x coordinates
    xy = Matrix{Float64}(reshape(1:N, N, 1))
    # another option: xy = (1.0:Float64(N))[:,:]

    # the weight matrix
    ev = ones(N - 1)               # subdiagonal entries
    W  = spdiagm(-1 => ev, 1 => ev)
    #W = spdiagm((ev, ev), (-1, 1)) # size(W) is (N, N)
    # MATLAB:
    # ev = ones(N,1);
    # W = spdiags([ev 0*ev ev], [-1 0 1], N, N);

    # create and return the GraphSig object
    return GraphSig(W, xy = xy, f = f, name = name)

end # of function gpath

"""
    function ggrid(Nx::Int, Ny::Int, connect::Symbol)

Generate a GraphSig object for a 2-D grid that is `Nx` by `Ny`.

### Input Arguments
* `Nx::Int`: the number of points in the x direction
* `Ny::Int`: the number of points in the y direction
* `connect::Symbol`: specifies the connectivity mode of the 2-D grid:
    choices are: :c4 (default), :c8, :full

### Output Argument
* `G::GraphSig`: the GraphSig ojbect representing the 2-D `Nx` by `Ny` grid.
"""
function ggrid(Nx::Int, Ny::Int, connect::Symbol = :c4)

    # the total number of nodes
    N = Nx * Ny

    # the xy coordinates
    # MATLAB:
    # xy = [(1:N)', repmat((1:Ny)',Nx,1)];
    # xy(:,1) = ceil(xy(:,1)/Ny);
    xy =hcat( ceil.(collect(1:N) / Ny), repeat(1:Ny, Nx) )

    if connect == :c4         # this is the default, i.e., l; r; u; d.
        # make 1-D weight matrices
        ex = ones(Nx - 1)
        Wx  = spdiagm(-1 => ex, 1 => ex)
        #Wx = spdiagm((ex, ex), (-1, 1))
        ey = ones(Ny - 1)
        Wy  = spdiagm(-1 => ey, 1 => ey)
        #Wy = spdiagm((ey, ey), (-1, 1))
        # form the 2-D weight matrix
        W = kron(sparse(1.0I, Nx, Nx), Wy) + kron(Wx, sparse(1.0I, Ny, Ny))
            # this keeps a sparse form
    else                # now for 8-connected or fully-connected cases
        W = sparse(pairwise(Euclidean(), xy')) # This is from Distances package.
        W[1:(N + 1):(N * N)] = 10 # set each diagonal entry to 10 to prevent blowup.
        W = W.^-1         # convert dist to affinity = inv. dist.
        if connect == :c8 # 8-connected case incl. diagonal connections
            W[find(W .< 0.7)] = 0
        elseif connect == :full # a complete graph case
            W[1:(N + 1):(N * N)] = 0 # set the diagonal to zero (i.e., no loops)
        else
            error("Cannot understand this connectivity: ", connect)
        end
    end

    return GraphSig(W, xy = xy, f = Matrix{Float64}((1:N)[:, :]), name = "$(Nx) by $(Ny) grid with connectivity mode $(connect)")

end # of function ggrid

"""
    function Adj2InvEuc(G::GraphSig)

Given a GraphSig object 'G' with a binary weight (adjacency) matrix, generate a GraphSig object 'GInvEuc' with an inverse Euclidean distance matrix

### Input Arguments
* `G::GraphSig`: a GraphSig object

### Output Argument
* `G::GraphSig`: a GraphSig object with inverse Euclidean weights
"""
function Adj2InvEuc(G::GraphSig)
    W = deepcopy(G.W)
    xy = deepcopy(G.xy)
    f = deepcopy(G.f)

    #find the edges
    (rows,cols) = findnz(W)

    #remove duplicate points
    for j = length(rows):-1:1
      if j<=length(rows) && norm(xy[rows[j],:]-xy[cols[j],:],2) < 10^3*eps()
        W = W[setdiff(1:end, rows[j]), :]
        W = W[:, setdiff(1:end, rows[j])]
        xy = xy[setdiff(1:end, rows[j]), :]
        f = f[setdiff(1:end, rows[j]),:]
        (rows,cols) = findnz(W)
      end
    end

    for j = 1:length(rows)
      W[rows[j],cols[j]] = 1/norm(xy[rows[j],:]-xy[cols[j],:],2)
    end

    return GraphSig(W, xy = xy, f = f, name = string(G.name, " (inverse Euclidean weight matrix)"), plotspecs = G.plotspecs)
end

"""
    function AddNoise(G::GraphSig; SNR::Float64=Float64(1.24), noisetype::String = "gaussian")

Add noise to the data of a GraphSig object

### Input Arguments
* `G::GraphSig`: a GraphSig object
* `SNR`: the SNR that the noisy signal should have
* `noisetype`: the type of noise: Gaussian (default) or Poisson

### Output Argument
* `G::GraphSig`: the GraphSig object with added noise
* `sigma`: the standard deviation of the noise
"""
function AddNoise(G::GraphSig; SNR::Float64=Float64(1.24), noisetype::String = "gaussian")
  f = deepcopy(G.f)

  if noisetype == "gaussian"
    # generate Gaussian noise
    noise = randn(size(f))

    # scale the noise to the desired SNR level
    sigma = norm(f,2)/10.0^(0.05*SNR)/norm(noise,2)
    noise = sigma*noise

    # generate the noisy signal
    f = f + noise

  elseif noisetype == "poisson"
    # generate Poisson noise
    (rows,cols) = size(f)
    noise = zeros(rows,cols)
    for i = 1:rows
      for j = 1:cols
        noise[i,j] = rand(Poisson(f[i,j])) - f[i,j]
      end
    end

    # scale the noise to the desired SNR level
    sigma = norm(f,2)/10.0^(0.05*SNR)/norm(noise,2)
    noise = sigma*noise

    # generate the noisy signal
    f = f + noise

  end

  return GraphSig(G.W, xy = G.xy, f = f, name = string("SNR = ",SNR," ", G.name), plotspecs = G.plotspecs)
end

end # of module GraphSignal


text:
include("wk.jl")
using DiffEqFlux
using Flux
using DifferentialEquations
using Random
using CSV
using DataFrames
using Optim
using FFTW
using IterTools: ncycle

#=
Learn waveforms from synthetic data, where we know an exact solution
exists
=#

#=
load real data
=#
function load_art(sample = true,
  loc = "/Users/daphne/Documents/mit/ThingsWithWK/samples.csv",
  start = 4170,
  delta = 280,
  fs = 360)

  out = CSV.read("/Users/daphne/Documents/mit/ThingsWithWK/samples.csv", DataFrame)
  col = out[!,"'ART'"][2:end]
  col_f = parse.(Float64, col)

  t = range(0.,delta/fs,length=delta)
  pressure_data =col_f[start:start+delta]

  if sample
    # list of ids at which to take pressure samples
    draw_list = randperm(length(pressure_data)-1)[1:n_points]
    draw_list = sort(draw_list)

    t = t[draw_list]
    pressure_data = pressure_data[draw_list]
  end

  t, pressure_data

end

#=
SETUP
pretty much all the params you will need to change are here
tspan = range over which to generate data
  perturbed from (0,1) so that timesteps don't fall exactly on
  parts of Q(t) which are not differentiable
u0 = initial condition
R = SVR
C = aortic compliance
R_a = aortic impedence
L = aortic inductance (blood inertia)

NOTE do not remove parameters from the list - models implemented
in wk.jl are designed to use the same parameter list.

=#
R = 20
C = 1.2
L = 0.002

R_a = 0.77

p = [R_a/60,C, R/60, L]

n_points = 30

t_p, pressure_data = load_art()

# get diastole (decreasing part of curve)
wh = t_p .> t_p[end]*2/5
t_p_dia = t_p[wh]
pressure_dia = pressure_data[wh]

wh = t_p .<= t_p[end]*2/5
t_p_sys = t_p[wh]
pressure_sys = pressure_data[wh]

u0 = pressure_sys[1]
tspan = (minimum(t_p_sys),maximum(t_p_sys))

t_synth = range(tspan[1],tspan[2],length=length(pressure_sys))
which_model = wk2

#=================================#

# solve ODE to get synthetic data
prob = ODEProblem(which_model,u0,tspan,p)
ode_data = Array(solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat=t_synth))

ode_data_dct = abs.(FFTW.rfft(ode_data,1)) #, FFTW.REDFT10))

# simple NN
dudt = FastChain((x, p) -> x,
                  FastDense(1, 32, elu),
                  FastDense(32, 16, elu),
                  FastDense(16, 8, elu),
                  FastDense(8, 1))

# initialize ODE to train for the synthetic data
n_ode = NeuralODE(dudt, tspan, saveat=t_synth)

n_ode_real = NeuralODE(dudt, tspan, Tsit5(), saveat=t_p_sys, reltol=1e-9, abstol=1e-9)

function predict_neuralode(p)
  reshape(Array(n_ode([u0], p)),:)
end

function predict_neuralode_real(p)
  reshape(Array(n_ode_real([u0], p)),:)
end

function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, (ode_data[1:length(pred)] .- pred[:]))
    return loss, pred
end

function loss_neuralode_real(p)
    pred = predict_neuralode_real(p)
    loss = sum(abs2, (pressure_sys[1:length(pred)] .- pred[:]))
    return loss, pred
end

function dct_loss_neuralode(p)
    pred = predict_neuralode(p)
    pred_dct = abs.(FFTW.rfft(pred[:],1))
    loss = sum(abs2, (ode_data_dct - pred_dct))/length(pred)
end

loss = []

iter = 0
callback = function (p, l, pred; doplot = true)
  global iter
  iter += 1

  display(l)
  push!(loss,l)
  if doplot
    # plot current prediction against data
    plt = scatter(t_synth, ode_data, label = "data")
    scatter!(plt, t_synth, pred[:], label = "prediction")
    xlabel!(plt," time (s)")
    ylabel!(plt, "pressure (mmHg)")
    display(plot(plt))
  end

  return false
end

iter = 0
callback_real = function (p, l, pred; doplot = true)
  global iter
  iter += 1

  display(l)
  push!(loss,l)
  if doplot
    # plot current prediction against data
    plt = scatter(t_p_sys, pressure_sys, label = "data")
    scatter!(plt, t_p_sys, pred[:], label = "prediction")
    xlabel!(plt," time (s)")
    ylabel!(plt, "pressure (mmHg)")
    display(plot(plt))
  end

  return false
end

k = 10
train_loader = Flux.Data.DataLoader((ode_data, t_synth); batchsize = k)
numEpochs = 1000

result_neuralode = DiffEqFlux.sciml_train(loss_neuralode,
                                          n_ode.p,
                                          #ADAM(0.05),
                                          Optim.KrylovTrustRegion(),
                                          #_data = ncycle(train_loader, numEpochs),
                                          cb = callback,
                                          maxiters = numEpochs,
                                          allow_f_increases = true)

final_neuralode = DiffEqFlux.sciml_train(loss_neuralode_real,
                                          result_neuralode.minimizer,
                                          #ADAM(0.05),
                                          Optim.KrylovTrustRegion(),
                                          #_data = ncycle(train_loader, numEpochs),
                                          cb = callback_real,
                                          maxiters = 200,
                                          allow_f_increases = true)

max_stars_repo_name: daphneschles/ThingsWithWK

max_issues_repo_name: daphneschles/ThingsWithWK

max_forks_repo_name: daphneschles/ThingsWithWK


text:
function L(member::Dict{Integer,Array{Integer}}, node::Dict{Integer,Array{Float64}})
	l = Dict{Integer, Float64}()
	for e = 1:length(member)
		l[e] = hypot(node[member[e][2]][1]-node[member[e][1]][1], node[member[e][2]][2]-node[member[e][1]][2])
	end
	return l
end

function LocalK(member::Dict{Integer,Array{Integer}}, node::Dict{Integer,Array{Float64}}, E::Dict{Integer, Float64}, ν::Dict{Integer, Float64}, A::Dict{Integer, Float64}, I:: Dict{Integer, Array{Float64}})
	K = Dict{Integer, Array{Float64}}()
	l = L(member,node)
	for e = 1:length(member)
		K[e] = zeros(Float64,12,12)
		K[e][1,1] = E[e]*A[e]/l[e]
		K[e][1,7] = -E[e]*A[e]/l[e]
		K[e][2,2] = E[e]*12*I[e][3]/l[e]^3
		K[e][2,6] = E[e]*6*I[e][3]/l[e]^2
		K[e][2,8] = -E[e]*12*I[e][3]/l[e]^3
		K[e][2,12] = E[e]*6*I[e][3]/l[e]^2
		K[e][3,3] = E[e]*12*I[e][2]/l[e]^3
		K[e][3,5] = -E[e]*6*I[e][2]/l[e]^2
		K[e][3,9] = -E[e]*12*I[e][2]/l[e]^3
		K[e][3,11] = -E[e]*6*I[e][2]/l[e]^2
		K[e][4,4] = E[e]*I[e][1]/2/(1+ν[e])/l[e]
		K[e][4,10] = -E[e]*I[e][1]/2/(1+ν[e])/l[e]
		K[e][5,3] = -E[e]*6*I[e][2]/l[e]^2
		K[e][5,5] = E[e]*4*I[e][2]/l[e]
		K[e][5,9] = E[e]*6*I[e][2]/l[e]^2
		K[e][5,11] = E[e]*2*I[e][2]/l[e]
		K[e][6,2] = E[e]*6*I[e][3]/l[e]^2
		K[e][6,6] = E[e]*4*I[e][3]/l[e]
		K[e][6,8] = -E[e]*6*I[e][3]/l[e]^2
		K[e][6,12] = E[e]*2*I[e][3]/l[e]
		K[e][7,1] = -E[e]*A[e]/l[e]
		K[e][7,7] = E[e]*A[e]/l[e]
		K[e][8,2] = -E[e]*12*I[e][3]/l[e]^3
		K[e][8,6] = -E[e]*6*I[e][3]/l[e]^2
		K[e][8,8] = E[e]*12*I[e][3]/l[e]^3
		K[e][8,12] = -E[e]*6*I[e][3]/l[e]^2
		K[e][9,3] = -E[e]*12*I[e][2]/l[e]^3
		K[e][9,5] = E[e]*6*I[e][2]/l[e]^2
		K[e][9,9] = E[e]*12*I[e][2]/l[e]^3
		K[e][9,11] = E[e]*6*I[e][2]/l[e]^2
		K[e][10,4] = -E[e]*I[e][1]/2/(1+ν[e])/l[e]
		K[e][10,10] = E[e]*I[e][1]/2/(1+ν[e])/l[e]
		K[e][11,3] = -E[e]*6*I[e][2]/l[e]^2
		K[e][11,5] = E[e]*2*I[e][2]/l[e]
		K[e][11,9] = E[e]*6*I[e][2]/l[e]^2
		K[e][11,11] = E[e]*4*I[e][2]/l[e]
		K[e][12,2] = E[e]*6*I[e][3]/l[e]^2
		K[e][12,6] = E[e]*2*I[e][3]/l[e]
		K[e][12,8] = -E[e]*6*I[e][3]/l[e]^2
		K[e][12,12] = E[e]*4*I[e][3]/l[e]
	end
	return K
end


text:
# using StaticArrays

abstract type Pulse end

abstract type BasicPulse <: Pulse end
abstract type PulseSequence end

"""
IdlePulse is a pulse where the amplitude and phase are zero
"""
struct IdlePulse{T} <: BasicPulse
    Ω::T
    ϕ::T
    T::T
    dt::T
    IdlePulse(T, dt) = new{Float64}(0.0, 0.0, T, dt)
end

"""
A πPulse is a pulse where the duration is chosen to produce a pi rotation
"""
struct PiPulse{T} <: BasicPulse
    Ω::T
    ϕ::T
    T::T
    dt::T
    PiPulse(Ω, ϕ, dt) = new{Float64}(Ω, ϕ, π/Ω, dt)
end

"""
A πHalfPulse is a pulse where the duration is chosen to produce a pi/2 rotation
"""
struct PiHalfPulse{T} <: BasicPulse
    Ω::T
    ϕ::T
    T::T
    dt::T
    PiHalfPulse(Ω, ϕ, dt) = new{Float64}(Ω, ϕ, π/Ω/2, dt)
end

"""
An ArbitraryPulse is a pulse where you control everything
"""
struct ArbitraryPulse{T} <: BasicPulse
    Ω::T
    ϕ::T
    T::T
    dt::T
    ArbitraryPulse(Ω, ϕ, T, dt) = new{Float64}(Ω, ϕ, T, dt)
end

"""
Piecewise constant pulse constructed from arrays of amplitude and phase
"""
struct PWPulse{T} <: Pulse
    Ω::Array{T,1}
    ϕ::Array{T,1}
    T::T
    TArray::Array{T,1}
    N::Integer
    dt::T
    # TODO we need to think about the internal time array here, step = dt instead of length?
    PWPulse(Ω, ϕ, T, N, dt) = new{Float64}(Ω, ϕ, T, collect(range(0, T, length = N)), N, dt)
end

function _pulse_t(p::BasicPulse, t)::Tuple{Float64,Float64}
    (p.Ω, p.ϕ)
    # p.Ω * exp(1im * p.ϕ)
end

function _pulse_t(p::PWPulse, t)::Tuple{Float64,Float64}
    idx = searchsortedlast(p.TArray, t)

    if idx == 0
        idx = 1
    end
    @inbounds (p.Ω[idx], p.ϕ[idx])
    # @inbounds p.Ω[idx] * exp(1im * p.ϕ[idx])
end

function get_length(p::BasicPulse)::Int64
    Int(floor(p.T/p.dt))
end

function get_length(p::PWPulse)::Int64
    p.N
end

function get_dt(p::Pulse)::Array{Float64,1}
    N = get_length(p)
    ones(N) .* p.dt
end

(c::BasicPulse)(t::Float64) = _pulse_t(c, t)
(c::PWPulse)(t::Float64) = _pulse_t(c, t)

"""
Sequence holds a list of pulses and other useful information to define a pulse sequence!
It's also possible to call (sequence(t)) to get back the value of the drive at that time
which makes visualising and simulation sequences much easier!
"""
struct Sequence
    pulses
    T_Array
    cum_T_Array
    duration
    function Sequence(pulseList)
        N = length(pulseList)
        TArray = zeros(N)
        for i=1:N
            TArray[i] = pulseList[i].T
        end
        cs = cumsum(TArray)
        prepend!(cs, 0.0)
        new(pulseList, TArray, cs ,sum(TArray))
    end
end

function _seq_t(c, t)::Tuple{Float64,Float64}
    ind = searchsortedlast(c.cum_T_Array, t)
    N = length(c.cum_T_Array)
    if ind >= N
        return (0.0, 0.0)
    elseif ind >= 2
        t = t - c.T_Array[ind - 1]
    end
    c.pulses[ind](t)
end

(c::Sequence)(t::Float64) = _seq_t(c, t)


text:
module Dets

export det

using ..Common.Generator

code_det_kernel(A::Matrix) = code_det_kernel(A, "subdet")
function code_det_kernel(A::Matrix, sdname::String)
    n, m = size(A)
    @assert m <= n
    @assert m >= 1

    code = {}
    N = 1<<n
#    subdets = {gensym("subdet$k") for k=1:(N-1)}
#    subdets = {symbol("$sdname$(bin(k))") for k=1:(N-1)}
    subdets = {symbol(string(sdname,
                             [(k&(1<<(j-1))!=0 ? j : "") for j=1:n ]...))
               for k=1:(N-1)}

    for i=1:n;   subdets[1<<(i-1)] = A[i,1];  end
    for mask = 1:(N-1)
        j = count_ones(mask)
        if j < 2 || j > m; continue; end

        terms = {}
        negate = false
        for i in n:-1:1
            bit = 1<<(i-1)
            if (mask & bit) == 0; continue; end

            term = A[i,j]
            if j > 1
                subdet = subdets[mask & ~bit]
                term = :($subdet*$term)
            end
            if negate;  term = :(-$term);  end
            push!(terms, term)
            negate = !negate
        end
        push!(code, :($(subdets[mask]) = +($(terms...))))
    end

    code, subdets
end

function incomplete_det_components(n::Int, subdets::Vector)
    ds = {}
    for k=1:n
        d = subdets[(1<<n)-1-(1<<(k-1))]
        if isodd(n-k);  d = :(-$d); end
        push!(ds, d)
    end
    ds
end

incomplete_det() = [1]

# ---- det(rows...), incomplete_det(rows...) ----
for NC = 1:8
    generators = [symbol("v$k") for k in 1:NC]
#    args = [:($generator::Generator{$NC}) for generator in generators]
    args = [:($generator) for generator in generators]

    A  = [ :($generator[$i]) for i=1:NC, generator in generators ]
    code, ds = code_det_kernel(A, "d")
    d = ds[end]

    @eval function det($(args...))
        $(code...)
        $d
    end

    if NC > 1
        code, ds = code_det_kernel(A[:,1:(end-1)], "d")
        ds = incomplete_det_components(NC, ds)
        iargs = args[1:(end-1)]
        @eval function incomplete_det($(iargs...))
            $(code...)
            [$(ds...)]
        end
    end
end

end # module


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 = "FlexExtract"
version = v"7.1.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://www.flexpart.eu/gitmob/flex_extract", "e0005c99ac81d12faa45a8ff799debbd592b0dc0"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd flex_extract
atomic_patch -p1 /workspace/srcdir/patches/flexextract-makefile.patch
cd Source/Fortran/
make -f makefile_local_gfortran
cp calc_etadot_fast.out $bindir/calc_etadot
install_license ../../LICENSE.md
"""

# platforms = [
#     Platform("x86_64", "linux"; libc = "glibc")
# ]

platforms = supported_platforms()
# Remove the unsupported platforms
filter!(!Sys.iswindows, platforms)
filter!(!Sys.isapple, platforms)
filter!(!Sys.isfreebsd, platforms)
filter!(p -> libc(p) != "musl", platforms)
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("calc_etadot", :calc_etadot),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("eccodes_jll"),
    BuildDependency("Emoslib_jll"),
]

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

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using SparseArrays
using LinearAlgebra

function euler_fluxes(r, ru, rv, rE)
    gamma = 7.0/5
    u = @. ru/r
    v = @. rv/r
    p = @. (gamma - 1.0)*(rE-(ru^2/r+rv^2/r)/2.0)

    Frx = ru
    Fry = rv
    Frux= @. r*u^2+p
    Fruy= @. ru*v
    Frvx= @. ru*v
    Frvy= @. r*v^2+p
    FrEx= @. u*(rE+p)
    FrEy= @. v*(rE+p)
    return Frx, Fry, Frux, Fruy, Frvx, Frvy, FrEx, FrEy
end

function compact_div(Fx, Fy, h)
    global m
    m=ceil(Int,1/h)
    #x = h * (1:m)
    #m = 4  # testing purpose
    #println("m=",m)
    LHS = SymTridiagonal(ones(m), ones(m-1)/4);
    #println(LHS)
    RHS = (3 / 4h) * Tridiagonal(-ones(m-1), zeros(m), ones(m-1));
    #println(RHS)
    RHS = sparse(RHS)
    #println(RHS)
    dFx = LHS \ (RHS * Fx)
    dFy = LHS \ (RHS * Fy)
    divF = dFx+dFy
    return divF
end

function compact_filter(u, alpha)
    a = 5/8+3*alpha/4
    b = alpha+1/2
    c = alpha/4-1/8
    global m
    LHS = SymTridiagonal(ones(m), alpha*ones(m-1));
    RHS=sparse([1:m;2:m;1:m-1;3:m;1:m-2], [1:m;1:m-1;2:m;1:m-2;3:m],
           [a*ones(m); b/2.0*ones(2*(m-1)); c/2.0*ones(2*(m-2))])
    u_filt = LHS \ (RHS * u)
    return u_filt
end

function euler_rhs(r, ru, rv, rE, h)
    Frx, Fry, Frux, Fruy, Frvx, Frvy, FrEx, FrEy = euler_fluxes(r, ru, rv, rE)
    #println("Frx=",Frx)
    #println("Fry=",Fry)
    fr  =-compact_div(Frx,  Fry,  h)
    fru =-compact_div(Frux, Fruy, h)
    frv =-compact_div(Frvx, Frvy, h)
    frE =-compact_div(FrEx, FrEy, h)
    return fr, fru, frv, frE
end

function ff(u,h)
    #println("size of u in ff=",size(u))
    #println("u in ff=",u)
    fr, fru, frv, frE = euler_rhs(u[:,1], u[:,2], u[:,3], u[:,4], h)
    return [fr fru frv frE]
end

function euler_rk4step(r, ru, rv, rE, h, k, alpha)
    u = [r ru rv rE]
    #println("u=",u)
    #println("size of u=",size(u))
    k1 = k * ff(u, h)
    k2 = k * ff(u + 0.5k1, h)
    k3 = k * ff(u + 0.5k2, h)
    k4 = k * ff(u + k3, h)
    u = u + (k1 + 2k2 + 2k3 + k4) / 6
    #println(u)
    u = compact_filter(u,alpha)
    return u[:,1], u[:,2], u[:,3], u[:,4]
end

if true
  ru = [1;2;3;4]
  rv = [5;6;7;8]
  r  = [2;3;4;2]
  rE = [1;2;2;2]
  println(euler_rk4step(r, ru, rv, rE, 0.25, 0.2, 0.3))
end

#f(x, y) = x * sqrt(y)
function rk4(f::Function, x0::Float64, y0::Float64, x1::Float64, n)
    vx = Array{Float64}(undef, n+1)
    vy = Array{Float64}(undef, n+1)
    vx[1] = x = x0
    vy[1] = y = y0
    h = (x1 - x0) / n
    for i in 1:n
        k1 = h * f(x, y)
        k2 = h * f(x + 0.5h, y + 0.5k1)
        k3 = h * f(x + 0.5h, y + 0.5k2)
        k4 = h * f(x + h, y + k3)
        vx[i + 1] = x = x0 + i * h
        vy[i + 1] = y = y + (k1 + 2k2 + 2k3 + k4) / 6
    end
    return vx, vy
end


text:
module TestDerivative
    using Base.Test
    import FiniteDiff

    negative_sin(x) = -sin(x)
    negative_cos(x) = -cos(x)
    square(x) = x * x
    double(x) = x + x
    sqrt_deriv(x) = 1 / (2 * sqrt(x))

    function check_error(y_true, y_approximate, err = 1 // 10_000)
        # TODO: Decide how to handle very small inputs.
        if y_true > eps(y_true)
            @test abs(y_true - y_approximate) / abs(y_true) < err
        end
    end

    function run_tests(n::Integer)
        output = Array(Float64, 1)

        funcs = (
            (sin, cos),
            (cos, negative_sin),
            (negative_sin, negative_cos),
            (square, double),
            (exp, exp),
            (log, inv),
            (sqrt, sqrt_deriv),
        )

        modes = (
            FiniteDiff.ForwardMode(),
            FiniteDiff.BackwardMode(),
            FiniteDiff.CentralMode(),
            FiniteDiff.ComplexMode(),
        )

        for (f, f′) in funcs
            for itr in 1:n
                x = 1.0 + itr

                for i in 1:length(modes)
                    m = modes[i]

                    y = FiniteDiff.derivative(f, x, m)
                    check_error(f′(x), y)

                    FiniteDiff.derivative!(output, f, x, m)
                    check_error(f′(x), output[1])

                    tmp = FiniteDiff.derivative(f, m, mutates=false)
                    y = tmp(x)
                    check_error(f′(x), y)

                    tmp! = FiniteDiff.derivative(f, m, mutates=true)
                    tmp!(output, x)
                    check_error(f′(x), output[1])
                end

                y = FiniteDiff.derivative(f, x)
                check_error(f′(x), y)
            end
        end
    end

    @testset "derivative tests" begin
        run_tests(100)
    end
end