text:
# Pie charts
#function PlotPieChart(label_ids, values::AbstractArray{<:Real}, args...)
#    return PlotPieChart(label_ids, Float64.(values), args...)
#end

function PlotPieChart(values::AbstractArray{T}, count::Integer, x, y, radius;
                      normalize::Bool=false, label_fmt::String="%.1f", angle0=90.0,
                      label_ids::Vector{String}=["" for _ in 1:length(values)]) where {T<:ImPlotData}
    return PlotPieChart(label_ids, values, count, x, y, radius, normalize, label_fmt,
                        angle0)
end

function PlotPieChart(values::AbstractArray{T}, count::Integer, x, y, radius;
                      normalize::Bool=false, label_fmt::String="%.1f", angle0=90.0,
                      label_ids::Vector{String}=["" for _ in 1:length(values)]) where {T<:Real}
    return PlotPieChart(label_ids, Float64.(values), count, x, y, radius, normalize,
                        label_fmt, angle0)
end


text:
"""
    inout_map_property(𝑃::PROPERTY,
                       partition::AbstractVector{<:AbstractVector{Int}},
                       blocks::AbstractVector{Int},
                       n::Int
                      )::PROPERTY where {PROPERTY<:Property}

Map a property to the dimensions of analyzed blocks.

### Input

- `𝑃`         -- property
- `partition` -- block partition; elements are start and end indices of a block
- `blocks`    -- list of all output block indices in the partition
- `n`         -- total number of input dimensions

### Output

A new property of reduced dimension, if necessary.

### Notes

The return type depends on the ambient dimension (`n`), the block indices in the
partition and on the type of property:

- The property `𝑃` is returned unchanged whenever `n` matches the dimensions
  in `blocks`.
- If the property is a `HalfSpace` (resp. `HPolyhedron`), the function `project`
  from `LazySets` returns a new property with a `HalfSpace` (resp. `HPolyhedron`)
  in the reduced dimensions, according to `blocks`.
- Otherwise, the dimensional reduction is implemented via a (lazy) `LinearMap`.
"""
function inout_map_property(𝑃::PROPERTY,
                            partition::AbstractVector{<:AbstractVector{Int}},
                            blocks::AbstractVector{Int},
                            n::Int
                           )::PROPERTY where {PROPERTY<:Property}

    # create a sorted list of all dimensions in `blocks` => available variables for projection
    proj = vcat(partition[blocks]...)

    # no change in the dimension => return the original property
    length(proj) == n && return 𝑃

    return inout_map_property_helper(𝑃, proj)
end

function inout_map_property_helper(𝑃::Conjunction, proj::Vector{Int})
    new_conjuncts = similar(𝑃.conjuncts)
    for (i, conjunct) in enumerate(𝑃.conjuncts)
        new_conjuncts[i] = inout_map_property_helper(conjunct, proj)
    end
    return Conjunction(new_conjuncts)
end

function inout_map_property_helper(𝑃::Disjunction, proj::Vector{Int})
    new_disjuncts = similar(𝑃.disjuncts)
    for (i, disjunct) in enumerate(𝑃.disjuncts)
        new_disjuncts[i] = inout_map_property_helper(disjunct, proj)
    end
    return Disjunction(new_disjuncts)
end

function inout_map_property_helper(𝑃::BadStatesProperty, proj::Vector{Int})
    return BadStatesProperty(project(𝑃.bad, proj))
end

function inout_map_property_helper(𝑃::SafeStatesProperty, proj::Vector{Int})
    return SafeStatesProperty(project(𝑃.safe, proj))
end

ext: jl


