text:
function UTCTAIBulletinAFilesLoader(arg0::JString)
    return UTCTAIBulletinAFilesLoader((JString,), arg0)
end

function UTCTAIBulletinAFilesLoader(arg0::JString, arg1::DataProvidersManager)
    return UTCTAIBulletinAFilesLoader((JString, DataProvidersManager), arg0, arg1)
end

function load_offsets(obj::UTCTAIBulletinAFilesLoader)
    return jcall(obj, "loadOffsets", List, ())
end


text:
function intersectingvertices(S1, S2; tol::Float64 = 1/10^10)
    # Dimension
    n = size(S1, 1)
    # Centroid and radii
	 c1, c2 = Circumsphere(S1)[2:n+1], Circumsphere(S2)[2:n+1]
  	 r1, r2 = Circumsphere(S1)[1], Circumsphere(S2)[1]

  	 # Orientation of simplices
  	 orientation_S1 = det([ones(1, n + 1); S1])
  	 orientation_S2 = det([ones(1, n + 1); S2])

  	 if abs(orientation_S1) < tol || abs(orientation_S2) < tol
        return Array{Float64, 2}(undef, 0, 0)
  	 end

  	# Set volume to zero initially. Change only if there is intersection
  	 IntVol = 0.0

    # -------------------------------------
    # Simplices intersect in some way
    # -------------------------------------

    # If the (distance between centroids)^2-(sum of radii)^2 < 0,
    # then the simplices intersect in some way.
    dist_difference::Float64 = (transpose(c1 - c2) * (c1 - c2) - (r1 + r2)^2)[1]

    if dist_difference < 0
        #@show "The simplices intersect in some way (possibly only along a boundary, with zero volume intersection)"
        # Find the number of points of each simplex contained within the
        # circumsphere of the other simplex
        vertices1InCircum2 = SomeVertexInCircumsphere(S1, r2, c2)
        vertices2InCircum1 = SomeVertexInCircumsphere(S2, r1, c1)

        if vertices1InCircum2 + vertices2InCircum1 >= 1
            βs1in2, βs2in1, ordered_vertices1, ordered_vertices2, numof1in2, numof2in1 =
            BarycentricCoordinates(S1,S2,orientation_S1,orientation_S2,tol)
            # Trivial intersections
            TriviallyContained = heaviside0([numof1in2 numof2in1] .- (n+1))
            IsSomeContained = sum(TriviallyContained, dims=2)[1]

            if IsSomeContained == 2.0 # The simplices coincide
                #@show "The simplices coincide"
                return copy(transpose(S1))
            elseif IsSomeContained == 1.0 # One simplex is contained in the other
                if TriviallyContained[1] == 1.0 # Simplex1 is contained in Simplex2
                    #@show "Simplex 1 is contained in simplex 2"
                    return copy(transpose(S1))
                else # Simplex2 is contained in Simplex1
                    #@show "Simplex 2 is contained in simplex 1"
                    return copy(transpose(S2))
                end
            else # No simplex contains the other
                #@show "No simplex contains the other"
                Ncomm, ordered_vertices1, ordered_vertices2 = SharedVertices(βs1in2,ordered_vertices1,ordered_vertices2,numof1in2,numof2in1)

                # Is there any shared face?
                if Ncomm == n
                    #@show "The simplices share a face"
                    IntVol = SharedFaceVolume(S2, βs1in2, ordered_vertices1, ordered_vertices2)
                    return transpose(SharedFaceVertices(S2, βs1in2, ordered_vertices1, ordered_vertices2))
                else # The simplices do not share a face.
                    #@show "The simplices do not share a face"
                    IntVert, ConvexExpIntVert = IntersectionOfBoundaries_NoStorage(S1,S2,βs1in2,βs2in1, ordered_vertices1, ordered_vertices2, numof1in2, numof2in1, Ncomm, tol)

                    if !isempty(IntVert)
                        IntVert,ConvexExpIntVert = PolytopeGeneratingVertices(S1,S2,IntVert,ConvexExpIntVert,βs1in2,βs2in1,ordered_vertices1,ordered_vertices2,numof1in2,numof2in1,Ncomm);
                        IntVol = VolumeComputation(IntVert, ConvexExpIntVert)

                        return IntVert

                    else
                        return Array{Float64, 2}(undef, 0, 0)
                    end
                end
            end
        else
            IntVert = Array{Float64, 2}(undef, 0, 0)
        end
    else
        IntVert = Array{Float64, 2}(undef, 0, 0)
    end

    return IntVert
end

ext: jl


text:
using CodecZlib: GzipCompressorStream
using Tar: Tar

const bindings_filename = "libp4est.jl"
const tarball_template = "P4est.vVERSION.tar.gz"

function create_bindings_artifact(command_line_arguments)
  # Iterate over all arguments to check if help is needed
  for arg in command_line_arguments
    if arg in ("-h", "-help", "--help")
      println("Usage: julia --project create_bindings_artifact.jl [-h|--help] BINDINGS_FILE VERSION")
      return
    end
  end

  # Abort if wrong number of arguments are given
  if length(command_line_arguments) != 2
    error("wrong number of command line arguments (use `-h` for more information)")
  end

  # Unpack arguments
  filepath, version_string = command_line_arguments

  # Check if input file exists
  if !isfile(filepath)
    error("bindings file `$filepath` does not exist")
  end

  # Parse version string
  version = VersionNumber(version_string)

  # Verify that it only contains MAJOR.MINOR.PATCH
  if version != VersionNumber(version.major, version.minor, version.patch)
    error("version number should only contain MAJOR.MINOR.PATCH, but found $version")
  end

  # OK, we got here, so everything *seems* to be OK to start creating the artifact archive

  # Determine canonical tarball filename
  tarball = replace(tarball_template, "VERSION" => string(version))

  # Create temporary directory, copy bindings file there, and create tarball
  mktempdir() do tmpdir
    cp(filepath, joinpath(tmpdir, bindings_filename))
    open(tarball, write=true) do tar_gz
      tar = GzipCompressorStream(tar_gz)
      Tar.create(tmpdir, tar)
      close(tar)
    end
  end

  # Print name of created file
  println(tarball)
end

create_bindings_artifact(ARGS)


text:
function repl_render(cur_player::AbstractPlayer)
  cur_string = []

  push!(cur_string, "\x1b[2J")

  push!(cur_string, crayon_dict["grey"])

  main_chunk = "\u258C" * " " ^ 20 * "\u2590"

  plain_row = " " ^ 29 * main_chunk * "\n"

  card_top_chunk = " " ^ 9 * "\u2584" ^ 10 * " " ^ 10

  card_top = card_top_chunk * main_chunk * reverse(card_top_chunk) * "\n"

  card_mid_chunk = " " ^ 9 * "\u258C" * " " ^ 8 * "\u2590" * " " ^ 9

  card_mid = card_mid_chunk * " " * main_chunk * " " * card_mid_chunk * "\n"

  card_bot_chunk = " " ^ 9 * "\u2580" ^ 10 * " " ^ 10

  card_bot = card_bot_chunk * main_chunk * reverse(card_bot_chunk) * "\n"

  bottom_row = " " ^ 29 * "\u2580" ^ 22 * "\n"

  push!(cur_string, plain_row)
  push!(cur_string, plain_row)

  for cur_index in 1:3
    push!(cur_string, card_top)
    push!(cur_string, card_mid)
    push!(cur_string, card_mid)
    push!(cur_string, card_bot)
    push!(cur_string, plain_row)
    push!(cur_string, plain_row)
  end

  push!(cur_string, bottom_row)

  push!(cur_string, "\x1b[s")

  push!(cur_string, "\x1b[7;11H")
  push!(cur_string, " level: ")

  push!(cur_string, "\x1b[8;11H")
  push!(cur_string, "   $(lpad(cur_player.level, 2, "0"))   ")

  push!(cur_string, "\x1b[13;11H")
  push!(cur_string, " score: ")

  push!(cur_string, "\x1b[14;11H")
  push!(cur_string, lpad(cur_player.score, 8, "0"))

  push!(cur_string, inv(crayon_dict["grey"]))

  push!(cur_string, "\x1b[u")
  print(cur_string...)

  raw_terminal() do
    while true
      tmp_input = readavailable(terminal.in_stream)

      if !cur_player.game.in_focus
        play(cur_player)
        continue
      end

      cur_player.state.is_playing || break

      if ( tmp_input == UInt8[0x1b] )
        cur_player.state.is_playing = false
        continue
      end

      ( tmp_input == UInt8[0x1b, 0x5b, 0x44] ) && slide_left(cur_player)
      ( tmp_input == UInt8[0x1b, 0x5b, 0x43] ) && slide_right(cur_player)
      ( tmp_input == UInt8[0x1b, 0x5b, 0x42] ) && soft_drop(cur_player)
      ( tmp_input == UInt8[0x1b, 0x5b, 0x41] ) && clockwise(cur_player)

      tmp_input = String(tmp_input)

      ( tmp_input == "c" ) && hold(cur_player)
      ( tmp_input == "z" ) && counter(cur_player)
      ( tmp_input == "p" ) && pause(cur_player)
      ( tmp_input == " " ) && hard_drop(cur_player)
    end
  end
end


text:
function draw(l::Line)
    x = collect(range(-1, stop=1, length=100))
    plot!(x, point_in_line(l, x))
end

function draw_lines(H...)
    r = 1.1
    c = palette(:Set1_8)

    x = collect(range(-r, stop=r, length=2))
    v = collect(range(0, stop=2*pi, length=100))
    dessin = plot(sin.(v), cos.(v), color=:gray) #plot the circle
    for (i, arg) in enumerate(H)
        for l in arg
            plot!(dessin, x, point_in_line(l, x), color=:thistle)
        end
    end
    plot!(dessin, xlims=(-r, r), ylims=(-r, r), legend=:none, aspect_ratio=:equal)
    return dessin
end

max_stars_repo_name: SimonCoste/LineIntersections

max_issues_repo_name: SimonCoste/LineIntersections

max_forks_repo_name: SimonCoste/LineIntersections


text:
@testset "uri_utils" begin
    @test string(uri_add_path("http://foo/bar", "baz", "bam")) ==
        "http://foo/bar/baz/bam"
end


text:
module GMM

#=
Gaussian Mixture Models (GMMs) in Julia
Written by: Jasper Ginn
Date: 12/6/20
=#

using Random, Plots, Logging, LinearAlgebra, ProgressMeter, Distributions
#logger = global_logger(SimpleLogger(stdout, Logging.Error)) # Change to Logging.Debug for detailed info

#= Expectation-maximization algorithm =#

#=
E-step (computing the probabilities of a point belonging to each cluster)
=#

