text:
using ZChop # For chopping small imaginary parts in ρ
using Distributed

"""
    (t, FI, QFI) = Eff_QFI_PD_pure(Nj, Ntraj, Tfinal, dt; kwargs... )

Evaluate the continuous-time FI and QFI of a final strong measurement for the
estimation of the frequency ω with continuous photo-detection monitoring
of each half-spin particle affected by noise at an angle θ,
with efficiency η = 1 using SSE (Stochastic Schrödinger equation).

The function returns a tuple `(t, FI, QFI)` containing the time vector and the
vectors containing the FI and average QFI

# Arguments

* `Nj`: number of spins
* `Ntraj`: number of trajectories for the SSE
* `Tfinal`: final time of evolution
* `dt`: timestep of the evolution
* `κ = 1`: the noise coupling
* `θ = 0`: noise angle (0 parallel, π/2 transverse)
* `ω = 0`: local value of the frequency
"""
function Eff_QFI_PD_pure(Nj::Int64,     # Number of spins
    Ntraj::Int64,                       # Number of trajectories
    Tfinal::Float64,                    # Final time
    dt::Float64;                        # Time step
    κ = 1.,                             # Noise coupling
    θ = 0.,                             # Noise angle
    ω = 0)                              # Frequency

    Ntime = Int(floor(Tfinal/dt)) # Number of timesteps

    dimJ = Int(2^Nj)   # Dimension of the corresponding Hilbert space

    # Define an array of noise channels (they are all sparse matrices)
    cj = sqrt(κ/2) *
        [(cos(θ) * σ_j(:z, j, Nj) + sin(θ) * σ_j(:x, j, Nj))
            for j = 1:Nj]

    # Hamiltonian of the Nj-atom system
    H = ω * σ(:z, Nj) / 2
    dH = σ(:z, Nj) / 2        # Derivative of H wrt the parameter ω

    dW() = sqrt(dt) * randn(Nj) # Function that returns a Wiener increment vector

    # Kraus-like operator, trajectory-independent part
    # see Eqs. (50-51)
    M0 = sparse(I - 1im * H * dt - 0.5 * dt * sum([c'*c for c in cj]) )
    M1 = sqrt(dt) * cj

    # Derivative of the Kraus-like operator wrt to ω
    dM0 = - 1im * dH * dt
    dM1 = 0

    # Probability of detecting a single photon (the probability of
    # detecting more than one in dt is negligible)
    # p_PD = η * sum(tr(ρ cj'*cj)) *dt,
    # but each noise operator gives 1/2 * κ * id when squared
    pPD = 0.5 * Nj * κ * dt;

    # Initial state of the dynamics
    ψ0 = ghz_state(Nj)

    t = (1 : Ntime) * dt

    # Run evolution for each trajectory, and build up the average
    # for FI and final strong measurement QFI
    # Uses Julia @distributed macro for simple parallelization
    result = @distributed (+) for ktraj = 1 : Ntraj

        ψ = copy(ψ0) # Assign initial state to each trajectory

        # Derivative of ρ wrt the parameter
        # Initial state does not depend on the paramter
        dψ = zero(ψ)
        ϕ = zero(ψ)

        # Vectors for the FI and QFI for the specific trajectory
        FisherT = zero(t)
        QFisherT = zero(t)

        for jt=1:Ntime
            # Has the photon been detected?
            if (rand() < real(pPD)) # Detected
                # Choose randomly which channel detected the photon
                # (with equal probabiltiy)
                ch = rand(1:Nj)

                new_ψ = M1[ch] * ψ

                norm_ψ = norm(new_ψ)

                ϕ = (dM1 * ψ + M1[ch] * ϕ) / norm_ψ
            else # Not detected
                new_ψ = M0 * ψ

                norm_ψ = norm(new_ψ)

                ϕ = (dM0 * ψ + M0 * ϕ) / norm_ψ
            end

            # Now we can renormalize ψ and evaluate its derivative wrt ω
            ψ = new_ψ / norm_ψ;
            dψ = ϕ - 0.5 * (ϕ' * ψ + ψ' * ϕ) * ψ;

            # We evaluate the classical FI for the continuous measurement
            FisherT[jt] = real((ψ' * ϕ + ϕ' * ψ)^2);

            # We evaluate the QFI for a final strong measurement done at time t
            QFisherT[jt] = 4 * real( dψ' * dψ + (dψ' * ψ)^2);
        end

        # Use the reduction feature of @distributed for
        # (at the end of each cicle, sum the result to result)
        hcat(FisherT, QFisherT)
    end

    return (t, result[:,1] / Ntraj, result[:,2] / Ntraj)
end


text:
@testset "Conjugates" begin
   k, a = wildanger_field(3, 13)

   @testset "Regulator" begin
     C_29 = Hecke.qAdicConj(k, 29)
     C_31 = Hecke.qAdicConj(k, 31)

     @test valuation(regulator(k, C_29)) == 2
     @test valuation(regulator(k, C_31)) == 2

     @test valuation(regulator(k, C_29, flat = true)) == 3
     @test valuation(regulator(k, C_31, flat = true)) == 2

     @test valuation(regulator(k, C_29, flat = false)) == 2
     @test valuation(regulator(k, C_31, flat = false)) == 2
   end

   @testset "Completion" begin

     L1, mL1 = completion(k, 37, 1)
     L2, mL2 = completion(k, 37, 3)

     @test degree(parent(mL1(a))) == 2
     @test degree(parent(mL2(a))) == 1

     lp = prime_decomposition(maximal_order(k), 37)

     @test valuation(a - preimage(mL1, mL1(a)), lp[1][1]) >= 10
     @test valuation(a - preimage(mL2, mL2(a)), lp[2][1]) >= 10
   end
end

ext: jl


text:
immutable Baseline{T} <: Model
    β₀::T
end

Base.show(io::IO, ::Type{Baseline}) = print(io, "Heuristic Baseline")

function Base.show(io::IO, m::Baseline)
    @printf(
        io,
        "Heuristic Baseline (β₀ = %.3f)",
        m.β₀,
    )
end

initial_parameters(::Type{Baseline}) = Float64[0.0]

canonical_parameters(m::Baseline) = (m.β₀, )

function Baseline(parameters::Vector)
    β₀ = parameters[1]
    return Baseline(β₀)
end

@inline function predict(m::Baseline, x2::Real, t2::Real, x1::Real, t1::Real)
    β₀ = m.β₀
    z = β₀
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::Baseline,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    gr = (1.0, )
    return gr
end


text:
"""
    DoubleState{ST<:DualValuedState} <: DualValuedComposedState

State composed by two identical states ST, used to represent states of a
density matrix.
"""
mutable struct DoubleState{ST<:FiniteBasisState} <: FiniteBasisState
    n::Int64 # number of sites
    σ_row::ST # left
    σ_col::ST # right

    space_dim::Int
end

## Constructors
DoubleState(sr::FiniteBasisState, sc::FiniteBasisState=deepcopy(sr)) =
    same_basis(sr, sc) ? DoubleState(nsites(sr),
                                     deepcopy(sr),
                                     deepcopy(sc),
                                     spacedimension(sr)^2) : error("not same basis")

# Specialized Constructors
DoubleState{ST}(n, i_σ=0) where ST = set!(DoubleState(ST(n, 0), ST(n, 0)), i_σ)

# Property Accessors
@inline spacedimension(state::DoubleState) = state.space_dim
@inline nsites(state::DoubleState) = 2*nsites(state.σ_row)
@inline local_dimension(state::DoubleState{ST}) where {ST} = local_dimension(ST)
@inline eltype(state::DoubleState) = eltype(row(state))

toint(state::DoubleState) = _toint(col(state), row(state)) #toint(state.σ_row, state.σ_col)
index(state::DoubleState) = toint(state) + 1 # was before
index_to_int(state::DoubleState, id) = (id -1)
flipped(a::DoubleState, b::DoubleState) = (flipped(a.σ_row, b.σ_row), flipped(a.σ_col, b.σ_col))

# custom accessor
row(v::DoubleState) = v.σ_row
col(v::DoubleState) = v.σ_col
@inline config(v::DoubleState) = (config(v.σ_row), config(v.σ_col))
#@inline config(v::DoubleState) = return config(v.σ_row), config(v.σ_col)

# checks
same_basis(v::DoubleState{ST}, v2::DoubleState{ST2}) where {ST, ST2} =
    ST==ST2 && same_basis(row(v), row(v2)) && same_basis(col(v), col(v2))

function flipat!(rng::AbstractRNG, v::DoubleState, i::Int)
    i > v.n ? flipat!(rng, v.σ_row, i-v.n) : flipat!(rng, v.σ_col, i)
end

function flipat_fast!(rng::AbstractRNG, v::DoubleState, i::Int)
    i > v.n ? flipat_fast!(rng, v.σ_row, i-v.n) : flipat_fast!(rng, v.σ_col, i)
end

function setat!(v::DoubleState, i::Int, val)
    i > v.n ? setat!(v.σ_row, i-v.n, val) : setat!(v.σ_col, i, val)
end

#=function apply!(state::DoubleState, changes::DoubleStateChanges)
    for (id, val)=row(changes)
        setat!(row(state), id, val)
    end
    for (id, val)=col(changes)
        setat!(col(state), id, val)
    end
    return state
end=#

set_index!(v::DoubleState, i::Integer) = set!(v, index_to_int(v, i))
function set!(v::DoubleState, i::Integer)
    row = div(i, spacedimension(v.σ_row)) #row = i>>(nsites(state.σ_row))
    col = i - row*spacedimension(v.σ_row)#col = i - (row<< nsites(state.σ_row))

    set!(v.σ_col, row) #i
    set!(v.σ_row, col) #j
    return v
end
set!(v::DoubleState, i_row, i_col) = (set!(row(v), i_row); set!(col(v), i_col);  v)

function rand!(rng::AbstractRNG, state::DoubleState)
    rand!(rng, state.σ_row)
    rand!(rng, state.σ_col)
    state
end

# -------------- Base.show extension for nice printing -------------- #
Base.show(io::IO, mm::MIME"text/plain", bs::DoubleState) =
    print(io, "DoubleState(",bs.n,") : \n\tRow: ", bs.σ_row,"\n\tCol: ", bs.σ_col)

Base.show(io::IO, bs::DoubleState) = print(io, "(",bs.σ_row,", ", bs.σ_col, ")")


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.

module MXNet

using Reexport

# we put everything in the namespace mx, because there are a lot of
# functions with the same names as built-in utilities like "zeros", etc.
export mx
module mx

using Base.Broadcast: Broadcasted, DefaultArrayStyle
using Libdl
using LinearAlgebra
using Markdown
using Printf
using Statistics
using Random

using Formatting
using MacroTools

# Functions from base that we can safely extend and that are defined by libmxnet.
import Base.Broadcast: broadcasted
import Base.Iterators: filter

###############################################################################
#  exports
###############################################################################

# exceptions.jl
export AbstractMXError,
       MXError

# symbolic-node.jl
export SymbolicNode,
       Variable,
       @var

# ndarray.jl
export NDArray,
       context,
       expand_dims,
       @inplace,
       # activation funcs
       σ,
       sigmoid,
       relu,
       softmax,
       log_softmax,
       # broadcast utils
       broadcast_to,
       broadcast_axis,
       broadcast_axes

# executor.jl
export Executor,
       bind,
       simple_bind,
       forward,
       backward

# context.jl
export Context,
       cpu,
       gpu,
       num_gpus

# model.jl
export AbstractModel,
       FeedForward,
       predict

# nn-factory.jl
export MLP

# metric.jl
export AbstractEvalMetric,
       ACE,
       Accuracy,
       MSE,
       MultiACE,
       MultiMetric,
       NMSE,
       SeqMetric

# kvstore.jl
export KVStore,
       init!,
       pull!,
       barrier,
       setoptimizer!,
       setupdater!

# initializer.jl
export AbstractInitializer,
       UniformInitializer,
       NormalInitializer,
       XavierInitializer

# optimizer.jl
export AbstractOptimizer,
       AdaDelta,
       AdaGrad,
       ADAM,
       AdaMax,
       Nadam,
       RMSProp,
       SGD,
       getupdater,
       normgrad!,
       update!

# io.jl
export AbstractDataProvider,
       AbstractDataBatch,
       DataBatch,
       ArrayDataProvider,
       ArrayDataBatch

# visualize.jl
export to_graphviz

###############################################################################
#  includes
###############################################################################

include("exceptions.jl")
include("base.jl")

include("runtime.jl")
include("context.jl")
include("util.jl")

include("ndarray.jl")

include("random.jl")
include("autograd.jl")

include("name.jl")
include("symbolic-node.jl")
include("executor.jl")

include("broadcast.jl")

include("metric.jl")
include("optimizer.jl")
include("initializer.jl")

include("io.jl")
include("kvstore.jl")

include("callback.jl")
include("model.jl")

include("visualize.jl")

include("nn-factory.jl")

include("deprecated.jl")

end # mx

@reexport using .mx

end # module MXNet

max_issues_repo_name: sheldon-pan/incubator-mxnet

max_forks_repo_name: sheldon-pan/incubator-mxnet


text:
# multigrid.jl : define multigrid cycle operations

"""
    MultigridMethod(A, sz, cycle_type)
    MultigridMethod(A, sz, cycle_type; kwargs...)
    MultigridMethod(f, sz, cycle_type)
    MultigridMethod(f, sz, cycle_type; kwargs...)

Geometric Multigrid method of type `cycle_type` for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh. A Galerkin approach is used to compose the coarse matrices, unless a function `f` is provided for direct discretization.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
f          : Function, function used for direct-discretization of the PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`
cycle_type : Multigrid cycle type, can be `V(ν₁,ν₂)` (V-cycle), `W(ν₁,ν₂)` (W-cycle) or `F(ν₀,ν₁,ν₂)` (FMG)

Options
=======
* max_iter : maximum number of iterations. Default value is `20` For FMG using the `F()`-cycle, this is set to 1.
* R_op     : restriction operator type, can be `Injection()` or `FullWeighting()` (default)
* P_op     : interpolation operator type, can be `Injection()`, `FullWeighting()` (default), or `Cubic()`
* ngrids   : total number of grids to use, default is `min.(⌊log₂(sz)⌋)`
* smoother : smoother, can be `GaussSeidel()` of `Jacobi()`
"""
MultigridMethod(A::Union{AbstractMatrix,Function}, sz::NTuple, cycle_type::MultigridCycle; max_iter::Int=20, R_op::TransferKind=FullWeighting(), P_op::TransferKind=FullWeighting(), ngrids=factor_twos.(sz), smoother::Smoother=GaussSeidel()) = MultigridIterable(coarsen(A, sz, R_op, P_op, ngrids), max_iter, cycle_type, smoother, Vector{Float64}(undef, 0))

"""
    V_cycle(A, sz)
    V_cycle(A, sz; kwargs...)

Geometric semicoarsened Multigrid V(2,1)-cycle for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`

For other options, see `MultigridMethod`.
"""
V_cycle(A::Union{AbstractMatrix,Function}, sz::NTuple; kwargs...) = MultigridMethod(A, sz, V(); kwargs...)

"""
    W_cycle(A, sz)
    W_cycle(A, sz; kwargs...)

Geometric semicoarsened Multigrid W(2,1)-cycle for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`

For other options, see `MultigridMethod`.
"""
W_cycle(A::Union{AbstractMatrix,Function}, sz::NTuple; kwargs...) = MultigridMethod(A, sz, W(); kwargs...)

"""
    F_cycle(A, sz)
    F_cycle(A, sz; kwargs...)

Geometric semicoarsened Full Multigrid for matrix `A`, that results from a discretization of a PDE on [0,1]^d using an `sz`-point mesh.

Inputs
======
A          : SparseMatrixCSC, sparse matrix with discretized PDE
sz         : NTuple, PDE grid size, e.g., `(n,m)`

For other options, see `MultigridMethod`.
"""
F_cycle(A::Union{AbstractMatrix,Function}, sz::NTuple; kwargs...) = MultigridMethod(A, sz, F(); max_iter=1, kwargs...)

# semicoarsening cycle
filter_at(R, s) = Base.Iterators.filter(i->sum(i.I)==s, R)

filter_fit(I1, Rall, R) = Base.Iterators.filter(i->i[2] ∈ Rall, enumerate(R))

child_iter(Rall, I1, I) = filter_fit(I1, Rall, Base.Iterators.reverse(filter_at(CartesianIndices(UnitRange.(Tuple(I-I1), Tuple(I))), sum(Tuple(I))-1)))

parent_iter(Rall, I1, I) = filter_fit(I1, Rall, filter_at(CartesianIndices(UnitRange.(Tuple(I), Tuple(I+I1))), sum(Tuple(I))+1))

grids_at_level(R::CartesianIndices{d}, s) where d = filter_at(R, s+d-1)

function μ_cycle!(grids::Array{G} where {G<:Grid}, μ::Int, ν₁::Int, ν₂::Int, grid_ptr::Int, smoother::Smoother)
    R = CartesianIndices(size(grids))
    I1, Iend = first(R), last(R)
    for I in grids_at_level(R, grid_ptr)
        smooth!(grids[I], ν₁, smoother)
    end
    if grid_ptr == sum(Tuple(Iend-I1)) + 1
        grids[Iend].x .= grids[Iend].A\grids[Iend].b # exact solve
    else
        for I in grids_at_level(R, grid_ptr+1)
            R_child = child_iter(R, I1, I)
            grids[I].b .= mean(map(i->grids[last(i)].R[first(i)]*residu(grids[last(i)]), R_child))
            fill!(grids[I].x, zero(eltype(grids[I].x)))
        end
        for i in 1:μ
            μ_cycle!(grids, μ, ν₁, ν₂, grid_ptr+1, smoother)
        end
        for I in grids_at_level(R, grid_ptr)
            R_parent = parent_iter(R, I1, I)
            # matrix-dependent prolongation
            λ = map(i->grids[I].A * high_freq_mode(first(i), grids[I].sz), R_parent)
            λ² = broadcast(i->broadcast(j->j^2, i), λ)
            ω = map(i->λ²[i]./sum(λ²), 1:length(λ)) # weight factors from [Naik, Van Rosendale]
            ip = map(i->grids[last(i)].P[first(i)]*grids[last(i)].x, R_parent)
            coarse_grid_correction!(grids, I, sum(map(i->ω[i].*ip[i], 1:length(ω))))
            smooth!(grids[I], ν₂, smoother)
        end
    end
end

high_freq_mode(dir,sz) = view(map(i->-iseven(i[dir])+isodd(i[dir]), CartesianIndices(sz.-1)), :)

function F_cycle!(grids::Array{G} where {G<:Grid}, ν₀::Int, ν₁::Int, ν₂::Int, grid_ptr::Int, smoother::Smoother)
    R = CartesianIndices(size(grids))
    I1, Iend = first(R), last(R)
    if grid_ptr == sum(Tuple(Iend-I1)) + 1
        fill!(grids[grid_ptr].x, zero(eltype(grids[grid_ptr].x)))
    else
        for I in grids_at_level(R, grid_ptr+1)
            R_child = child_iter(R, I1, I)
            grids[I].b .= mean(map(i->grids[last(i)].R[first(i)]*grids[last(i)].b, R_child))
        end
        F_cycle!(grids, ν₀, ν₁, ν₂, grid_ptr+1, smoother)
        for I in grids_at_level(R, grid_ptr)
            R_parent = parent_iter(R, I1, I)
            # matrix-dependent prolongation
            λ = map(i->grids[I].A * high_freq_mode(first(i), grids[I].sz), R_parent)
            λ² = broadcast(i->broadcast(j->j^2, i), λ)
            ω = map(i->λ²[i]./sum(λ²),1:length(λ)) # weight factors from [Naik, Van Rosendale]
            ip = map(i->P̃(first(i), Cubic(), grids[last(i)].sz...) * grids[last(i)].x, R_parent)
            grids[I].x .= sum(map(i->ω[i].*ip[i], 1:length(ω)))
        end
    end
    for i in 1:ν₀
        μ_cycle!(grids, 1, ν₁, ν₂, grid_ptr, smoother)
    end
end


text:
function quick_sort!(arr::Vector{T},l::Int=1,r::Int=length(arr))where T
    if l>=r
        return
    end
    mid=arr[(l+r)>>1]
    l_pos=l
    r_pos=r
    while l_pos<=r_pos
        while arr[l_pos]<mid l_pos+=1 end
        while arr[r_pos]>mid r_pos-=1 end
        if l_pos<=r_pos
            t=arr[l_pos]
            arr[l_pos]=arr[r_pos]
            arr[r_pos]=t
            l_pos+=1
            r_pos-=1
        end
    end
    if l<r_pos quick_sort!(arr,l,r_pos) end
    if l_pos<r quick_sort!(arr,l_pos,r) end
end

max_stars_repo_name: frankschmitt/Julia

max_issues_repo_name: frankschmitt/Julia

max_forks_repo_name: frankschmitt/Julia


text:
using MyExampleInJulia
using Documenter

makedocs(;
    modules=[MyExampleInJulia],
    authors="Alberto Cuoci, Politecnico di Milano",
    repo="https://github.com/acuoci/MyExampleInJulia.jl/blob/{commit}{path}#L{line}",
    sitename="MyExampleInJulia.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://acuoci.github.io/MyExampleInJulia.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/acuoci/MyExampleInJulia.jl",
)


text:
using JLD2, LinearAlgebra, Statistics

# prefix = "fluxernathy_tracers_restarted_j1_k1_averages.jld2"
# prefix = "fluxernathy_nh_ns_averages.jld2"
# prefix = "fluxernathy_averages.jld2"
# prefix = "fluxernathy_tracers_restarted_j1_k10_averages.jld2"
# prefix = "fluxernathy_tracers_restarted_j1_k20_averages.jld2"
prefix = "fluxernathy_tracers_restarted_j128_k1_averages.jld2"
# prefix = "check_h200.jld2"
include("utils.jl")
jl_file = jldopen(prefix, "r+")

time_index = 22
b = get_field("b", time_index, jl_file)
y,z = get_grid(b, jl_file; ghost = 3)

y̅ = avg(y)
z̅ = avg(z)
bz_p = ∂z(b, z)
by_p = ∂y(b, y)

bz = avgy(bz_p)
by = avgz(by_p)

v = get_field("v", time_index, jl_file)
w = get_field("w", time_index, jl_file)

# tracer dependent
cs   = []
cys  = []
czs  = []
vcps = []
wcps = []

∇c∇ᵖbs  = []  # isopycnal tracer gradients ∇c ⋅ ∇ᵖb
∇c∇bs  = []  # diapycnal tracer gradients ∇c ⋅ ∇b
u⃗c∇ᵖbs = [] # isopycnal eddy flux u⃗'c' ⋅ ∇ᵖb
u⃗c∇bs = [] # diapycnal eddy flux u⃗'c' ⋅ ∇b

for tracer_string in ["c1", "c2", "c3", "c4"]
    local c = get_field(tracer_string, time_index, jl_file)
    local vc = get_field("v"*tracer_string, time_index, jl_file)
    local wc = get_field("w"*tracer_string, time_index, jl_file)

    local vcp = avgy(avgz( avgy(vc) - avgy(v) .* c ))
    local wcp = avgy(avgz( avgz(wc) - avgz(w) .* c ))

    local cz = avgy(∂z(c, z))
    local cy = avgz(∂y(c, y))

    local vcp_avg = avgz(avgy(vcp))
    local wcp_avg = avgz(avgy(wcp))

    local ∇c∇b  = cy .* by + cz .* bz
    local ∇c∇ᵖb = cy .* bz - cz .* by

    local u⃗c∇b  = vcp .* by + wcp .* bz
    local u⃗c∇ᵖb = vcp .* bz - wcp .* by

    push!(cs, c)
    push!(cys, cy)
    push!(czs, cz)
    push!(vcps, vcp_avg)
    push!(wcps, wcp_avg)

    push!(∇c∇bs, ∇c∇b)
    push!(∇c∇ᵖbs, ∇c∇ᵖb)
    push!(u⃗c∇bs, u⃗c∇b)
    push!(u⃗c∇ᵖbs, u⃗c∇ᵖb)
end

ti = 1 # tracer index
mask = zeros(size(∇c∇bs[1]))
mask[40:255-40, 10:30] .= 1.0
mask[100:150, :] .= 0.0
# mask[:,:] .= 1.0

∇c∇b  = ∇c∇bs[ti]  .* mask
∇c∇ᵖb = ∇c∇ᵖbs[ti] .* mask

u⃗c∇b  = u⃗c∇bs[ti] .* mask
u⃗c∇ᵖb = u⃗c∇ᵖbs[ti] .* mask

A = [mean(∇c∇b .* ∇c∇b) mean(∇c∇b .* ∇c∇ᵖb); mean(∇c∇b .* ∇c∇ᵖb) mean(∇c∇ᵖb .* ∇c∇ᵖb)]
b⃗ =  -1.0 .* [mean(∇c∇b .* u⃗c∇b) ; mean(∇c∇ᵖb .* u⃗c∇b)]
κ¹¹, κ¹² = A \ b⃗

b⃗ = -1.0 .* [mean(∇c∇b .* u⃗c∇ᵖb) ; mean(∇c∇ᵖb .* u⃗c∇ᵖb)]
κ²¹, κ²² = A \ b⃗

K = [κ¹¹ κ¹²; κ²¹ κ²²]
S = (K + K') * 0.5
A = (K - K') * 0.5

println("The eigenvalues of K are ", eigvals(K))

Ld = sqrt(mean( (u⃗c∇b .+  κ¹¹ .*  ∇c∇b .+  κ¹² .* ∇c∇ᵖb) .^2 )) / sqrt( mean( u⃗c∇b .^2 ) )
Li = sqrt(mean( (u⃗c∇ᵖb .+  κ²¹ .*  ∇c∇b .+  κ²² .* ∇c∇ᵖb) .^2 )) / sqrt( mean(u⃗c∇ᵖb .^2 ) )

lq, uq = [0.05, 0.95]
residuald = (u⃗c∇b .+  κ¹¹ .*  ∇c∇b .+  κ¹² .* ∇c∇ᵖb) ./ sqrt( mean( u⃗c∇b .^2 ) )
residuali = (u⃗c∇ᵖb .+  κ²¹ .*  ∇c∇b .+  κ²² .* ∇c∇ᵖb) ./ sqrt( mean(u⃗c∇ᵖb .^2 ) )
climsi = quantile.(Ref(residuali[:]), [lq, uq])
fig, ax, hm = heatmap(y̅, z̅, residuali, colormap = :balance, colorrange = climsi, interpolate = true)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using DataDeps
using Base.Test

@testset "Downloading and unpacking" begin
  @test DataDeps.add("https://www.cs.toronto.edu/~kriz/cifar-10-binary.tar.gz", "test") == 0
end

@testset "List" begin
  @test DataDeps.list()=="test\n"
end

x=0; y=0
@testset "Segregate" begin
  @test DataDeps.setseries("data_batch_*.bin", "series", 1, 5)==0
  global x,y
  x,y=DataDeps.traindata("test")
end

println("X: ")
println(x)
println("Y: ")
println(y)


text:
module Solvers

using Senuk.PrintUtils.print_iter_head
using Senuk.PrintUtils.print_iter
using Senuk.PrintUtils.print_blank_line

function gd(δf :: Function, X :: Matrix{Float64}, y :: Vector{Float64},
            step_size :: Float64, max_iter :: Int64, verbose :: Bool)
  w = zeros(size(X)[2])
  print_steps = 50

  print_iter_head()
  for i in 1:max_iter
    cost, Δcost = δf(w, X, y)
    if verbose && i % print_steps == 0
      print_iter(i, cost)
    end
    w = w - step_size*Δcost
  end
  print_blank_line()
  return w
end

end

max_stars_repo_name: upul/Senuk

max_issues_repo_name: upul/Senuk

max_forks_repo_name: upul/Senuk


text:
using Nemo, GibbsTypePriors, JLD, DataFramesMeta

include("common_functions_for_tests.jl")

grid_n = [collect(10:100:4000); 4000]
grid_s = GibbsTypePriors.RR.(0.1:0.2:0.9)

accuracy_Cnk_sigma = expand_grid(grid_n, grid_s) |>
    df -> rename(df, [:n, :s]) |>
    df -> vcat(@transform(df, k = :n), @transform(df, k = repeat([1], inner = length(:n)))) |>
    df -> @transform(df, Accuracy = GibbsTypePriors.Cnk_rec.(:n, :k, :s) |> x -> accuracy_bits.(x))
save("test/graphical_tests/saves_for_graphical_tests/accuracy_Cnk_rec_sigma.jld", "accuracy_Cnk_sigma", accuracy_Cnk_sigma)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
__precompile__()
module DiscontinuousGalerkin
  using Polynomials
  using FastGaussQuadrature
  using LsqFit
  using Parameters, Reexport
  using DiffEqBase
  @reexport using OrdinaryDiffEq
  using RecipesBase

  import DiffEqBase: solve, @def
  import Base: show

  abstract type AbstractConservationLawProblem{islinear,isstochastic,MeshType} <: DEProblem end
  abstract type AbstractFEAlgorithm <: DEAlgorithm end
  abstract type AbstractDGLimiter end
  abstract type AbstractFVMesh end

  include("mesh.jl")
  include("problem.jl")
  include("riemann_solvers.jl")
  include("basis.jl")
  include("solve.jl")
  include("solution.jl")
  include("plots.jl")
  include("limiter.jl")

  export Uniform1DFVMesh
  export legendre_basis
  export cell_faces
  export cell_centers, cell_volume, cell_indices, numcells
  export ConservationLawsProblem
  export advection_num_flux, rusanov_euler_num_flux, glf_num_flux
  export legendre_basis, PolynomialBasis, poly_jacobi, poly_legendre, reference_to_interval
  export DGLimiter, Linear_MUSCL_Limiter, WENO_Limiter
  export solve, DiscontinuousGalerkinScheme
  export fluxρ, myblock, apply_boundary
end

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
using ZipFile

function getAhornModDirs()
    if !get(config, "load_plugins_ahorn", true)
        return String[]
    end

    modsPath = joinpath(storageDirectory, "plugins")
    targetFolders = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if isdir(joinpath(modsPath, fn))
                push!(targetFolders, joinpath(modsPath, fn))
            end
        end
    end

    return targetFolders
end

function getAhornModZips()
    if !get(config, "load_plugins_ahorn_zip", true)
        return String[]
    end

    modsPath = joinpath(storageDirectory, "plugins")
    targetFolders = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if isfile(joinpath(modsPath, fn)) && hasExt(fn, ".zip")
                push!(targetZips, joinpath(modsPath, fn))
            end
        end
    end

    return targetFolders
end

function getCelesteModDirs()
    if !get(config, "load_plugins_celeste", true)
        return String[]
    end

    celesteDir = get(config, "celeste_dir", "")
    modsPath = joinpath(celesteDir, "Mods")

    targetFolders = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if isdir(joinpath(modsPath, fn))
                push!(targetFolders, joinpath(modsPath, fn))
            end
        end
    end

    return targetFolders
end

function getModRoot(fn::String)
    celesteDir = get(config, "celeste_dir", "")
    modsPath = normpath(joinpath(celesteDir, "Mods"))

    path = normpath(fn)

    while true
        parent = dirname(path)

        if parent == modsPath
            if isdir(path)
                return true, path

            else
                return true, parent
            end
        end

        if path == parent
            return false, ""
        end

        path = parent
    end
end

function getCelesteModZips()
    if !get(config, "load_plugins_celeste_zip", true)
        return String[]
    end

    celesteDir = get(config, "celeste_dir", "")
    modsPath = joinpath(celesteDir, "Mods")

    targetZips = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if isfile(joinpath(modsPath, fn)) && hasExt(fn, ".zip")
                push!(targetZips, joinpath(modsPath, fn))
            end
        end
    end

    return targetZips
end

function findExternalSprite(resource::String, atlas::String="Gameplay")
    targetFolders = getCelesteModDirs()
    targetZips = getCelesteModZips()
    atlasPath = joinpath("Graphics", "Atlases", atlas)

    for target in targetFolders
        fn = joinpath(target, atlasPath, resource * ".png")
        if isfile(fn)
            return fn
        end
    end

    # ZipFile uses unix paths
    gameplayUnixPath = replace(atlasPath, "\\" => "/")

    for target in targetZips
        fn = gameplayUnixPath * "/" * resource * ".png"
        zipfh = ZipFile.Reader(target)

        for file in zipfh.files
            if file.name == fn
                return target
            end
        end

        close(zipfh)
    end
end

const warnedFilenames = Set{String}()

function warnBadExtensions(fixable::Array{String, 1}, nonFixable::Array{String, 1}, correctExt::String)
    res = false

    fixable = String[fn for fn in fixable if !(fn in warnedFilenames)]
    nonFixable = String[fn for fn in nonFixable if !(fn in warnedFilenames)]

    if !isempty(fixable)
        fixableCount = length(fixable)
        plural = fixableCount == 1 ? "" : "s"
        dialogText = "You have $fixableCount image$plural with bad extension$plural.\n" *
            "These are not loadable by Celeste (Everest).\n" *
            "Do you want to rename them?\n\n" *
            join(fixable, "\n")
        confirmed = ask_dialog(dialogText, Ahorn.window)
        res = confirmed

        if confirmed
            for fn in fixable
                path, ext = splitext(fn)
                out = path * correctExt
                tmp = out * ".tmp"

                # Workaround for Windows...
                mv(fn, tmp)
                mv(tmp, out)
            end

        else
            push!(warnedFilenames, fixable...)
        end
    end

    if !isempty(nonFixable)
        nonFixableCount = length(nonFixable)
        plural = nonFixableCount == 1 ? "" : "s"
        dialogText = "You have $nonFixableCount file$plural with bad extension$plural that cannot be automatically fixed.\n" *
            "These are not loadable by Celeste (Everest) or Ahorn.\n" *
            "Please rename the file$plural manually with the proper extension '$correctExt' or contact the mod maker.\n\n" *
            join(nonFixable, "\n")
        info_dialog(dialogText, Ahorn.window)

        push!(warnedFilenames, nonFixable...)
    end

    return res
end

const spriteZipCache = Dict{String, Any}()

# Resource prefix filters uneeded files, uses unix paths
function cacheZipContent(filename::String, resourcePrefix::String="Graphics/Atlases/Gameplay/")
    debug.log("Caching zip content for $filename", "ZIP_CACHING_VERBOSE")

    spriteZipCache[filename] = Dict{String, Any}()

    # Always uses Unix path
    zipfh = ZipFile.Reader(filename)

    for file in zipfh.files
        if startswith(file.name, resourcePrefix)
            parts = split(file.name, "/")
            atlas = parts[3]
            resource = join(parts[3:end], "/")

            if hasExt(file.name, ".png")
                spriteZipCache[filename][resource] = Cairo.read_from_png(file)

            elseif hasExt(file.name, ".meta.yaml")
                try
                    spriteZipCache[filename][resource] = YAML.load(String(read(file)))

                catch
                    # Bad yaml file
                end
            end
        end
    end

    close(zipfh)
end

# Cleanup and destroy surfaces
function uncacheZipContent(filename::String)
    debug.log("Uncaching zip content for $filename", "ZIP_CACHING_VERBOSE")

    fileCache = get(spriteZipCache, filename, nothing)

    if fileCache !== nothing
        for (resource, data) in fileCache
            if hasExt(resource, ".png")
                deleteSurface(data)
            end
        end
    end
end

function findExternalSpritesInZip(filename::String, atlasesPath::String, atlasesUnixPath::String, badExt::Array{String, 1}, allowRetry::Bool=true)
    res = Tuple{String, String, String}[]
    fileCache = get(spriteZipCache, filename, nothing)

    if fileCache !== nothing
        for (resource, data) in fileCache
            if hasExt(resource, ".png")
                if !hasExt(resource, ".png", false)
                    # Case sensitive check for warnings
                    push!(badExt, resource)
                end

                parts = split(resource, "/")
                atlas = parts[1]
                resourceNoAtlas = join(parts[2:end], "/")

                push!(res, (resourceNoAtlas, filename, atlas))
            end
        end

    else
        if allowRetry
            cacheZipContent(filename)

            return findExternalSpritesInZip(filename, atlasesPath, atlasesUnixPath, badExt, false)

        else
            return res
        end
    end

    return res
end

function findExternalSpritesInFolder(filename::String, atlasesPath::String, badExt::Array{String, 1}, targetAtlases::Array{String, 1}=["Gameplay"])
    res = Tuple{String, String, String}[]
    atlasPath = joinpath(filename, atlasesPath)

    if !isdir(atlasPath)
        return res
    end

    for atlas in readdir(atlasPath)
        if !(atlas in targetAtlases)
            continue
        end

        path = joinpath(atlasPath, atlas)

        if isdir(path)
            for (root, dir, files) in walkdir(path)
                for file in files
                    if hasExt(file, ".png")
                        if !hasExt(file, ".png", false)
                            # Case sensitive check for warnings
                            push!(badExt, joinpath(root, file))
                        end

                        rawpath = joinpath(root, file)
                        push!(res, (relpath(rawpath, path), rawpath, atlas))
                    end
                end
            end
        end
    end

    return res
end

function findAllExternalSprites()
    warnOnBadExts = get(config, "prompt_bad_resource_exts", true)
    res = Tuple{String, String, String}[]

    folderFileBadExt = String[]
    zipFileBadExt = String[]

    targetFolders = getCelesteModDirs()
    targetZips = getCelesteModZips()
    atlasesPath = joinpath("Graphics", "Atlases")

    # ZipFile uses unix paths
    atlasesUnixPath = replace(atlasesPath, "\\" => "/")

    for target in targetFolders
        append!(res, findExternalSpritesInFolder(target, atlasesPath, folderFileBadExt))
    end

    for target in targetZips
        append!(res, findExternalSpritesInZip(target, atlasesPath, atlasesUnixPath, zipFileBadExt))
    end

    if warnOnBadExts
        confirmed = warnBadExtensions(folderFileBadExt, zipFileBadExt, ".png")

        if confirmed
            return findAllExternalSprites()
        end
    end

    return res
end

function findChangedExternalSprites()
    warnOnBadExts = get(config, "prompt_bad_resource_exts", true)
    res = Tuple{String, String, String}[]

    changedFilenames = FileWatcher.processWatchEvents(FileWatcher.basePath)

    folderFileBadExt = String[]
    zipFileBadExt = String[]

    # ZipFile uses unix paths
    atlasesPath = joinpath("Graphics", "Atlases")
    atlasesUnixPath = replace(atlasesPath, "\\" => "/")

    for filename in changedFilenames
        if hasExt(filename, ".zip")
            append!(res, findExternalSpritesInZip(filename, atlasesPath, atlasesUnixPath, zipFileBadExt))

        elseif hasExt(filename, ".png")
            hasModRoot, modRoot = getModRoot(filename)

            if hasModRoot
                relative = relpath(filename, modRoot)
                parts = splitpath(relative)

                if parts[1] == "Graphics" && parts[2] == "Atlases"
                    push!(res, (joinpath(parts[4:end]...), filename, parts[3]))
                end
            end
        end
    end

    if warnOnBadExts
        confirmed = warnBadExtensions(folderFileBadExt, zipFileBadExt, ".png")

        if confirmed
            return findChangedExternalSprites()
        end
    end

    return res
end

function loadLangdata()
    targetFolders = joinpath.(getCelesteModDirs(), "Ahorn")
    targetZips = getCelesteModZips()

    lang = get(config, "language", "en_gb")
    internalFilename = abspath("../lang/$lang.lang")

    res = parseLangfile(String(read(internalFilename)))

    for folder in targetFolders
        fn = joinpath(folder, "lang", "$(lang).lang")
        if ispath(fn)
            content = String(read(fn))
            parseLangfile(content, init=res)
        end
    end

    for fn in targetZips
        zipfh = ZipFile.Reader(fn)

        for file in zipfh.files
            name = file.name
            path, fnext = splitext(name)

            if path == "Ahorn/lang/$(lang)" && hasExt(name, ".lang")
                content = String(read(file))
                parseLangfile(content, init=res)
            end
        end

        close(zipfh)
    end

    return res
end

function findExternalModules(args::String...)
    res = String[]

    targetFolders = vcat(
        joinpath.(getCelesteModDirs(), "Ahorn"),
        getAhornModDirs()
    )

    for folder in targetFolders
        path = joinpath(folder, args...)
        if isdir(path)
            for file in readdir(path)
                if hasExt(file, ".jl")
                    push!(res, joinpath(path, file))
                end
            end
        end
    end

    return res
end

function loadExternalModules!(loadedModules::Dict{String, Module}, loadedNames::Array{String, 1}, args::String...)
    # ZipFile uses unix paths
    path = joinpath("Ahorn", args...)
    path = replace(path, "\\" => "/")

    targets = vcat(
        getCelesteModZips(),
        getAhornModZips()
    )

    for target in targets
        zipfh = ZipFile.Reader(target)

        for file in zipfh.files
            name = file.name

            if startswith(name, path) && file.uncompressedsize > 0
                # Add .from_zip to the filename
                # Prevents hotswaping from trying to access invalid file
                fakeFn = name * ".from_zip"

                try
                    content = String(read(file))
                    loadedModules[fakeFn] = eval(Meta.parse(content))

                    if !(fakeFn in loadedNames)
                        push!(loadedNames, fakeFn)
                    end

                catch e
                    println(Base.stderr, "! Failed to load \"$name\" from \"$target\"")
                    println(Base.stderr, e)
                end
            end
        end

        close(zipfh)
    end
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
@testset "Testing quadsums.jl" begin
    rough_kwargs = Dict(:atol => 1e-10, :rtol => 1e-10)

    @testset "test simple var sum" begin
        beta = .95
        A = 1.
        C = 0.
        H = 1.
        x0 = 1.

        val = var_quadratic_sum(A, C, H, beta, x0)

        @test isapprox(val, 20.0; rough_kwargs...)
    end

    @testset "test identity var sum" begin
        beta = .95
        A = eye(3)
        C = zeros(3, 3)
        H = eye(3)
        x0 = ones(3)

        val = var_quadratic_sum(A, C, H, beta, x0)

        @test isapprox(val, 60.0; rough_kwargs...)
    end

end  # facts


text:
@doc raw"""
    dominoes(imgname, colr, constrnt)

Approximate a grayscale image with complete sets of double-nine dominoes.

# Resources
Watch Bob Bosch of Oberlin College talk at Google about domino mosaics, and his interpretation of this art form as an integer program.
https://www.youtube.com/watch?v=g3CiVrN-BnY

IP modeling details can be found in this article.
http://www.optimization-online.org/DB_FILE/2003/09/722.pdf

# Examples
```julia-repl
julia> include("ModuleDominoes.jl")
julia> using .Dominoes
julia> createMosaic("Frankenstein", "white", "constrained")
```

# Author
Kalyani Nagaraj
March 2021
"""
function dominoes(imgname, colr, constrnt)

   ## READ DATA
   g = readdlm("../../Data/Compressed_Images/" * imgname * ".txt", ' ', Int)

   if cmp(colr, "black") == 0 # Flip the scale, i.e., map 0-9 to 9-0,
      g = 9 .- g;                 # if final image is to be generated
   end                           # from white tiles with black dots.

   R = size(g,1)         # Number of rows of squares in the domino image
   C = size(g,2)         # Number of columns of squares in the domino image
   s = Int(R*C/110)      # Number of sets of dominoes

   ## MODEL THE OPTIMIZATION PROBLEM
   # I. Define a model and the optimizer it will use
   dominoes = Model(GLPK.Optimizer)
   # set_optimizer_attribute(dominoes, "msg_lev", GLPK.GLP_MSG_ON)

   # II. Create binary variables
   # x_mnoij = 1 if domino (m,n), where m<=n,  has oriented o starting in cell
   # (i,j)
   # number of variables = 55*4*R*C
   @variable(dominoes, x[m=0:9, n=m:9, o=1:4, i=1:R, j=1:C], Bin);
   # undo_relax = relax_integrality(dominoes); # Returns a function that can be
                                               # called without any arguments to
                                               # restore the original model.
                                               # The behavior of this function is
                                               # undefined if additional changes are
                                               # made to the affected variables in
                                               # the meantime.

   # III. Define the objective function
   @objective(dominoes, Min, sum(( (m-g[i,j])^2 + (n-g[i+1,j])^2 ) * x[m,n,1,i,j]  for m=0:9, n=m:9, i=1:R-1, j=1:C)
                           + sum(( (m-g[i+1,j])^2 + (n-g[i,j])^2 ) * x[m,n,2,i,j]  for m=0:9, n=m:9, i=1:R-1, j=1:C)
                           + sum(( (m-g[i,j])^2 + (n-g[i,j+1])^2 ) * x[m,n,3,i,j]  for m=0:9, n=m:9, i=1:R, j=1:C-1)
                           + sum(( (m-g[i,j+1])^2 + (n-g[i,j])^2 ) * x[m,n,4,i,j]  for m=0:9, n=m:9, i=1:R, j=1:C-1))

   # IV. Now, the constraints:
   # A. Type-1 constraint: Domino (m,n) can be used exactly s times, if using exactly s complete sets of double-none dominoes.
   if constrnt == "constrained"
      for m = 0:9, n = m:9
         # 55 type-1 constraints
         @constraint(dominoes, sum(x[m,n,o,i,j] for o=1:4, i=1:R, j=1:C) == s)
      end
   elseif constrnt == "unconstrained"
      # a total of 55s domino tiles are used
      @constraint(dominoes, sum(x[m,n,o,i,j] for m=0:9, n=m:9, o=1:4, i=1:R, j=1:C) == 55*s)
   end

   # B. Type-2 constraints: Each square in the image is covered by one domino
   # either placed horizontally or vertically

   # All but first row and first column:
   for i = 2:R, j = 2:C # a total of 110s type-2 constraints
      @constraint(dominoes, sum(x[m,m,1,i,j] for m=0:9) +
                            sum(x[m,m,1,i-1,j] for m=0:9) +
                            sum(x[m,m,3,i,j] for m=0:9) +
                            sum(x[m,m,3,i,j-1] for m=0:9) +

                            sum(x[m,n,1,i,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,2,i,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,1,i-1,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,2,i-1,j] for m=0:9, n=m+1:9) +

                            sum(x[m,n,3,i,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,4,i,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,3,i,j-1] for m=0:9, n=m+1:9) +
                            sum(x[m,n,4,i,j-1] for m=0:9, n=m+1:9) == 1)
   end

   # top row cells (i=1), all except top left cell
   for j = 2:C
      # C-1 constraints
      @constraint(dominoes, sum(x[m,m,1,1,j] for m=0:9) +
                            sum(x[m,m,3,1,j] for m=0:9) +
                            sum(x[m,m,3,1,j-1] for m=0:9) +

                            sum(x[m,n,1,1,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,2,1,j] for m=0:9, n=m+1:9) +

                            sum(x[m,n,3,1,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,4,1,j] for m=0:9, n=m+1:9) +
                            sum(x[m,n,3,1,j-1] for m=0:9, n=m+1:9) +
                            sum(x[m,n,4,1,j-1] for m=0:9, n=m+1:9) == 1)
   end

   # first column (j=1), all except top left cell
   for i = 2:R
      # R-1 constraints
      @constraint(dominoes, sum(x[m,m,1,i,1] for m=0:9) +
                            sum(x[m,m,1,i-1,1] for m=0:9) +
                            sum(x[m,m,3,i,1] for m=0:9) +

                            sum(x[m,n,1,i,1] for m=0:9, n=m+1:9) +
                            sum(x[m,n,2,i,1] for m=0:9, n=m+1:9) +
                            sum(x[m,n,1,i-1,1] for m=0:9, n=m+1:9) +
                            sum(x[m,n,2,i-1,1] for m=0:9, n=m+1:9) +

                            sum(x[m,n,3,i,1] for m=0:9, n=m+1:9) +
                            sum(x[m,n,4,i,1] for m=0:9, n=m+1:9) == 1)
   end

   # 1 constraint corresponding to top left cell
   @constraint(dominoes, sum(x[m,m,1,1,1] for m=0:9) +
                         sum(x[m,m,3,1,1] for m=0:9) +

                         sum(x[m,n,1,1,1] for m=0:9, n=m+1:9) +
                         sum(x[m,n,2,1,1] for m=0:9, n=m+1:9) +

                         sum(x[m,n,3,1,1] for m=0:9, n=m+1:9) +
                         sum(x[m,n,4,1,1] for m=0:9, n=m+1:9) == 1)
   # print(dominoes)

   f = open("../../Data/Optimization_Output/Model_files/" * imgname * "_" * colr * "_" * constrnt * ".txt", "w")
   println("This optimization model has ", num_variables(dominoes), " variables, and ", 55 + R*C, " constraints.")

   print(f, dominoes)
   close(f)

   ## SOLVE THE OPTIMIZATION PROBLEM
   # V. And finally, solve problem
   optimize!(dominoes)
   if primal_status(dominoes) != MOI.FEASIBLE_POINT
      print("No feasible solution found!\n")
   else

      println("Feasible solution found after ", solve_time(dominoes), " seconds!")
      x_val = value.(x)
      println("Optimal objective value = ", objective_value(dominoes));
      println("Average difference in cell value = ", sqrt(objective_value(dominoes)/(110*s)));

      display(x_val)
      ## SAVE SOLUTION
      # Create a matrix to store the solution
      sol = zeros(Int,R,C)  # 9x9 matrix of integers
      for m=0:9, n=m:9, o=1:4, i=1:R, j=1:C
         # Integer programs are solved as a series of linear programs
         # so the values might not be precisely 0 and 1. We can just
         # round them to the nearest integer to make it easier
         if round(Int,x_val[m,n,o,i,j]) == 1
            if o == 1
               sol[i,j] = m
               sol[i+1,j] = n
            elseif o == 2
               sol[i,j] = n
               sol[i+1,j] = m
            elseif o == 3
               sol[i,j] = m
               sol[i,j+1] = n
            elseif o == 4
               sol[i,j] = n
               sol[i,j+1] = m
            end
         end
      end

      open("../../Data/Optimization_Output/Images/" * imgname * "_" * colr * "_" * constrnt * ".txt", "w") do io
          writedlm(io, sol, ' ')
      end

      opt_x_ind = [[m, n, o, i, j] for i=1:R for j=1:C for m=0:9 for n=m:9 for o=1:4 if round(Int, x_val[m,n,o,i,j]) == 1]
      opt_x_ind = reshape( collect(Iterators.flatten(opt_x_ind)), (5, :) )'
      open("../../Data/Optimization_Output/Indices/" * imgname * "_" * colr * "_" * constrnt * ".txt", "w") do io
         writedlm(io, [R C], ' ')
         writedlm(io, opt_x_ind, ' ')
      end

   end

end

max_stars_repo_name: kalyaninagaraj/Domino-Mosaics

max_issues_repo_name: kalyaninagaraj/Domino-Mosaics

max_forks_repo_name: kalyaninagaraj/Domino-Mosaics


text:
@doc raw"""
```
display_3D(img::Array{RGB, 3}, layer::Int=1, suppress_feedback::Bool=false)
```

Function for displaying 3D images in a user-friendly way using Makie.

# Arguments

The function arguments are described in more detail below.

##  `img`

An `AbstractArray` containing the 3d image to be viewed.

##  `layer`

An `Int` which controls the layer that is displayed initially. If not defined
then the first layer will be shown.

##  `suppress_feedback`

A `Bool` keyword argument which controls if console feedback is to be returned.

# Controls

Right mouse button: move image.
Scroll: zoom.
Left arrow: change layer up.
Right arrow: change layer down.

# Example

Create and view synthetic 3D image.
```julia
using Images, PictureSegmentation

#create image of cone with multiple colours
img = Array{RGB{Float64}, 3}(undef, 100, 100, 100)
for i in CartesianIndices(img)
    if (i[1] - 50)^2 + (i[2] - 50)^2 < (i[3] / 4)^2
        img[i] = RGB(i[3] / 100, 0, 1 - (i[3] / 100))
    else
        img[i] = RGB(0,0,0)
    end
end

#display
display_3D(img)
```
"""
function display_3D(img::AbstractArray, layer::Int=1, suppress_feedback::Bool=false)
    scene = Scene()
    image!(scene,img[:,:,layer], show_axis = false)
    display(scene)
    #change layer
    on(scene.events.keyboardbuttons) do buttons
        if ispressed(scene, Keyboard.right)
            if layer != length(img[1,1,:])
                layer += 1
                image!(scene,img[:,:,layer], show_axis = false)
                display(scene)
                if suppress_feedback == false
                    println("layer: ", layer)
                end
            end
        end
        if ispressed(scene, Keyboard.left)
            if layer != 1
                layer -= 1
                image!(scene,img[:,:,layer], show_axis = false)
                display(scene)
                if suppress_feedback == false
                    println("layer: ", layer)
                end
            end
        end
        return
    end
end

@doc raw"""
```
img = convert_labels(labels::Array{Int})
img = convert_labels(labels::Array{Int}, img::AbstractArray)
```

Converts a label array of type `Int` into an RGB image with each region having a
unique colour. If the original image is passed in then the region will be assigned
the average colour of its pixels.

# Arguments

The function arguments are described in more detail below.

#  `labels`

An array of `Int` containing the assigned label for each pixel.

##  `img`

An `AbstractArray` containing the original image.
"""
function convert_labels(labels::Array{Int})
    img = zeros(RGB,axes(labels))
    regions=maximum(labels)
    for i in CartesianIndices(img)
        cVal = π*labels[i]/regions
        img[i] = RGB(sin(cVal)/2+0.5, sin(cVal+π/2)/2+0.5, sin(cVal-π/2)/2+0.5)
    end
    return img
end

function convert_labels(labels::Array{Int}, img::AbstractArray)
    img2 = zeros(RGB{Float64},axes(labels))
    regions = maximum(labels)
    colours = zeros(RGB{Float64},regions)
    counts = zeros(Int,regions)
    for i in CartesianIndices(img)
        if labels[i] != 0
            colours[labels[i]] = colours[labels[i]] + img[i]
            counts[labels[i]] += 1
        end
    end
    for i in eachindex(colours)
        colours[i] = colours[i]/counts[i]
    end
    for i in CartesianIndices(img2)
        if labels[i] != 0
            img2[i] = colours[labels[i]]
        else
            img2[i] = RGB(0.0, 0.0, 0.0)
        end
    end
    return img2
end

@doc raw"""
```
Contours = make_contour(label::Array{Int})
```

Creates a binary mask for the contours of the segmentation of an image using
an array of labels.

# Arguments

The function arguments are described in more detail below.

#  `labels`

An array of `Int` containing each pixel's assigned label.
"""
function make_contour(label::Array{Int})
    contour_mask = zeros(Int, size(label))
    r, c = size(label)
    for i in CartesianIndices(label)
        rₘ, cₘ = i.I
        r₀ = max(rₘ-1, 1)
        r₁ = min(rₘ+1, r)
        c₀ = max(cₘ-1, 1)
        c₁ = min(cₘ+1, c)
        if label[i] != label[r₀, cₘ]
            contour_mask[i] = 1
        elseif label[i] != label[r₁, cₘ]
            contour_mask[i] = 1
        elseif label[i] != label[rₘ, c₀]
            contour_mask[i] = 1
        elseif label[i] != label[rₘ, c₁]
            contour_mask[i] = 1
        end
    end
    return contour_mask
end


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

## thread/task locking abstraction

@noinline function concurrency_violation()
    # can be useful for debugging
    #try; error(); catch; ccall(:jlbacktrace, Cvoid, ()); end
    throw(ConcurrencyViolationError("lock must be held"))
end

"""
    AbstractLock

Abstract supertype describing types that
implement the synchronization primitives:
[`lock`](@ref), [`trylock`](@ref), [`unlock`](@ref), and [`islocked`](@ref).
"""
abstract type AbstractLock end
function lock end
function unlock end
function trylock end
function islocked end
unlockall(l::AbstractLock) = unlock(l) # internal function for implementing `wait`
relockall(l::AbstractLock, token::Nothing) = lock(l) # internal function for implementing `wait`
assert_havelock(l::AbstractLock, tid::Integer) =
    (islocked(l) && tid == Threads.threadid()) ? nothing : concurrency_violation()
assert_havelock(l::AbstractLock, tid::Task) =
    (islocked(l) && tid === current_task()) ? nothing : concurrency_violation()
assert_havelock(l::AbstractLock, tid::Nothing) = concurrency_violation()

"""
    AlwaysLockedST

This struct does not implement a real lock, but instead
pretends to be always locked on the original thread it was allocated on,
and simply ignores all other interactions.
It also does not synchronize tasks; for that use a real lock such as [`RecursiveLock`](@ref).
This can be used in the place of a real lock to, instead, simply and cheaply assert
that the operation is only occurring on a single cooperatively-scheduled thread.
It is thus functionally equivalent to allocating a real, recursive, task-unaware lock
immediately calling `lock` on it, and then never calling a matching `unlock`,
except that calling `lock` from another thread will throw a concurrency violation exception.
"""
struct AlwaysLockedST <: AbstractLock
    ownertid::Int16
    AlwaysLockedST() = new(Threads.threadid())
end
assert_havelock(l::AlwaysLockedST) = assert_havelock(l, l.ownertid)
lock(l::AlwaysLockedST) = assert_havelock(l)
unlock(l::AlwaysLockedST) = assert_havelock(l)
trylock(l::AlwaysLockedST) = l.ownertid == Threads.threadid()
islocked(::AlwaysLockedST) = true

## condition variables

"""
    GenericCondition

Abstract implementation of a condition object
for synchronizing tasks objects with a given lock.
"""
struct GenericCondition{L<:AbstractLock}
    waitq::InvasiveLinkedList{Task}
    lock::L

    GenericCondition{L}() where {L<:AbstractLock} = new{L}(InvasiveLinkedList{Task}(), L())
    GenericCondition{L}(l::L) where {L<:AbstractLock} = new{L}(InvasiveLinkedList{Task}(), l)
    GenericCondition(l::AbstractLock) = new{typeof(l)}(InvasiveLinkedList{Task}(), l)
end

assert_havelock(c::GenericCondition) = assert_havelock(c.lock)
lock(c::GenericCondition) = lock(c.lock)
unlock(c::GenericCondition) = unlock(c.lock)
trylock(c::GenericCondition) = trylock(c.lock)
islocked(c::GenericCondition) = islocked(c.lock)

lock(f, c::GenericCondition) = lock(f, c.lock)

# have waiter wait for c
function _wait2(c::GenericCondition, waiter::Task)
    ct = current_task()
    assert_havelock(c)
    push!(c.waitq, waiter)
    # since _wait2 is similar to schedule, we should observe the sticky bit now
    if waiter.sticky && Threads.threadid(waiter) == 0
        # Issue #41324
        # t.sticky && tid == 0 is a task that needs to be co-scheduled with
        # the parent task. If the parent (current_task) is not sticky we must
        # set it to be sticky.
        # XXX: Ideally we would be able to unset this
        ct.sticky = true
        tid = Threads.threadid()
        ccall(:jl_set_task_tid, Cvoid, (Any, Cint), waiter, tid-1)
    end
    return
end

"""
    wait([x])

Block the current task until some event occurs, depending on the type of the argument:

* [`Channel`](@ref): Wait for a value to be appended to the channel.
* [`Condition`](@ref): Wait for [`notify`](@ref) on a condition and return the `val`
  parameter passed to `notify`.
* `Process`: Wait for a process or process chain to exit. The `exitcode` field of a process
  can be used to determine success or failure.
* [`Task`](@ref): Wait for a `Task` to finish. If the task fails with an exception, a
  `TaskFailedException` (which wraps the failed task) is thrown.
* [`RawFD`](@ref): Wait for changes on a file descriptor (see the `FileWatching` package).

If no argument is passed, the task blocks for an undefined period. A task can only be
restarted by an explicit call to [`schedule`](@ref) or [`yieldto`](@ref).

Often `wait` is called within a `while` loop to ensure a waited-for condition is met before
proceeding.
"""
function wait(c::GenericCondition)
    ct = current_task()
    _wait2(c, ct)
    token = unlockall(c.lock)
    try
        return wait()
    catch
        ct.queue === nothing || list_deletefirst!(ct.queue, ct)
        rethrow()
    finally
        relockall(c.lock, token)
    end
end

"""
    notify(condition, val=nothing; all=true, error=false)

Wake up tasks waiting for a condition, passing them `val`. If `all` is `true` (the default),
all waiting tasks are woken, otherwise only one is. If `error` is `true`, the passed value
is raised as an exception in the woken tasks.

Return the count of tasks woken up. Return 0 if no tasks are waiting on `condition`.
"""
notify(c::GenericCondition, @nospecialize(arg = nothing); all=true, error=false) = notify(c, arg, all, error)
function notify(c::GenericCondition, @nospecialize(arg), all, error)
    assert_havelock(c)
    cnt = 0
    while !isempty(c.waitq)
        t = popfirst!(c.waitq)
        schedule(t, arg, error=error)
        cnt += 1
        all || break
    end
    return cnt
end

notify_error(c::GenericCondition, err) = notify(c, err, true, true)

n_waiters(c::GenericCondition) = length(c.waitq)

"""
    isempty(condition)

Return `true` if no tasks are waiting on the condition, `false` otherwise.
"""
isempty(c::GenericCondition) = isempty(c.waitq)

# default (Julia v1.0) is currently single-threaded
# (although it uses MT-safe versions, when possible)
"""
    Condition()

Create an edge-triggered event source that tasks can wait for. Tasks that call [`wait`](@ref) on a
`Condition` are suspended and queued. Tasks are woken up when [`notify`](@ref) is later called on
the `Condition`. Edge triggering means that only tasks waiting at the time [`notify`](@ref) is
called can be woken up. For level-triggered notifications, you must keep extra state to keep
track of whether a notification has happened. The [`Channel`](@ref) and [`Threads.Event`](@ref) types do
this, and can be used for level-triggered events.

This object is NOT thread-safe. See [`Threads.Condition`](@ref) for a thread-safe version.
"""
const Condition = GenericCondition{AlwaysLockedST}

lock(c::GenericCondition{AlwaysLockedST}) =
    throw(ArgumentError("`Condition` is not thread-safe. Please use `Threads.Condition` instead for multi-threaded code."))
unlock(c::GenericCondition{AlwaysLockedST}) =
    throw(ArgumentError("`Condition` is not thread-safe. Please use `Threads.Condition` instead for multi-threaded code."))

max_stars_repo_name: pepijndevos/julia

max_issues_repo_name: pepijndevos/julia

max_forks_repo_name: pepijndevos/julia


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

module REPLCompletions

export completions, shell_completions, bslash_completions, completion_text

using Base.Meta
using Base: propertynames, something

abstract type Completion end

struct KeywordCompletion <: Completion
    keyword::String
end

struct PathCompletion <: Completion
    path::String
end

struct ModuleCompletion <: Completion
    parent::Module
    mod::String
end

struct PackageCompletion <: Completion
    package::String
end

struct PropertyCompletion <: Completion
    value
    property::Symbol
end

struct FieldCompletion <: Completion
    typ::DataType
    field::Symbol
end

struct MethodCompletion <: Completion
    func
    input_types::Type
    method::Method
end

struct BslashCompletion <: Completion
    bslash::String
end

struct ShellCompletion <: Completion
    text::String
end

struct DictCompletion <: Completion
    dict::AbstractDict
    key::String
end

# interface definition
function Base.getproperty(c::Completion, name::Symbol)
    if name === :keyword
        return getfield(c, :keyword)::String
    elseif name === :path
        return getfield(c, :path)::String
    elseif name === :parent
        return getfield(c, :parent)::Module
    elseif name === :mod
        return getfield(c, :mod)::String
    elseif name === :package
        return getfield(c, :package)::String
    elseif name === :property
        return getfield(c, :property)::Symbol
    elseif name === :field
        return getfield(c, :field)::Symbol
    elseif name === :method
        return getfield(c, :method)::Method
    elseif name === :bslash
        return getfield(c, :bslash)::String
    elseif name === :text
        return getfield(c, :text)::String
    elseif name === :key
        return getfield(c, :key)::String
    end
    return getfield(c, name)
end

_completion_text(c::KeywordCompletion) = c.keyword
_completion_text(c::PathCompletion) = c.path
_completion_text(c::ModuleCompletion) = c.mod
_completion_text(c::PackageCompletion) = c.package
_completion_text(c::PropertyCompletion) = string(c.property)
_completion_text(c::FieldCompletion) = string(c.field)
_completion_text(c::MethodCompletion) = sprint(io -> show(io, c.method))
_completion_text(c::BslashCompletion) = c.bslash
_completion_text(c::ShellCompletion) = c.text
_completion_text(c::DictCompletion) = c.key

completion_text(c) = _completion_text(c)::String

const Completions = Tuple{Vector{Completion}, UnitRange{Int}, Bool}

function completes_global(x, name)
    return startswith(x, name) && !('#' in x)
end

function appendmacro!(syms, macros, needle, endchar)
    for s in macros
        if endswith(s, needle)
            from = nextind(s, firstindex(s))
            to = prevind(s, sizeof(s)-sizeof(needle)+1)
            push!(syms, s[from:to]*endchar)
        end
    end
end

function filtered_mod_names(ffunc::Function, mod::Module, name::AbstractString, all::Bool = false, imported::Bool = false)
    ssyms = names(mod, all = all, imported = imported)
    filter!(ffunc, ssyms)
    syms = String[string(s) for s in ssyms]
    macros =  filter(x -> startswith(x, "@" * name), syms)
    appendmacro!(syms, macros, "_str", "\"")
    appendmacro!(syms, macros, "_cmd", "`")
    filter!(x->completes_global(x, name), syms)
    return [ModuleCompletion(mod, sym) for sym in syms]
end

# REPL Symbol Completions
function complete_symbol(sym::String, ffunc, context_module::Module=Main)
    mod = context_module
    name = sym

    lookup_module = true
    t = Union{}
    val = nothing
    if something(findlast(in(non_identifier_chars), sym), 0) < something(findlast(isequal('.'), sym), 0)
        # Find module
        lookup_name, name = rsplit(sym, ".", limit=2)

        ex = Meta.parse(lookup_name, raise=false, depwarn=false)

        b, found = get_value(ex, context_module)
        if found
            val = b
            if isa(b, Module)
                mod = b
                lookup_module = true
            elseif Base.isstructtype(typeof(b))
                lookup_module = false
                t = typeof(b)
            end
        else # If the value is not found using get_value, the expression contain an advanced expression
            lookup_module = false
            t, found = get_type(ex, context_module)
        end
        found || return Completion[]
    end

    suggestions = Completion[]
    if lookup_module
        # We will exclude the results that the user does not want, as well
        # as excluding Main.Main.Main, etc., because that's most likely not what
        # the user wants
        p = let mod=mod, modname=nameof(mod)
            s->(!Base.isdeprecated(mod, s) && s != modname && ffunc(mod, s)::Bool)
        end
        # Looking for a binding in a module
        if mod == context_module
            # Also look in modules we got through `using`
            mods = ccall(:jl_module_usings, Any, (Any,), context_module)::Vector
            for m in mods
                append!(suggestions, filtered_mod_names(p, m::Module, name))
            end
            append!(suggestions, filtered_mod_names(p, mod, name, true, true))
        else
            append!(suggestions, filtered_mod_names(p, mod, name, true, false))
        end
    elseif val !== nothing # looking for a property of an instance
        for property in propertynames(val, false)
            # TODO: support integer arguments (#36872)
            if property isa Symbol && startswith(string(property), name)
                push!(suggestions, PropertyCompletion(val, property))
            end
        end
    else
        # Looking for a member of a type
        if t isa DataType && t != Any
            # Check for cases like Type{typeof(+)}
            if t isa DataType && t.name === Base._TYPE_NAME
                t = typeof(t.parameters[1])
            end
            # Only look for fields if this is a concrete type
            if isconcretetype(t)
                fields = fieldnames(t)
                for field in fields
                    s = string(field)
                    if startswith(s, name)
                        push!(suggestions, FieldCompletion(t, field))
                    end
                end
            end
        end
    end
    suggestions
end

const sorted_keywords = [
    "abstract type", "baremodule", "begin", "break", "catch", "ccall",
    "const", "continue", "do", "else", "elseif", "end", "export", "false",
    "finally", "for", "function", "global", "if", "import",
    "let", "local", "macro", "module", "mutable struct",
    "primitive type", "quote", "return", "struct",
    "true", "try", "using", "while"]

function complete_keyword(s::Union{String,SubString{String}})
    r = searchsorted(sorted_keywords, s)
    i = first(r)
    n = length(sorted_keywords)
    while i <= n && startswith(sorted_keywords[i],s)
        r = first(r):i
        i += 1
    end
    Completion[KeywordCompletion(kw) for kw in sorted_keywords[r]]
end

function complete_path(path::AbstractString, pos::Int; use_envpath=false, shell_escape=false)
    if Base.Sys.isunix() && occursin(r"^~(?:/|$)", path)
        # if the path is just "~", don't consider the expanded username as a prefix
        if path == "~"
            dir, prefix = homedir(), ""
        else
            dir, prefix = splitdir(homedir() * path[2:end])
        end
    else
        dir, prefix = splitdir(path)
    end
    local files
    try
        if isempty(dir)
            files = readdir()
        elseif isdir(dir)
            files = readdir(dir)
        else
            return Completion[], 0:-1, false
        end
    catch
        return Completion[], 0:-1, false
    end

    matches = Set{String}()
    for file in files
        if startswith(file, prefix)
            id = try isdir(joinpath(dir, file)) catch; false end
            # joinpath is not used because windows needs to complete with double-backslash
            push!(matches, id ? file * (@static Sys.iswindows() ? "\\\\" : "/") : file)
        end
    end

    if use_envpath && length(dir) == 0
        # Look for files in PATH as well
        local pathdirs = split(ENV["PATH"], @static Sys.iswindows() ? ";" : ":")

        for pathdir in pathdirs
            local actualpath
            try
                actualpath = realpath(pathdir)
            catch
                # Bash doesn't expect every folder in PATH to exist, so neither shall we
                continue
            end

            if actualpath != pathdir && in(actualpath,pathdirs)
                # Remove paths which (after resolving links) are in the env path twice.
                # Many distros eg. point /bin to /usr/bin but have both in the env path.
                continue
            end

            local filesinpath
            try
                filesinpath = readdir(pathdir)
            catch e
                # Bash allows dirs in PATH that can't be read, so we should as well.
                if isa(e, Base.IOError) || isa(e, Base.ArgumentError)
                    continue
                else
                    # We only handle IOError and ArgumentError here
                    rethrow()
                end
            end

            for file in filesinpath
                # In a perfect world, we would filter on whether the file is executable
                # here, or even on whether the current user can execute the file in question.
                if startswith(file, prefix) && isfile(joinpath(pathdir, file))
                    push!(matches, file)
                end
            end
        end
    end

    matchList = Completion[PathCompletion(shell_escape ? replace(s, r"\s" => s"\\\0") : s) for s in matches]
    startpos = pos - lastindex(prefix) + 1 - count(isequal(' '), prefix)
    # The pos - lastindex(prefix) + 1 is correct due to `lastindex(prefix)-lastindex(prefix)==0`,
    # hence we need to add one to get the first index. This is also correct when considering
    # pos, because pos is the `lastindex` a larger string which `endswith(path)==true`.
    return matchList, startpos:pos, !isempty(matchList)
end

function complete_expanduser(path::AbstractString, r)
    expanded = expanduser(path)
    return Completion[PathCompletion(expanded)], r, path != expanded
end

# Determines whether method_complete should be tried. It should only be done if
# the string endswiths ',' or '(' when disregarding whitespace_chars
function should_method_complete(s::AbstractString)
    method_complete = false
    for c in reverse(s)
        if c in [',', '(']
            method_complete = true
            break
        elseif !(c in whitespace_chars)
            method_complete = false
            break
        end
    end
    method_complete
end

# Returns a range that includes the method name in front of the first non
# closed start brace from the end of the string.
function find_start_brace(s::AbstractString; c_start='(', c_end=')')
    braces = 0
    r = reverse(s)
    i = firstindex(r)
    in_single_quotes = false
    in_double_quotes = false
    in_back_ticks = false
    while i <= ncodeunits(r)
        c, i = iterate(r, i)
        if !in_single_quotes && !in_double_quotes && !in_back_ticks
            if c == c_start
                braces += 1
            elseif c == c_end
                braces -= 1
            elseif c == '\''
                in_single_quotes = true
            elseif c == '"'
                in_double_quotes = true
            elseif c == '`'
                in_back_ticks = true
            end
        else
            if !in_back_ticks && !in_double_quotes &&
                c == '\'' && i <= ncodeunits(r) && iterate(r, i)[1] != '\\'
                in_single_quotes = !in_single_quotes
            elseif !in_back_ticks && !in_single_quotes &&
                c == '"' && i <= ncodeunits(r) && iterate(r, i)[1] != '\\'
                in_double_quotes = !in_double_quotes
            elseif !in_single_quotes && !in_double_quotes &&
                c == '`' && i <= ncodeunits(r) && iterate(r, i)[1] != '\\'
                in_back_ticks = !in_back_ticks
            end
        end
        braces == 1 && break
    end
    braces != 1 && return 0:-1, -1
    method_name_end = reverseind(s, i)
    startind = nextind(s, something(findprev(in(non_identifier_chars), s, method_name_end), 0))::Int
    return (startind:lastindex(s), method_name_end)
end

# Returns the value in a expression if sym is defined in current namespace fn.
# This method is used to iterate to the value of a expression like:
# :(REPL.REPLCompletions.whitespace_chars) a `dump` of this expression
# will show it consist of Expr, QuoteNode's and Symbol's which all needs to
# be handled differently to iterate down to get the value of whitespace_chars.
function get_value(sym::Expr, fn)
    sym.head !== :. && return (nothing, false)
    for ex in sym.args
        fn, found = get_value(ex, fn)
        !found && return (nothing, false)
    end
    return (fn, true)
end
get_value(sym::Symbol, fn) = isdefined(fn, sym) ? (getfield(fn, sym), true) : (nothing, false)
get_value(sym::QuoteNode, fn) = isdefined(fn, sym.value) ? (getfield(fn, sym.value), true) : (nothing, false)
get_value(sym, fn) = (sym, true)

# Return the value of a getfield call expression
function get_value_getfield(ex::Expr, fn)
    # Example :((top(getfield))(Base,:max))
    val, found = get_value_getfield(ex.args[2],fn) #Look up Base in Main and returns the module
    (found && length(ex.args) >= 3) || return (nothing, false)
    return get_value_getfield(ex.args[3], val) #Look up max in Base and returns the function if found.
end
get_value_getfield(sym, fn) = get_value(sym, fn)

# Determines the return type with Base.return_types of a function call using the type information of the arguments.
function get_type_call(expr::Expr)
    f_name = expr.args[1]
    # The if statement should find the f function. How f is found depends on how f is referenced
    if isa(f_name, GlobalRef) && isconst(f_name.mod,f_name.name) && isdefined(f_name.mod,f_name.name)
        ft = typeof(eval(f_name))
        found = true
    else
        ft, found = get_type(f_name, Main)
    end
    found || return (Any, false) # If the function f is not found return Any.
    args = Any[]
    for ex in expr.args[2:end] # Find the type of the function arguments
        typ, found = get_type(ex, Main)
        found ? push!(args, typ) : push!(args, Any)
    end
    # use _methods_by_ftype as the function is supplied as a type
    world = Base.get_world_counter()
    matches = Base._methods_by_ftype(Tuple{ft, args...}, -1, world)
    length(matches) == 1 || return (Any, false)
    match = first(matches)
    # Typeinference
    interp = Core.Compiler.NativeInterpreter()
    return_type = Core.Compiler.typeinf_type(interp, match.method, match.spec_types, match.sparams)
    return_type === nothing && return (Any, false)
    return (return_type, true)
end

# Returns the return type. example: get_type(:(Base.strip("", ' ')), Main) returns (String, true)
function try_get_type(sym::Expr, fn::Module)
    val, found = get_value(sym, fn)
    found && return Core.Typeof(val), found
    if sym.head === :call
        # getfield call is special cased as the evaluation of getfield provides good type information,
        # is inexpensive and it is also performed in the complete_symbol function.
        a1 = sym.args[1]
        if isa(a1,GlobalRef) && isconst(a1.mod,a1.name) && isdefined(a1.mod,a1.name) &&
            eval(a1) === Core.getfield
            val, found = get_value_getfield(sym, Main)
            return found ? Core.Typeof(val) : Any, found
        end
        return get_type_call(sym)
    elseif sym.head === :thunk
        thk = sym.args[1]
        rt = ccall(:jl_infer_thunk, Any, (Any, Any), thk::Core.CodeInfo, fn)
        rt !== Any && return (rt, true)
    elseif sym.head === :ref
        # some simple cases of `expand`
        return try_get_type(Expr(:call, GlobalRef(Base, :getindex), sym.args...), fn)
    elseif sym.head === :.  && sym.args[2] isa QuoteNode # second check catches broadcasting
        return try_get_type(Expr(:call, GlobalRef(Core, :getfield), sym.args...), fn)
    end
    return (Any, false)
end

try_get_type(other, fn::Module) = get_type(other, fn)

function get_type(sym::Expr, fn::Module)
    # try to analyze nests of calls. if this fails, try using the expanded form.
    val, found = try_get_type(sym, fn)
    found && return val, found
    return try_get_type(Meta.lower(fn, sym), fn)
end

function get_type(sym, fn::Module)
    val, found = get_value(sym, fn)
    return found ? Core.Typeof(val) : Any, found
end

# Method completion on function call expression that look like :(max(1))
function complete_methods(ex_org::Expr, context_module::Module=Main)
    args_ex = Any[]
    func, found = get_value(ex_org.args[1], context_module)::Tuple{Any,Bool}
    !found && return Completion[]

    funargs = ex_org.args[2:end]
    # handle broadcasting, but only handle number of arguments instead of
    # argument types
    if ex_org.head === :. && ex_org.args[2] isa Expr
        for _ in (ex_org.args[2]::Expr).args
            push!(args_ex, Any)
        end
    else
        for ex in funargs
            val, found = get_type(ex, context_module)
            push!(args_ex, val)
        end
    end

    out = Completion[]
    t_in = Tuple{Core.Typeof(func), args_ex...} # Input types
    na = length(args_ex)+1
    ml = methods(func)
    for method in ml
        ms = method.sig

        # Check if the method's type signature intersects the input types
        if typeintersect(Base.rewrap_unionall(Tuple{(Base.unwrap_unionall(ms)::DataType).parameters[1 : min(na, end)]...}, ms), t_in) !== Union{}
            push!(out, MethodCompletion(func, t_in, method))
        end
    end
    return out
end

include("latex_symbols.jl")
include("emoji_symbols.jl")

const non_identifier_chars = [" \t\n\r\"\\'`\$><=:;|&{}()[],+-*/?%^~"...]
const whitespace_chars = [" \t\n\r"...]
# "\"'`"... is added to whitespace_chars as non of the bslash_completions
# characters contain any of these characters. It prohibits the
# bslash_completions function to try and complete on escaped characters in strings
const bslash_separators = [whitespace_chars..., "\"'`"...]

# Aux function to detect whether we're right after a
# using or import keyword
function afterusing(string::String, startpos::Int)
    (isempty(string) || startpos == 0) && return false
    str = string[1:prevind(string,startpos)]
    isempty(str) && return false
    rstr = reverse(str)
    r = findfirst(r"\s(gnisu|tropmi)\b", rstr)
    r === nothing && return false
    fr = reverseind(str, last(r))
    return occursin(r"^\b(using|import)\s*((\w+[.])*\w+\s*,\s*)*$", str[fr:end])
end

function bslash_completions(string::String, pos::Int)
    slashpos = something(findprev(isequal('\\'), string, pos), 0)
    if (something(findprev(in(bslash_separators), string, pos), 0) < slashpos &&
        !(1 < slashpos && (string[prevind(string, slashpos)]=='\\')))
        # latex / emoji symbol substitution
        s = string[slashpos:pos]
        latex = get(latex_symbols, s, "")
        if !isempty(latex) # complete an exact match
            return (true, (Completion[BslashCompletion(latex)], slashpos:pos, true))
        end
        emoji = get(emoji_symbols, s, "")
        if !isempty(emoji)
            return (true, (Completion[BslashCompletion(emoji)], slashpos:pos, true))
        end
        # return possible matches; these cannot be mixed with regular
        # Julian completions as only latex / emoji symbols contain the leading \
        if startswith(s, "\\:") # emoji
            namelist = Iterators.filter(k -> startswith(k, s), keys(emoji_symbols))
        else # latex
            namelist = Iterators.filter(k -> startswith(k, s), keys(latex_symbols))
        end
        return (true, (Completion[BslashCompletion(name) for name in sort!(collect(namelist))], slashpos:pos, true))
    end
    return (false, (Completion[], 0:-1, false))
end

function dict_identifier_key(str::String, tag::Symbol, context_module::Module = Main)
    if tag === :string
        str_close = str*"\""
    elseif tag === :cmd
        str_close = str*"`"
    else
        str_close = str
    end

    frange, end_of_identifier = find_start_brace(str_close, c_start='[', c_end=']')
    isempty(frange) && return (nothing, nothing, nothing)
    obj = context_module
    for name in split(str[frange[1]:end_of_identifier], '.')
        Base.isidentifier(name) || return (nothing, nothing, nothing)
        sym = Symbol(name)
        isdefined(obj, sym) || return (nothing, nothing, nothing)
        obj = getfield(obj, sym)
    end
    (isa(obj, AbstractDict) && length(obj)::Int < 1_000_000) || return (nothing, nothing, nothing)
    begin_of_key = something(findnext(!isspace, str, nextind(str, end_of_identifier) + 1), # +1 for [
                             lastindex(str)+1)
    return (obj::AbstractDict, str[begin_of_key:end], begin_of_key)
end

# This needs to be a separate non-inlined function, see #19441
@noinline function find_dict_matches(identifier::AbstractDict, partial_key)
    matches = String[]
    for key in keys(identifier)
        rkey = repr(key)
        startswith(rkey,partial_key) && push!(matches,rkey)
    end
    return matches
end

function project_deps_get_completion_candidates(pkgstarts::String, project_file::String)
    loading_candidates = String[]
    d = Base.parsed_toml(project_file)
    pkg = get(d, "name", nothing)::Union{String, Nothing}
    if pkg !== nothing && startswith(pkg, pkgstarts)
        push!(loading_candidates, pkg)
    end
    deps = get(d, "deps", nothing)::Union{Dict{String, Any}, Nothing}
    if deps !== nothing
        for (pkg, _) in deps
            startswith(pkg, pkgstarts) && push!(loading_candidates, pkg)
        end
    end
    return Completion[PackageCompletion(name) for name in loading_candidates]
end

function completions(string::String, pos::Int, context_module::Module=Main)
    # First parse everything up to the current position
    partial = string[1:pos]
    inc_tag = Base.incomplete_tag(Meta.parse(partial, raise=false, depwarn=false))

    # if completing a key in a Dict
    identifier, partial_key, loc = dict_identifier_key(partial, inc_tag, context_module)
    if identifier !== nothing
        matches = find_dict_matches(identifier, partial_key)
        length(matches)==1 && (lastindex(string) <= pos || string[nextind(string,pos)] != ']') && (matches[1]*=']')
        length(matches)>0 && return Completion[DictCompletion(identifier, match) for match in sort!(matches)], loc::Int:pos, true
    end

    # otherwise...
    if inc_tag in [:cmd, :string]
        m = match(r"[\t\n\r\"`><=*?|]| (?!\\)", reverse(partial))
        startpos = nextind(partial, reverseind(partial, m.offset))
        r = startpos:pos

        expanded = complete_expanduser(replace(string[r], r"\\ " => " "), r)
        expanded[3] && return expanded  # If user expansion available, return it

        paths, r, success = complete_path(replace(string[r], r"\\ " => " "), pos)

        if inc_tag === :string &&
           length(paths) == 1 &&  # Only close if there's a single choice,
           !isdir(expanduser(replace(string[startpos:prevind(string, first(r))] * paths[1].path,
                                     r"\\ " => " "))) &&  # except if it's a directory
           (lastindex(string) <= pos ||
            string[nextind(string,pos)] != '"')  # or there's already a " at the cursor.
            paths[1] = PathCompletion(paths[1].path * "\"")
        end

        #Latex symbols can be completed for strings
        (success || inc_tag==:cmd) && return sort!(paths, by=p->p.path), r, success
    end

    ok, ret = bslash_completions(string, pos)
    ok && return ret

    # Make sure that only bslash_completions is working on strings
    inc_tag==:string && return Completion[], 0:-1, false
    if inc_tag === :other && should_method_complete(partial)
        frange, method_name_end = find_start_brace(partial)
        # strip preceding ! operator
        s = replace(partial[frange], r"\!+([^=\(]+)" => s"\1")
        ex = Meta.parse(s * ")", raise=false, depwarn=false)

        if isa(ex, Expr)
            if ex.head === :call
                return complete_methods(ex, context_module), first(frange):method_name_end, false
            elseif ex.head === :. && ex.args[2] isa Expr && (ex.args[2]::Expr).head === :tuple
                return complete_methods(ex, context_module), first(frange):(method_name_end - 1), false
            end
        end
    elseif inc_tag === :comment
        return Completion[], 0:-1, false
    end

    dotpos = something(findprev(isequal('.'), string, pos), 0)
    startpos = nextind(string, something(findprev(in(non_identifier_chars), string, pos), 0))
    # strip preceding ! operator
    if (m = match(r"^\!+", string[startpos:pos])) !== nothing
        startpos += length(m.match)
    end

    ffunc = (mod,x)->true
    suggestions = Completion[]
    comp_keywords = true
    if afterusing(string, startpos)
        # We're right after using or import. Let's look only for packages
        # and modules we can reach from here

        # If there's no dot, we're in toplevel, so we should
        # also search for packages
        s = string[startpos:pos]
        if dotpos <= startpos
            for dir in Base.load_path()
                if basename(dir) in Base.project_names && isfile(dir)
                    append!(suggestions, project_deps_get_completion_candidates(s, dir))
                end
                isdir(dir) || continue
                for pname in readdir(dir)
                    if pname[1] != '.' && pname != "METADATA" &&
                        pname != "REQUIRE" && startswith(pname, s)
                        # Valid file paths are
                        #   <Mod>.jl
                        #   <Mod>/src/<Mod>.jl
                        #   <Mod>.jl/src/<Mod>.jl
                        if isfile(joinpath(dir, pname))
                            endswith(pname, ".jl") && push!(suggestions,
                                                            PackageCompletion(pname[1:prevind(pname, end-2)]))
                        else
                            mod_name = if endswith(pname, ".jl")
                                pname[1:prevind(pname, end-2)]
                            else
                                pname
                            end
                            if isfile(joinpath(dir, pname, "src",
                                               "$mod_name.jl"))
                                push!(suggestions, PackageCompletion(mod_name))
                            end
                        end
                    end
                end
            end
        end
        ffunc = (mod,x)->(Base.isbindingresolved(mod, x) && isdefined(mod, x) && isa(getfield(mod, x), Module))
        comp_keywords = false
    end
    startpos == 0 && (pos = -1)
    dotpos < startpos && (dotpos = startpos - 1)
    s = string[startpos:pos]
    comp_keywords && append!(suggestions, complete_keyword(s))
    # The case where dot and start pos is equal could look like: "(""*"").d","". or  CompletionFoo.test_y_array[1].y
    # This case can be handled by finding the beginning of the expression. This is done below.
    if dotpos == startpos
        i = prevind(string, startpos)
        while 0 < i
            c = string[i]
            if c in [')', ']']
                if c==')'
                    c_start='('; c_end=')'
                elseif c==']'
                    c_start='['; c_end=']'
                end
                frange, end_of_identifier = find_start_brace(string[1:prevind(string, i)], c_start=c_start, c_end=c_end)
                startpos = first(frange)
                i = prevind(string, startpos)
            elseif c in ('\'', '\"', '\`')
                s = "$c$c"*string[startpos:pos]
                break
            else
                break
            end
            s = string[startpos:pos]
        end
    end
    append!(suggestions, complete_symbol(s, ffunc, context_module))
    return sort!(unique(suggestions), by=completion_text), (dotpos+1):pos, true
end

function shell_completions(string, pos)
    # First parse everything up to the current position
    scs = string[1:pos]
    local args, last_parse
    try
        args, last_parse = Base.shell_parse(scs, true)::Tuple{Expr,UnitRange{Int}}
    catch
        return Completion[], 0:-1, false
    end
    ex = args.args[end]::Expr
    # Now look at the last thing we parsed
    isempty(ex.args) && return Completion[], 0:-1, false
    arg = ex.args[end]
    if all(s -> isa(s, AbstractString), ex.args)
        arg = arg::AbstractString
        # Treat this as a path

        # As Base.shell_parse throws away trailing spaces (unless they are escaped),
        # we need to special case here.
        # If the last char was a space, but shell_parse ignored it search on "".
        ignore_last_word = arg != " " && scs[end] == ' '
        prefix = ignore_last_word ? "" : join(ex.args)

        # Also try looking into the env path if the user wants to complete the first argument
        use_envpath = !ignore_last_word && length(args.args) < 2

        return complete_path(prefix, pos, use_envpath=use_envpath, shell_escape=true)
    elseif isexpr(arg, :incomplete) || isexpr(arg, :error)
        partial = scs[last_parse]
        ret, range = completions(partial, lastindex(partial))
        range = range .+ (first(last_parse) - 1)
        return ret, range, true
    end
    return Completion[], 0:-1, false
end

end # module

max_stars_repo_name: shikhargoswami/julia

max_issues_repo_name: PallHaraldsson/julia

max_forks_repo_name: PallHaraldsson/julia