text:
@testset "RoutingData" begin
    g = MetaDiGraph() # Directed graph.
    add_vertex!(g, :name, "A")
    add_vertex!(g, :name, "B")
    add_vertex!(g, :name, "C")
    add_vertex!(g, :name, "D")
    add_edge!(g, 1, 2, :capacity, 1.0)
    add_edge!(g, 1, 3, :capacity, 1.0)
    add_edge!(g, 2, 4, :capacity, 1.0)
    add_edge!(g, 3, 4, :capacity, 1.0)

    ug = MetaGraph() # Undirected graph (otherwise same as above).
    add_vertex!(ug, :name, "A")
    add_vertex!(ug, :name, "B")
    add_vertex!(ug, :name, "C")
    add_vertex!(ug, :name, "D")
    add_edge!(ug, 1, 2, :capacity, 1.0)
    add_edge!(ug, 1, 3, :capacity, 1.0)
    add_edge!(ug, 2, 4, :capacity, 1.0)
    add_edge!(ug, 3, 4, :capacity, 1.0)

    g_unreachable = MetaDiGraph()
    add_vertex!(g_unreachable, :name, "A")
    add_vertex!(g_unreachable, :name, "B")
    add_vertex!(g_unreachable, :name, "C")
    add_vertex!(g_unreachable, :name, "D")
    add_edge!(g_unreachable, 1, 2, :capacity, 1.0)
    add_edge!(g_unreachable, 2, 4, :capacity, 1.0)
    # Node C has no edge.

    k = SimpleDiGraph(4)
    add_edge!(k, 1, 4)

    uk = SimpleGraph(4)
    add_edge!(k, 1, 4)

    k_unreachable = SimpleDiGraph(4)
    add_edge!(k_unreachable, 1, 3)

    @testset "Helpers" begin
        @testset "remove_unreachable_nodes!" begin
            g_ = copy(g_unreachable)
            @test remove_unreachable_nodes!(g_) == 1
            @test ne(g_) == 2
            @test nv(g_) == 3
        end

        @testset "remove_unsatisfiable_demands!" begin
            g_ = copy(g_unreachable)
            k_ = copy(k_unreachable)
            @test remove_unsatisfiable_demands!(g_, k_) == 1
            @test ne(g_) == 2
            @test nv(g_) == 4
            @test ne(k_) == 0
        end
    end

    @testset "Parameters" begin
        @testset "Undirected graph" begin
            mt = ModelType(Load(), MinimumTotal(), PathFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            @test_throws ErrorException RoutingData(copy(ug), copy(k), ECOS.Optimizer, mt)
            @test_throws ErrorException RoutingData(copy(g), copy(uk), ECOS.Optimizer, mt)
            @test_throws ErrorException RoutingData(copy(ug), copy(uk), ECOS.Optimizer, mt)
        end

        @testset "Nonexistent folder" begin
            mt = ModelType(Load(), MinimumTotal(), PathFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            @test_throws ErrorException RoutingData(copy(g), copy(k), ECOS.Optimizer, mt, output_folder="/42/84/bullshit/")
        end

        @testset "Unsatisfiable or unreachable parts" begin
            mt = ModelType(Load(), MinimumTotal(), PathFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())

            tmp = AbstractString[]
            lm = (x) -> push!(tmp, x)
            RoutingData(copy(g_unreachable), copy(k), ECOS.Optimizer, mt, verbose=true, logmessage=lm)
            @test length(tmp) == 2
            @test "Removed 1 vertices because they were not reachable." in tmp
            @test "Removed 1 demands because they were not routable." in tmp
        end
    end

    @testset "RoutingData" begin
        set = [("Flow-based", FlowFormulation()), ("Path-based", PathFormulation())]

        @testset "$name" for (name, type) in set
            mt = ModelType(Load(), MinimumTotal(), type, false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            rd = RoutingData(copy(g), copy(k), ECOS.Optimizer, mt)

            @test ne(rd.g) == ne(g)
            @test nv(rd.g) == nv(g)
            @test ne(rd.k) == ne(k)
            @test nv(rd.k) == nv(k)
            @test rd.solver == ECOS.Optimizer
            @test rd.name == ""
            @test rd.time_precompute_ms > 0.0
            @test rd.model_type == mt
            @test rd.sub_model_type == mt
            @test length(rd.traffic_matrix) == 0 # TODO: is this right?

            @test length(edges(rd)) == 4
            @test length(vertices(rd)) == 4
            @test length(demands(rd)) == 1
            @test n_nodes(rd) == 4
            @test n_edges(rd) == 4
            @test n_demands(rd) == 1
            @test n_paths(rd) == (type == PathFormulation() ? 2 : 0)

            if type == PathFormulation()
                @test [Edge(1, 2), Edge(2, 4)] in rd.paths_edges
                @test [Edge(1, 3), Edge(3, 4)] in rd.paths_edges
                @test length(rd.demand_to_path_ids[Edge(1, 4)]) == 2
                @test rd.path_id_to_demand[1] == Edge(1, 4)
                @test rd.path_id_to_demand[2] == Edge(1, 4)
            else
                @test type == FlowFormulation()
                @test length(rd.paths_edges) == 0
                @test length(rd.demand_to_path_ids) == 0
                @test length(rd.path_id_to_demand) == 0
            end

            # TODO: what about the submodel type? Testing that it is essentially the same as the main object, also try to generate a different one.
        end

        @testset "Warning when generating too many paths" begin
            # Flow formulation.
            mt = ModelType(Load(), MinimumTotal(), FlowFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            @test_nowarn RoutingData(g, k, ECOS.Optimizer, mt, npaths=1)

            # Without column generation: warning.
            mt = ModelType(Load(), MinimumTotal(), PathFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            wrnmsg = "The limit of paths to generate (1) has been reached; the result will not be " *
                  "guaranteed to be optimum. Switch to a column-generation formulation for exact results."
            @test_logs (:warn, wrnmsg) RoutingData(g, k, ECOS.Optimizer, mt, npaths=1)

            # With column generation: no warning.
            mt = ModelType(Load(), MinimumTotal(), PathFormulation(), true, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            @test_nowarn RoutingData(g, k, ECOS.Optimizer, mt, npaths=1)
        end

        @testset "find_path_ids_with_edge" begin
            # Flow formulation: throw errors.
            mt = ModelType(Load(), MinimumTotal(), FlowFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            rd = RoutingData(copy(g), copy(k), ECOS.Optimizer, mt)
            @test_throws ErrorException find_path_ids_with_edge(rd, Edge(1, 2))

            # Path formulation.
            mt = ModelType(Load(), MinimumTotal(), PathFormulation(), true, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
            rd = RoutingData(copy(g), copy(k), ECOS.Optimizer, mt)

            p_id = find_path_ids_with_edge(rd, Edge(1, 2))
            @test length(p_id) == 1
            @test Edge(1, 2) in rd.paths_edges[p_id[1]]

            p_id = find_path_ids_with_edge(rd, Edge(1, 5))
            @test length(p_id) == 0
        end
    end

    @testset "Objective functions: support_min/max" begin
        @test ! supports_min(TestEdgeWise())
        @test ! supports_max(TestEdgeWise())

        @test supports_min(Load())
        @test supports_max(Load())

        @test supports_min(KleinrockLoad())
        @test ! supports_max(KleinrockLoad())

        @test supports_min(KleinrockLoad(true))
        @test supports_max(KleinrockLoad(true))

        @test supports_min(FortzThorupLoad())
        @test ! supports_max(FortzThorupLoad())

        @test supports_min(AlphaFairness(0.0))
        @test supports_max(AlphaFairness(0.0))

        @test ! supports_min(AlphaFairness(0.5))
        @test supports_max(AlphaFairness(0.5))
    end
end


text:
### Question 1

# I make this global variable `const` to ensure code stability
const codes = Dict('U'=>(1,0), 'D'=>(-1,0), 'R'=>(0,1), 'L'=>(0,-1))

function parsecode(code)
    op = code[1]
    dist = parse(Int, code[2:end])
    codes[op], dist   # captures the global variable `codes`
end

function evolve(line)
    ret = Tuple{Int, Int}[]
    pos = (0, 0)
    for code in line
        dir, dist = parsecode(code)
        for i in 1:dist
            pos = pos .+ dir
            push!(ret, pos)
        end
    end
    ret
end

manhattan(x) = sum(abs.(x))

function runlines(l1, l2)
    crosses = intersect(evolve(l1), evolve(l2))
    findmin(manhattan.(crosses))[1]
end

getlines(x::Vector) = split.(x, Ref(','))
# test case

using Test
test1 = """
R75,D30,R83,U83,L12,D49,R71,U7,L72
U62,R66,U55,R34,D71,R55,D58,R83
"""
l1, l2 = getlines(split(test1, '\n'))
@test runlines(l1, l2) == 159

test2 = """
R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51
U98,R91,D20,R16,D67,R40,U7,R15,U6,R7
"""
l1, l2 = getlines(split(test2, '\n'))
@test runlines(l1, l2) == 135

# Final result
l1, l2 = getlines(readlines("data/Dec3_data.txt"))
runlines(l1, l2)

### Question 2
function runlines2(l1, l2)
    line1, line2 = evolve(l1), evolve(l2)
    crosses = intersect(line1, line2)
    dist = [findfirst(==(x), line1) + findfirst(==(x), line2) for x in crosses]
    findmin(dist)[1]
end

runlines2(l1, l2)


text:
using BandedMatrices, LinearAlgebra, UnicodePlots
using GPUifyLoops, StaticArrays
using Test

using CUDAapi # this will NEVER fail
if has_cuda()
  try
    using CUDAnative, CuArrays
    CuArrays.allowscalar(false)
  catch ex
    # something is wrong with the user's set-up (or there's a bug in CuArrays)
    @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())
  end
end

function forward!(b, LU::AbstractArray{T,N}, ::Val{Nq}, ::Val{Nfields},
                  ::Val{Ne_vert}, ::Val{Ne_horz},
                  ::Val{element_bandwidth}) where {T, N, Nq, Nfields,
                                                   Ne_vert, Ne_horz,
                                                   element_bandwidth}
  FT = eltype(b)
  n = Nfields * Nq * Ne_vert
  p = q = element_bandwidth * Nfields * Nq

  l_b = MArray{Tuple{p+1}, FT}(undef)

  @inbounds @loop for h in (1:Ne_horz; blockIdx().x)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        @unroll for v = 1:element_bandwidth
          @unroll for k = 1:Nq
            @unroll for f = 1:Nfields
              ii  = f + (k - 1) * Nfields + (v - 1) * Nfields * Nq
              l_b[ii] =  Ne_vert ≥ v ? b[i, j, k, f, v, h] : zero(FT)
            end
          end
        end
        l_b[p+1] = Ne_vert > element_bandwidth ?
            b[i, j, 1, 1, element_bandwidth+1, h] : zero(FT)

        for v = 1:Ne_vert
          @unroll for k = 1:Nq
            @unroll for f = 1:Nfields
              jj = f + (k - 1) * Nfields + (v - 1) * Nfields * Nq

              @unroll for ii = 2:p+1
                Lii = N == 2 ? LU[ii+q, jj] : LU[i, j, ii+q, jj, h]
                l_b[ii] -= Lii * l_b[1]
              end

              b[i, j, k, f, v, h] = l_b[1]

              @unroll for ii = 1:p
                l_b[ii] = l_b[ii + 1]
              end

              # idx = jj + p
              # fi = (idx % Nfields) + 1
              # idx2 = idx ÷ Nfields
              # ki = (idx2 % Nq) + 1
              # vi = (idx2 ÷ Nq) + 1

              if jj + p < n
                (idx, fi) = fldmod1(jj + p + 1, Nfields)
                (vi, ki) = fldmod1(idx, Nq)

                l_b[p + 1] = b[i, j, ki, fi, vi, h]
              end
            end
          end
        end
      end
    end
  end

  nothing
end

function backward!(b, LU::AbstractArray{T, N}, ::Val{Nq}, ::Val{Nfields},
                   ::Val{Ne_vert}, ::Val{Ne_horz},
                   ::Val{element_bandwidth}) where {T, N, Nq, Nfields,
                                                    Ne_vert, Ne_horz,
                                                    element_bandwidth}
  FT = eltype(b)
  n = Nfields * Nq * Ne_vert
  q = Nfields * Nq * element_bandwidth

  l_b = MArray{Tuple{q+1}, FT}(undef)

  @inbounds @loop for h in (1:Ne_horz; blockIdx().x)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        @unroll for v = Ne_vert:-1:(Ne_vert - element_bandwidth + 1)
          @unroll for k = Nq:-1:1
            @unroll for f = Nfields:-1:1
              ii  = f + (k - 1) * Nfields + (Ne_vert - v) * Nfields * Nq
              l_b[ii+1] =  b[i, j, k, f, v, h]
            end
          end
        end
        l_b[1] = Ne_vert - element_bandwidth > 0 ?
            b[i, j, Nq, Nfields, Ne_vert - element_bandwidth, h] : zero(FT)

        for v = Ne_vert:-1:1
          @unroll for k = Nq:-1:1
            @unroll for f = Nfields:-1:1
              jj = f + (k - 1) * Nfields + (v - 1) * Nfields * Nq

              l_b[q + 1] /= N == 2 ? LU[q + 1, jj] : LU[i, j, q + 1, jj, h]

              @unroll for ii = 1:q
                Uii = N == 2 ? LU[ii, jj] : LU[i, j, ii, jj, h]
                l_b[ii] -= Uii * l_b[q + 1]
              end

              b[i, j, k, f, v, h] = l_b[q + 1]

              @unroll for ii = q:-1:1
                l_b[ii+1] = l_b[ii]
              end

              if jj - q  > 1
                (idx, fi) = fldmod1(jj - q - 1, Nfields)
                (vi, ki) = fldmod1(idx, Nq)

                l_b[1] = b[i, j, ki, fi, vi, h]
              end
            end
          end
        end
      end
    end
  end

  nothing
end

function band_lu!(A, ::Val{Nq}, ::Val{Nfields}, ::Val{Ne_vert},
                  ::Val{Ne_horz},
                  ::Val{element_bandwidth}) where {Nq, Nfields, Ne_vert,
                                                   Ne_horz, element_bandwidth}
  FT = eltype(A)
  n = Nfields * Nq * Ne_vert
  p = q = Nfields * Nq * element_bandwidth

  @inbounds @loop for h in (1:Ne_horz; blockIdx().x)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        for v = 1:Ne_vert
          for k = 1:Nq
            for f = 1:Nfields
              kk = f + (k - 1) * Nfields + (v - 1) * Nfields * Nq

              Aq = A[i, j, q + 1, kk, h]
              for ii = 1:p
                A[i, j, q + ii + 1, kk, h] /= Aq
              end

              for jj = 1:q
                if jj + kk ≤ n
                  Ajj = A[i, j, q - jj + 1, jj + kk, h]
                  for ii = 1:p
                    A[i, j, q + ii - jj + 1, jj + kk, h] -=
                        A[i, j, q + ii + 1, kk, h] * Ajj
                  end
                end
              end
            end
          end
        end
      end
    end
  end
end

function band_forward!(b, L, p)
  n = size(b, 1)
  for j = 1:n, i = (j+1):min(j+p, n)
    b[i] -= L[i, j]*b[j]
  end
end

function band_backward!(b, U, q)
  n = size(b, 1)
  for j = n:-1:1
    b[j] /= U[j, j]
    for i = max(1,j-q):(j-1)
      b[i] -= U[i, j]*b[j]
    end
  end
end

function extract_banded!(B::AbstractMatrix, A::AbstractMatrix, p, q)
  @assert p ≥ 0
  @assert q ≥ 0

  m, n = size(A)

  for j = 1:n, i = max(1, j - q):min(j + p, n)
    B[q + i - j + 1, j] = A[i, j]
  end

  return B
end

function extract_banded(A::AbstractMatrix, p, q)
  B = similar(A, p + q + 1, size(A, 2))
  fill!(B, 0)
  extract_banded!(B, A, p, q)
end

let

  T = Float64
  A = T[2 1 0
       1 2 1
       0 1 2]
  F = lu(A, Val(false))

  m = n = size(A, 1)
  p = q = 1

  b = T[4;8;8]

  A\b

  bb = copy(b)
  band_forward!(bb, F.L, p)
  band_backward!(bb, F.U, q)

end

let
  Nq = 2
  Nfields = 2
  Ne_vert = 10
  Ne_horz = 2
  EB = 2

  FT = Float64
  m = n = Nq * Nfields * Ne_vert
  p = q = Nq * Nfields * EB

  A = Array(brand(FT, m, n, p, q)) + 10I
  F = lu(A, Val(false))

  A = extract_banded(A, p, q)
  L = extract_banded(F.L, p, 0)
  U = extract_banded(F.U, 0, q)

  A = repeat(reshape(repeat(A, inner=(Nq*Nq, 1)), Nq, Nq, p+q+1, n), outer=(1,1,1,1,Ne_horz))
  L = repeat(reshape(repeat(L, inner=(Nq*Nq, 1)), Nq, Nq, p+1,   n), outer=(1,1,1,1,Ne_horz))
  U = repeat(reshape(repeat(U, inner=(Nq*Nq, 1)), Nq, Nq, q+1,   n), outer=(1,1,1,1,Ne_horz))

  threads = (Nq, Nq)
  blocks = Ne_horz
  @launch(CPU(), threads=threads, blocks=blocks, band_lu!(A, Val(Nq), Val(Nfields), Val(Ne_vert), Val(Ne_horz), Val(EB)))

  LL = A[:, :, q+1:end, :, :]
  LL[:, :, 1,:, :] .= one(eltype(LL))
  UU = A[:, :, 1:q+1, :, :]

  @test U ≈ UU
  @test L ≈ LL

end

let

  N = 5
  Nq = N + 1

  # H-S
  # Nfields = 5
  # Ne_vert = 8
  # Ne_horz = 6*6*6

  # Dycoms
  Nfields = 5
  Ne_vert = 60
  Ne_horz = 22*22
  EB = 1

  FT = Float64
  m = n = Nq * Nfields * Ne_vert
  p = q = Nq * Nfields * EB

  A = Array(brand(FT, m, n, p, q)) + 10I
  F = lu(A, Val(false))

  # display(spy(A))
  # display(spy(F.L))
  # display(spy(F.U))

  A = extract_banded(A, p, q)
  L = extract_banded(F.L, p, 0)
  U = extract_banded(F.U, 0, q)

  b = rand(FT, Nq, Nq, Nq, Nfields, Ne_vert, Ne_horz)
  borig = copy(b)
  x = similar(b)
  fill!(x, NaN)

  # This is the operation we want to reproduce with a GPUifyLoops kernel
  perm = (4, 3, 5, 1, 2, 6)
  xp = reshape(PermutedDimsArray(x, perm), Nfields*Nq*Ne_vert, Nq*Nq*Ne_horz)
  bp = reshape(PermutedDimsArray(b, perm), Nfields*Nq*Ne_vert, Nq*Nq*Ne_horz)
  xp .= F \ bp

  # Move data to the GPU
  DefaultArray = has_cuda() ? CuArray : Array
  device = has_cuda() ? CUDA() : CPU()

  # d_L = repeat(reshape(repeat(DefaultArray(L), inner=(Nq*Nq, 1)), Nq, Nq, p+1, n), outer=(1,1,1,1,Ne_horz))
  # d_U = repeat(reshape(repeat(DefaultArray(U), inner=(Nq*Nq, 1)), Nq, Nq, q+1, n), outer=(1,1,1,1,Ne_horz))

  d_A = DefaultArray(repeat(reshape(repeat(A, inner=(Nq*Nq, 1)), Nq, Nq, p+q+1, n), outer=(1,1,1,1,Ne_horz)))

  threads = (Nq, Nq)
  blocks = Ne_horz
  @launch(device, threads=threads, blocks=blocks, band_lu!(d_A, Val(Nq), Val(Nfields), Val(Ne_vert), Val(Ne_horz), Val(EB)))
  d_LU = d_A

  # d_L = d_A[:, :, q+1:end, :, :]
  # d_L[:, :, 1,:, :] .= one(eltype(d_L))
  # d_U = d_A[:, :, 1:q+1, :, :]

  d_b = DefaultArray(b)

  threads = (Nq, Nq)
  blocks = Ne_horz

  LU = Array(d_LU)
  @launch(CPU(), threads=threads, blocks=blocks, forward!(b, LU, Val(Nq), Val(Nfields), Val(Ne_vert), Val(Ne_horz), Val(EB)))
  @launch(CPU(), threads=threads, blocks=blocks, backward!(b, LU, Val(Nq), Val(Nfields), Val(Ne_vert), Val(Ne_horz), Val(EB)))

  @test x ≈ Array(b)

  @launch(device, threads=threads, blocks=blocks, forward!(d_b, d_LU, Val(Nq), Val(Nfields), Val(Ne_vert), Val(Ne_horz), Val(EB)))
  @launch(device, threads=threads, blocks=blocks, backward!(d_b, d_LU, Val(Nq), Val(Nfields), Val(Ne_vert), Val(Ne_horz), Val(EB)))

  @test x ≈ Array(d_b)
end


text:
# 'cstrt_tamvakis_ngon' provides the vertices coordinates (a,b) of the
# Tamvakis n-gon
# Please, cite...
#   N. K. Tamvakis. On the perimeter and the area of the convex polygon of
#   a given diameter. Bull. Greek Math. Soc, 28: 115-132, 1987.
function cstrt_tamvakis_ngon(n)
if n>=3
    a = zeros(n-1); b = zeros(n-1);
    if mod(n,3) == 0
        a[Int(n/3)] = 1/2; a[Int(2*n/3)] = -a[Int(n/3)];
        b[Int(n/3)] = sqrt(3)/2; b[Int(2*n/3)] = b[Int(n/3)];
        if n>3
            for j = 1:Int((n-3)/3)
                a[Int(n/3+j)] = cos(pi/3+j*pi/n);
                b[Int(n/3+j)] = sin(pi/3+j*pi/n);
                a[Int(2*n/3+j)] = a[Int(n/3)] + cos(pi+j*pi/n);
                b[Int(2*n/3+j)] = b[Int(n/3)] + sin(pi+j*pi/n);
                a[Int(n/3-j)] = -a[Int(2*n/3+j)];
                b[Int(n/3-j)] = b[Int(2*n/3+j)];
            end
        end
    elseif mod(n,3) == 1
        a[Int((n-1)/3)] = 1/2; a[Int((2*n+1)/3)] = -a[Int((n-1)/3)];
        b[Int((n-1)/3)] = sqrt(3)/2; b[Int((2*n+1)/3)] = b[Int((n-1)/3)];
        for j = 1:Int((n-1)/3)
            a[Int((n-1)/3+j)] = cos(pi/3+j*pi/(n+2));
            b[Int((n-1)/3+j)] = sin(pi/3+j*pi/(n+2));
        end
        if n>4
            for j = 1:Int((n-4)/3)
                a[Int((2*n+1)/3+j)] = a[Int((n-1)/3)] + cos(pi+j*pi/(n-1));
                b[Int((2*n+1)/3+j)] = b[Int((n-1)/3)] + sin(pi+j*pi/(n-1));
                a[Int((n-1)/3-j)] = -a[Int((2*n+1)/3+j)];
                b[Int((n-1)/3-j)] = b[Int((2*n+1)/3+j)];
            end
        end
    else
        a[Int((n+1)/3)] = 1/2; a[Int((2*n-1)/3)] = -a[Int((n+1)/3)];
        b[Int((n+1)/3)] = sqrt(3)/2; b[Int((2*n-1)/3)] = b[Int((n+1)/3)];
        for j = 1:Int((n-2)/3)
            a[Int((2*n-1)/3+j)] = a[Int((n+1)/3)] + cos(pi+j*pi/(n+1));
            b[Int((2*n-1)/3+j)] = b[Int((n+1)/3)] + sin(pi+j*pi/(n+1));
            a[Int((n+1)/3-j)] = -a[Int((2*n-1)/3+j)];
            b[Int((n+1)/3-j)] = b[Int((2*n-1)/3+j)];
        end
        if n>5
            for j = 1:Int((n-5)/3)
                a[Int((n+1)/3+j)] = cos(pi/3+j*pi/(n-2));
                b[Int((n+1)/3+j)] = sin(pi/3+j*pi/(n-2));
            end
        end
    end
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
# Step 1
struct Fib
    n::Int
end

# Step 2
function Base.iterate(fib::Fib, state=(i=1, vals=(1, 1)))
    state.i > fib.n && return nothing

    if state.i <= 2
        return 1, (i=state.i + 1, vals=(1, 1))
    end

    aₙ₋₂, aₙ₋₁ = state.vals
    aₙ = aₙ₋₂ + aₙ₋₁

    aₙ, (i=state.i + 1, vals=(aₙ₋₁, aₙ))
end

# Step 3
Base.length(fib::Fib) = fib.n

# Step 4
Base.eltype(::Type{Fib}) = Int


text:
# Test for gaussjacobi().
tol = 1e-14
# Test a small n:
n = 42
a = -.1
b = .3
x, w = gaussjacobi(n, a, b)
@test all(length(x) == n) && all(length(w) == n)
@test abs(x[37] - 0.912883347814032) < tol
@test abs(w[37] - 0.046661910947553) < tol

@test dot(w,exp.(x)) ≈ 2.7568520361985516

# Test a larger n (using ASY)
a = -.7
b = 1.3
n = 251
x, w = gaussjacobi(n, a, b)
@test all(length(x) == n) && all(length(w) == n)
@test abs(x[37] + 0.893103435898983) < tol
@test abs(w[37] - 1.962534523788093e-04) < tol

@test dot(w,exp.(x)) ≈ 16.722957039404044

# Test n = 1:
a = 1.0
b = 2.0
x, w = gaussjacobi(1, a, b)
@test abs(x[1] - (b - a) / (a + b + 2)) < tol
@test abs(w[1] - 2^(a + b + 1) * beta(a + 1, b + 1)) < tol

@test dot(w,fill(1,length(x))) ≈ 1.3333333333333333

x, w = gaussjacobi(1013, .9, -.1)
@test abs(x[2] + 0.999986012231899) < tol
@test abs(x[13] + 0.999225722939832) < tol
@test abs(w[2] - 9.314674169892358e-05) < tol
@test abs(w[13] - 4.654651764553262e-04) < tol

@test dot(w,exp.(x)) ≈ 1.6915068974063106

x, w = gaussjacobi(10013, .9, -.1)
@test abs(x[2]  +0.999999856605293) < tol
@test abs(x[13]  + 0.999992061552711) < tol
@test abs(w[2] - 1.509654630405615e-06) < tol
@test abs(w[13] - 7.548275262993863e-06) < tol

@test dot(w,exp.(x)) ≈ 1.6915068974063106

# tests bug where row vectors were returned
@test isa(x,Vector{Float64})
@test isa(w,Vector{Float64})

# test last alpha and beta parameters:
x, w = gaussjacobi(100, 19., 21.)
@test abs(x[87] - 0.832211446176040) < tol
@test abs(w[50] - 0.026363584978877305) < tol # new test without condition sum(w)==1

# test for small alpha and beta:
x, w = gaussjacobi(10000, .1, .2)
@test abs(x[1] - -0.999999963363548) < tol
@test abs(w[500] - 2.183393039546711e-05) < tol

x,w=gaussjacobi(20,11.,0.)
# tests bug where row vectors were returned
@test isa(x,Vector{Float64})
@test isa(w,Vector{Float64})

# test JacobiRec against a precomputed rule in BigFloat
x,w = FastGaussQuadrature.JacobiRec(10, big(0), big(0))
@test abs(x[6] - big(1488743389816312108848260011297199846175648594206916957079892535159036173556674)/big(10)^79) < 1e-70
@test abs(w[4] - big(2692667193099963550912269215694693528597599384608837958005632762421534323191819)/big(10)^79) < 1e-70

x,w = FastGaussQuadrature.JacobiRec(10, big(2//10), big(-1//30))
@test abs(x[7] - big(4146701176053244724446765533149623814305606864483746748725151214196203624388371)/big(10)^79) < 1e-70
@test abs(w[3] - big(2482452398859023537636458227096017466104571146686194593251519240011574719667464)/big(10)^79) < 1e-70


text:
# TODO Make a PR to Distributions.jl

using Distributions
import Distributions.@check_args
################################################################
function Distributions.Uniform(::Type{T}, a, b) where {T <: Real}
    return Uniform(T(a), T(b))
end
################################################################
struct Uniform2{T}
    a::T
    b::T

    # Uniform2{T}(a::T, b::T) constructor
    function Uniform2{T}(a::T, b::T; check_args=true) where {T <: Real}
        check_args && @check_args(Uniform2, a < b)
        return new{T}(a, b)
    end

    # Abstract a,b - Uniform2{T}(a, b) constructor
    function Uniform2{T}(a::Real, b::Real; check_args=true) where {T <: Real}
        check_args && @check_args(Uniform2, a < b)
        return new{T}(T(a), T(b))
    end

    # Uniform2{T}(a::T, b::T) constructor
    function Uniform2{T}(a::T, b::T) where {T <:Complex}
        new{T}(a, b)
    end

    # Abstract a,b - Uniform2{T}(a, b) constructor
    function Uniform2{T}(a, b) where {T <: Complex}
        return new{T}(T(a), T(b))
    end
end
# Real

# zeros() like constructor (Uniform2(T, a::T, b::T))
function Uniform2(::Type{T}, a::T, b::T; check_args=true) where {T <: Real}
    return Uniform2{T}(a, b, check_args = check_args)
end

# Abstract a,b - zeros() like constructor (Uniform2(T, a, b))
function Uniform2(::Type{T}, a, b; check_args=true) where {T <: Real}
    return Uniform2{T}(T(a), T(b), check_args = check_args)
end

# No type specified constructor:
function Uniform2(a::Float64, b::Float64; check_args=true)
    return Uniform2{Float64}(a, b, check_args = check_args)
end

# Abstract a,b - no type specified constructor:
function Uniform2(a, b; check_args=true)
    return Uniform2{Float64}(Float64(a), Float64(b), check_args = check_args)
end

# Complex

# zeros() like constructor (Uniform2(T, a::T, b::T))
function Uniform2(::Type{T}, a::T, b::T) where {T <: Complex}
    return Uniform2{T}(a, b)
end

# Abstract a,b - zeros() like constructor (Uniform2(T, a, b))
function Uniform2(::Type{T}, a, b) where {T <: Complex}
    return Uniform2{T}(T(a), T(b))
end

# No type specified constructor:
function Uniform2(a::ComplexF64, b::ComplexF64)
    return Uniform2{ComplexF64}(a, b)
end

Base.rand(d::Uniform2{T}, dims::Integer...) where {T} =   d.a .+ (d.b - d.a) .* Base.rand(T, dims)


text:
function read_csv(csv_path::String)
    readstr = open(f->read(f, String), csv_path, "r")
    rows = split(readstr, '\n')
    columns = split(rows[1], ',')
    values = [[] for col in columns]
    columns = [Symbol(c) for c in Array{String}(columns)]
    deleteat!(rows, 1)
    for row in rows
        data = (split(row, ','))
        for col in 1:length(columns)
            push!(values[col], data[col])
        end
    end
    return(columns, values)
end

function read_types(columns)
    types = []
    for col in columns
        try
            parse(Bool, col[1])
            parse(Bool, col[2])
            push!(types, Bool)
        catch
        try
            parse(Int64, col[1])
            parse(Int64, col[2])
            push!(types, Float64)
        catch
            try
                parse(Float64, col[1])
                parse(Float64, col[2])
                push!(types, Int64)
            catch
                push!(types, String)
            end
            end
        end
    end
    newcolumns = []
    for (n, col) in enumerate(columns)
        if types[n] != String
            push!(newcolumns, [parse(types[n], val) for val in col])

        else
            push!(newcolumns, col)
        end

    end
    return(types, newcolumns)
end


text:
export sample, gibbs_sampling, gibbs_sampling2, gibbs_sampling_query, get_sho_samples

export GMSampler, Gibbs

export SamplingRegime, S_regime, T_regime, M_regime

using StatsBase
using LinearAlgebra

abstract type GMSampler end

struct Gibbs <: GMSampler end

abstract type SamplingRegime end

# Time-Series and Histogramming
struct T_regime <: SamplingRegime end

# Multiple Restarts and Histogramming
struct M_regime <: SamplingRegime end

# Time-Series without Histogramming
struct S_regime <: SamplingRegime end

function int_to_spin(int_representation::Int, spin_number::Int)
    spin = 2*digits(int_representation, base=2, pad=spin_number) .- 1
    return spin
end

function weigh_proba(int_representation::Int, adj::Array{T,2}, prior::Array{T,1}) where T <: Real
    spin_number = size(adj,1)
    spins = int_to_spin(int_representation, spin_number)
    return exp(((0.5) * spins' * adj * spins + prior' * spins)[1])
end

bool_to_spin(bool::Int) = 2*bool-1

function weigh_proba(int_representation::Int, adj::Array{T,2}, prior::Array{T,1}, spins::Array{Int,1}) where T <: Real
    digits!(spins, int_representation, base=2)
    spins .= bool_to_spin.(spins)
    return exp((0.5 * spins' * adj * spins + prior' * spins)[1])
end

# assumes second order
function sample_generation_ising(gm::FactorGraph{T}, samples_per_bin::Integer, bins::Int) where T <: Real
    @assert bins >= 1

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector = copy(transpose(diag(adjacency_matrix)))[1,:]

    items   = [i for i in 0:(config_number-1)]
    assignment_tmp = [0 for i in 1:spin_number] # pre allocate assignment memory
    weights = [weigh_proba(i, adjacency_matrix, prior_vector, assignment_tmp) for i in (0:config_number-1)]

    raw_sample = StatsBase.sample(items, StatsBase.Weights(weights), samples_per_bin*bins, ordered=false)
    raw_sample_bins = reshape(raw_sample, bins, samples_per_bin)

    spin_samples = []
    for b in 1:bins
        # Find number of samples of one configuration
        raw_binning = countmap(raw_sample_bins[b,:])

        spin_sample = [ vcat(raw_binning[i], int_to_spin(i, spin_number)) for i in keys(raw_binning)]
        push!(spin_samples, hcat(spin_sample...)')
    end
    return spin_samples
end

function weigh_proba(int_representation::Int, gm::FactorGraph{T}, spins::Array{Int,1}) where T <: Real
    digits!(spins, int_representation, base=2)
    spins .= bool_to_spin.(spins)
    evaluation = sum( weight*prod(spins[i] for i in term) for (term, weight) in gm)
    return exp(evaluation)
end

function sample_generation(gm::FactorGraph{T}, samples_per_bin::Integer, bins::Int) where T <: Real
    @assert bins >= 1
    #info("use general sample model")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    items   = [i for i in 0:(config_number-1)]
    assignment_tmp = [0 for i in 1:spin_number] # pre allocate assignment memory
    weights = [weigh_proba(i, gm, assignment_tmp) for i in (0:config_number-1)]

    raw_sample = StatsBase.sample(items, StatsBase.Weights(weights), samples_per_bin*bins, ordered=false)
    raw_sample_bins = reshape(raw_sample, bins, samples_per_bin)

    spin_samples = []
    for b in 1:bins
        raw_binning = countmap(raw_sample_bins[b,:])
        spin_sample = [ vcat(raw_binning[i], int_to_spin(i, spin_number)) for i in keys(raw_binning)]
        push!(spin_samples, hcat(spin_sample...)')
    end
    return spin_samples
end

## Different methods to create histograms given samples from Glauber dynamics
function hist_countmap(samples_pairs_T::Array, spin_number::Integer, config_number::Integer)
    # Works well for smaller number of spins
    spin_pairs_configs = vec(collect(Iterators.product(1:spin_number,0:(config_number-1),0:(config_number-1))))
    dict_spin_pairs_configs = Dict(spin_pairs_configs[i+1] => i for i = 0:(length(spin_pairs_configs)-1))
    dict_spin_pairs_strings_configs = Dict(dict_spin_pairs_configs[spin_pairs_configs[i+1]] => vcat(collect(spin_pairs_configs[i+1])[1],2*digits(collect(spin_pairs_configs[i+1])[2],2,spin_number) - 1, 2*digits(collect(spin_pairs_configs[i+1])[3],2,spin_number) - 1) for i = 0:(length(spin_pairs_configs)-1))

    samples_pairs_configs_T = Array{Int64}(size(samples_pairs_T,1))
    for i = 1:size(samples_pairs_T,1)
        samples_pairs_configs_T[i] = dict_spin_pairs_configs[tuple(samples_pairs_T[i,1],samples_pairs_T[i,2],samples_pairs_T[i,3])]
    end

    raw_binning_T = countmap(samples_pairs_configs_T)
    samples_T = [ vcat(raw_binning_T[i], dict_spin_pairs_strings_configs[i]) for i in keys(raw_binning_T)]
    samples_T = hcat(samples_T...)'

    return samples_T
end

function countfreq(A)
    #=
    Ref: https://stackoverflow.com/questions/39101839/the-number-of-occurences-of-elements-in-a-vector-julia
    =#
    d = Dict{Array{eltype(A),1}, Int}()
    for i in 1:size(A,1)
        a = A[i,:]
        if a in keys(d)
            d[a] += 1
        else
            d[a] = 1
        end
    end
    return d
end

function hist_countfreq(samples_pairs_T::Array, spin_number::Integer, config_number::Integer)
    # Dictionary over observed configurations with corr frequency
    dict_obs_configs = countfreq(samples_pairs_T)

    samples_T = Array{Int64,2}(undef,length(dict_obs_configs),2+2*spin_number)
    keys_dict_obs_configs = [key for key in keys(dict_obs_configs)]
    for i = 1:length(keys_dict_obs_configs)
        obs_config = keys_dict_obs_configs[i]
        samples_T[i,:] = vcat(dict_obs_configs[obs_config],obs_config[1],2*digits(obs_config[2],base=2,pad=spin_number) .- 1,2*digits(obs_config[3],base=2,pad=spin_number) .- 1)'
    end

    return samples_T
end

function hist_countfreq2(samples_pairs_T::Array, spin_number::Integer, config_number::Integer)
    # Dictionary over observed configurations with corr frequency
    dict_obs_configs = countfreq(samples_pairs_T)

    samples_T = Array{Int64,2}(undef,length(dict_obs_configs),1+spin_number)
    keys_dict_obs_configs = [key for key in keys(dict_obs_configs)]
    for i = 1:length(keys_dict_obs_configs)
        obs_config = keys_dict_obs_configs[i]
        samples_T[i,:] = vcat(dict_obs_configs[obs_config],obs_config)'
    end

    return samples_T
end

function hist_countfreq_glauber_dynamics(samples_pairs_T::Array, spin_number::Integer, config_number::Integer)
    # Dictionary over observed configurations with corr frequency
    dict_obs_configs = countfreq(samples_pairs_T)

    samples_T = Array{Int64,2}(undef,length(dict_obs_configs),2+2*spin_number)
    keys_dict_obs_configs = [key for key in keys(dict_obs_configs)]
    for i = 1:length(keys_dict_obs_configs)
        obs_config = keys_dict_obs_configs[i]
        samples_T[i,:] = vcat(dict_obs_configs[obs_config],obs_config)'
    end

    return samples_T
end

# Create a dictionary out of a histogram
function dictionary_samples_hist(samples_T::Array)
    d = Dict{Array{eltype(samples_T),1}, Int}()
    for i in 1:size(samples_T,1)
        a = samples_T[i,2:end]
        if a in keys(d)
            d[a] += samples_T[i,1]
        else
            d[a] = samples_T[i,1]
        end
    end
    return d
end

# Add two different histograms (note that there may be new rows)
function add_histograms(samples_T1::Array, samples_T2::Array)
    d1 = dictionary_samples_hist(samples_T1)
    d2 = dictionary_samples_hist(samples_T2)

    for a in keys(d2)
        if a in keys(d1)
            d1[a] += d2[a]
        else
            d1[a] = d2[a]
        end
    end

    samples_T = Array{Int64,2}(undef, length(d1), size(samples_T1,2))
    keys_dict_obs_configs = [key for key in keys(d1)]
    for i = 1:length(keys_dict_obs_configs)
        obs_config = keys_dict_obs_configs[i]
        samples_T[i,:] = vcat(d1[obs_config],obs_config)'
    end

    return samples_T
end

# Convert samples to dictionary over nodes
function get_sho_samples(samples_T::Array)
    num_spins = (size(samples_T,2) - 2)
    num_spins = Int(num_spins/2)

    # Create a dictionary over samples_T
    dict_samples_sho = Dict{Int, Array{Int64, 2}}()

    # Node number is the key of the dictionary
    # Should include those which are flips and those which are not
    # For the cases where there are flips, we know node i is being updated because this is single-site dynamics and we see a change
    # Where there are no flips, we don't know which node is being updated so include all possibilities

    # For each dict[key], array should be of form (n_samples, FLAG_flip, \sigma^{(t)}, \sigma^{(t+1)})

    # First get all those samples where there is no flip
    ind_samples_T_no_flip = []
    ind_samples_T_flip = []

    for ind_config = 1:size(samples_T,1)
        if dot(samples_T[ind_config,3:(2+num_spins)],samples_T[ind_config,(3+num_spins):end]) == num_spins
            push!(ind_samples_T_no_flip,ind_config)
        else
            push!(ind_samples_T_flip,ind_config)
        end
    end

    samples_T_flip = samples_T[ind_samples_T_flip,:]
    samples_T_no_flip = samples_T[ind_samples_T_no_flip,:]

    samples_T_no_flip[:,2] .= 0
    # Now create the dictionary
    for current_spin = 1:num_spins
        samples_current_spin = samples_T_flip[findall(isequal(current_spin),samples_T_flip[:,2]),:]
        samples_current_spin[:,2] .= 1
        dict_samples_sho[current_spin] = vcat(samples_current_spin,samples_T_no_flip)
    end

    return dict_samples_sho
end

function gibbs_sampling_ising(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::S_regime) where T <: Real
    @info("using Glauber dynamics v1 to generate samples")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Allocate memory for matrix of samples (time-series)
    #spin_samples = Array{Int8, 2}(num_samples, spin_number)
    spin_samples = Array{Int64, 1}(undef, num_samples)

    # Allocate memory for node chosen for updating
    node_selected_samples = Array{Int8,1}(undef, num_samples-1)

    spin_configs = []
    for i = 0:(config_number - 1)
        spin_tmp = 2*digits(i, base = 2, pad = spin_number).-1
        push!(spin_configs, string(spin_tmp))
    end
    dict_spin_configs = Dict(spin_configs[i+1] => i for i = 0:(config_number-1))

    # Start Gibbs sampling

    # generate a random state initially
    sigma = [-1 for i in 1:spin_number]

    sigma_new = deepcopy(sigma)
    #spin_samples[1,:] = sigma
    spin_samples[1] = dict_spin_configs[string(sigma)]

    for ind_step = 1:(num_samples-1)
        i = rand(1:spin_number)
        #i = ((ind_step - 1) % spin_number) + 1

        # Proposed state
        sigma_new = deepcopy(sigma)
        sigma_new[i] = 1

        # Calculate the probability of accepting new state
        num_prob = exp( 2.0*(dot(adj_matrix[i,:],sigma) + prior_vector[i])[1] )
        denom_prob = 1.0 + num_prob
        acceptance_prob = num_prob / denom_prob

        # Accept new state or not
        r = rand()
        if r < acceptance_prob
            sigma = deepcopy(sigma_new)
        else
            sigma_new[i] = -1
            sigma = deepcopy(sigma_new)
        end

        # Update samples
        #spin_samples[ind_step+1,:] = sigma
        spin_samples[ind_step+1] = dict_spin_configs[string(sigma)]

        # Update node selected
        node_selected_samples[ind_step] = copy(i)
    end

    return spin_samples, node_selected_samples
end

function gibbs_sampling_ising(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::T_regime) where T <: Real
    @info("using Glauber dynamics v1 to generate T-regime samples")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Allocate memory for matrix of samples (time-series)
    spin_samples = Array{Int64, 1}(undef, num_samples)

    # Allocate memory for node chosen for updating
    node_selected_samples = Array{Int8,1}(undef, num_samples-1)

    spin_configs = []
    for i = 0:(config_number - 1)
        spin_tmp = 2*digits(i, base = 2, pad = spin_number).-1
        push!(spin_configs, string(spin_tmp))
    end
    dict_spin_configs = Dict(spin_configs[i+1] => i for i = 0:(config_number-1))

    # Start Gibbs sampling

    # generate a random state initially
    sigma_conf = rand(0:(config_number - 1))
    sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1

    sigma_new = deepcopy(sigma)
    #spin_samples[1,:] = sigma
    spin_samples[1] = dict_spin_configs[string(sigma)]

    for ind_step = 1:(num_samples-1)
        i = rand(1:spin_number)
        #i = ((ind_step - 1) % spin_number) + 1

        # Proposed state
        sigma_new = deepcopy(sigma)
        sigma_new[i] = 1

        # Calculate the probability of accepting new state
        num_prob = exp( 2.0*(dot(adj_matrix[i,:],sigma) + prior_vector[i])[1] )
        denom_prob = 1.0 + num_prob
        acceptance_prob = num_prob / denom_prob

        # Accept new state or not
        r = rand()
        if r < acceptance_prob
            sigma = deepcopy(sigma_new)
        else
            sigma_new[i] = -1
            sigma = deepcopy(sigma_new)
        end

        # Update samples
        spin_samples[ind_step+1] = dict_spin_configs[string(sigma)]

        # Update node selected
        node_selected_samples[ind_step] = copy(i)
    end

    ## T Samples
    samples_pairs_T = hcat(node_selected_samples,spin_samples[1:end-1],spin_samples[2:end])
    samples_T = hist_countfreq(samples_pairs_T, spin_number, config_number)

    ## Mixing the T samples
    raw_binning = countmap(spin_samples)
    samples_mixed = [ vcat(raw_binning[i], 2*digits(i, base=2, pad=spin_number).-1) for i in keys(raw_binning)]
    samples_mixed = hcat(samples_mixed...)'

    return samples_T, samples_mixed
end

function gibbs_sampling_ising(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::M_regime) where T <: Real
    @info("using Glauber dynamics v1 to generate M-regime samples")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Allocate memory for matrix of samples
    # Arranged as [node selected, \sigma^{(t)}, \sigma^{(t+1)}]
    samples_pairs_T = Array{Int64, 2}(undef, num_samples, 3)

    spin_configs = []
    for i = 0:(config_number - 1)
        spin_tmp = 2*digits(i, base=2, pad=spin_number).-1
        push!(spin_configs, string(spin_tmp))
    end
    dict_spin_configs = Dict(spin_configs[i+1] => i for i = 0:(config_number-1))

    # Start Gibbs Sampling for Multiple Restarts
    for ind_step = 1:num_samples
        # generate a random state initially
        sigma_conf = rand(0:(config_number - 1))
        sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1

        # Random generation of site to be changed
        i = rand(1:spin_number)
        #i = ((ind_step - 1) % spin_number) + 1

        # Proposed state
        sigma_new = deepcopy(sigma)
        sigma_new[i] = 1

        # Calculate the probability of accepting new state
        num_prob = exp( 2.0*(dot(adj_matrix[i,:],sigma) + prior_vector[i])[1] )
        denom_prob = 1.0 + num_prob
        acceptance_prob = num_prob / denom_prob

        # Accept new state or not
        r = rand()
        if r < acceptance_prob
            sigma = deepcopy(sigma_new)
        else
            sigma_new[i] = -1
            sigma = deepcopy(sigma_new)
        end

        # Update Samples Array
        samples_pairs_T[ind_step,1] = copy(i)
        samples_pairs_T[ind_step,2] = deepcopy(sigma_conf)
        samples_pairs_T[ind_step,3] = dict_spin_configs[string(sigma)]
    end

    ## M Samples
    samples_T = hist_countfreq(samples_pairs_T, spin_number, config_number)

    ## Mixing the M samples
    spin_samples = samples_pairs_T[:,3]
    raw_binning = countmap(spin_samples)
    samples_mixed = [ vcat(raw_binning[i], 2*digits(i, base=2, pad=spin_number).-1) for i in keys(raw_binning)]
    samples_mixed = hcat(samples_mixed...)'

    return samples_T, samples_mixed
end

function gibbs_sampling_ising2(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::S_regime) where T <: Real
    @info("using Glauber dynamics v2 to generate samples")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    assignment_tmp = [0 for i in 1:spin_number] # pre allocate assignment memory

    # Allocate memory for matrix of samples (time-series)
    spin_samples = Array{Int8, 2}(undef, num_samples, spin_number)

    # Allocate memory for node chosen for updating
    node_selected_samples = Array{Int64,1}(undef, num_samples-1)

    # Do not save dictionary over spin_configs as of length config_number

    # Start Gibbs sampling

    # generate a random state initially
    sigma = [-1 for i in 1:spin_number]

    sigma_new = copy(sigma)
    spin_samples[1,:] = sigma

    for ind_step = 1:(num_samples-1)
        i = ((ind_step - 1) % spin_number) + 1

        # Proposed state - flip ith state
        sigma_new = copy(sigma)
        sigma_new[i] = -1*sigma_new[i]

        # prob_sigma = prob_cond_ising(sigma, i, adjacency_matrix, prior_vector)
        # prob_sigma_new = prob_cond_ising(sigma_new, i, adjacency_matrix, prior_vector)

        # Calculate the probability of accepting new state which is also prob_sigma_new/prob_sigma
        acceptance_prob = exp( -2*sigma[i]*(adj_matrix[i,:]'*sigma + prior_vector[i])[1] )

        # Accept new state or not
        transitionProbability = min(1, acceptance_prob)
        if rand() < transitionProbability
            sigma = copy(sigma_new)
        end

        # Update samples
        spin_samples[ind_step+1,:] = sigma

        # Update node selected
        node_selected_samples[ind_step] = copy(i)
    end

    samples = hist_countfreq2(spin_samples, spin_number, config_number)

    return samples, node_selected_samples
end

function gibbs_sampling_ising2(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::T_regime) where T <: Real
    @info("using Glauber dynamics v2 to generate T-regime samples")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    assignment_tmp = [0 for i in 1:spin_number] # pre allocate assignment memory

    # Allocate memory for matrix of samples (time-series)
    spin_samples = Array{Int8, 2}(undef, num_samples, spin_number)

    # Allocate memory for node chosen for updating
    node_selected_samples = Array{Int8,1}(undef, num_samples-1)

    # Not creating dictionary over spin configs as of length config_number

    # Start Gibbs sampling

    # generate a random state initially
    sigma_conf = rand(0:(config_number - 1))
    sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1

    sigma_new = copy(sigma)
    spin_samples[1,:] = copy(sigma)

    for ind_step = 1:(num_samples-1)
        i = ((ind_step - 1) % spin_number) + 1

        # Proposed state - flip ith state
        sigma_new = copy(sigma)
        sigma_new[i] = -1*sigma_new[i]

        # prob_sigma = prob_cond_ising(sigma, i, adjacency_matrix, prior_vector)
        # prob_sigma_new = prob_cond_ising(sigma_new, i, adjacency_matrix, prior_vector)

        # Calculate the probability of accepting new state which is also prob_sigma_new/prob_sigma
        acceptance_prob = exp( -2*sigma[i]*(adj_matrix[i,:]'*sigma + prior_vector[i])[1] )

        # Accept new state or not
        transitionProbability = min(1, acceptance_prob)
        if rand() < transitionProbability
            sigma = copy(sigma_new)
        end

        # Update samples
        spin_samples[ind_step+1,:] = copy(sigma)

        # Update node selected
        node_selected_samples[ind_step] = copy(i)
    end

    # Histogramming for the T-regime
    samples_pairs_T = hcat(node_selected_samples,spin_samples[1:end-1,:],spin_samples[2:end,:])
    samples_T = hist_countfreq_glauber_dynamics(samples_pairs_T, spin_number, config_number)

    # Mixing the T samples
    samples_mixed = hist_countfreq2(spin_samples, spin_number, config_number)

    return samples_T, samples_mixed
end

function gibbs_sampling_ising2(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::M_regime) where T <: Real
    @info("using Glauber dynamics v2 to generate M-regime samples")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Allocate memory for matrix of samples
    # Arranged as [node selected, \sigma^{(t)}, \sigma^{(t+1)}]
    samples_pairs_T = Array{Int64, 2}(undef, num_samples, 1 + 2*spin_number)

    # Not creating dictionary over possible configurations as of large size

    # Start Gibbs Sampling for Multiple Restarts
    for ind_step = 1:num_samples
        # generate a random state initially
        sigma_conf = rand(0:(config_number - 1))
        sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1

        # Random generation of site to be changed
        i = rand(1:spin_number)
        #i = ((ind_step - 1) % spin_number) + 1

        # Proposed state - flip ith state
        sigma_old = deepcopy(sigma)
        sigma_new = deepcopy(sigma)
        sigma_new[i] = -1*sigma_new[i]

        # prob_sigma = prob_cond_ising(sigma, i, adjacency_matrix, prior_vector)
        # prob_sigma_new = prob_cond_ising(sigma_new, i, adjacency_matrix, prior_vector)

        # Calculate the probability of accepting new state which is also prob_sigma_new/prob_sigma
        acceptance_prob = exp( -2*sigma[i]*(adj_matrix[i,:]'*sigma + prior_vector[i])[1] )

        # Accept new state or not
        transitionProbability = min(1, acceptance_prob)
        if rand() < transitionProbability
            sigma = copy(sigma_new)
        end

        # Update Samples Array
        samples_pairs_T[ind_step,1] = copy(i)
        samples_pairs_T[ind_step,2:(spin_number+1)] = deepcopy(sigma_old)
        samples_pairs_T[ind_step,(spin_number+2):end] = deepcopy(sigma)
    end

    # Histogramming for the M-regime
    samples_T = hist_countfreq_glauber_dynamics(samples_pairs_T, spin_number, config_number)

    # Mixing the M samples
    spin_samples = samples_pairs_T[:,(spin_number+2):end]
    samples_mixed = hist_countfreq2(spin_samples, spin_number, config_number)

    return samples_T, samples_mixed
end

function gibbs_sampling_ising2_binning(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::T_regime) where T <: Real
    @info("using Glauber dynamics v2 to generate T-regime samples with binning")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Not creating dictionary over spin configs as of length config_number

    # Start Gibbs sampling

    # create bins of samples to look at as the overall number of samples is too large
    max_samples_per_bin = 10000000
    min_bins = div(num_samples,max_samples_per_bin)
    n_bins = copy(min_bins)
    rem_samples = rem(num_samples,max_samples_per_bin)
    if rem_samples > 0
        n_bins += 1
    end
    samples_bin_array = [i <= min_bins ? max_samples_per_bin : rem_samples for i=1:n_bins]

    # generate a random state initially
    sigma_conf = rand(0:(config_number - 1))
    sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1
    sigma_new = copy(sigma)

    # Some temporary structures
    samples_T = []
    samples_T_temp = []
    samples_mixed = []
    samples_mixed_temp = []

    # One bin at a time -- make sure to reuse the last computed sigma from prev bin to initialize new bin
    for ind_bin = 1:n_bins
        num_samples_bin = samples_bin_array[ind_bin]
        # Allocate memory for matrix of samples (time-series)
        spin_samples = Array{Int8, 2}(undef,num_samples_bin, spin_number)
        spin_samples[1,:] = copy(sigma)

        # Allocate memory for node chosen for updating
        node_selected_samples = Array{Int8,1}(undef,num_samples_bin-1)

        for ind_step = 1:(num_samples_bin-1)
            i = ((ind_step - 1) % spin_number) + 1

            # Proposed state - flip ith state
            sigma_new = copy(sigma)
            sigma_new[i] = -1*sigma_new[i]

            # prob_sigma = prob_cond_ising(sigma, i, adjacency_matrix, prior_vector)
            # prob_sigma_new = prob_cond_ising(sigma_new, i, adjacency_matrix, prior_vector)

            # Calculate the probability of accepting new state which is also prob_sigma_new/prob_sigma
            acceptance_prob = exp( -2*sigma[i]*(adj_matrix[i,:]'*sigma + prior_vector[i])[1] )

            # Accept new state or not
            transitionProbability = min(1, acceptance_prob)
            if rand() < transitionProbability
                sigma = copy(sigma_new)
            end

            # Update samples
            spin_samples[ind_step+1,:] = copy(sigma)

            # Update node selected
            node_selected_samples[ind_step] = copy(i)
        end

        # Histogramming for the T-regime
        samples_pairs_T = hcat(node_selected_samples,spin_samples[1:end-1,:],spin_samples[2:end,:])
        if ind_bin == 1
            samples_T = hist_countfreq_glauber_dynamics(samples_pairs_T, spin_number, config_number)
        else
            samples_T_temp = hist_countfreq_glauber_dynamics(samples_pairs_T, spin_number, config_number)
            samples_T = add_histograms(samples_T_temp, samples_T)
        end

        # Mixing the T samples
        if ind_bin == 1
            samples_mixed = hist_countfreq2(spin_samples, spin_number, config_number)
        else
            samples_mixed_temp = hist_countfreq2(spin_samples, spin_number, config_number)
            samples_mixed = add_histograms(samples_mixed_temp, samples_mixed)
        end
    end

    return samples_T, samples_mixed
end

function gibbs_sampling_ising2_binning(gm::FactorGraph{T}, num_samples::Integer, sampling_regime::M_regime) where T <: Real
    @info("using Glauber dynamics v2 to generate M-regime samples with binning")

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Not creating dictionary over spin configs as of length config_number

    # Start Gibbs sampling

    # create bins of samples to look at as the overall number of samples is too large
    max_samples_per_bin = 10000000
    min_bins = div(num_samples,max_samples_per_bin)
    n_bins = copy(min_bins)
    rem_samples = rem(num_samples,max_samples_per_bin)
    if rem_samples > 0
        n_bins += 1
    end
    samples_bin_array = [i <= min_bins ? max_samples_per_bin : rem_samples for i=1:n_bins]

    # generate a random state initially
    sigma_conf = rand(0:(config_number - 1))
    sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1
    sigma_new = copy(sigma)

    # Some temporary structures
    samples_T = []
    samples_T_temp = []
    samples_mixed = []
    samples_mixed_temp = []

    # One bin at a time -- make sure to reuse the last computed sigma from prev bin to initialize new bin
    for ind_bin = 1:n_bins
        num_samples_bin = samples_bin_array[ind_bin]

        # Allocation memory for samples
        samples_pairs_T = Array{Int64, 2}(undef,num_samples_bin, 1 + 2*spin_number)

        # Allocate memory for node chosen for updating
        node_selected_samples = Array{Int8,1}(undef,num_samples_bin)

        # Start Gibbs Sampling for Multiple Restarts
        for ind_step = 1:num_samples_bin
            # generate a random state initially
            sigma_conf = rand(0:(config_number - 1))
            sigma = 2*digits(sigma_conf, base=2, pad=spin_number)-1

            # Random generation of site to be changed
            i = rand(1:spin_number)
            #i = ((ind_step - 1) % spin_number) + 1

            # Proposed state - flip ith state
            sigma_old = deepcopy(sigma)
            sigma_new = deepcopy(sigma)
            sigma_new[i] = -1*sigma_new[i]

            # prob_sigma = prob_cond_ising(sigma, i, adjacency_matrix, prior_vector)
            # prob_sigma_new = prob_cond_ising(sigma_new, i, adjacency_matrix, prior_vector)

            # Calculate the probability of accepting new state which is also prob_sigma_new/prob_sigma
            acceptance_prob = exp( -2*sigma[i]*(adj_matrix[i,:]'*sigma + prior_vector[i])[1] )

            # Accept new state or not
            transitionProbability = min(1, acceptance_prob)
            if rand() < transitionProbability
                sigma = copy(sigma_new)
            end

            # Update Samples Array
            samples_pairs_T[ind_step,1] = copy(i)
            samples_pairs_T[ind_step,2:(spin_number+1)] = deepcopy(sigma_old)
            samples_pairs_T[ind_step,(spin_number+2):end] = deepcopy(sigma)
        end

        # Histogramming for the M-regime
        if ind_bin == 1
            samples_T = hist_countfreq_glauber_dynamics(samples_pairs_T, spin_number, config_number)
        else
            samples_T_temp = hist_countfreq_glauber_dynamics(samples_pairs_T, spin_number, config_number)
            samples_T = add_histograms(samples_T_temp, samples_T)
        end

        # Mixing the T samples
        spin_samples = samples_pairs_T[:,(spin_number+2):end]
        if ind_bin == 1
            samples_mixed = hist_countfreq2(spin_samples, spin_number, config_number)
        else
            samples_mixed_temp = hist_countfreq2(spin_samples, spin_number, config_number)
            samples_mixed = add_histograms(samples_mixed_temp, samples_mixed)
        end
    end

    return samples_T, samples_mixed
end

sample(gm::FactorGraph{T}, number_sample::Integer) where T <: Real = sample(gm, number_sample, 1, Gibbs())[1]
sample(gm::FactorGraph{T}, number_sample::Integer, replicates::Integer) where T <: Real = sample(gm, number_sample, replicates, Gibbs())

function sample(gm::FactorGraph{T}, number_sample::Integer, replicates::Integer, sampler::Gibbs) where T <: Real
    if gm.alphabet != :spin
        error("sampling is only supported for spin FactorGraphs, given alphabet $(gm.alphabet)")
    end

    if gm.order <= 2
        samples = sample_generation_ising(gm, number_sample, replicates)
    else
        samples = sample_generation(gm, number_sample, replicates)
    end

    return samples
end

# Function to generate samples according to Glauber dynamics
gibbs_sampling(gm::FactorGraph{T}, number_sample::Integer) where T <: Real = gibbs_sampling(gm, number_sample, S_regime())
gibbs_sampling2(gm::FactorGraph{T}, number_sample::Integer) where T <: Real = gibbs_sampling2(gm, number_sample, S_regime())

function gibbs_sampling(gm::FactorGraph{T}, number_sample::Integer, sampling_regime::SamplingRegime) where T <: Real
    if gm.alphabet != :spin
        error("sampling is only supported for spin FactorGraphs, given alphabet $(gm.alphabet)")
    end

    if gm.order <= 2
        samples = gibbs_sampling_ising(gm, number_sample, sampling_regime)
    else
        error("not supported yet")
    end

    return samples
end

function gibbs_sampling2(gm::FactorGraph{T}, number_sample::Integer, sampling_regime::SamplingRegime) where T <: Real
    max_num_samples = 10000000

    if gm.alphabet != :spin
        error("sampling is only supported for spin FactorGraphs, given alphabet $(gm.alphabet)")
    end

    if gm.order <= 2 && number_sample <= max_num_samples
        samples = gibbs_sampling_ising2(gm, number_sample, sampling_regime)
    elseif gm.order <= 2 && number_sample > max_num_samples
        samples = gibbs_sampling_ising2_binning(gm, number_sample, sampling_regime)
    else
        error("not supported yet")
    end

    return samples
end

# Function to generate samples according to Glauber dynamics and given a query distribution for σ0
function gibbs_sampling_query(gm::FactorGraph{T}, num_samples::Integer, X_σ::Array, q::Array, sampling_regime::SamplingRegime, FLAG_verbose=true) where T <: Real
    #=
    Inputs:
    X_σ is the alphabet of interest (array of integers to be converted to binary strings)
    q is the query distribution
    =#
    if FLAG_verbose
        @info("using Glauber dynamics query sampling v1 to generate M-regime samples")
    end

    # Error messages
    if gm.alphabet != :spin
        error("sampling is only supported for spin FactorGraphs, given alphabet $(gm.alphabet)")
    end

    if gm.order > 2
        error("not supported yet")
    end

    spin_number   = gm.varible_count
    config_number = 2^spin_number

    adjacency_matrix = convert(Array{T,2}, gm)
    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
    adj_matrix = adjacency_matrix - diagm(0 => diag(adjacency_matrix))

    # Allocate memory for matrix of samples
    # Arranged as [node selected, \sigma^{(t)}, \sigma^{(t+1)}]
    samples_pairs_T = Array{Int64, 2}(undef, num_samples, 3)

    spin_configs = []
    for i = 0:(config_number - 1)
        spin_tmp = 2*digits(i, base=2, pad=spin_number).-1
        push!(spin_configs, string(spin_tmp))
    end
    dict_spin_configs = Dict(spin_configs[i+1] => i for i = 0:(config_number-1))

    # Start Gibbs Sampling for Multiple Restarts
    for ind_step = 1:num_samples
        # generate a random state from X_σ according to distrn q
        sigma_conf = StatsBase.sample(X_σ, StatsBase.weights(q))
        sigma = 2*digits(sigma_conf, base=2, pad=spin_number).-1

        # Random generation of site to be changed
        i = rand(1:spin_number)
        #i = ((ind_step - 1) % spin_number) + 1

        # Proposed state
        sigma_new = deepcopy(sigma)
        sigma_new[i] = 1

        # Calculate the probability of accepting new state
        num_prob = exp( 2.0*(dot(adj_matrix[i,:],sigma) + prior_vector[i])[1] )
        denom_prob = 1.0 + num_prob
        acceptance_prob = num_prob / denom_prob

        # Accept new state or not
        r = rand()
        if r < acceptance_prob
            sigma = deepcopy(sigma_new)
        else
            sigma_new[i] = -1
            sigma = deepcopy(sigma_new)
        end

        # Update Samples Array
        samples_pairs_T[ind_step,1] = copy(i)
        samples_pairs_T[ind_step,2] = deepcopy(sigma_conf)
        samples_pairs_T[ind_step,3] = dict_spin_configs[string(sigma)]
    end

    ## M Samples
    samples_T = hist_countfreq(samples_pairs_T, spin_number, config_number)

    ## Mixing the M samples
    spin_samples = samples_pairs_T[:,3]
    raw_binning = countmap(spin_samples)
    samples_mixed = [ vcat(raw_binning[i], 2*digits(i, base=2, pad=spin_number).-1) for i in keys(raw_binning)]
    samples_mixed = hcat(samples_mixed...)'

    return samples_T, samples_mixed
end


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int, m::Int = 1_000_000_007)::Int
    result = ones(Int, n)
    for prime ∈ prime_numbers(n)
        past_super, last_super, new_super = 0, 0, 0
        last_hyper, new_hyper = 0, 0
        inv_mod = invmod(prime - 1, m)
        for i ∈ prime:n
            number, j = i, 0
            while number % prime == 0
                j += 1
                number ÷= prime
            end
            new_super, new_hyper = last_super + j, last_hyper + i * j
            past_super = last_super = last_super + past_super
            result[i - 1] *= ((powermod(prime, last_hyper - last_super + 1, m) - 1) * inv_mod) % m
            result[i - 1] %= m
            last_super, last_hyper = new_super, new_hyper
        end
        result[n] *= ((powermod(prime, new_hyper - new_super - past_super + 1, m) - 1) * inv_mod) % m
        result[n] %= m
    end
    return sum(result) % m
end

compute(5)

compute(10)

compute(100)

compute(20_000)

@benchmark compute(20_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Test
# using RCall

# @testset "DataConvenience.jl" begin
#     for i in 1:100
#         # Write your own tests here.
#         x = rand(100, 5)
#         y = rand(100, 5)
#
#         @rput x
#         @rput y
#         R"""
#         res = cancor(x,y)$cor[1]
#         """
#         @rget res
#         @test res ≈  canonicalcor(x,y)
#     end
# end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
using GeoInterfaceRecipes
using GeoInterface
using Plots
using Test

abstract type MyAbstractGeom{N} end
# Implement interface
struct MyPoint{N} <: MyAbstractGeom{N} end
struct MyCurve{N} <: MyAbstractGeom{N} end
struct MyPolygon{N} <: MyAbstractGeom{N} end
struct MyMultiPoint{N} <: MyAbstractGeom{N} end
struct MyMultiCurve{N} <: MyAbstractGeom{N} end
struct MyMultiPolygon{N} <: MyAbstractGeom{N} end
struct MyCollection{N} <: MyAbstractGeom{N} end

GeoInterfaceRecipes.@enable_geo_plots MyAbstractGeom

GeoInterface.isgeometry(::MyAbstractGeom) = true
GeoInterface.is3d(::GeoInterface.AbstractGeometryTrait, ::MyAbstractGeom{N}) where {N} = N == 3
GeoInterface.ncoord(::GeoInterface.AbstractGeometryTrait, geom::MyAbstractGeom{N}) where {N} = N
GeoInterface.coordnames(::GeoInterface.AbstractGeometryTrait, ::MyAbstractGeom{2}) = (:X, :Y)
GeoInterface.coordnames(::GeoInterface.AbstractGeometryTrait, ::MyAbstractGeom{3}) = (:X, :Y, :Z)

GeoInterface.geomtrait(::MyPoint) = GeoInterface.PointTrait()
GeoInterface.getcoord(::GeoInterface.PointTrait, geom::MyPoint{2}, i::Integer) = (rand(1:10), rand(11:20))[i]
GeoInterface.getcoord(::GeoInterface.PointTrait, geom::MyPoint{3}, i::Integer) = (rand(1:10), rand(11:20), rand(21:30))[i]

GeoInterface.geomtrait(::MyCurve) = GeoInterface.LineStringTrait()
GeoInterface.ngeom(::GeoInterface.LineStringTrait, geom::MyCurve) = 3
GeoInterface.getgeom(::GeoInterface.LineStringTrait, geom::MyCurve{N}, i) where {N} = MyPoint{N}()
GeoInterface.convert(::Type{MyCurve}, ::GeoInterface.LineStringTrait, geom) = geom

GeoInterface.geomtrait(::MyPolygon) = GeoInterface.PolygonTrait()
GeoInterface.ngeom(::GeoInterface.PolygonTrait, geom::MyPolygon) = 2
GeoInterface.getgeom(::GeoInterface.PolygonTrait, geom::MyPolygon{N}, i) where {N} = MyCurve{N}()

GeoInterface.geomtrait(::MyMultiPolygon) = GeoInterface.MultiPolygonTrait()
GeoInterface.ngeom(::GeoInterface.MultiPolygonTrait, geom::MyMultiPolygon) = 2
GeoInterface.getgeom(::GeoInterface.MultiPolygonTrait, geom::MyMultiPolygon{N}, i) where {N} = MyPolygon{N}()

GeoInterface.geomtrait(::MyMultiPoint) = GeoInterface.MultiPointTrait()
GeoInterface.ngeom(::GeoInterface.MultiPointTrait, geom::MyMultiPoint) = 10
GeoInterface.getgeom(::GeoInterface.MultiPointTrait, geom::MyMultiPoint{N}, i) where {N} = MyPoint{N}()

GeoInterface.geomtrait(geom::MyCollection) = GeoInterface.GeometryCollectionTrait()
GeoInterface.ncoord(::GeoInterface.GeometryCollectionTrait, geom::MyCollection{N}) where {N} = N
GeoInterface.ngeom(::GeoInterface.GeometryCollectionTrait, geom::MyCollection) = 4
GeoInterface.getgeom(::GeoInterface.GeometryCollectionTrait, geom::MyCollection{N}, i) where {N} = MyMultiPolygon{N}()

@testset "plot" begin
    # We just check if they actually run
    # 2d
    plot(MyPoint{2}())
    plot(MyCurve{2}())
    plot(MyMultiPoint{2}())
    plot(MyPolygon{2}())
    plot(MyMultiPolygon{2}())
    plot(MyCollection{2}())
    # 3d
    plot(MyPoint{3}())
    plot(MyCurve{3}())
    plot(MyMultiPoint{3}())
    plot(MyPolygon{3}())
    plot(MyMultiPolygon{3}())
    plot(MyCollection{3}())
end


text:
using Zygote
using FiniteDifferences
using KernelFunctions
using Distances
using Test

x = rand(5)
y = rand(5)

@testset "Testing Zygote adjoints" begin
gzeucl = first(Zygote.gradient(xy->evaluate(Euclidean(),xy[1],xy[2]),[x,y]))
gzsqeucl =  first(Zygote.gradient(xy->evaluate(SqEuclidean(),xy[1],xy[2]),[x,y]))
gzdotprod = first(Zygote.gradient(xy->evaluate(KernelFunctions.DotProduct(),xy[1],xy[2]),[x,y]))

FDM = central_fdm(5,1)

gfeucl = collect(first(FiniteDifferences.grad(FDM,xy->evaluate(Euclidean(),xy[1],xy[2]),(x,y))))
gfsqeucl = collect(first(FiniteDifferences.grad(FDM,xy->evaluate(SqEuclidean(),xy[1],xy[2]),(x,y))))
gfdotprod =collect(first(FiniteDifferences.grad(FDM,xy->evaluate(KernelFunctions.DotProduct(),xy[1],xy[2]),(x,y))))

@test all(gzeucl .≈ gfeucl)
@test all(gzsqeucl .≈ gfsqeucl)
@test all(gzdotprod .≈ gfdotprod)

end


text:
function equiv(x::BoolVec, y::BoolVec)::Bool
    if length(exprs1) == length(exprs2)
        for e1, e2 in zip(exprs1, exprs2)
            if !equiv(e1, e2)
                return false
            end
        end
        true
    else
        false
    end
end

equiv(x::True, y::True) = true
equiv(x::False, y::False) = true
equiv(x::Var, y::Var) = x.name == y.name
equiv(x::And, y::And) = equiv(x.exprs, y.exprs)
equiv(x::Or, y::Or)   = equiv(x.exprs, y.exprs)
equiv(x::Not, y::Not)   = equiv(x.expr, y.expr)
equiv(x::BoolExpr, y::BoolExpr) = false


text:
const ScreenID = UInt8
const ZIndex = Int
# ID, Area, clear, is visible, background color
const ScreenArea = Tuple{ScreenID, Node{IRect2D}, Node{Bool}, Node{Bool}, Node{RGBAf0}}

abstract type GLScreen <: AbstractScreen end

mutable struct Screen <: GLScreen
    glscreen::GLFW.Window
    framebuffer::GLFramebuffer
    rendertask::RefValue{Task}
    screen2scene::Dict{WeakRef, ScreenID}
    screens::Vector{ScreenArea}
    renderlist::Vector{Tuple{ZIndex, ScreenID, RenderObject}}
    cache::Dict{UInt64, RenderObject}
    cache2plot::Dict{UInt16, AbstractPlot}
    framecache::Tuple{Matrix{RGB{N0f8}}, Matrix{RGB{N0f8}}}
    function Screen(
            glscreen::GLFW.Window,
            framebuffer::GLFramebuffer,
            rendertask::RefValue{Task},
            screen2scene::Dict{WeakRef, ScreenID},
            screens::Vector{ScreenArea},
            renderlist::Vector{Tuple{ZIndex, ScreenID, RenderObject}},
            cache::Dict{UInt64, RenderObject},
            cache2plot::Dict{UInt16, AbstractPlot},
        )
        s = size(framebuffer)
        obj = new(
            glscreen, framebuffer, rendertask, screen2scene,
            screens, renderlist, cache, cache2plot,
            (Matrix{RGB{N0f8}}(undef, s), Matrix{RGB{N0f8}}(undef, reverse(s)))
        )
        finalizer(obj) do obj
            # save_print("Freeing screen")
            empty!.((obj.renderlist, obj.screens, obj.cache, obj.screen2scene, obj.cache2plot))
            return
        end
        obj
    end
end
GeometryTypes.widths(x::Screen) = size(x.framebuffer.color)

Base.wait(x::Screen) = isassigned(x.rendertask) && wait(x.rendertask[])
Base.wait(scene::Scene) = wait(global_gl_screen()) # TODO per scene screen

function insertplots!(screen::GLScreen, scene::Scene)
    for elem in scene.plots
        insert!(screen, scene, elem)
    end
    foreach(s-> insertplots!(screen, s), scene.children)
end

function Base.empty!(screen::GLScreen)
    empty!(screen.renderlist)
    empty!(screen.screen2scene)
    empty!(screen.screens)
end

function destroy!(screen::Screen)
    empty!(screen)
    empty!(screen.cache)
    empty!(screen.cache2plot)
    destroy!(screen.glscreen)
end

function Base.resize!(window::GLFW.Window, resolution...)
    if isopen(window)
        retina_scale = retina_scaling_factor(window)
        w, h = resolution ./ retina_scale
        GLFW.SetWindowSize(window, round(Int, w), round(Int, h))
    end
end

function Base.resize!(screen::Screen, w, h)
    nw = to_native(screen)
    resize!(nw, w, h)
    fb = screen.framebuffer
    resize!(fb, (w, h))
end

function Base.display(screen::Screen, scene::Scene)
    empty!(screen)
    resize!(screen, size(scene)...)
    GLFW.PollEvents() # let the size change go through (TODO is this necessary?)
    register_callbacks(scene, to_native(screen))
    insertplots!(screen, scene)
    AbstractPlotting.update!(scene)
    return
end

function to_jl_layout!(A, B)
    ind1, ind2 = axes(A)
    n = first(ind2) + last(ind2)
    for i in ind1
        @simd for j in ind2
            @inbounds B[n-j, i] = ImageCore.clamp01nan(A[i, j])
        end
    end
    return B
end

function fast_color_data!(dest::Array{RGB{N0f8}, 2}, source::Texture{T, 2}) where T
    GLAbstraction.bind(source)
    glPixelStorei(GL_PACK_ALIGNMENT, 1)
    glGetTexImage(source.texturetype, 0, GL_RGB, GL_UNSIGNED_BYTE, dest)
    GLAbstraction.bind(source, 0)
    nothing
end

function AbstractPlotting.colorbuffer(screen::Screen)
    if isopen(screen)
        GLFW.PollEvents()
        render_frame(screen) # let it render
        GLFW.SwapBuffers(to_native(screen))
        glFinish() # block until opengl is done rendering
        ctex = screen.framebuffer.color
        if size(ctex) != size(screen.framecache[1])
            s = size(ctex)
            screen.framecache = (Matrix{RGB{N0f8}}(undef, s), Matrix{RGB{N0f8}}(undef, reverse(s)))
        end
        fast_color_data!(screen.framecache[1], ctex)
        to_jl_layout!(screen.framecache...)
        return screen.framecache[2]
    else
        error("Screen not open!")
    end
end

Base.isopen(x::Screen) = isopen(x.glscreen)
function Base.push!(screen::GLScreen, scene::Scene, robj)
    filter!(screen.screen2scene) do (k, v)
        k.value != nothing
    end
    screenid = get!(screen.screen2scene, WeakRef(scene)) do
        id = length(screen.screens) + 1
        bg = lift(to_color, scene.theme[:backgroundcolor])
        clear = lift(identity, scene.theme[:clear])
        visible = lift(identity, scene.theme[:visible])
        push!(screen.screens, (id, scene.px_area, clear, visible, bg))
        id
    end
    push!(screen.renderlist, (0, screenid, robj))
    return robj
end

to_native(x::Screen) = x.glscreen
const gl_screens = GLFW.Window[]

"""
OpenGL shares all data containers between shared contexts, but not vertexarrays -.-
So to share a robjs between a context, we need to rewrap the vertexarray into a new one for that
specific context.
"""
function rewrap(robj::RenderObject{Pre}) where Pre
    RenderObject{Pre}(
        robj.main,
        robj.uniforms,
        GLVertexArray(robj.vertexarray),
        robj.prerenderfunction,
        robj.postrenderfunction,
        robj.boundingbox,
    )
end

const _global_gl_screen = Ref{Screen}()

# will get overloaded later
function renderloop end

# TODO a global is not very nice, but it's the simplest way right now to swap out
# the rendering loop
const opengl_renderloop = Ref{Function}(renderloop)

function Screen(; resolution = (10, 10), visible = false, kw_args...)
    if !isempty(gl_screens)
        for elem in gl_screens
            isopen(elem) && destroy!(elem)
        end
        empty!(gl_screens)
    end

    window = GLFW.Window(
        name = "Makie", resolution = (10, 10), # 10, because smaller sizes seem to error on some platforms
        windowhints = [
            (GLFW.SAMPLES,      0),
            (GLFW.DEPTH_BITS,   0),

            # SETTING THE ALPHA BIT IS REALLY IMPORTANT ON OSX, SINCE IT WILL JUST KEEP SHOWING A BLACK SCREEN
            # WITHOUT ANY ERROR -.-
            (GLFW.ALPHA_BITS,   8),
            (GLFW.RED_BITS,     8),
            (GLFW.GREEN_BITS,   8),
            (GLFW.BLUE_BITS,    8),

            (GLFW.STENCIL_BITS, 0),
            (GLFW.AUX_BUFFERS,  0)
        ],
        visible = false,
        kw_args...
    )
    GLFW.SetWindowIcon(window , AbstractPlotting.icon())

    # tell GLAbstraction that we created a new context.
    # This is important for resource tracking, and only needed for the first context
    GLAbstraction.switch_context!(window)
    GLAbstraction.empty_shader_cache!()
    push!(gl_screens, window)

    GLFW.SwapInterval(0)

    # Retina screens on osx have a different scaling!
    retina_scale = retina_scaling_factor(window)
    resolution = round.(Int, retina_scale .* resolution)
    # Set the resolution for real now!
    GLFW.SetWindowSize(window, resolution...)
    fb = GLFramebuffer(Int.(resolution))

    screen = Screen(
        window, fb,
        RefValue{Task}(),
        Dict{WeakRef, ScreenID}(),
        ScreenArea[],
        Tuple{ZIndex, ScreenID, RenderObject}[],
        Dict{UInt64, RenderObject}(),
        Dict{UInt16, AbstractPlot}(),
    )
    if visible
        GLFW.ShowWindow(window)
    else
        GLFW.HideWindow(window)
    end
    screen.rendertask[] = @async((opengl_renderloop[])(screen))
    screen
end

function global_gl_screen()
    screen = if isassigned(_global_gl_screen) && isopen(_global_gl_screen[])
        _global_gl_screen[]
    else
        _global_gl_screen[] = Screen()
        _global_gl_screen[]
    end
    GLFW.set_visibility!(to_native(screen), AbstractPlotting.use_display[])
    screen
end

# TODO per scene screen
getscreen(scene) = global_gl_screen()

function pick_native(scene::SceneLike, xy::VectorTypes{2}, sid = Base.RefValue{SelectionID{UInt16}}())
    screen = getscreen(scene)
    screen == nothing && return SelectionID{Int}(0, 0)
    window_size = widths(screen)
    fb = screen.framebuffer
    buff = fb.objectid
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1])
    glReadBuffer(GL_COLOR_ATTACHMENT1)
    x, y = Int.(floor.(xy))
    w, h = window_size
    if x > 0 && y > 0 && x <= w && y <= h
        glReadPixels(x, y, 1, 1, buff.format, buff.pixeltype, sid)
        return convert(SelectionID{Int}, sid[])
    end
    return SelectionID{Int}(0, 0)
end

pick(scene::SceneLike, xy...) = pick(scene, Float64.(xy))

function pick(scene::SceneLike, xy::VectorTypes{2})
    sid = pick_native(scene, xy)
    screen = getscreen(scene)
    if screen != nothing && haskey(screen.cache2plot, sid.id)
        plot = screen.cache2plot[sid.id]
        return (plot, sid.index)
    end
    return (nothing, 0)
end

# TODO does this actually needs to be a global?
const _mouse_selection_id = Base.RefValue{SelectionID{UInt16}}()
function mouse_selection_native(scene::SceneLike)
    function query_mouse(buff, w, h)
        xy = events(scene).mouseposition[]
        x, y = Int.(floor.(xy))
        if x > 0 && y > 0 && x <= w && y <= h
            glReadPixels(x, y, 1, 1, buff.format, buff.pixeltype, _mouse_selection_id)
        end
        return
    end
    if !(query_mouse in selection_queries)
        push!(selection_queries, query_mouse)
    end
    convert(SelectionID{Int}, _mouse_selection_id[])
end

function mouse_selection(scene::SceneLike)
    sid = mouse_selection_native(scene)
    screen = getscreen(scene)
    if screen != nothing && haskey(screen.cache2plot, sid.id)
        plot = screen.cache2plot[sid.id]
        return (plot, sid.index)
    end
    return (nothing, 0)
end

function mouseover(scene::SceneLike, plots::AbstractPlot...)
    p, idx = mouse_selection(scene)
    p in flatten_plots(plots)
end

function flatten_plots(x::Atomic, plots = AbstractPlot[])
    push!(plots, x)
    plots
end

function flatten_plots(x::Combined, plots = AbstractPlot[])
    for elem in x.plots
        flatten_plots(elem, plots)
    end
    plots
end

function flatten_plots(array, plots = AbstractPlot[])
    for elem in array
        flatten_plots(elem, plots)
    end
    plots
end

function onpick(f, scene::SceneLike, plots::AbstractPlot...)
    fplots = flatten_plots(plots)
    map_once(events(scene).mouseposition) do mp
        p, idx = mouse_selection(scene)
        (p in fplots) && f(idx)
        return
    end
end

function pick(screen::Screen, rect::IRect2D)
    window_size = widths(screen)
    buff = screen.framebuffer.objectid
    sid = zeros(SelectionID{UInt16}, widths(rect)...)
    glReadBuffer(GL_COLOR_ATTACHMENT1)
    x, y = minimum(rect)
    rw, rh = widths(rect)
    w, h = window_size
    if x > 0 && y > 0 && x <= w && y <= h
        glReadPixels(x, y, rw, rh, buff.format, buff.pixeltype, sid)
        return map(unique(vec(SelectionID{Int}.(sid)))) do sid
            screen.cache2plot[sid.id], Int(sid.index)
        end
    end
    return SelectionID{Int}[]
end


text:
# Various color scales.

# Weighted mean of some number of colors within the same space.
#
# Args:
#  cs: Colors.
#  ws: Weights of the same length as cs.
#
# Returns:
#   A weighted mean color of type T.
#
function weighted_color_mean{S <: Number}(
        cs::AbstractArray{LAB{Float64},1}, ws::AbstractArray{S,1})
    l = 0.0
    a = 0.0
    b = 0.0
    sumws = sum(ws)
    for (c, w) in zip(cs, ws)
        w /= sumws
        l += w * c.l
        a += w * c.a
        b += w * c.b
    end
    LAB(l, a, b)
end

# Discrete scales
# ---------------

# Generate colors in the LCHab (LCHuv, resp.) colorspace by using a fixed
# luminance and chroma, and varying the hue.
#
# Args:
#   l: luminance
#   c: chroma
#   h0: start hue
#   n: number of colors
#
function lab_rainbow(l, c, h0, n)
    ColorValue[LCHab(l, c, h0 + 360.0 * (i - 1) / n) for i in 1:n]
end

function luv_rainbow(l, c, h0, n)
    ColorValue[LCHuv(l, c, h0 + 360.0 * (i - 1) / n) for i in 1:n]
end

# Helpful for Experimenting
function plot_color_scale{T <: ColorValue}(colors::Vector{T})
    println(colors)
    canvas(UnitBox(length(colors), 1)) <<
            (compose([rectangle(i-1, 0, 1, 1) << fill(c)
                      for (i, c) in enumerate(colors)]...) << stroke(nothing))
end

# Continuous scales
# -----------------

# Generate a gradient between n >= 2, colors.

# Then functions return functions suitable for ContinuousColorScales.
function lab_gradient(cs::ColorValue...)
    if length(cs) < 2
        error("Two or more colors are needed for gradients")
    end

    cs_lab = [convert(LAB, c) for c in cs]
    n = length(cs_lab)
    function f(p::Float64)
        @assert 0.0 <= p <= 1.0
        i = 1 + min(n - 2, max(0, int(floor(p*(n-1)))))
        w = p*(n-1) + 1 - i
        weighted_color_mean([cs_lab[i], cs_lab[i+1]], [1.0 - w, w])
    end
    f
end

function lchabmix(c0_, c1_, r, power)
    c0 = convert(LCHab, c0_)
    c1 = convert(LCHab, c1_)
    w = r^power
    lspan = c1.l - c0.l
    cspan = c1.c - c0.c
    hspan = c1.h - c0.h
    return LCHab(c0.l + w*lspan, c0.c + w*cspan, c0.h + w*hspan)
end