function E_step(X::Array{Float64}, π::Array{Float64}, μ::Array{Float64}, Σ::Array{Float64})::Array{Float64}
    #=
    Perform the E-step of the expectation-maximization (EM) algorithm for a GMM
    :param X: matrix of dimensions (n x m) that contains the input data
    :param π: array of dimensions (k x 1) that contains the mixing proportions
    :param μ: array of dimensions (k x m) that contains the mixture component means
    :param Σ: tensor of dimensions (k x m x m) that contains the mixture component covariance matrices
    :return: matrix Γ of dimensions (n x k) that contains the probabilities of clusters for the objects
    :details: During the E-step, we solve for the posterior distribution of p(t_i =k | x_i, θ).
        Here, θ is a vector containing all parameters of the GMM.

        1. Given some variational distribution q(t) and parameter vector θ, we can define

           L(θ, q) = KL( q(t) || p(t | X, θ) ). The variational distribution that minimizes
           the KL distance is p(t | X, θ) = Γ.

        2. We can compute gamma as:

            (1) Γ = p(t = k | X, θ) = p(X | t = k) p(t = k) / [Σ p(X | t = k) p(t = k)]

           Where:

            (2) p(X | t = k) = N(X | μ_k, Σ_k)
            (3) p(t = k) = π_k

           Hence, we have:

            (4) Γ = [π_k * N(X | μ_k, Σ_k)] / [∑^K π_j * N(X | μ_j, Σ_j)]
    =#
    # Number of subjects / variables
    N, M = size(X)
    # Number of clusters
    K = size(π)[1]
    @debug "Dimensions are N=$N, M=$M, K=$K ..."
    # Init gamma
    Γ = zeros((N, K))
    # For each cluster, compute the E-step
    for k ∈ 1:K
        # Get mixing proportion, μ_k, Σ_k
        π_k = π[k]
        μ_k = μ[k, :]
        Σ_k = Σ[:, :, k]
        # Emit log
        @debug "Performing E-step for cluster $k with following parameters ..." π_k μ_k Σ_k
        # Draw pdf from multivariate normal
        # Multiply the likelihood times the prior
        posterior = MvNormal(μ_k, Σ_k |>
                                    Hermitian |>
                                    Matrix) |>
                    x -> pdf(x, X') |>
                    x -> π_k .* x
        # Add to gamma
        Γ[:, k] = posterior
    end;
    # Normalize to make proper distribution
    Γ ./= sum(Γ, dims=2)
    # Return
    return Γ
end;

#=
M-step (computing the values of the Gaussian parameters)
=#

function M_step(X::Array{Float64}, Γ::Array{Float64})::Tuple{Array{Float64}, Array{Float64}, Array{Float64}}
    #=
    Perform the M-step of the EM algorithm
    :param X: matrix of dimensions (n x m) that contains the input data
    :param Γ: matrix Γ of dimensions (n x k) that contains the probabilities of clusters for the objects
    :return: a 3-tuple containing (I) (k x m) mean vector μ, (II) (k x 1) mixing proportions vector π, and
             (III) (k x m x m) covariance tensor Σ
    :details: We want to maximize π_k, μ_k and Σ_k given the formula:

        (1) Q(X, θ) = ∑^N ∑^K p(t=k | X, θ) ln[p(x | t=k, θ) * p(t=k)] =
                      ∑^N ∑^K γ_{ik} ln[N(X; μ_k, Σ_k)] + ln[π_k]

    This optimization problem uses the probabilities computed during the E-step (γ_nk)
    and is subject to the following constraints:

        (a) ∑^K π_k = 1 --> g(π_k) = ∑^K π_k - 1
        (b) π_k ≥ 0 ∀π_k

    Constraint (b) can be ignored in this problem for the sake of simplicity, which leaves:

        (2) Z(X, θ) = Q(X, θ) - λ[g(π_k)] =
            ∑^N ∑^K {γ_{ik} (ln[N(X; μ_k, Σ_k)] + ln[π_k])} - λ[∑^K π_k - 1]

    And we want:

        (3) arg max_{π_k, μ_k, Σ_k} Z(X, θ)

    - Mixing components

    To obtain the derivative of the mixing proportions π_k, we solve:

        (4) ∂Z/∂π_k ∑^N {γ_nk * ln[π_k]} - λ(π_k - 1) =
            ∑^N ∂Z/∂π_k {γ_nk * ln[π_k]} - ∂Z/∂π_k λ(π_k - 1) =
            ∑^N γ_nk / π_k - λ

    We set this to 0 and solve for π_k

        (5) ∑^N γ_nk / π_k - λ = 0 -->
            1/π_k ∑^N γ_nk = λ -->
            π_k = ∑^N γ_nk / λ

    Notice that, if we multiply both sides with ∑^K, we have:

        (c) ∑^K π_k = 1
        (d) ∑^N ∑^k γ_nk = N (because sum across probabilities for each person equals 1,
                              and summing this across subjects yields N)
        (c) 1 = N / λ --> λ = N

    Hence, if we substitute these results in (5), we get:

        (6) π_k = ∑^N γ_nk / N

    - Means

    The means for each cluster are given by:

        μ_k = ∑^N γ_nk x_n / ∑^N γ_nk

    - Covariances

    The covariances for each cluster are given by:

        Σ_k = ∑^N γ_nk (x_n - μ_k) (x_n - μ_k)' / ∑^N γ_nk
    =#
    # Record dimensions
    N, M = size(X)
    _, K = size(Γ)
    @debug "Dimensions are N=$N, M=$M, K=$K ..."
    # Allocate memory
    μ = zeros((K, M))
    Σ = zeros((M, M, K))
    μ_k = zeros((1, M))
    Σ_k = zeros((M, M))
    # Compute the column sum of Γ
    Γ_colsum = sum(Γ, dims=1)
    # Compute new mixing probabilities
    π = Γ_colsum / N
    # For each cluster, compute values
    for k ∈ 1:K
        @debug "Performing M-step for cluster $k ..."
        # Compute means
        μ_k = sum(Γ[:,k] .* X, dims=1) ./ Γ_colsum[k]
        @debug "Computed M-step for μ ..."
        # Compute sigma
        X_centered = X .- μ_k
        Σ_k = (Γ[:,k] .* X_centered)' * X_centered ./ Γ_colsum[k]
        @debug "Computed M-step for Σ ..."
        # Emit
        @debug "Obtained the following parameter estimates on the M-step for cluster $k ..." round.(π[k], digits=2) round.(μ_k, digits=2) round.(Σ_k, digits=2)
        # Store
        μ[k,:] = μ_k
        Σ[:,:,k] = Σ_k
    end;
    # Return
    return π', μ, Σ
end;

#=
Loss function
=#

function L(X::Array{Float64}, π::Array{Float64}, μ::Array{Float64}, Σ::Array{Float64}, Γ::Array{Float64})::Float64
    #=
    Compute the loss value given the current parameters
    :param X: matrix of dimensions (n x m) that contains the input data
    :param π: array of dimensions (k x 1) that contains the mixing proportions
    :param μ: array of dimensions (k x m) that contains the mixture component means
    :param Σ: tensor of dimensions (k x m x m) that contains the mixture component covariance matrices
    :param Γ: matrix of dimensions (n x k) that contains the probabilities of clusters for the objects
    :return: scalar, Loss value
    :details: Compute the loss for the likelihood of the data

        (1) log[p(X)] = ∑^n log[∑^k π_k * N(x_n | μ_k, Σ_k)]

    =#
    # Get dimensions
    N, M = size(X)
    _, K = size(Γ)
    @debug "Dimensions are N=$N, M=$M, K=$K ..."
    # Allocate
    LL = zeros(N)
    # For each cluster
    for k ∈ 1:K
        @debug "Computing log-likelihood for cluster $k ..."
        # Retrieve params for cluster k
        μ_k = μ[k, :]
        Σ_k = Σ[:, :, k]
        π_k = π[k]
        # Compute likelihood
        LL_k = MvNormal(μ_k, Σ_k |>
                            Hermitian |>
                            Matrix) |>
                x -> pdf(x, X') |>
                x -> π_k .* x
        # Add to running total
        LL += LL_k
    end;
    # Apply logarithm and sum
    LL = sum(log.(LL))
    @debug "Log-likelihood is $(round(LL, digits=3)) ..."
    # Return
    return LL
end;

#=
Initialization of parameters
=#

function init(M::Int64, K::Int64)::Tuple{Array{Float64}, Array{Float64}, Array{Float64}}
    #=
    Initialize parameter vectors/matrices π_k, μ_k and Σ_k
    :param N: Number of subjects
    :param M: Number of variables
    :param K: Number of clusters
    :return: tuple containing π0, μ0 and Σ0
    =#
    # Allocate
    π0 = zeros((K,))
    μ0 = rand(Float64, (K, M))
    Σ0 = zeros((M, M, K))
    @debug "Randomly initialized μ0 for all clusters..." μ0
    # For each cluster, init parameters
    for k ∈ 1:K
        # Draw π_k:
            # 1) from uniform distribution with a = 0 and b = sum(π0) if k < K
            # 2) As the remainder of 1 - sum(π)
        π0[k] = k < K ? Uniform(0, 1 - sum(π0)) |> rand : 1 - sum(π0)
        # Initialize Sigma
        # NB: must be symmetric + positive semidefinite
        Σ0[:, :, k] = Matrix(I, (M, M))
        # Assert that covariance matrix is symmetric
        @debug "Initialized π0 and Σ0 on cluster $k ..." π0[k] Σ0[:, :, k]
        @assert Σ0[:, :, k] == Σ0[:, :, k]' "Initial covariance matrix on cluster $k is not symmetric ..."
        # Assert that covariance matrix is positive semidefinite
        @assert (eigen(Σ0[:, :, k]).values .> 0) |> all "Initial covariance matrix on cluster $k not positive semidefinite ..."
    end;
    @debug "π0 sums to $(sum(π0)) ..."
    # Assert that π0 sums to 1
    @assert sum(π0) == 1 "Initial mixing proportions do not sum exactly to 1 ..."
    # Return
    return π0, μ0, Σ0
end;

#=
Training loop for GMM using EM algorithm
=#

function train_GMM(X::Array{Float64}, K::Int64, ϵ::Float64=1e-3, maxiter::Int64=100, epochs::Int64=10)::Tuple{Array{Float64}, Array{Float64}, Array{Float64}, Float64, Array{Float64}}
    #=
    Implement the training loop for the GMM using EM
    :param X: matrix of dimensions (n x m) that contains the input data
    :param K: number of clusters
    :param ϵ: delta value between current and last loss. If delta > ϵ the program will exit
    :param maxiter: maximum number of iterations in an epoch
    :param epochs: number of times we run the EM algorithm using different starting values
    :return: tuple containing:
        1. π_best (k x 1)
        2. μ_best (k x m)
        3. Σ_best (m x m x k)
        4. loss_best
        5. loss_history (maxiter x epochs)
    =#
    # Get dimensions
    N, M = size(X)
    # Allocate
    π_best = zeros((K,))
    μ_best = zeros((K, M))
    Σ_best = zeros((K, M, M))
    Γ = zeros((N, K))
    Σ = zeros((M, M, K))
    μ = zeros((K, 1))
    loss_δ = 0.0
    # Best los
    loss_best = NaN
    # Loss history
    LH = zeros((maxiter, epochs))
    # Run epochs
    @showprogress "Fitting GMM ..." for epoch ∈ 1:epochs
        @debug "Started EM algorithm for epoch $epoch ..."
        loss_previous = NaN
        loss_current = 0.0
        # Initialize means etc.
        π, μ, Σ = init(M, K)
        # Run iterations
        for iter ∈ 1:maxiter
            @debug "Started iteration $iter on epoch $epoch ..."
            # Perform E-step
            try
                Γ = E_step(X, π, μ, Σ)
            catch e
                @error "Encountered $e error while computing E-step ..."
                break
            end;
            @debug "Finished E-step ..."
            # Perform M-step
            π, μ, Σ = M_step(X, Γ)
            @debug "Finished M-step ..."
            # Compute loss
            try
                loss_current = L(X, π, μ, Σ, Γ)
            catch e
                @error "Encountered $e error while computing loss ..."
                break
            end;
            @debug "Computed log-likelihood ..."
            # If best loss NaN
            if isnan(loss_best)
                LH[iter, epoch] = loss_best = loss_previous = loss_current
                π_best, μ_best, Σ_best = π, μ, Σ
                continue
            end;
            # If previous loss NaN
            if isnan(loss_previous)
                LH[iter, epoch] = loss_previous = loss_current
                continue
            end;
            # Assert that current loss is not smaller than previous loss
            @assert loss_current > loss_previous "Loss is decreasing. There must be a bug in your code somewhere \\_(:/)_/ ..."
            # Store parameters if loss current beats loss best
            if loss_current > loss_best
                loss_best = loss_current
                π_best, μ_best, Σ_best = π, μ, Σ
            end;
            # Check delta in losses
            loss_δ = -1 * ((loss_current - loss_previous) / loss_previous)
            if loss_δ ≤ ϵ
                @debug "Stopping condition reached. Exiting epoch $epoch ..."
                break
            end;
            # Set current loss as previous loss and add to history
            LH[iter, epoch] = loss_previous = loss_current
            @debug "Finished iteration $iter on epoch $epoch ..."
        end;
        @debug "Finished epoch $epoch ..."
    end;
    @debug "Finished program. Exiting now ..."
    # Return best values
    return π_best, μ_best, Σ_best, loss_best, LH
end;

#=
Wrapper functions / utility functions
=#

"""
    clust(X, K, ϵ=1.e-3, maxiter=100, epochs=10)::Tuple{Array{Float64}, Array{Float64}}

Fit a Gaussian Mixture Model (GMM) on the input data X using EM estimation.

# Arguments
- `X::Array{Float64}`:
- `K::Int64`: Number of clusters
- `ϵ::Float64=1e-3`: Value that governs the stopping condition. If the delta loss is smaller than this value, then quit the program.
- `maxiter::Int64=100`: Maximum number of iterations the algorithm will execute in each epoch.
- `epochs::Int64=10`: Number of epochs the algorithm will run.

# Returns
A tuple containing
- Cluster labels. This is an array of dimensions (N x 1)
- Loss history for each epoch. This is an array of dimensions (maxiter x epochs)

# Examples
```julia-repl

```
"""
function clust(X::Array{Float64}, K::Int64; ϵ::Float64=1e-3, maxiter::Int64=100, epochs::Int64=10)::Tuple{Tuple{Array{Float64}, Array{Float64}, Array{Float64}, Float64},Array{Float64}, Array{Float64}}
    # Record dimensions
    N, M = size(X)
    # Run algorithm
    π1, μ1, Σ1, L1, history = train_GMM(X, K, ϵ, maxiter, epochs)
    # Compute labels using best values
    lbls = [i for i in mapslices(argmax, E_step(X, π1, μ1, Σ1), dims=2)[:]]
    # Return history and labels
    return (π1, μ1, Σ1, L1), lbls, history
end;

function plot_history(history::Array{Float64})
    # Plot loss
    maxiter, epochs = size(history)
    hc = history[:,1];
    plot(hc[hc .< 0], legend = false)
    # For each epoch, plot
    for epoch ∈ 1:(size(history)[2]-1)
        hc = history[:,epoch]
        plot!(hc[hc .< 0])
    end;
    hc = history[:,end];
    p = plot!(hc[hc .< 0], title = "Change in loss for $epochs epochs")
    return p
end;

"""
    simulate_gmm(K::Int64, N::Int64, μ::Array{Float64}, Σ::Array{Float64})::Tuple{Array{Float64}, Array{Float64}}

Simulate a dataset that is a mixture of K multivariate Gaussians

# Arguments
- `K::Int64`: Number of clusters
- `N::Array{Int64}`: Array of dimensions (1 x K) containing the number of data points for each cluster.
- `μ::Array{Float64}`: Array of dimensions (K x M) containing the M means for each cluster.
- `Σ::Array{Float64}`: Array of dimensions (M x M X K) containing the K M X M (co)variances. These matrices must be symmetric and positive semi-definite.

# Returns
An array of dimensions (∑N x M) containing the data points.

# Examples
```julia-repl
julia>
```
"""
function simulate(K::Int64, N::Array{Int64}, μ::Array{Float64}, Σ::Array{Float64})::Tuple{Array{Float64}, Array{Float64}}
    # Assert dimensions
    @assert K == size(Σ)[3] "You must pass a covariance matrix for each separate cluster ..."
    @assert K == size(μ)[2] "You must pass a mean array for each separate cluster ..."
    @assert K == size(N)[2] "You must pass a sample size (N) for each separate cluster ..."
    # Assert that dim(μ) and dim(Σ) are equal in the number of columns
    @assert size(Σ)[2] == size(μ)[1]
    # Dimensions
    N_total = sum(N)
    M = size(Σ)[1]
    # Open up result matrix
    X = zeros((N_total, M))
    lbls = zeros((N_total,1))
    # For each cluster, generate data
    start_idx = 1
    for k ∈ 1:K
        # Dims
        N_k = N[k]
        μ_k = μ[:,k]
        Σ_k = Σ[:,:,k]
        # Create end idx
        end_idx = k == 1 ? N_k : start_idx + N_k - 1
        # Draw points
        X[start_idx:end_idx, :] = MvNormal(μ_k, Σ_k) |>
            x -> rand(x, N_k)';
        # Labels
        lbls[start_idx:end_idx, :] = fill(k, N_k)
        # Overwrite start values
        if k < K
            start_idx = N_k + start_idx
        end;
    end;
    # Cat X and labels
    Xtmp = hcat(X, lbls);
    # Shuffle
    Xtmp = Xtmp[shuffle(1:end), :];
    # Return X and labels
    return Xtmp[:, 1:M], Xtmp[:, end]
end;

## End module
end;


text:
function network_advanced(T, V, cfg)::Matrix{T}

    # Get data
    data = get_network_data(T, V, cfg)

    # Get flags
    flags = get_network_flags(cfg)

    # Compute advanced data
    advanced_data = compute_advanced_data(data, flags, cfg)

    # Send to main kernel
    v , _ = advanced_kernel(advanced_data, flags, cfg)

    v
end

function compute_advanced_data(data::NetworkData{T,V},
                            flags, cfg)::AdvancedProblem{T,V} where {T,V}

    i,j,v = data.coords

    idx = findfirst(x -> x < 1, i)
    idx != nothing && throw("Indices no good")
    idx = findfirst(x -> x < 1, j)
    idx != nothing && throw("Indices no good")

    m = max(i[end], j[end])
    A = sparse(i,j,v,m,m)
    A = A + A'

    cc = connected_components(SimpleWeightedGraph(A))

    t = @elapsed G = laplacian(A)
    csinfo("Time taken to construct graph laplacian = $t", cfg["suppress_messages"] in TRUELIST)

    nodemap, polymap = Matrix{V}(undef,0,0), Matrix{V}(undef,0,0)
    cellmap = Matrix{T}(undef,0,0)

    solver = get_solver(cfg)

    sources, grounds, finite_grounds =
                get_sources_and_grounds(data, flags, G, nodemap)

    source_map = Matrix{eltype(A)}(undef,0,0)
    AdvancedProblem(G, cc, nodemap, polymap, RasterMeta(),
                 sources, grounds, source_map, finite_grounds, V(-1), V(0), cellmap, solver)

end


text:
"""
    optimal(s::General, d::KullbackLeibler, m::DRO)

Solves the DRO model `m` with Kullback-Leibler divergence using the Ipopt solver.
"""
function optimal(s::General, d::KullbackLeibler, m::DRO)
    stats.optimizer = "Ipopt"

    model = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@NLconstraint(model, sum(p[i]*log(p[i]/m.q[i]) for i in eachindex(p)) <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::Burg, m::DRO)

Solves the DRO model `m` with Burg entropy using the Ipopt solver.
"""
function optimal(s::General, d::Burg, m::DRO)
    stats.optimizer = "Ipopt"

    model = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@NLconstraint(model, sum(m.q[i]*log(m.q[i]/p[i]) for i in eachindex(p)) <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::Hellinger, m::DRO)

Solves the DRO model `m` with Hellinger distance using the Ipopt solver.
"""
function optimal(s::General, d::Hellinger, m::DRO)
    stats.optimizer = "Ipopt"

    model = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@NLconstraint(model, sum((sqrt(p[i]) - sqrt(m.q[i]))^2 for i in eachindex(p)) <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::ChiSquare, m::DRO)

Solves the DRO model `m` with χ²-distance using the Ipopt solver.
"""
function optimal(s::General, d::ChiSquare, m::DRO)
    stats.optimizer = "Ipopt"

    model = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@NLconstraint(model, sum(((p[i] - m.q[i])^2)/p[i] for i in eachindex(p)) <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::ModifiedChiSquare, m::DRO)

Solves the DRO model `m` with modified χ²-distance using the Ipopt solver.
"""
function optimal(s::General, d::ModifiedChiSquare, m::DRO)
    stats.optimizer = "Ipopt"

    model = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@NLconstraint(model, sum(((p[i] - m.q[i])^2)/m.q[i] for i in eachindex(p)) <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::Linf, m::DRO)

Solves the DRO model `m` with l-infinity using the CPLEX solver.
"""
function optimal(s::General, d::Linf, m::DRO)
    stats.optimizer = "CPLEX"

    model = JuMP.Model(JuMP.with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)
    JuMP.@variable(model, p_max >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, p_max .>= p - m.q)
    JuMP.@constraint(model, p_max .>= - p + m.q)
    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@constraint(model, p_max <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::Lone, m::DRO)

Solves the DRO model `m` with l-1 using the CPLEX solver.
"""
function optimal(s::General, d::Lone, m::DRO)
    stats.optimizer = "CPLEX"

    model = JuMP.Model(JuMP.with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)
    JuMP.@variable(model, p_abs[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)

    JuMP.@constraint(model, p_abs .>= p - m.q)
    JuMP.@constraint(model, p_abs .>= - p + m.q)
    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@constraint(model, sum(p_abs) <= m.ε)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, d::Ltwo, m::DRO)

Solves the DRO model `m` with l-2 using the CPLEX solver.
"""
function optimal(s::General, d::Ltwo, m::DRO)
    stats.optimizer = "CPLEX"

    model = JuMP.Model(JuMP.with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    JuMP.@variable(model, p[1:length(m.q)] >= 0)

    JuMP.@objective(model, Max, m.c'*p)
    JuMP.@constraint(model, sum(p) == 1)
    JuMP.@constraint(model, sum((p[i] - m.q[i])^2 for i in eachindex(p)) <= m.ε^2)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

"""
    optimal(s::General, m::Simplex)

Solves the Simplex model using the CPLEX solver.
"""
function optimal(s::General, m::Simplex)
    stats.optimizer = "CPLEX"

    model = JuMP.Model(JuMP.with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    JuMP.@variable(model, m.lb[i] <= p[i = 1:length(m.q)] <= m.ub[i])

    JuMP.@objective(model, Min, sum(0.5 * (p[i] - m.q[i])^2 for i in eachindex(p)))
    JuMP.@constraint(model, sum(p) == 1)

    JuMP.optimize!(model)
    return JuMP.value.(p)
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

###############
# Prepare FMU #
###############

cd(dirname(@__FILE__))
pathToFMU = joinpath(pwd(), "../model/IO.fmu")

myFMU = fmiLoad(pathToFMU)

# create arrays for in place getter
vR = zeros(Real, 3)
vI = zeros(Integer, 3)
vB = [false, false, false]
vS = ["test","test"]

c1 = fmiInstantiate!(myFMU; loggingOn=true)

#########################
# Testing Setter/Getter #
#########################

fmiEnterInitializationMode(c1)
@test fmiSetReal(c1, "p_real", 5.0) == 0
@test fmiSetReal(c1, ["u_real", "p_real"], [7.0, 8.0]) == 0
@test fmiSetInteger(c1, "p_integer", 6) == 0
@test fmiSetInteger(c1, ["u_integer"], [4]) == 0
@test fmiSetBoolean(c1, "p_boolean", false) == 0
@test fmiSetBoolean(c1, ["u_boolean"], [true]) == 0
@test fmiSetString(c1, "p_string", "New String") == 0
fmiGetReal!(c1, ["p_real", "y_real", "u_real"], vR)
@test vR == [8.0, 7.0, 7.0]
p = fmiGetReal(c1, "p_real")
@test p == 8.0
fmiGetInteger!(c1, ["u_integer", "y_integer", "p_integer"], vI)
@test vI == [4, 4, 6]
i = fmiGetInteger(c1, "y_integer")
@test i == 4
fmiGetBoolean!(c1, ["u_boolean", "y_boolean", "p_boolean"], vB)
@test vB == [true, true, false]
b = fmiGetBoolean(c1, "p_boolean")
@test b == false
fmiGetString!(c1, ["p_string", "p_string"], vS)
@test vS == ["New String", "New String"]
string = fmiGetString(c1, "p_string")
@test string == "New String"
fmiExitInitializationMode(c1)

############
# Clean up #
############

fmiReset(c1)
fmiTerminate(c1)
fmiUnload(myFMU)


text:
function DQN_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* maximum(value(π, 𝒟[:sp]), dims=1)
end

function DQN(;π::DiscreteNetwork,
              N::Int,
              ΔN=4,
              π_explore=ϵGreedyPolicy(LinearDecaySchedule(1., 0.1, floor(Int, N/2)), π.outputs),
              c_opt::NamedTuple=(;),
              log::NamedTuple=(;),
              c_loss=td_loss(),
              target_fn=DQN_target,
              prefix="",
              kwargs...)

     OffPolicySolver(;agent=PolicyParams(π=π, π_explore=π_explore, π⁻=deepcopy(π)),
                      log=LoggerParams(;dir="log/dqn", log...),
                      N=N,
                      ΔN=ΔN,
                      c_opt = TrainingParams(;loss=c_loss, name=string(prefix, "critic_"), epochs=ΔN, c_opt...),
                      target_fn=target_fn,
                      kwargs...)
end


text:
user = ENV["USER"]
host = readlines("/etc/hostname")[1]
kernel = split(readlines("/proc/version")[1])[3]
term = ENV["TERM"]
shell = ENV["SHELL"]
tasks = length(filter(x -> isa(tryparse(Int, x), Int), readdir("/proc")))

mem = readlines("/proc/meminfo")[1:3]
total, _, avail = map(x -> fld(parse(Int, split(x)[2]), 1000), mem)

uptime = parse(Int, split(readlines("/proc/uptime")[1], '.')[1])
d = trunc(Int, (uptime / 60 / 60 / 24))
h = trunc(Int, (uptime / 60 / 60 % 24))
m = trunc(Int, (uptime / 60 % 60))

print("$user@$host
kernel\t$kernel
term\t$term
shell\t$shell
tasks\t$tasks
mem\t$(avail)m / $(total)m
uptime\t$(d)d $(h)h $(m)m")

max_stars_repo_name: aosync/fetch

max_issues_repo_name: aosync/fetch

max_forks_repo_name: aosync/fetch


text:
function maxmin(t, cone, a, z)
    # Solves (ASP) as a part of Algorithm 3 (Step 4)
    cd(string(root,"Results//ex",id,"//opts",t))
    n = length(a[1,:])
    c = length(a[:,1])
    iteration = length(cone[:])
    m2 = Model(solver=BaronSolver(CplexLibName="/panfs/roc/groups/10/qizh/rishabh/Mangi/CPLEX_Studio128/cplex/bin/x86-64_linux/libcplex1280.so", PrLevel = 0, MaxTime = 100))
    # m2 = Model(solver=BaronSolver(PrLevel = 0, MaxTime = 100))

    M = 4000
    epsilon = 1

    @variable(m2, 160 >= y[1:n] >= -5)
    @variable(m2, d[1:2*n] >= 0)
    @variable(m2,  1 >= gamma[1:iteration, 1:c] >= 0)
    @variable(m2, 200 >= t[1:n] >= 0)
    @variable(m2, gamma2[1:c] >= 0)

    @objective(m2, Min, -sum(t[i] for i = 1:n))

    @NLconstraint(m2, sum(y[i]*d[i] for i = 1:n) - sum(y[i]*d[i+n] for i = 1:n) == sum(t[i] for i = 1:n))

    for i = 1:n
        @constraint(m2, -t[i] - sum(gamma2[j]*a[j,i]*z[j] for j = 1:c) <= -y[i])
    end

    for i = 1:n
        @constraint(m2, -t[i] + sum(gamma2[j]*a[j,i]*z[j] for j = 1:c) <= y[i])
    end

    for i = 1:n
        @constraint(m2, d[i]+d[n+i] <= 1)
    end

    for j = 1:c
        if z[j] >= 0.1
            @constraint(m2, sum(a[j,i]*d[i] for i = 1:n) - sum(a[j,i]*d[i+n] for i = 1:n) <= 0)
        end
    end

    @constraint(m2, CC[i = 1:iteration, j = 1:n], y[j] == sum(gamma[i,k]*cone[i][k,j] for k = 1:length(cone[i][:,1])))

    status = solve(m2; suppress_warnings = true)
    if status != :Infeasible
        return -1*getobjectivevalue(m2), -1*getobjbound(m2)
    else
        return 0, 0
    end
end

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP


text:
using CUDArt

_getbytes(x::DataType,d) = sizeof(x)

_getbytes(x::NTuple,d)=sum(map(y->_getbytes(y,d), x))

_getbytes(x::AbstractCudaArray,d)=(haskey(d,x) ? 0 : (d[x]=1; length(x) * sizeof(eltype(x))))

function _getbytes(x::DenseArray,d)
    haskey(d,x) && return 0; d[x]=1
    total = sizeof(x)
    if !isbits(eltype(x))
        for i = 1:length(x)
            isize = _getbytes(x[i],d)
            # @show (typeof(x), i, isize)
            total += isize
        end
    end
    # @show (typeof(x), total)
    # @show (eltype(x), size(x), total)
    return total
end

function _getbytes(x,d)
    total = sizeof(x)
    isbits(x) && return total
    haskey(d,x) && return 0; d[x]=1
    fieldNames = typeof(x).names
    if fieldNames != ()
        for fieldName in fieldNames
            isdefined(x, fieldName) || continue
            f = x.(fieldName)
            fieldBytes = _getbytes(f,d)
            fieldSize = (isa(f, Union(AbstractArray,CudaArray)) ? (eltype(f),size(f)) : ())
            # @show (typeof(x), fieldName, fieldSize, fieldBytes)
            total += fieldBytes
        end
    end
    return total
end

getbytes(x)=_getbytes(x, ObjectIdDict())


text:
using FinalExample
using Test

@testset "FinalExample.jl" begin
    @test my_fu(2, 2) == 10
    @test my_fu(2, 3) == 13
    @test my_fu(2, 3) == 12
    # Write your own tests here.
end


text:
include("functions_encoding.jl")

struct FastaRecord
    id::Array{UInt8, 1}
    description::Array{UInt8, 1}
    seq::Array{UInt8, 1}
end

struct fasta_record
    id::String
    description::String
    seq::String
end

# get_index_of_diff_sites returns an array of integers which are the sites
# at which two vectors of bit-level encoded nucleotides differ
# NB nuc_bit_V_A should be the reference, as it stands, because we ignore N/-/?
# in nuc_bit_V_B
function get_index_of_diff_sites(nuc_bit_V_A, nuc_bit_V_B)

    A = Array{Int, 1}()

    for i in 1:length(nuc_bit_V_A)

        @inbounds x = nuc_bit_V_A[i]
        @inbounds y = nuc_bit_V_B[i]

        if x != y && y < 240
            push!(A, i)
        end
    end

    return A
end

# get_list_of_differences applies get_index_of_diff_sites between nuc_bit_V_A
# and every sequence stored in nuc_bit_array
function get_list_of_differences(nuc_bit_array, nuc_bit_V_A)

    L = Array{Array{Int, 1}, 1}(undef, size(nuc_bit_array, 2))

    for i in 1:length(L)
        L[i] = get_index_of_diff_sites(nuc_bit_V_A, view(nuc_bit_array, :, i))
    end

    return L
end

function readlinebyte(s::IO)::Vector{UInt8}
    line = readuntil(s, 0x0a, keep=true)
    i = length(line)
    if i == 0 || line[i] != 0x0a
        return line
    elseif i < 2 || line[i-1] != 0x0d
        return line[1:length(line) - 1]
    else
        return line[1:length(line) - 2]
    end
end

function readuntilspace(V::Vector{UInt8})::Vector{UInt8}
    N = Vector{UInt8}()
    for c in V
        if c == 0x20
            return N
        end
        push!(N, c)
    end
    return N
end

function get_alignment_dimensions(filepath::AbstractString)

    n = 0 # number of sequences
    l = 0 # length of sequence
    first = true

    open(filepath, "r") do io
        while !eof(io)

            if n > 1
                first = false
            end

            line = readlinebyte(io)

            if line[1] == 0x3e # '>'
                n+=1
            end

            if first && line[1] != 0x3e # '>'
                l+=length(line)
            end
        end

    end

    return l, n
end

function get_fasta_descriptions(filepath::AbstractString)

    A = Array{String,1}()

    open(filepath, "r") do io
        while !eof(io)
            l = readline(io)
            if string(l[1]) == ">"
                description = lstrip(l, '>')
                id = split(description)[1]
                push!(A, description)
            end
        end
    close(io)
    end

    return A
end

function read_fasta_alignment_bytes(filename, channel::Channel)
    open(filename, "r") do io

        firstHeader = true
        seqbuffer = IOBuffer()
        description = Vector{UInt8}
        id = Vector{UInt8}

        while !eof(io)

            L = readlinebyte(io)

            if firstHeader
                description = L[2:length(L)]
                id = readuntilspace(description)
                firstHeader = false

            elseif L[1] == 0x3e # '>'

                FR = FastaRecord(id, description, take!(seqbuffer))
                put!(channel, FR)

                description = L[2:length(L)]
                id = readuntilspace(description)

            else
                write(seqbuffer, L)
            end
        end

        FR = FastaRecord(id, description, take!(seqbuffer))
        put!(channel, FR)

    end
end

function read_fasta_alignment(filepath::AbstractString, channel::Channel)

    open(filepath, "r") do io

        n = 0

        seq_buffer = String("")
        id = String("")
        description = String("")

        while !eof(io)

            l = readline(io)
            n += 1

            if n == 1
                if string(l[1]) != ">"
                    e = error("badly formed fasta file")
                    throw(e)
                end

                description = lstrip(l, '>')
                id = split(description)[1]

                first = false

            elseif eof(io)

                seq_buffer = seq_buffer * uppercase(l)

                fr = fasta_record(id, description, seq_buffer)
                put!(channel, fr)

                break

            elseif string(l[1]) == ">"

                fr = fasta_record(id, description, seq_buffer)
                put!(channel, fr)

                description = lstrip(l, '>')
                id = split(description)[1]
                seq_buffer = String("")

            else

                seq_buffer = seq_buffer * uppercase(l)

            end
        end

    close(io)
    end
end

function populate_byte_array_get_names(filepath::AbstractString)
    byte_2_byte = make_byte_array()

    alignment_dim = get_alignment_dimensions(filepath)
    height = alignment_dim[1]
    width = alignment_dim[2]

    A = Array{UInt8,2}(undef, height, width)

    ids = Array{String,1}(undef, width)

    ch = Channel{FastaRecord}((channel_arg) -> read_fasta_alignment_bytes(filepath, channel_arg))

    j = 1
    @inbounds for record in ch
        ids[j] = String(record.id)
        i = 1
        for nuc in record.seq
            if byte_2_byte[nuc] == 0
                e = error("unknown nucleotide in alignment: ", Char(nuc))
                throw(e)
            end
            A[i, j] = byte_2_byte[nuc]
            # A[i, j] = byte_dict[nuc]
            i+=1
        end

        j+=1
    end

    return A, ids
end

function get_seq_from_1D_byte_array(byte_V)
    nuc_dict = make_nuc_dict()

    char_V = Array{Char, 1}(undef, length(byte_V))

    for i in 1:size(byte_V, 1)
        char_V[i] = nuc_dict[byte_V[i]]
    end

    return join(char_V)
end

# function score_alignment_column(nuc_bit_sequence)
#     score_dict = make_score_dict()
#
#     score = 0::Int64
#     for nuc in nuc_bit_sequence
#         score+=score_dict[nuc]
#     end
#     return score
# end
#
# # TODO: parallelise this (should be easy - but needs testing)
# function score_alignment(nuc_bit_array)
#
#     A = fill(0::Int64, size(nuc_bit_array, 2))
#
#     for i in 1:size(nuc_bit_array, 2)
#         A[i] = score_alignment_column(view(nuc_bit_array, :, i))
#     end
#
#     return A
# end

function score_alignment_column2(nuc_bit_sequence)
    score_array = make_score_array()

    score = 0::Int64
    for nuc in nuc_bit_sequence
        score+=score_array[nuc]
    end
    return score
end

function score_alignment2(nuc_bit_array)

    A = fill(0::Int64, size(nuc_bit_array, 2))

    Threads.@threads for i in 1:size(nuc_bit_array, 2)
        A[i] = score_alignment_column2(view(nuc_bit_array, :, i))
    end

    return A
end

function consensus(nuc_bit_array)
    score = Vector{Int}(undef, 256)
    for i in 1:length(score)
        score[i] = 5
    end
    score[136] = 1 # => 'A'
    score[72] = 2 # => 'G'
    score[40] = 3 # => 'C'
    score[24] = 4 # => 'T'

    consensus_array = zeros(Int, size(nuc_bit_array, 1), 5)
    for j in 1:size(nuc_bit_array, 2)
        for i in 1:size(nuc_bit_array, 1)
            nuc = nuc_bit_array[i, j]
            consensus_array[i, score[nuc]] += 1
        end
    end

    backtranslate = [0x88, 0x48, 0x28, 0x18]
    consensus = Vector{UInt8}(undef, size(nuc_bit_array, 1))
    for i in 1:size(consensus_array, 1)
        maxval = consensus_array[i,1]
        maxidx = 1
        for j in 2:4
            if consensus_array[i,j] > maxval
                maxval = consensus_array[i,j]
                maxidx = j
            end
        end
        consensus[i] = backtranslate[maxidx]
    end

    return consensus
end

max_stars_repo_name: snake-flu/julialign

max_issues_repo_name: snake-flu/julialign

max_forks_repo_name: snake-flu/julialign


text:
using ArgParse
using BSON
using MarkovModels
using Distributed

# BUT specific package to work with the SGE cluster
# See https://github.com/BUTSpeechFIT/BUTSGEManager
using BUTSGEManager

function getargs()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--add-sil", "-s"
            help = "add silence add the beginning/end of each utterance"
            action = :store_true
        "--jobs", "-j"
            help = "number of jobs to run in parallel"
            arg_type = Int
            default = 4
        "--jobs-args", "-a"
            help = "arguments to pass to the parallel engine"
            arg_type = String
            default = ""
        "trans"
            help = "list of word level transcription"
            arg_type = String
            required = true
        "lexicon"
            help = "mapping word -> pronunciation"
            arg_type = String
            required = true
        "hmms"
            help = "unit hmms"
            arg_type = String
            required = true
        "outdir"
            help = "output directory where will be stored the alignment graphs"
            arg_type = String
            required = true
    end
    args = parse_args(s)
end

const args = getargs()
const trans = args["trans"]
const lexiconfile = args["lexicon"]
const hmmsfile = args["hmms"]
const outdir = args["outdir"]

@info "loading the pronunciations"
const pronuns = Dict{String, Vector{FSM}}()
open(lexiconfile, "r") do f
    for line in eachline(f)
        tokens = split(line)
        word = tokens[1]
        pronun = tokens[2:end]
        ps = get(pronuns, word, FSM[])
        push!(ps, LinearFSM(pronun))
        pronuns[word] = ps
    end
end

@info "compiling pronunciation graphs"
const lexicon = Dict{String, FSM}()
for (k, v) in pronuns
    lexicon[k] = union(v...) |> minimize!
end

# Make sure the directory exists
run(`mkdir -p $outdir`)

@info "initializing the workers"
addprocs(SGEManager(args["jobs"]), args = args["jobs-args"],
         exeflags = "--project=$(Base.active_project())")
@everywhere using BSON
@everywhere using MarkovModels
@everywhere const outdir = $outdir
@everywhere const lexicon = $lexicon
@everywhere const hmmsfile = $hmmsfile
@everywhere hmms = BSON.load(hmmsfile)[:hmms]

@info "compiling alignment graphs"
@sync @distributed for line in readlines(trans)
    tokens = split(line)
    uttid = tokens[1]
    words = tokens[2:end]
    if args["add-sil"]
        insert!(words, 1, "sil")
        push!(words, "sil")
    end

    ali = compose!(compose!(LinearFSM(words), lexicon), hmms) |> removenilstates!
    bson(joinpath(outdir, uttid * ".bson"), Dict(:ali => ali))
end

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
using FuzzyFinder
using Test

@testset "FuzzyFinder.jl" begin
	global const f = FuzzyFinder
	global const collections = ["migrations.jl",
            "random_migrations.jl",
            "random_admin_log.jl",
            "api_user.doc",
            "user_group.doc",
            "users.txt",
            "accounts.txt",
            "123.jl",
            "test123test.jl"
            ]
	global const cased_collections = ["MIGRATIONS.jl",
            "random_MiGRations.jl",
            "random_admin_log.jl",
            "migrations.doc",
            "user_group.doc",
            "users.txt",
            ]
	global const dict_collections =[
							  Dict("name"=> "migrations.jl"),
							  Dict("name"=> "random_migrations.jl"),
							  Dict("name"=> "random_admin_log.jl"),
							  Dict("name"=> "api_user.doc"),
							  Dict("name"=> "user_group.doc"),
							  Dict("name"=> "users.txt"),
							  Dict("name"=> "accounts.txt"),
							  ]
	result = f.fuzzyfinder("txt", collections)
	expected = ["users.txt", "accounts.txt"]
	@test collect(result) == expected

	result = f.fuzzyfinder("miGr", cased_collections)
	expected = ["MIGRATIONS.jl", "migrations.doc", "random_MiGRations.jl"]
	@test collect(result) == expected

	collection_arr = ["fuuz", "fuz", "fufuz"]
	result = f.fuzzyfinder("fuz", collection_arr)
	expected = ["fuz", "fuuz", "fufuz"]
	@test collect(result) == expected

	result = f.fuzzyfinder(".txt", collections)
	expected = ["users.txt", "accounts.txt"]
	@test collect(Set(result)) == expected

	result = f.fuzzyfinder("rami", collections)
	expected = ["random_migrations.jl", "random_admin_log.jl"]
	@test collect(result) == expected

	result = f.fuzzyfinder("mi", collections)
	expected = ["migrations.jl", "random_migrations.jl", "random_admin_log.jl"]
	@test collect(result) == expected

	result = f.fuzzyfinder("user", collections)
	expected = ["user_group.doc", "users.txt", "api_user.doc"]
	@test collect(result) == expected

	result = f.fuzzyfinder(123, collections)
	expected = ["123.jl", "test123test.jl"]
	@test collect(result) == expected

	result = f.fuzzyfinder("user", dict_collections, fetcher = x -> get(x, "name", ""))
	expected = [Dict("name"=> "user_group.doc"), Dict("name" => "users.txt"), Dict("name" => "api_user.doc")]
	@test collect(result) == expected

	collection = ["zzfuz", "nnfuz", "aafuz", "ttfuz", "wow!", "julia"]
	result = f.fuzzyfinder("fuz", collection, sorted=false)
	expected = ["zzfuz", "nnfuz", "aafuz", "ttfuz"]
	@test collect(result) == expected

end


text:
import StringUtils

function lines_to_bigints(a::Array{String,1})
    result::Array{BigInt,1} = []
    for s in a
        big = parse(BigInt, chomp(s))
        push!(result, big)
    end
    return result
end

function main()
    if length(ARGS) > 0
        filename = ARGS[1]
        a = readlines(filename)
        m = lines_to_bigints(a)
        s = string(sum(m))
        @printf("The first 10 digits are %s\n", s[1:10])
    else
        error("this program needs a single numeric argument")
    end
end

main()


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/MbedTLS_jll.jl

baremodule MbedTLS_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libmbedcrypto, libmbedtls, libmbedx509

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libmbedcrypto_handle = C_NULL
libmbedcrypto_path = ""
libmbedtls_handle = C_NULL
libmbedtls_path = ""
libmbedx509_handle = C_NULL
libmbedx509_path = ""

if Sys.iswindows()
    const libmbedcrypto = "libmbedcrypto.dll"
    const libmbedtls = "libmbedtls.dll"
    const libmbedx509 = "libmbedx509.dll"
elseif Sys.isapple()
    const libmbedcrypto = "@rpath/libmbedcrypto.5.dylib"
    const libmbedtls = "@rpath/libmbedtls.13.dylib"
    const libmbedx509 = "@rpath/libmbedx509.1.dylib"
else
    const libmbedcrypto = "libmbedcrypto.so.5"
    const libmbedtls = "libmbedtls.so.13"
    const libmbedx509 = "libmbedx509.so.1"
end

function __init__()
    global libmbedcrypto_handle = dlopen(libmbedcrypto)
    global libmbedcrypto_path = dlpath(libmbedcrypto_handle)
    global libmbedtls_handle = dlopen(libmbedtls)
    global libmbedtls_path = dlpath(libmbedtls_handle)
    global libmbedx509_handle = dlopen(libmbedx509)
    global libmbedx509_path = dlpath(libmbedx509_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = dirname(libmbedtls_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libmbedcrypto_path() =libmbedcrypto_path
get_libmbedtls_path() = libmbedtls_path
get_libmbedx509_path() = libmbedx509_path

end  # module MbedTLS_jll

max_stars_repo_name: ranocha/julia

max_issues_repo_name: ranocha/julia

max_forks_repo_name: ranocha/julia


text:
@testset "754.reach-a-number.jl" begin
    @test reach_number(3) == 2
    @test reach_number(2) == 3
end


text:
module Solution
export multipleof3regex

multipleof3regex = r"^(0?(1(01*0)*1)*)+$"

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export Ensemble, DataEnsemble, AbstractEnsemble, length, output,
    ActivationEnsemble

Base.length(ens::AbstractEnsemble)         = Base.length(ens.neurons)
Base.start(ens::AbstractEnsemble)          = Base.start(ens.neurons)
Base.done(ens::AbstractEnsemble, i::Int64) = Base.done(ens.neurons, i)
Base.next(ens::AbstractEnsemble, i::Int64) = Base.next(ens.neurons, i)
Base.ndims(ens::AbstractEnsemble)          = Base.ndims(ens.neurons)
Base.size(ens::AbstractEnsemble)           = Base.size(ens.neurons)
Base.size(ens::AbstractEnsemble, i::Int64) = Base.size(ens.neurons, i)
Base.getindex(ens::AbstractEnsemble, a...) = Base.getindex(ens.neurons, a...)

function Base.length(conn::Connection)
    conn.size
end

function Ensemble{T,N}(net::Net, name::Symbol, neurons::Array{T,N},
                       params::Vector=Param[]; phase::Phase=TrainTest)
    ens = Ensemble{T,N}(name, neurons, [], [:value, :∇, :inputs, :∇inputs],
                        Dict(), params, phase, 1)
    add_ensemble(net, ens)
    ens
end

function AbstractString(ensemble::AbstractEnsemble)
    return "$(ensemble.name) $(size(ensemble.neurons))"
end

function init_inputs(ensemble::AbstractEnsemble, net::Net)
    batch_size = net.batch_size
    neuron = ensemble.neurons[1]
    shape = size(ensemble)
    for name in filter((x) -> x in [:inputs, :∇inputs], fieldnames(neuron))
        field = getfield(neuron, name)
        typ = typeof(field)
        key = symbol(ensemble.name,name)
        for t = 1:net.time_steps
            target = name == :inputs ? :value : :∇
            for (source_index, connection) in enumerate(ensemble.connections)
                key = symbol(ensemble.name,name,source_index)
                source_target = symbol(connection.source.name, target)
                if connection.source.net_subgroup != ensemble.net_subgroup
                    set_buffer(net, source_target, zeros(Float32, size(connection.source)..., batch_size))
                end
                source = get(net.buffers[1][t],
                             source_target,
                             nothing)
                if source == nothing
                    @assert !connection.recurrent
                    sink = Array(eltype(typ), 0, shape..., batch_size)
                    set_buffer(net, key, sink, t)
                elseif all(connection.is_dim_fixed)
                    connection.copy = false
                    if connection.recurrent && t == 1
                        set_buffer(net, key, zeros(eltype(typ), connection.size, batch_size), t)
                    elseif connection.recurrent && t > 1
                        source = get_buffer(net, source_target, t - 1)
                        set_buffer(net, key, reshape(source, (connection.size, batch_size)), t)
                    else
                        set_buffer(net, key, reshape(source, (connection.size, batch_size)), t)
                    end
                elseif connection.is_one_to_one
                    connection.copy = false
                    if connection.recurrent && t == 1
                        set_buffer(net, key, zeros(eltype(typ), size(ensemble)..., batch_size), t)
                    elseif connection.recurrent && t > 1
                        source = get_buffer(net, source_target, t - 1)
                        set_buffer(net, key, source, t)
                    else
                        set_buffer(net, key, source, t)
                    end
                else
                    _shape = [shape...]
                    if LATTE_DISABLE_TILING && length(_shape) >= 3
                        num_tiles = batch_size
                    else
                        _shape[2] = TILE_SIZE * 2
                        num_tiles = num_threads
                    end
                    _shape = _shape[!connection.is_dim_fixed]
                    sink = zeros(eltype(typ), connection.size, _shape..., num_tiles)
                    set_buffer(net, key, sink, t)
                end
            end
        end
    end
end

function init_params(ensemble::AbstractEnsemble, net::Net)
    if :params in fieldnames(ensemble)
        for param in ensemble.params
            param.value = get_buffer(net, param.name)
            param.gradient = get_buffer(net, param.gradient_name)
            param.hist = zeros(param.value)
            set_buffer(net, param.hist_name, param.hist)
            @latte_mpi param.request = @eval ccall((:init_request, $libComm), Cint, ())
        end
    end
end

@doc """
Initialize `ensemble` of neuron type `T`
Allocate a buffer for each field in T

If the neurons are connected to a contigious region in another ensemble, map the
inputs array using ArrayViews and remove the need for a CopyTask
"""
@inbounds function init(ensemble::AbstractEnsemble, net::Net)
    batch_size = net.batch_size
    neuron = ensemble.neurons[1]
    shape = size(ensemble)
    for name in fieldnames(neuron)
        # log_info("Initializing field $(name) for ensemble $(ensemble.name)")
        # @time begin
        field = getfield(neuron, name)
        typ = typeof(field)
        key = symbol(ensemble.name,name)
        if name in [:inputs, :∇inputs]
            # skip
        elseif name in [:value, :∇]
            # :value and :∇ should allocate space for every item in a batch
            init_buffer(net, key, (shape..., batch_size))
        elseif isa(field, Batch)
            arr = Array(typeof(field.init), shape..., batch_size)
            fill!(arr, field.init)
            set_buffer(net, key, arr)
            push!(ensemble.batch_fields, name)
        elseif !isbits(typ)
            uniform_across_dim = [true for _ in 1:length(shape)]
            first = getfield(ensemble.neurons[ones(Int, length(shape))...], name)
            for d in 1:length(shape)
                for i in 1:shape[d]
                    idx = [x == d ? i : 1 for x in 1:length(shape)]
                    if first !== getfield(ensemble.neurons[idx...], name)
                        uniform_across_dim[d] = false
                        break
                    end
                end
            end
            ensemble.arg_dim_info[key] = uniform_across_dim
            _shape = []
            iter = []
            for i in 1:length(shape)
                if !uniform_across_dim[i]
                    push!(_shape, shape[i])
                    push!(iter, shape[i])
                else
                    push!(iter, 1)
                end
            end
            elem = getfield(ensemble.neurons[ones(Int, length(shape))...], name)
            if !LOSSY_GRADIENTS && contains(string(name), "∇")
                _shape = (_shape..., num_threads)
            end
            buf = Array(eltype(typ), size(elem)..., _shape...)
            for index in CartesianRange(CartesianIndex(tuple(ones(Int, length(iter))...)),
                                        CartesianIndex(tuple(iter...)))
                _index = []
                for i in 1:length(shape)
                    if !uniform_across_dim[i]
                        push!(_index, index.I[i])
                    end
                end
                if !LOSSY_GRADIENTS && contains(string(name), "∇")
                    for i in 1:num_threads
                        @inbounds buf[:, _index..., i] = copy(getfield(ensemble.neurons[index], name))
                    end
                else
                    @inbounds buf[:, _index...] = getfield(ensemble.neurons[index], name)
                end
            end
            set_buffer(net, key, buf)
        else
            # Neuron specific fields, only one instance per neuron instead of
            # per batch
            arr = map((elem) -> getfield(elem, name), ensemble.neurons)
            for t = 1:net.time_steps
                set_buffer(net, key, arr, t)
            end
        end
        # end  # @time begin
    end
end

"""
Get the forward function for a neuron in `ens`
"""
function get_forward(ens::AbstractEnsemble)
    neuron = ens.neurons[1]
    return forward(neuron)
end

"""
Get the backward function for a neuron in `ens`
"""
function get_backward(ens::AbstractEnsemble)
    neuron = ens.neurons[1]
    return backward(neuron)
end

function one_to_one(ndims)
    idx = [symbol(:i_, d) for d in 1:ndims]
    fn = Expr(:(->), Expr(:tuple, idx...),
                     Expr(:block, Expr(:tuple, idx...)))
    return eval(fn)
end

type ActivationEnsemble{T <: Neuron, N} <: AbstractEnsemble
    name         :: Symbol
    neurons      :: Array{T,N}
    connections  :: Vector{Connection}
    batch_fields :: Vector{Symbol}
    arg_dim_info :: Dict{Symbol, Vector{Bool}}
    params       :: Vector{Param}
    phase        :: Phase
    net_subgroup :: Cint
end

function ActivationEnsemble{T,N}(net::Net, name::Symbol, neurons::Array{T,N},
                                 input_ensemble::AbstractEnsemble; phase::Phase=TrainTest)
    params = Param[]
    ens = ActivationEnsemble{T,N}(name, neurons, [],
        [:value, :∇, :inputs, :∇inputs], Dict(), params, phase, 1)
    add_ensemble(net, ens)
    mapping = one_to_one(ndims(input_ensemble))
    add_connections(net, input_ensemble, ens, mapping)
    ens
end

function init_inputs(ensemble::ActivationEnsemble, net::Net)
    for conn in ensemble.connections
        # Activation ensembles occur in place so they should be
        # in the same subgroup as its input
        @assert conn.source.net_subgroup == ensemble.net_subgroup
    end
end

@inbounds function init(ensemble::ActivationEnsemble, net::Net)
    batch_size = net.batch_size
    neuron = ensemble.neurons[1]
    shape = size(ensemble)
    for name in fieldnames(neuron)
        field = getfield(neuron, name)
        typ = typeof(field)
        key = symbol(ensemble.name,name)
        @assert length(ensemble.connections) == 1
        if name in [:inputs, :∇inputs]
            target = name == :inputs ? :value : :∇
            for (source_index, connection) in enumerate(ensemble.connections)
                connection.copy = false
                key = symbol(ensemble.name,name,source_index)
                for t = 1:net.time_steps
                    source = get_buffer(net, symbol(connection.source.name, target), t)
                    set_buffer(net, key, reshape(source, (1, size(source)...)), t)
                    key = symbol(ensemble.name,target)
                    set_buffer(net, key, source, t)
                end
            end
        elseif name in [:value, :∇]
            # Skip
        elseif isa(field, Batch)
            arr = Array(typeof(field.init), shape..., batch_size)
            fill!(arr, field.init)
            set_buffer(net, key, arr)
            push!(ensemble.batch_fields, name)
        else
            arr = map((elem) -> getfield(elem, name), ensemble.neurons)
            for t = 1:net.time_steps
                set_buffer(net, key, arr, t)
            end
        end
    end
end

macro output(ensemble_type, fn)
    return Expr(:escape, quote
        function $(fn.args[1].args[1])(ens::$ensemble_type)
            return $(Expr(:quote, fn))
        end
    end)
end

function init_inputs(ensemble::NormalizationEnsemble, net::Net)
    for connection in ensemble.connections
        if connection.source.net_subgroup != ensemble.net_subgroup
            key = symbol(connection.source.name, :value)
            set_buffer(net, key, zeros(Float32, size(connection.source)..., net.batch_size))
        end
    end
end


text:
### This code was taken and modified from jl.
# jl is licensed under the MIT License:

# > Copyright (c) 2012: John Myles White, Tim Holy, and other contributors.
# > Copyright (c) 2016: Patrick Kofod Mogensen, John Myles White, Tim Holy,
# >                     and other contributors.
# > Copyright (c) 2017: Patrick Kofod Mogensen, Asbjørn Nilsen Riseth,
# >                     John Myles White, Tim Holy, and other contributors.
# >
# > Permission is hereby granted, free of charge, to any person
# > obtaining a copy of this software and associated documentation files
# > (the "Software"), to deal in the Software without restriction,
# > including without limitation the rights to use, copy, modify, merge,
# > publish, distribute, sublicense, and/or sell copies of the Software,
# > and to permit persons to whom the Software is furnished to do so,
# > subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be
# > included in all copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# > NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# > BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# > ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# > CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.

struct LightOptions{T} #add iterations?
    x_tol::T
    f_tol::T
    g_tol::T
end
LightOptions() = LightOptions{Float64}(0.0,0.0,1e-8)

# mutable struct BFGSState{Tx, Tm, T,G} <: AbstractOptimizerState
#     x::Tx
#     x_previous::Tx
#     g_previous::G
#     f_x_previous::T
#     dx::Tx
#     dg::Tx
#     u::Tx
#     invH::Tm
#     s::Tx
#     @add_linesearch_fields()
# end
function initial_state!(state, method::BFGS{M}, d, initial_x::AbstractArray{T}) where {T,M}
    n = length(initial_x)
    copyto!(state.x, initial_x)
    retract!(M, initial_x)
    value_gradient!!(d, initial_x)
    project_tangent!(M, gradient(d), initial_x)
    copyto!(state.g_previous, gradient(d))
    method.initial_invH((state.invH, initial_x))
    state.alpha = one(T)#/2
    nothing
end
# function initial_state!(state, method::BFGS, d, initial_x::AbstractArray{T}) where T
#     n = length(initial_x)
#     copyto!(state.x, initial_x)
#     retract!(method.manifold, initial_x)
#     value_gradient!!(d, initial_x)
#     project_tangent!(method.manifold, gradient(d), initial_x)
#     copyto!(state.g_previous, gradient(d))
#     method.initial_invH(state.invH, initial_x)
#     state.alpha = one(T)
#     nothing
# end

function uninitialized_state(initial_x::SizedVector{P,T}) where {T,P}
    BFGSState(similar(initial_x), # Maintain current state in state.x
        similar(initial_x), # Maintain previous state in state.x_previous
        similar(initial_x), # Store previous gradient in state.g_previous
        T(NaN), # Store previous f in state.f_x_previous
        similar(initial_x), # Store changes in position in state.dx
        similar(initial_x), # Store changes in gradient in state.dg
        similar(initial_x), # Buffer stored in state.u
        SymmetricMatrix{T,P}(),
        #RecursiveMatrix{T,P,P}(), # Store current invH in state.invH
        similar(initial_x), # Store current search direction in state.s
        # T(NaN),            # Keep track of previous descent value ⟨∇f(x_{k-1}), s_{k-1}⟩
        similar(initial_x), # Buffer of x for line search in state.x_ls
        one(T))
end
function initial_convergence(d, state, method::AbstractOptimizer, initial_x, options)
    gradient!(d, initial_x)
    # norm(gradient(d), Inf) < options.g_tol
    SIMDArrays.maximum_abs(gradient(d)) < options.g_tol
end

function optimize_heavy(d::D, initial_x::Tx, method::M,
                  options::LightOptions = LightOptions(),
                  state = initial_state(method, options, d, initial_x)) where {D<:DifferentiableObject, M<:AbstractOptimizer, Tx <: AbstractArray}
    # @show state

    f_limit_reached, g_limit_reached, h_limit_reached = false, false, false
    x_converged, f_converged, f_increased = false, false, false

    g_converged = initial_convergence(d, state, method, initial_x, options)
    # @show g_converged
    converged = g_converged

    # prepare iteration counter (used to make "initial state" trace entry)
    iteration = 0
    # println("\n")
    while !converged && iteration < 60#options.iterations #uncomment if you add back iterations field to LightOptions
        iteration += 1
        # @show gradient(d)
        update_state!(d, state, method) && break # it returns true if it's forced by something in update! to stop (eg dx_dg == 0.0 in BFGS, or linesearch errors)
        # @show state.x_previous, state.g_previous
        update_g!(d, state, method) # TODO: Should this be `update_fg!`?
        x_converged, f_converged,
        g_converged, converged, f_increased = assess_convergence(state, d, options)
        !converged && update_h!(d, state, method) # only relevant if not converged
        # debug() && # @show state.f_x_previous, value(d)
        # println("\n")
    end # while
    # println("\n")
    # @show iteration, x_converged, f_converged, g_converged, converged
    f_incr_pick = disallow_f_increases(options) && f_increased
    f_incr_pick ? ( state.x_previous, state.f_x_previous ) : ( state.x, value(d) )
end
disallow_f_increases(options) = false
# disallow_f_increases(options::Options) = options.allow_f_increases

function update_g!(d, state, method::FirstOrderOptimizer{M}) where M# Union{FirstOrderOptimizer, Newton}
    # Update the function value and gradient
    value_gradient!(d, state.x)
    # if M <: FirstOrderOptimizer #only for methods that support manifold optimization
    project_tangent!(method.manifold, gradient(d), state.x)
    # end
end


text:
factorial5 = n -> reduce(*, [1:n;])
println(factorial5(5))

function get_pi(n)
  s = 0
  for (x, y) in zip(rand(n), rand(n))
    if x^2 + y^2 <= 1
      s += 1
    end
  end
  return s / n * 4
end

get_pi2 = n -> sum([x^2 + y^2 <= 1 ? 1 : 0 for (x, y) in zip(rand(n), rand(n))])/n*4
get_pi3 = n -> sum([x^2 + y^2 <= 1 ? 1 : 0 for x in rand(n), y in rand(n)])/n*4

pis = [get_pi2(10000) for i in range(1, 100)]
println(mean(pis))


text:
export matrix_algebra

################################################################################
#
#  Basic field access
#
################################################################################

degree(A::AlgMat) = A.degree

dim(A::AlgMat) = A.dim

base_ring(A::AlgMat{T, S}) where {T, S} = A.base_ring::parent_type(T)

coefficient_ring(A::AlgMat) = A.coefficient_ring

basis(A::AlgMat) = A.basis

has_one(A::AlgMat) = true

elem_type(A::AlgMat{T, S}) where { T, S } = AlgMatElem{T, AlgMat{T, S}, S}

elem_type(::Type{AlgMat{T, S}}) where { T, S } = AlgMatElem{T, AlgMat{T, S}, S}

order_type(::AlgMat{fmpq, S}) where { S } = AlgAssAbsOrd{AlgMat{fmpq, S}, elem_type(AlgMat{fmpq, S})}
order_type(::Type{AlgMat{fmpq, S}}) where { S } = AlgAssAbsOrd{AlgMat{fmpq, S}, elem_type(AlgMat{fmpq, S})}

order_type(::AlgMat{T, S}) where { T <: NumFieldElem, S } = AlgAssRelOrd{T, fractional_ideal_type(order_type(parent_type(T))), AlgMat{T, S}}
order_type(::Type{AlgMat{T, S}}) where { T <: NumFieldElem, S } = AlgAssRelOrd{T, fractional_ideal_type(order_type(parent_type(T))), AlgMat{T, S}}

matrix_algebra_type(K::Field) = matrix_algebra_type(typeof(K))

matrix_algebra_type(::Type{T}) where {T <: Field} = AlgMat{elem_type(T), dense_matrix_type(elem_type(T))}

# Returns the dimension d of the coefficient_ring of A, so that dim(A) = degree(A)^2 + d.
function dim_of_coefficient_ring(A::AlgMat)
  if base_ring(A) == coefficient_ring(A)
    return 1
  end
  @assert _base_ring(coefficient_ring(A)) == base_ring(A)
  return dim(coefficient_ring(A))
end

################################################################################
#
#  Commutativity
#
################################################################################

iscommutative_known(A::AlgMat) = (A.iscommutative != 0)

@doc Markdown.doc"""
    iscommutative(A::AlgMat) -> Bool

Returns `true` if $A$ is a commutative ring and `false` otherwise.
"""
function iscommutative(A::AlgMat)
  if iscommutative_known(A)
    return A.iscommutative == 1
  end
  dcr = dim_of_coefficient_ring(A)
  if dim(A) == degree(A)^2*dcr
    A.iscommutative = 2
    return false
  end

  mt = multiplication_table(A, copy = false)
  for i = 1:dim(A)
    for j = i + 1:dim(A)
      if mt[i, j, :] != mt[j, i, :]
        A.iscommutative = 2
        return false
      end
    end
  end
  A.iscommutative = 1
  return true
end

################################################################################
#
#  Basis matrix
#
################################################################################

function assure_has_basis_matrix(A::AlgMat)
  if isdefined(A, :basis_matrix)
    return nothing
  end

  d2 = degree(A)^2
  if coefficient_ring(A) == base_ring(A)
    M = zero_matrix(base_ring(A), dim(A), d2)
    for i = 1:dim(A)
      N = matrix(A[i])
      @assert length(N) == d2
      for (j, n) in enumerate(N)
        M[i, j] = n
      end
    end
    A.basis_matrix = M
  else
    dcr = dim_of_coefficient_ring(A)
    M = zero_matrix(base_ring(A), dim(A), d2*dcr)
    for i = 1:dim(A)
      N = matrix(A[i])
      for (j, n) in enumerate(N)
        jj = (j - 1)*dcr
        for k = 1:dcr
          M[i, jj + k] = coefficients(n, copy = false)[k]
        end
      end
    end
    A.basis_matrix = M
  end
  return nothing
end

function basis_matrix(A::AlgMat; copy::Bool = true)
  assure_has_basis_matrix(A)
  if copy
    return deepcopy(A.basis_matrix)
  else
    return A.basis_matrix
  end
end

################################################################################
#
#  Multiplication Table
#
################################################################################

function assure_has_multiplication_table(A::AlgMat{T, S}) where { T, S }
  if isdefined(A, :mult_table)
    return nothing
  end

  B = basis(A)
  d = dim(A)
  mt = Array{T, 3}(undef, d, d, d)
  for i in 1:d
    for j in 1:d
      mt[i, j, :] = coefficients(B[i]*B[j], copy = false)
    end
  end
  A.mult_table = mt
  return nothing
end

@doc Markdown.doc"""
    multiplication_table(A::AlgMat; copy::Bool = true) -> Array{RingElem, 3}

Given an algebra $A$ this function returns the multiplication table of $A$:
If the function returns $M$ and the basis of $A$ is $e_1,\dots, e_n$ then
it holds $e_i \cdot e_j = \sum_k M[i, j, k] \cdot e_k$.
"""
function multiplication_table(A::AlgMat; copy::Bool = true)
  assure_has_multiplication_table(A)
  if copy
    return deepcopy(A.mult_table)
  else
    return A.mult_table
  end
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    matrix_algebra(R::Ring, n::Int) -> AlgMat

Returns $\mathrm{Mat}_n(R)$.
"""
function matrix_algebra(R::Ring, n::Int)
  A = AlgMat{elem_type(R), dense_matrix_type(elem_type(R))}(R)
  n2 = n^2
  A.dim = n2
  A.degree = n
  B = Vector{elem_type(A)}(undef, n2)
  for i = 1:n
    ni = n*(i - 1)
    for j = 1:n
      M = zero_matrix(R, n, n)
      M[j, i] = one(R)
      B[ni + j] = A(M)
    end
  end
  A.basis = B
  A.one = identity_matrix(R, n)
  A.canonical_basis = 1
  A.issimple = 1
  A.issemisimple = 1
  return A
end

# Constructs Mat_n(S) as an R-algebra
@doc Markdown.doc"""
    matrix_algebra(R::Ring, S::Ring, n::Int) -> AlgMat

Returns $\mathrm{Mat}_n(S)$ as an $R$-algebra.
It is assumed that $S$ is an $R$-algebra.
"""
function matrix_algebra(R::Ring, S::Ring, n::Int)
  A = AlgMat{elem_type(R), dense_matrix_type(elem_type(S))}(R, S)
  n2 = n^2
  A.dim = n2*dim(S)
  A.degree = n
  B = Vector{elem_type(A)}(undef, dim(A))
  for k = 1:dim(S)
    n2k = n2*(k - 1)
    for i = 1:n
      ni = n2k + n*(i - 1)
      for j = 1:n
        M = zero_matrix(S, n, n)
        M[j, i] = S[k]
        B[ni + j] = A(M)
      end
    end
  end
  A.basis = B
  A.one = identity_matrix(S, n)
  A.canonical_basis = 2
  return A
end

@doc Markdown.doc"""
    matrix_algebra(R::Ring, gens::Vector{<: MatElem}; isbasis::Bool = false)
      -> AlgMat

Returns the matrix algebra over $R$ generated by the matrices in `gens`.
If `isbasis` is `true`, it is assumed that the given matrices are an $R$-basis
of this algebra, i. e. that the spanned $R$-module is closed under
multiplication.
"""
function matrix_algebra(R::Ring, gens::Vector{<:MatElem}; isbasis::Bool = false)
  @assert length(gens) > 0
  A = AlgMat{elem_type(R), dense_matrix_type(elem_type(R))}(R)
  A.degree = nrows(gens[1])
  A.one = identity_matrix(R, degree(A))
  if isbasis
    A.dim = length(gens)
    bas = Vector{elem_type(A)}(undef, dim(A))
    for i = 1:dim(A)
      bas[i] = A(gens[i])
    end
    A.basis = bas
    return A
  end
  A.gens = A.(gens)

  d = degree(A)
  d2 = degree(A)^2
  span = deepcopy(gens)
  push!(span, identity_matrix(R, d))
  M = zero_matrix(R, max(d2, length(span)), d2) # the maximal possible dimension is d^2
  pivot_rows = zeros(Int, d2)
  new_elements = Set{Int}()
  cur_rank = 0
  for i = 1:length(span)
    cur_rank == d2 ? break : nothing
    new_elt = _add_row_to_rref!(M, reshape(collect(span[i]), :), pivot_rows, cur_rank + 1)
    if new_elt
      push!(new_elements, i)
      cur_rank += 1
    end
  end

  # Build all possible products
  while !isempty(new_elements)
    cur_rank == d2 ? break : nothing
    i = pop!(new_elements)
    b = span[i]

    n = length(span)
    for r = 1:n
      s = b*span[r]
      for l = 1:n
        t = span[l]*s
        new_elt = _add_row_to_rref!(M, reshape(collect(t), :), pivot_rows, cur_rank + 1)
        if !new_elt
          continue
        end
        push!(span, t)
        cur_rank += 1
        push!(new_elements, length(span))
        cur_rank == d2 ? break : nothing
      end
      cur_rank == d2 ? break : nothing
    end
  end

  A.dim = cur_rank
  A.basis_matrix = sub(M, 1:cur_rank, 1:d2)
  bas = Vector{elem_type(A)}(undef, dim(A))
  for i = 1:dim(A)
    N = zero_matrix(R, degree(A), degree(A))
    for j = 1:d
      jd = (j - 1)*d
      for k = 1:d
        N[k, j] = basis_matrix(A, copy = false)[i, jd + k]
      end
    end
    bas[i] = A(N)
  end
  A.basis = bas

  return A
end

@doc Markdown.doc"""
    matrix_algebra(R::Ring, S::Ring, gens::Vector{<: MatElem};
                   isbasis::Bool = false)
      -> AlgMat

Returns the matrix algebra over $S$ generated by the matrices in `gens` as
an algebra with base ring $R$.
It is assumed $S$ is an algebra over $R$.
If `isbasis` is `true`, it is assumed that the given matrices are an $R$-basis
of this algebra, i. e. that the spanned $R$-module is closed under
multiplication.
"""
function matrix_algebra(R::AbstractAlgebra.NCRing, S::AbstractAlgebra.NCRing, gens::Vector{<:MatElem}; isbasis::Bool = false)
  @assert length(gens) > 0
  A = AlgMat{elem_type(R), dense_matrix_type(elem_type(S))}(R, S)
  A.degree = nrows(gens[1])
  A.one = identity_matrix(S, degree(A))
  if isbasis
    A.dim = length(gens)
    bas = Vector{elem_type(A)}(undef, dim(A))
    for i = 1:dim(A)
      bas[i] = A(gens[i])
    end
    A.basis = bas

    return A
  end

  d = degree(A)
  d2 = degree(A)^2
  span = deepcopy(gens)
  push!(span, identity_matrix(S, d))
  dcr = dim(S)
  max_dim = d2*dcr
  M = zero_matrix(R, max(length(gens), max_dim), max_dim)
  pivot_rows = zeros(Int, max_dim)
  new_elements = Set{Int}()
  cur_rank = 0
  v = Vector{elem_type(R)}(undef, max_dim)
  for i = 1:length(span)
    cur_rank == max_dim ? break : nothing
    @assert length(span[i]) == d2
    for (j, s) in enumerate(span[i])
      jj = (j - 1)*dcr
      for k = 1:dcr
        v[jj + k] = coefficients(s, copy = false)[k]
      end
    end
    new_elt = _add_row_to_rref!(M, v, pivot_rows, cur_rank + 1)
    if new_elt
      push!(new_elements, i)
      cur_rank += 1
    end
  end

  # Build all possible products
  while !isempty(new_elements)
    cur_rank == max_dim ? break : nothing
    i = pop!(new_elements)
    b = span[i]

    n = length(span)
    for r = 1:n
      s = b*span[r]
      for l = 1:n
        t = span[l]*s
        @assert length(t) == d2
        for (j, s) in enumerate(t)
          jj = (j - 1)*dcr
          for k = 1:dcr
            v[jj + k] = coefficients(s, copy = false)[k]
          end
        end
        new_elt = _add_row_to_rref!(M, v, pivot_rows, cur_rank + 1)
        if !new_elt
          continue
        end
        push!(span, t)
        cur_rank += 1
        push!(new_elements, length(span))
        cur_rank == max_dim ? break : nothing
      end
      cur_rank == max_dim ? break : nothing
    end
  end

  A.dim = cur_rank
  A.basis_matrix = sub(M, 1:cur_rank, 1:max_dim)
  bas = Vector{elem_type(A)}(undef, dim(A))
  for i = 1:dim(A)
    N = zero_matrix(S, degree(A), degree(A))
    for j = 1:d
      jd = (j - 1)*d
      for k = 1:d
        jkd = (jd + k - 1)*dcr
        t = Vector{elem_type(_base_ring(S))}(undef, dcr)
        for l = 1:dcr
          t[l] = basis_matrix(A, copy = false)[i, jkd + l]
        end
        N[k, j] = S(t)
      end
    end
    bas[i] = A(N)
  end
  A.basis = bas

  return A
end

###############################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::AlgMat)
  print(io, "Matrix algebra of dimension ", dim(A), " over ", base_ring(A))
end

################################################################################
#
#  Deepcopy
#
################################################################################

function deepcopy_internal(A::AlgMat{T, S}, dict::IdDict) where { T, S }
  B = AlgMat{T, S}(base_ring(A))
  for x in fieldnames(typeof(A))
    if x != :base_ring x != :coefficient_ring && isdefined(A, x)
      setfield!(B, x, deepcopy_internal(getfield(A, x), dict))
    end
  end
  B.base_ring = A.base_ring
  B.coefficient_ring = A.coefficient_ring
  return B
end

################################################################################
#
#  Inclusion of matrices
#
################################################################################

function _check_matrix_in_algebra(M::S, A::AlgMat{T, S}, short::Type{Val{U}} = Val{false}) where {S, T, U}
  if nrows(M) != degree(A) || ncols(M) != degree(A)
    if short == Val{true}
      return false
    end
    return false, zeros(base_ring(A), dim(A))
  end

  d2 = degree(A)^2
  B = basis_matrix(A, copy = false)
  if coefficient_ring(A) == base_ring(A)
    t = zero_matrix(base_ring(A), 1, d2)
    @assert length(M) == d2
    for (i, m) in enumerate(M)
      t[1, i] = m
    end
  else
    dcr = dim_of_coefficient_ring(A)
    t = zero_matrix(base_ring(A), 1, d2*dcr)
    @assert length(M) == d2
    for (i, m) in enumerate(M)
      ii = (i - 1)*dcr
      for j = 1:dcr
        t[1, ii + j] = coefficients(m, copy = false)[j]
      end
    end
  end
  b, N = can_solve_with_solution(B, t, side = :left)
  if short == Val{true}
    return b
  end
  s = elem_type(base_ring(A))[ N[1, i] for i = 1:dim(A) ]
  return b, s
end

################################################################################
#
#  Center
#
################################################################################

@doc Markdown.doc"""
    center(A::AlgMat) -> AlgAss, AbsAlgAssMor

Returns the center $C$ of $A$ and the inclusion $C \to A$.
"""
function center(A::AlgMat{T, S}) where {T, S}
  if isdefined(A, :center)
    return A.center::Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}
  end

  # Unlike for AlgAss, we should cache the centre even if A is commutative
  # since it is of a different type, so A !== center(A)[1].
  # Otherwise center(A)[1] !== center(A)[1] which is really annoying.
  B, mB = AlgAss(A)
  C, mC = center(B)
  mD = compose_and_squash(mB, mC)
  A.center = C, mD
  return C, mD
end

################################################################################
#
#  Conversion to AlgAss
#
################################################################################

function AlgAss(A::AlgMat{T, S}) where {T, S}
  K = base_ring(A)
  B = AlgAss(K, multiplication_table(A))
  B.issimple = A.issimple
  B.issemisimple = A.issemisimple
  AtoB = hom(A, B, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
  if isdefined(A, :center)
    Z, ZtoA = center(A)
    B.center = (Z, compose_and_squash(AtoB, ZtoA))
  end
  if isdefined(A, :decomposition)
    dec = Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(B))}[]
    for (C, CtoA) in A.decomposition
      CtoB = compose_and_squash(AtoB, CtoA)
      push!(dec, (C, CtoB))
    end
    B.decomposition = dec
  end
  if isdefined(A, :maps_to_numberfields)
    fields_and_maps = Tuple{AnticNumberField, AbsAlgAssToNfAbsMor{typeof(B), elem_type(B)}}[]
    for (K, AtoK) in A.maps_to_numberfields
      BtoK = AbsAlgAssToNfAbsMor(B, K, AtoK.mat, AtoK.imat)
      push!(fields_and_maps, (K, BtoK))
    end
    B.maps_to_numberfields = fields_and_maps
  end
  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
end

################################################################################
#
#  Canonical basis
#
################################################################################

# Checks whether A[(j - 1)*n + i] == E_ij, where E_ij = (e_kl)_kl with e_kl = 1 if i =k and j = l and e_kl = 0 otherwise.
function iscanonical(A::AlgMat)
  if A.canonical_basis != 0
    return A.canonical_basis == 1
  end

  if coefficient_ring(A) != base_ring(A)
    A.canonical_basis = 2
    return false
  end

  n = degree(A)
  if dim(A) != n^2
    A.canonical_basis = 2
    return false
  end

  for j = 1:n
    nj = (j - 1)*n
    for i = 1:n
      m = matrix(A[nj + i], copy = false)
      for k = 1:n
        for l = 1:n
          if k == i && l == j
            if !isone(m[k, l])
              A.canonical_basis = 2
              return false
            end
          else
            if !iszero(m[k, l])
              A.canonical_basis = 2
              return false
            end
          end
        end
      end
    end
  end

  A.canonical_basis = 1
  return true
end

################################################################################
#
#  Opposite algebra
#
################################################################################

function opposite_algebra(A::AlgMat)
  B, BtoA = AlgAss(A)
  O, BtoO = opposite_algebra(B)
  return O, compose_and_squash(BtoO, inv(BtoA))
end


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

function pdhg_restart_x_y(
    problem::StandardLinearProgram,
    exitcriterion::ExitCriterion;
    γ=1.0, L=10, restartfreq=Inf,
    io=nothing)

    @info("Running pdhg_restart_x_y...")
    @info("γ = $(γ)")
    @info("L = $(L)")
    @info("restartfreq = $(restartfreq)")
    if !isnothing(io)
        flush(io)
    end

    A_T, b, c = problem.A_T, problem.b, problem.c
    prox = problem.prox

    d, n = size(A_T)
    x0 = zeros(d)
    y0 = zeros(n)
    m = 1

    ##### Start of pdhg_restart

    # Log initial measure
    starttime = time()
    results = Results()
    init_fvaluegap, init_metricLP = compute_fvaluegap_metricLP(x0, y0, problem)
    logresult!(results, 1, 0.0, init_fvaluegap, init_metricLP)

    outer_k = 1
    exitflag = false

    while !exitflag
        # Init of PDHG
        τ, σ =  1.0 / (γ * L), 1.0 * γ / L
        x_bar = deepcopy(x0)
        x = deepcopy(x0)
        y = deepcopy(y0)
        x_tilde = zero(x0)
        y_tilde = zero(y0)

        k = 1
        restartflag = false
        while !exitflag && !restartflag

            y[:] = y[:] + σ * ((x_bar' * A_T)' - b)
            x_pre = deepcopy(x)
            x[:] = x[:]  - τ * (A_T * y[:] + c[:])
            x = prox(x, τ)
            x_bar[:] = 2 * x[:] -  x_pre[:]  # Correct term from original paper
            # x_bar[:] = x_bar[:] + (x[:] -  x_pre[:])  # Wrong term from CB
            x_tilde[:] = x_tilde[:] + x[:]
            y_tilde[:] = y_tilde[:] + y[:]

            # Logging and checking exit condition
            # set restartflag when reached some measure
            if outer_k % exitcriterion.loggingfreq == 0
                x_out = x_tilde / k
                y_out = y_tilde / k

                # Progress measures
                fvaluegap, metricLP = compute_fvaluegap_metricLP(x_out, y_out, problem)

                elapsedtime = time() - starttime
                @info "elapsedtime: $elapsedtime"
                @info "outer_k: $(outer_k), fvaluegap: $(fvaluegap), metricLP: $(metricLP)"
                elapsedtime = time() - starttime
                logresult!(results, outer_k, elapsedtime, fvaluegap, metricLP)
                if !isnothing(io)
                    flush(io)
                end

                exitflag = checkexitcondition(exitcriterion, outer_k, elapsedtime, metricLP)
                if exitflag
                    break
                end

                if k >= restartfreq * m || (restartfreq == Inf && metricLP <= 0.5 * init_metricLP)
                    @info "<===== RESTARTING"
                    @info "k ÷ m: $(k ÷ m)"
                    @info "elapsedtime: $elapsedtime"
                    @info "outer_k: $(outer_k), fvaluegap: $(fvaluegap), metricLP: $(metricLP)"
                    if !isnothing(io)
                        flush(io)
                    end

                    x0, y0 = deepcopy(x_out), deepcopy(y_out)
                    init_fvaluegap = fvaluegap
                    init_metricLP = metricLP
                    restartflag = true
                    break
                end
            end

            k += 1
            outer_k += 1
        end
    end

    return results
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP