text:
"""
```
koopman_smoother(y, T, R, C, Q, Z, D, E, s_0, P_0, s_pred, P_pred;
    Nt0 = 0)

koopman_smoother(regime_indices, y, Ts, Rs, Cs, Qs, Zs, Ds, Es,
    s_0, P_0, s_pred, P_pred; Nt0 = 0)
```

This is a Kalman smoothing program based on S.J. Koopman's \"Disturbance
Smoother for State Space Models\" (Biometrika, 1993), as specified in
Durbin and Koopman's \"A Simple and Efficient Simulation Smoother for
State Space Time Series Analysis\" (Biometrika, 2002).

Unlike other Kalman smoothing programs, there is no need to invert
singular matrices using the Moore-Penrose pseudoinverse (`pinv`), which
should lead to efficiency gains and fewer inversion problems. Also, the
states vector and the corresponding matrices do not need to be augmented
to include the shock innovations. Instead they are saved automatically
in the `ϵ_smth` matrix.

The state space is given by:

```
s_{t+1} = C + T*s_t + R*ϵ_t    (transition equation)
y_t     = D + Z*s_t + u_t      (measurement equation)

ϵ_t ∼ N(0, Q)
u_t ∼ N(0, E)
Cov(ϵ_t, u_t) = 0
```

### Inputs

- `y`: `Ny` x `Nt` matrix containing data `y_1, ... , y_T`
- `s_0`: `Ns` x 1 initial state vector
- `P_0`: `Ns` x `Ns` initial state covariance matrix
- `s_pred`: `Ns` x `Nt` matrix of one-step predicted state vectors `s_{t|t-1}`
  (from the Kalman filter)
- `P_pred`: `Ns` x `Ns` x `Nt` array of mean squared errors `P_{t|t-1}` of
  predicted state vectors

**Method 1 only:** state-space system matrices `T`, `R`, `C`, `Q`, `Z`,
`D`, `E`. See `?kalman_filter`

**Method 2 only:** `regime_indices` and system matrices for each regime `Ts`,
`Rs`, `Cs`, `Qs`, `Zs`, `Ds`, `Es`. See `?kalman_filter`

where:

- `Nt`: number of periods for which we have data
- `Ns`: number of states
- `Ne`: number of shocks
- `Ny`: number of observables

### Keyword Arguments

- `Nt0`: if greater than 0, the returned smoothed states and
  shocks matrices will be shorter by that number of columns (taken from the
  beginning)

### Outputs

- `s_smth`: `Ns` x `Nt` matrix of smoothed states `s_{t|T}`
- `ϵ_smth`: `Ne` x `Nt` matrix of smoothed shocks `ϵ_{t|T}`
"""
function koopman_smoother(y::AbstractMatrix,
    T::Matrix{S}, R::Matrix{S}, C::Vector{S}, Q::Matrix{S},
    Z::Matrix{S}, D::Vector{S}, E::Matrix{S},
    s_0::Vector{S}, P_0::Matrix{S}, s_pred::Matrix{S}, P_pred::Array{S, 3};
    Nt0::Int = 0) where {S<:AbstractFloat}

    Nt = size(y, 2)
    koopman_smoother(UnitRange{Int}[1:Nt], y, Matrix{S}[T], Matrix{S}[R], Vector{S}[C],
        Matrix{S}[Q], Matrix{S}[Z], Vector{S}[D], Matrix{S}[E],
        s_0, P_0, s_pred, P_pred; Nt0 = Nt0)
end

function koopman_smoother(regime_indices::Vector{UnitRange{Int}}, y::AbstractMatrix,
    Ts::Vector{Matrix{S}}, Rs::Vector{Matrix{S}}, Cs::Vector{Vector{S}}, Qs::Vector{Matrix{S}},
    Zs::Vector{Matrix{S}}, Ds::Vector{Vector{S}}, Es::Vector{Matrix{S}},
    s_0::Vector{S}, P_0::Matrix{S}, s_pred::Matrix{S}, P_pred::Array{S, 3};
    Nt0::Int = 0) where {S<:AbstractFloat}

    # Dimensions
    n_regimes = length(regime_indices)
    Nt = size(y, 2)     # number of periods of data
    Ns = size(Ts[1], 1) # number of states
    Ne = size(Rs[1], 2) # number of shocks

    # Call disturbance smoother
    s_dist, _ = koopman_disturbance_smoother(regime_indices, y, Ts, Rs, Qs,
                                             Zs, Ds, Es, s_pred, P_pred; Nt0 = 0)

    # Initialize outputs
    s_smth = zeros(S, Ns, Nt)
    ϵ_smth = zeros(S, Ne, Nt)

    s_t = zeros(S, Ns) # initialize dummy value s.t. s_t is in scope

    for i = 1:n_regimes
        # Get state-space system matrices for this regime
        T, R, Q = Ts[i], Rs[i], Qs[i]

        for t in regime_indices[i]
            r_t = s_dist[:, t]
            s_t = if t == 1
                s_0 + P_0*r_t
            else
                T*s_t + R*Q*R'*r_t
            end

            s_smth[:, t] = s_t
            ϵ_smth[:, t] = Q*R'*r_t
        end
    end

    # Trim the presample if needed
    if Nt0 > 0
        insample = Nt0+1:Nt

        s_smth = s_smth[:, insample]
        ϵ_smth = ϵ_smth[:, insample]
    end

    return s_smth, ϵ_smth
end

"""
```
koopman_disturbance_smoother(y, T, R, Q, Z, D, E, s_pred, P_pred;
    Nt0 = 0)

koopman_smoother(regime_indices, y, Ts, Rs, Qs, Zs, Ds, Es,
    s_pred, P_pred; Nt0 = 0)
```

This disturbance smoother is intended for use with the state smoother
`koopman_smoother` from S.J. Koopman's \"Disturbance Smoother for State Space
Models\" (Biometrika, 1993), as specified in Durbin and Koopman's \"A Simple and
Efficient Simulation Smoother for State Space Time Series Analysis\"
(Biometrika, 2002).

The state space is given by:

```
s_{t+1} = C + T*s_t + R*ϵ_t    (transition equation)
y_t     = D + Z*s_t + u_t      (measurement equation)

ϵ_t ∼ N(0, Q)
u_t ∼ N(0, E)
Cov(ϵ_t, u_t) = 0
```

### Inputs

- `y`: `Ny` x `Nt` matrix containing data `y_1, ... , y_T`
- `s_pred`: `Ns` x `Nt` matrix of one-step predicted state vectors `s_{t|t-1}`
  (from the Kalman filter)
- `P_pred`: `Ns` x `Ns` x `Nt` array of mean squared errors `P_{t|t-1}` of
  predicted state vectors

**Method 1 only:** state-space system matrices `T`, `R`, `Q`, `Z`, `D`, `E`. See
`?kalman_filter`

**Method 2 only:** `regime_indices` and system matrices for each regime `Ts`,
`Rs`, `Qs`, `Zs`, `Ds`, `Es`. See `?kalman_filter`

where:

- `Nt`: number of periods for which we have data
- `Ns`: number of states
- `Ne`: number of shocks
- `Ny`: number of observables

### Keyword Arguments

- `Nt0`: if greater than 0, the returned smoothed disturbances and shocks
  matrices will be shorter by that number of columns (taken from the beginning)

### Outputs

- `s_dist`: `Ns` x `Nt` matrix of transition equation disturbances `r_t`
- `y_dist`: `Ny` x `Nt` matrix of measurement equation disturbances `e_t`
"""
function koopman_disturbance_smoother(y::AbstractArray,
    T::Matrix{S}, R::Matrix{S}, Q::Matrix{S},
    Z::Matrix{S}, D::Vector{S}, E::Matrix{S},
    s_pred::Matrix{S}, P_pred::Array{S, 3}; Nt0::Int = 0) where {S<:AbstractFloat}

    Nt = size(y, 2)
    koopman_disturbance_smoother(UnitRange{Int}[1:Nt], y, Matrix{S}[T], Matrix{S}[R],
        Matrix{S}[Q], Matrix{S}[Z], Vector{S}[D], Matrix{S}[E],
        s_pred, P_pred; Nt0 = Nt0)
end

function koopman_disturbance_smoother(regime_indices::Vector{UnitRange{Int}}, y::AbstractArray,
    Ts::Vector{Matrix{S}}, Rs::Vector{Matrix{S}}, Qs::Vector{Matrix{S}},
    Zs::Vector{Matrix{S}}, Ds::Vector{Vector{S}}, Es::Vector{Matrix{S}},
    s_pred::Matrix{S}, P_pred::Array{S, 3}; Nt0::Int = 0) where {S<:AbstractFloat}

    # Dimensions
    Nt = size(y ,    2) # number of periods of data
    Ns = size(Ts[1], 1) # number of states
    Ne = size(Rs[1], 2) # number of shocks
    Ny = size(Zs[1], 1) # number of observables

    # Initialize outputs
    s_dist = zeros(S, Ns, Nt)
    y_dist = zeros(S, Ny, Nt)

    r_t = zeros(S, Ns) # r_0 = 0

    for i = length(regime_indices):-1:1
        # Get state-space system matrices for this regime
        T, R, Q = Ts[i], Rs[i], Qs[i]
        Z, D, E = Zs[i], Ds[i], Es[i]

        for t in reverse(regime_indices[i])
            # Keep rows of measurement equation corresponding to nonmissing observables
            nonmissing = .!map(x -> ismissing(x) ? true : isnan(x), y[:, t]) #.!ismissing.(y[:, t])
            y_t = y[nonmissing, t]
            Z_t = Z[nonmissing, :]
            D_t = D[nonmissing]
            E_t = E[nonmissing, nonmissing]

            s_pred_t = s_pred[:, t]            # s_{t|t-1}
            P_pred_t = P_pred[:, :, t]         # P_{t|t-1} = Var s_{t|t-1}

            y_pred_t = Z_t*s_pred_t + D_t      # y_{t|t-1} = Z*s_{t|t-1} + D
            V_pred_t = Z_t*P_pred_t*Z_t' + E_t # V_{t|t-1} = Var y_{t|t-1} = Z*P_{t|t-1}*Z' + E
            dy = y_t - y_pred_t                # dy = y_t - y_{t|t-1} = prediction error
            K = T*P_pred_t*Z_t'/V_pred_t       # K_t = T*P_{t|t-1}'Z'/V_{t|t-1} = Kalman gain

            e_t = V_pred_t\dy - K'*r_t         # e_t     = (1/V_{t|t-1})dy - K_t'*r_t
            r_t = Z_t'*e_t + T'*r_t            # r_{t-1} = Z'*e_t + T'*r_t

            s_dist[:,      t] = r_t
            y_dist[nonmissing, t] = e_t
        end # of loop backward through this regime's periods
    end # of loop backward through regimes

    # Trim the presample if needed
    if Nt0 > 0
        insample = Nt0+1:Nt
        s_dist = s_dist[:, insample]
        y_dist = y_dist[:, insample]
    end

    return s_dist, y_dist
end


text:
#!/usr/bin/env julia
# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy
using RainFARM
using ArgParse
using Printf

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--radius", "-r"
            help = "Smoothing radius (in grid units)"
            arg_type = Float64
            default = 0.
        "--lapse", "-l"
            help = "Lapse rate [K/km]"
            arg_type =  Float64
            default = 6.5
        "--varname", "-v"
            help = "Input variable name (in orofile)"
            arg_type = AbstractString
            default = ""
        "--coord", "--box", "-b"
            help = "Coordinates of box in which to downscale [lon1 lon2 lat1 lat2]"
 	    nargs = 4
            arg_type =  Union{Float64, Nothing}
	    default =  [nothing, nothing, nothing, nothing]
        "--orocoarse", "-o"
            help = "Coarse orography of original data"
            arg_type = AbstractString
            default = ""
        "orofile"
            help = "The input file to use for orography"
            arg_type = AbstractString
            required = true
        "infile"
            help = "The file to downscale"
            arg_type = AbstractString
            required = true
        "outfile"
            help = "The output file name"
            arg_type = AbstractString
            required = true
    end

    s.description = "Downscale temperature with a given lapse rate"
    s.version = "0.1"
    s.add_version = true

    return parse_args(s)
end

args = parse_commandline()
lapse = abs(args["lapse"])
radius = args["radius"]
filein = args["infile"]
fileout = args["outfile"]
fileoro = args["orofile"]
fileorocoarse = args["orocoarse"]
varname = args["varname"]
coord = args["coord"]

lon1 = coord[1]
lon2 = coord[2]
lat1 = coord[3]
lat2 = coord[4]

(tin0, lonl0, latl0) = read_netcdf2d(filein, varname)
if length(size(lonl0))>1
    dxl = max(lonl0[2,1]-lonl0[1,1], lonl0[1,2]-lonl0[1,1])
    if lon1 == nothing
        lon1 = lonl0[1,1]
        lon2 = lonl0[end,1]
        lat1 = latl0[1,1]
        lat2 = latl0[1,end]
    end
else
    dxl = lonl0[2]-lonl0[1]
    if lon1 == nothing
        lon1 = lonl0[1]
        lon2 = lonl0[end]
        lat1 = latl0[1]
        lat2 = latl0[end]
    end
end

println("Downscaling in box: ", lon1, "/", lon2, "/", lat1, "/", lat2)

println("dx=", dxl)
if radius==0
    radius = dxl/2
end

#Add buffer
lonl1 = lon1-radius
lonl2 = lon2+radius
latl1 = lat1-radius
latl2 = lat2+radius
println("box + buffer=", lonl1, "/", lonl2, "/", latl1, "/", latl2)

rr = round.(Int32,rand()*1000000)

run(`cdo -s sellonlatbox,$lonl1,$lonl2,$latl1,$latl2 $fileoro orocut$rr.nc`)
(oro, lonl2, latl2, oroname) = read_netcdf2d("orocut$rr.nc", "")
if length(size(lonl2))>1
    dxf = max(lonl2[2,1]-lonl2[1,1], lonl2[1,2]-lonl2[1,1])
else
    dxf = lonl2[2]-lonl2[1]
end

println("Remapping input data ...")
run(`cdo -s -b F32 remapnn,orocut$rr.nc $filein input_nn$rr.nc`)
(tin, lonl, latl, varname) = read_netcdf2d("input_nn$rr.nc", varname)
oro = float(oro) # convert to float
(nx, ny, nt) = size(tin)

nf2 = div(radius, abs(dxf))
println("Smoothing radius = ", radius, " = ", nf2, " pixel")

println("Preparing correction ...")
if fileorocoarse==""
    oros = smooth(oro, nf2)
else
    run(`cdo -s -b F32 remapnn,orocut$rr.nc $fileorocoarse orocut_coarse$rr.nc`)
    (oroc, lonl, latl, orocname) = read_netcdf2d("orocut_coarse$rr.nc", "")
    oros = smooth(oroc, nf2)
    run(`rm orocut_coarse$rr.nc`)
end
#println("oro=",mean(oro)," oros=",mean(oros))
oro = -(oro-oros)*lapse/1000.

write_netcdf2d("orocorr_temp$rr.nc", oro, lonl, latl, oroname, "orocut$rr.nc")
run(`cdo -s sellonlatbox,$lon1,$lon2,$lat1,$lat2 orocorr_temp$rr.nc orocorr$rr.nc `)

println("Downscaling ...")
for i = 1:nt
    println("t=",i)
    tins = smooth(tin[:,:,i], nf2)
    tin[:,:,i] = tins+oro
#println("tin=",mean(tin[:,:,i])," tins=",mean(tins))
end

write_netcdf2d("out_temp$rr.nc", tin, lonl, latl, varname, "input_nn$rr.nc")
run(`cdo -s sellonlatbox,$lon1,$lon2,$lat1,$lat2 out_temp$rr.nc $fileout `)
run(`rm orocut$rr.nc input_nn$rr.nc out_temp$rr.nc orocorr_temp$rr.nc orocorr$rr.nc`)

ext: jl


text:
abstract type AbstractNeuronSelector end

"""
    MaxActivationSelector()

Neuron selector that picks the output neuron with the highest activation.
"""
struct MaxActivationSelector <: AbstractNeuronSelector end
function (::MaxActivationSelector)(out::AbstractArray{T,N}) where {T,N}
    N < 2 && throw(BATCHDIM_MISSING)
    return Vector{CartesianIndex{N}}([argmax(out; dims=1:(N - 1))...])
end

"""
    IndexSelector(index)

Neuron selector that picks the output neuron at the given index.
"""
struct IndexSelector{I} <: AbstractNeuronSelector
    index::I
end
function (s::IndexSelector{<:Integer})(out::AbstractArray{T,N}) where {T,N}
    N < 2 && throw(BATCHDIM_MISSING)
    return CartesianIndex{N}.(s.index, 1:size(out, N))
end
function (s::IndexSelector{I})(out::AbstractArray{T,N}) where {I,T,N}
    N < 2 && throw(BATCHDIM_MISSING)
    return CartesianIndex{N}.(s.index..., 1:size(out, N))
end

"""
    AugmentationSelector(index)

Neuron selector that passes through an augmented neuron selection.
"""
struct AugmentationSelector{I} <: AbstractNeuronSelector
    indices::I
end
(s::AugmentationSelector)(out) = s.indices


text:
"""
	write_mccs(of, MCCs::AbstractArray)

Write MCCs to file.
"""
function write_mccs(of, MCCs::AbstractArray, mode="w")
	open(of, mode) do w
		for (i,m) in enumerate(MCCs)
			for x in m[1:end-1]
				write(w, x * ",")
			end
			if i < length(MCCs)
				write(w, m[end] * "\n")
			else
				write(w, m[end])
			end
		end
	end

	return nothing
end

function read_mccs(file)
	map(eachline(file)) do m
       String.(split(m, ','))
    end
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
module InterfaceTesting
using Base.Test

export partial_method_exists,
    test_iterator_interface,
    test_index_interface,
    test_abstractarray_interface

"""
Checks if a method exist.
The first parameter is the function.
The second parameter is a list of types, or wild cards.
THe `:` is used as a wildcard.
Wild cards match any type (Not just `Any`)

"""
function partial_method_exists(f, types)
    meths = methods(f)
    for mm in meths
        @show mm
        sig_types =  mm.sig.types[2:end]
        length(types) != length(sig_types)  && continue

        for (spec, sig) in zip(types, sig_types)
            spec==(:) && continue
            spec <: sig || @goto next_meth
        end
        return true
        @label next_meth
    end
    return false
end

"""
Checks that all the method expected of an Interator are defined.
This is the most formal and common of informal interfaced.
It is practically useless to not implement it properly.
"""
function test_iterator_interface(itertype; testset_type=Test.DefaultTestSet)
    @testset testset_type "$(itertype)" begin
        #Basic Methods
        @test method_exists(start, (itertype,))
        @test partial_method_exists(next, (itertype,:))
        @test partial_method_exists(done, (itertype,:))

        #Iterators size

        @test which(Base.iteratorsize, (itertype,)) == which(Base.iteratorsize, (Any,))#Method for Instance of type should always be the fall back -- never set directly

        if Base.iteratorsize(itertype)==Base.HasShape()
            @test method_exists(size, (itertype,))
            @test method_exists(length, (itertype,))
        elseif Base.iteratorsize(itertype)==Base.HasLength()
            @test method_exists(length, (itertype,))
        else
            @test Base.iteratorsize(itertype) in (Base.SizeUnknown(), Base.IsInfinite())
        end

        #Iterator Eltype
        if Base.iteratoreltype(itertype) == Base.HasEltype()
            @test method_exists(eltype, (itertype,))
        end

    end
end

"""
Checks Indexing Interface implemented.
This is the least formal of informal interfaces.
There a plenty of useful partial implementatons of it.
Eg a read-only things
"""
function test_index_interface(kind; testset_type=Test.DefaultTestSet)
    @testset testset_type "$(kind)" begin
        if method_exists(getindex, (kind, Vararg))
            @test partial_method_exists(size, (kind, :)) #otherwise `end` won't work in multidimentional index expressions
        else #Assume 1D indexing only has been defined
            @test partial_method_exists(getindex, (kind,:))
        end
        @test partial_method_exists(setindex!, (kind,:,:)) #Got to have something. 1D will do. Might have more, kinda complex
        @test method_exists(endof, (kind,)) #If not defined then `end` won't work

    end
end

"""
Checks AbstractArray Interface implemented.
This is required for all subtypes of AbstractArray.
Filling this will mean automatically furfilling `test_index_interface`
"""
function test_abstractarray_interface(kind; testset_type=Test.DefaultTestSet)
    @testset testset_type "$(kind)" begin
        @test method_exists(size, (kind,))
        @test method_exists(getindex, (kind, Vararg{Int}))
        @test partial_method_exists(setindex!, (kind, :, :))
    end
end

end # module


text:
### A Pluto.jl notebook ###
# v0.14.6

using Markdown
using InteractiveUtils

# ╔═╡ 0f6cd296-6cb8-4a9d-a21f-02b689dfe4e9
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 646a9fc2-3a73-4f8d-a952-2cdcf0e6287a
begin
	# Packages
	using LightGraphs
	using GraphPlot
	using Clustering
	using SparseArrays
	using LinearAlgebra
	using Plots
end

# ╔═╡ 7cd3af86-f459-4e4b-8e35-c29bc6966eed
begin
	# Define sizes
	using Random
	m=[200,100,100]
	n=[100,200,100]
	density=[0.5,0.7,0.2]
	Pts=Array{Any}(undef,3)
	Random.seed!(421)
	for i=1:3
	    # Generate sparse random
	    Pts[i]=sprand(m[i],n[i],density[i])
	end
	B=blockdiag(Pts[1],Pts[2],Pts[3])
end

# ╔═╡ 9ab98160-3e37-44c3-b219-a3e1bf95fd68
using Arpack

# ╔═╡ 89299cb7-a4eb-48f3-b3f1-ab2b5c6cfa2e
md"""
# Spectral Partitioning of Bipartite Graphs

Typical example of bipartite graph is a graph obtained from a collection of documents presented as a _term $\times$ document_  matrix.

__Prerequisites__

The reader should be familiar with k-means algorithm and spectral graph partitioning theory and algorithms.

__Competences__

The reader should be able to apply spectral partitioning of bipartite graphs to data clustering problems.

__Credits.__ The notebook was initially derived from M.Sc. Thesis of Ivančica Mirošević.
"""

# ╔═╡ a92d0e70-479d-444e-b06d-e47d2e385b10
md"""
## Definitions

__Undirected bipartite graph__ $G$ is a triplet $G=(T,D,E)$, where $T=\{t_{1},\cdots ,t_{m}\}$ and $D=\{d_{1},...,d_{n}\}$ are two sets of vertices and $E=\{(t_{i},d_{j}):t_{i}\in R,d_{j}\in D\}$, is a set of edges.

 $G$ is __weighted__ if there is weight $\omega(e)$ associated with each edge $e\in E$.

For example, $D$ is a set of documents, $T$ is a set of terms (words) and  edge $e=(t_{i},d_{j})$ exists if document $d_{j}$ contains term $t_{i}$. Weight $\omega(e)$ can be number of appearances of the term $t_i$ in the document $d_j$.

A __term-by-document-matrix__ is a matrix $A\in\mathbb{R}^{m\times n}$ with $A_{ij}=\omega((t_i,d_j))$.
"""

# ╔═╡ fddcc2a7-f718-4900-a135-764eb36496a7
md"""
## Facts

1. The weight matrix of $G$ is $W=\begin{bmatrix}0 & A \\ A^{T} & 0 \end{bmatrix}$.

2. The Laplacian matrix of $G$ is

$$L=\begin{bmatrix} \Delta_{1} & -A \\ -A^{T} & \Delta_{2}\end{bmatrix},$$

where $\Delta_1$ and $\Delta_2$ are diagonal matrices with elements

$\Delta_{1,ii}=\sum\limits_{j=1}^n A_{ij},\quad i=1,\ldots,m,$

$\Delta_{2,jj}=\sum\limits_{i=1}^m A_{ij},\quad j=1,\ldots,n.$

3. The normalized Laplacian matrix of $G$ is

$$L_n=\begin{bmatrix}
I & -\Delta_{1}^{-\frac{1}{2}}A\Delta_{2}^{-\frac{1}{2}} \\
-\Delta_{2}^{-\frac{1}{2}}A^T\Delta_{1}^{-\frac{1}{2}} & I
\end{bmatrix} \equiv
\begin{bmatrix} I & -A_n \\ -A_n^T & I \end{bmatrix}.$$
"""

# ╔═╡ b3b44560-203a-11eb-2782-170e4c12627a
md"""
4. Let $\lambda$ be an eigenvalue of $L_n$ with an eigenvector $w=\begin{bmatrix} u \\ v\end{bmatrix}$, where $u\in \mathbb{R}^{m}$ $v\in\mathbb{R}^{n}$. Then $L_n w=\lambda w$ implies $A_n v =(1-\lambda)u$ and $A_n^T u=(1-\lambda)v$. Vice versa, if $(u,\sigma,v)$ is a singular triplet of $A_n$, then $1-\sigma$ is an eigenvalue of $L_n$ with (non-unit) eigenvector $w=\begin{bmatrix} u \\ v\end{bmatrix}$.

5. The second largest singular value of $A_n$ corresponds to the second smallest eigenvalue of $L_n$, and computing the former is numerically more stable.
"""

# ╔═╡ 1be90f30-203b-11eb-34da-c3116523ec3a
md"

## Algorithms

### Bipartitioning algorithm

1. For given $A$ compute $A_{n}$.
2. Compute singular vectors of $A_{n}$, $u^{[2]}$ and $v^{[2]}$, which correspond to the second largest singular value, $\sigma_2(A_n)$.
3. Assign the partitions $T=\{T_1,T_2\}$ and $D=\{D_1,D_2\}$ according to the signs of $u^{[2]}$ and $v^{[2]}$. The pair $(T,D)$ is now partitioned as $\{(T_1,D_1),(T_2,D_2)\}$.
"

# ╔═╡ 5211ce2e-203b-11eb-0b6f-55e809177fb9
md"
### Recursive bipartitioning algorithm

1. Compute the bipartition $\pi=\{(T_1,D_1),(T_2,D_2)\}$ of $(T,D)$. Set the counter $c=2$.

2. While $c<k$ repeat

   - compute bipartitions of each of the subpartitions of $(T,D)$,
   - among all $(c+1)$-subpartitions, choose the one with the smallest $\mathop{\mathrm{pcut}}(\pi_{c+1})$ or $\mathop{\mathrm{ncut}}(\pi_{c+1})$, respectively.
   - Set $c=c+1$

3. Stop
"

# ╔═╡ 38333a20-203c-11eb-2b28-59e07b3e50ee
md"
### Multipartitioning algorithm
1. For given $A$ compute $A_{n}$.
2. Compute $k$ left and right singular vectors, $u^{[1]},\ldots,u^{[k]}$ and $v^{[1]},\ldots,v^{[k]}$, which correspond to $k$ largest singular values $\sigma_1\geq \cdots \geq \sigma_k$ of $A_n$.
3. Partition the rows of matrices $\Delta_{1}^{-\frac{1}{2}}\begin{bmatrix} u^{[1]} & \ldots & u^{[k]}\end{bmatrix}$ and $\Delta_{2}^{-\frac{1}{2}}\begin{bmatrix} v^{[1]} & \ldots & v^{[k]}\end{bmatrix}$ with the k-means algorithm.
"

# ╔═╡ 1bc20329-3a6e-435d-a25c-8cce79fb7d9d
md"
## Examples

### Small term-by- document matrix
"

# ╔═╡ b9d0830a-0a33-40c1-a08a-49538cf7a879
begin
	# Make a nicer spy function
	import Plots.spy
	spy(A)=heatmap(A,yflip=true,color=:RdBu,aspectratio=1,clim=(-1,1.0))
end

# ╔═╡ f6492dfe-faac-407f-9a0e-5dcf321a4f2d
begin
	# Sources, targets, and weights
	dn=[6,6,7,6,7,7]
	tn=[1,2,2,3,4,5]
	wn=[3,1,3,2,2,3]
	[dn tn wn]
end

# ╔═╡ 1a8b6f7b-4098-4cd2-adad-8ee788b47333
mynames=["Term 1";"Term 2";"Term 3";"Term 4";"Term 5";"Doc 1";"Doc 2"]

# ╔═╡ 4079aa8e-2bb5-4d51-bca8-04d0997db2bf
begin
	G=Graph(7)
	for i=1:length(dn)
	    add_edge!(G,tn[i],dn[i])
	end
	gplot(G, nodelabel=mynames, edgelabel=wn)
end

# ╔═╡ b6d41ab1-24f2-417d-b604-3e12384e65f8
begin
	# We define some functions
	function WeightMatrix(src::Array,dst::Array,weights::Array)
	    n=nv(G)
	    sparse([src;dst],[dst;src],[weights;weights],n,n)
	end

	Laplacian(W::AbstractMatrix)=spdiagm(0=>vec(sum(W,dims=2)))-W

	function NormalizedLaplacian(L::AbstractMatrix)
	    D=1.0./sqrt.(diag(L))
	    n=length(D)
	    [L[i,j]*(D[i]*D[j]) for i=1:n, j=1:n]
	end
end

# ╔═╡ d5d8ee5d-9cce-45fc-a31b-b2ed8237cc86
W=WeightMatrix(tn,dn,wn)

# ╔═╡ 86fd4f5a-af57-490d-a17f-2339a9414057
Matrix(W)

# ╔═╡ 1f53b5cf-ab6a-4107-b53b-9ce96afcdb7f
begin
	L=Laplacian(W)
	Matrix(L)
end

# ╔═╡ 63403411-cca4-4d71-8a76-8dea8638f69f
Lₙ=NormalizedLaplacian(L)

# ╔═╡ 3cf471bc-d316-4de9-af63-71db9b035c40
begin
	A=W[1:5,6:7]
	Δ₁=sqrt.(sum(A,dims=2))
	Δ₂=sqrt.(sum(A,dims=1))
	Aₙ=[A[i,j]/(Δ₁[i]*Δ₂[j]) for i=1:size(A,1), j=1:size(A,2)]
end

# ╔═╡ ac0b5eda-299c-4ae2-ae78-bda92b378f19
# The partitioning - explain the results!
U,σ,V=svd(Aₙ)

# ╔═╡ 9bf6bbca-f96c-4b38-a4fd-5188ec29a2a2
U[:,2]

# ╔═╡ aa51688c-21a0-40f4-a6d6-01d441bfd4b7
V[:,2]

# ╔═╡ 5f27f1a2-80de-41d1-b1e2-da3a99493217
md"
### Sets of points
"

# ╔═╡ e1149b05-d594-4d75-ae56-1ce18bd0f60d
spy(Matrix(B))

# ╔═╡ b95d4070-c127-4db0-987b-5fc17421b667
# The structure of singular vectors reflects the block
S,rest₀=svds(B,nsv=3);

# ╔═╡ 57d472c1-548d-4678-b884-85f863771aa7
# S is a structure
S.S

# ╔═╡ 4806378b-63cb-445a-8c71-583458fce846
begin
	# Plot the first three left singular vectors
	k=size(B,1)
	x=collect(1:k)
	scatter(x,S.U[:,1],title="Left Singular Vectors",label="U[:,1]")
	scatter!(x,S.U[:,2],label="U[:,2]",legend=:topleft)
	scatter!(x,S.U[:,3],label="U[:,3]")
end

# ╔═╡ a08c99a3-9c10-47be-9dc6-987d92a85087
begin
	# Plot the first three right singular vectors
	scatter(x,S.Vt[1,:],title="Right Singular Vectors",label="V[:,1]")
	scatter!(x,S.Vt[2,:],label="V[:,2]")
	scatter!(x,S.Vt[3,:],label="V[:,3]")
end

# ╔═╡ 6888e8a9-f306-4490-9147-2460fbda4889
begin
	# Add random noise
	noise=sprand(k,k,0.3)
	C=B+noise
	spy(Matrix(C))
end

# ╔═╡ 3bd5ad15-31cd-4efd-b4e3-7cf7ede3edae
begin
	# Apply random permutation to rows and columns of C
	D=C[randperm(k),randperm(k)]
	spy(Matrix(D))
end

# ╔═╡ e4487f35-5de0-492e-96f1-69318d7d8e8b
md"""
__Question.__ Given D, can we recover C?

__Answer.__ Yes (with spectral partitioning)!
"""

# ╔═╡ 81d99f24-82b9-4216-8d17-6a18317beb77
Sₙ,rest=svds(D,nsv=3);

# ╔═╡ 8f925ff7-90b4-4fe4-a869-86dec693c0f4
# K-means on rows of U
outU=kmeans(Matrix(transpose(Sₙ.U)),3)

# ╔═╡ d446724a-97d1-403f-baaf-d0619203f351
# K-means on Vt
outV=kmeans(Sₙ.Vt,3)

# ╔═╡ 7572c241-dbcc-4094-97e3-ee3172255294
sortperm(outU.assignments)

# ╔═╡ 22e90c11-ba4e-4073-9f78-f3b7beb24e69
begin
	# RECOVERY of B
	E=D[sortperm(outU.assignments),sortperm(outV.assignments)]
	spy(Matrix(E))
end

# ╔═╡ Cell order:
# ╟─0f6cd296-6cb8-4a9d-a21f-02b689dfe4e9
# ╟─89299cb7-a4eb-48f3-b3f1-ab2b5c6cfa2e
# ╟─a92d0e70-479d-444e-b06d-e47d2e385b10
# ╟─fddcc2a7-f718-4900-a135-764eb36496a7
# ╟─b3b44560-203a-11eb-2782-170e4c12627a
# ╟─1be90f30-203b-11eb-34da-c3116523ec3a
# ╟─5211ce2e-203b-11eb-0b6f-55e809177fb9
# ╟─38333a20-203c-11eb-2b28-59e07b3e50ee
# ╟─1bc20329-3a6e-435d-a25c-8cce79fb7d9d
# ╠═646a9fc2-3a73-4f8d-a952-2cdcf0e6287a
# ╠═b9d0830a-0a33-40c1-a08a-49538cf7a879
# ╠═b6d41ab1-24f2-417d-b604-3e12384e65f8
# ╠═f6492dfe-faac-407f-9a0e-5dcf321a4f2d
# ╠═1a8b6f7b-4098-4cd2-adad-8ee788b47333
# ╠═4079aa8e-2bb5-4d51-bca8-04d0997db2bf
# ╠═d5d8ee5d-9cce-45fc-a31b-b2ed8237cc86
# ╠═86fd4f5a-af57-490d-a17f-2339a9414057
# ╠═1f53b5cf-ab6a-4107-b53b-9ce96afcdb7f
# ╠═63403411-cca4-4d71-8a76-8dea8638f69f
# ╠═3cf471bc-d316-4de9-af63-71db9b035c40
# ╠═ac0b5eda-299c-4ae2-ae78-bda92b378f19
# ╠═9bf6bbca-f96c-4b38-a4fd-5188ec29a2a2
# ╠═aa51688c-21a0-40f4-a6d6-01d441bfd4b7
# ╟─5f27f1a2-80de-41d1-b1e2-da3a99493217
# ╠═7cd3af86-f459-4e4b-8e35-c29bc6966eed
# ╠═e1149b05-d594-4d75-ae56-1ce18bd0f60d
# ╠═9ab98160-3e37-44c3-b219-a3e1bf95fd68
# ╠═b95d4070-c127-4db0-987b-5fc17421b667
# ╠═57d472c1-548d-4678-b884-85f863771aa7
# ╠═4806378b-63cb-445a-8c71-583458fce846
# ╠═a08c99a3-9c10-47be-9dc6-987d92a85087
# ╠═6888e8a9-f306-4490-9147-2460fbda4889
# ╠═3bd5ad15-31cd-4efd-b4e3-7cf7ede3edae
# ╟─e4487f35-5de0-492e-96f1-69318d7d8e8b
# ╠═81d99f24-82b9-4216-8d17-6a18317beb77
# ╠═8f925ff7-90b4-4fe4-a869-86dec693c0f4
# ╠═d446724a-97d1-403f-baaf-d0619203f351
# ╠═7572c241-dbcc-4094-97e3-ee3172255294
# ╠═22e90c11-ba4e-4073-9f78-f3b7beb24e69

max_stars_repo_path: src/Module C - Applications/L12 Spectral Partitioning of Bipartite Graphs.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
using TightBindingLattice

function make_kagome_lattice(size_matrix ::AbstractMatrix{<:Integer}; compute_symmetry::Bool=false)
    latticevectors = [1 -0.5; 0 0.5*sqrt(3.0)];
    unitcell = make_unitcell(latticevectors, SiteType=String)
    addsite!(unitcell, "A", carte2fract(unitcell, [0.5, 0.0]))
    addsite!(unitcell, "B", carte2fract(unitcell, [0.25, 0.25*sqrt(3.0)]))
    addsite!(unitcell, "C", carte2fract(unitcell, [0.5+0.25, 0.25*sqrt(3.0)]))

    nnbondtypes = [
        ([ 0, 0], "A", [ 0, 0], "B", 1),
        ([ 0, 0], "A", [ 0, 0], "C", 1),
        ([ 0, 0], "B", [ 0, 0], "C", 1),
        ([ 1, 1], "A", [ 1, 0], "B",-1),
        ([ 1, 0], "B", [ 0, 0], "C",-1),
        ([ 0, 0], "C", [ 1, 1], "A",-1),
    ]

    nnnbondtypes = [
        ([ 0, 0], "A", [ 1, 0], "B", 1), # ◁
        ([ 1, 0], "B", [ 0,-1], "C", 1),
        ([ 0,-1], "C", [ 0, 0], "A", 1),
        ([ 0, 0], "C", [ 1, 0], "A",-1),
        ([ 1, 0], "A", [ 0,-1], "B",-1), # ▷
        ([ 0,-1], "B", [ 0, 0], "C",-1),
    ]

    lattice = make_lattice(unitcell, size_matrix)
    orthocube = lattice.orthocube
    supercell = lattice.supercell
    tsym = TranslationSymmetry(lattice)
    psym = project(PointSymmetryDatabase.get(25), [1 0 0; 0 1 0])

    nnbonds = []
    nnnbonds = []

    for r in lattice.bravais_coordinates
        for (rowvec, roworb, colvec, colorb, bondsign) in nnbondtypes
            R_row, r_row = orthocube.wrap(r .+ rowvec)
            R_col, r_col = orthocube.wrap(r .+ colvec)
            roworb_super = (roworb, r_row)
            colorb_super = (colorb, r_col)
            irow = get(supercell.siteindices, roworb_super, -1)
            icol = get(supercell.siteindices, colorb_super, -1)
            push!(nnbonds, ((irow, icol), R_col-R_row, bondsign))
        end
        for (rowvec, roworb, colvec, colorb, bondsign) in nnnbondtypes
            R_row, r_row = orthocube.wrap(r .+ rowvec)
            R_col, r_col = orthocube.wrap(r .+ colvec)
            roworb_super = (roworb, r_row)
            colorb_super = (colorb, r_col)
            irow = get(supercell.siteindices, roworb_super, -1)
            icol = get(supercell.siteindices, colorb_super, -1)
            push!(nnnbonds, ((irow, icol), R_col-R_row, bondsign))
        end
    end

    nn_triangles = []
    for r in lattice.bravais_coordinates
      triangle = []
      for (rowvec, roworb, colvec, colorb, bondsign) in nnbondtypes[1:3]
        R_row, r_row = orthocube.wrap(r .+ rowvec)
        R_col, r_col = orthocube.wrap(r .+ colvec)
        roworb_super = (roworb, r_row)
        colorb_super = (colorb, r_col)
        irow = get(supercell.siteindices, roworb_super, -1)
        icol = get(supercell.siteindices, colorb_super, -1)
        push!(triangle, ((irow, icol), R_col-R_row))
      end
      push!(nn_triangles, (triangle, 1))

      triangle = []
      for (rowvec, roworb, colvec, colorb, bondsign) in nnbondtypes[4:6]
        R_row, r_row = orthocube.wrap(r .+ rowvec)
        R_col, r_col = orthocube.wrap(r .+ colvec)
        roworb_super = (roworb, r_row)
        colorb_super = (colorb, r_col)
        irow = get(supercell.siteindices, roworb_super, -1)
        icol = get(supercell.siteindices, colorb_super, -1)
        push!(triangle, ((irow, icol), R_col-R_row))
      end
      push!(nn_triangles, (triangle, -1))
    end

    return (unitcell=unitcell,
            lattice=lattice,
            translation_symmetry=tsym,
            point_symmetry=psym,
            nearest_neighbor_bonds=nnbonds,
            next_nearest_neighbor_bonds=nnnbonds,
            nearest_neighbor_triangles=nn_triangles)
end


text:
using Plots, DifferentialEquations
gr()

function f′(du, u, p, t)
    du[1] = p.input(t) - p.α*u[1]
    du[2] = -p.β * u[2]
    du[3] = p.β*(p.input(t) - u[3])
end

spike_condition(u,t,integrator) = u[1]-integrator.p.δ

function spike!(integrator)
    integrator.u[1] = 0.0
    integrator.u[2] += 1.0
    push!(integrator.p.spiketimes, integrator.t)
end

u₀ = [0.0,1.0,1.0]
tspan = (0.0,10.0)
p = (
    input=t-> 0.5*(sin(t)+1) + 0.3*(cos(0.7+1.5*t)+1),
    α = 0.5,
    β = 0.5,
    δ = 0.2,
    spiketimes = Float64[0.0]
)

cb = ContinuousCallback(spike_condition, spike!, nothing)
sol = solve(ODEProblem(f′, u₀, tspan, p); callback=cb)

function make_analytical_solution(spike_times; α=1.0, δ=0.2, basis=:sample_and_hold)
    return function analytical_reconstruction(t)
        i = searchsortedfirst(spike_times,t)
        return if i > length(spike_times)
            0.0
        else
            t₁ = spike_times[i]
            t₀ = i<=1 ? -Inf : spike_times[i-1]
            dt = t₁-t₀
            norm = √((1.0-exp(-2.0α*dt))/(2.0α))
            f = δ/norm # fourier coefficient

            if basis == :sample_and_hold
                1.0/√(dt) * f
            else#if basis == :truncated_exponential
                exp(-α*(t₁-t))/norm * f
            end
        end
    end
end

ŝ = make_analytical_solution(p.spiketimes; α=p.α, δ=p.δ, basis=:sample_and_hold)
mean_dt = (tspan[2]-tspan[1])/length(p.spiketimes)
mean_scale = 2*p.δ*p.β/(1+exp(-p.α*mean_dt))

begin
p1=vline(p.spiketimes, linewidth=2, color="#babdb6ff", label="spikes ", xlims=(0,10.5), ylims=(0,1.7))
plot!(p.input, tspan..., linewidth=2, color="#729fcfff", label="input ")
plot!(t->sol(t,idxs=3), tspan..., linewidth=2, color="#fcaf3eff", linestyle=:dash, label="filtered input  ")
plot!(ŝ, tspan..., linewidth=2, color="#729fcfff", linestyle=:dash, label="acausal rec.  ")
plot!(t->sol(t,idxs=2)*mean_scale+exp(-p.β*t)*(1-mean_scale), tspan..., linewidth=2, color="#fcaf3eff", label="causal LTI rec.  ")
plot!(xlabel="time", legend=false)

idx = searchsortedlast(p.spiketimes, 5.0)
(t₀,t₁,t₂) = p.spiketimes[idx:idx+2]
xlims = (t₀-0.1, t₂+0.1)

p2 = vline(p.spiketimes, linewidth=2, color="#babdb6ff", xlims=xlims, ylims=(0,1.0), label="spikes ", legend=false)
plot!(p.input, xlims..., linewidth=2, color="#729fcfff", label="input ")
plot!(t->sol(t,idxs=3), tspan..., linewidth=2, color="#fcaf3eff", linestyle=:dash, label="filtered input ")
plot!(t->sol(t,idxs=1), xlims..., linewidth=2, color="#ad7fa8ff", label="state ")
hline!([p.δ], color="#ef2929ff", linewidth=2, linestyle=:dash, label="threshold ")
plot!(t->sol(t,idxs=2)*mean_scale+exp(-p.β*t)*(1-mean_scale), xlims..., linewidth=2, color="#fcaf3eff", label="causal LTI rec.")
plot!(t->sol(t,idxs=2)*mean_scale+exp(-p.β*t)*(1-mean_scale)-exp(-p.β*(t-t₁))*(t>t₁)*mean_scale, t₁, xlims[end], linestyle=:dot, linewidth=2, color="#fcaf3eff", label="")

p3 = plot(showaxis=false, framestyle=:none, showgrid=false)
plot!([],[], linewidth=2, color="#babdb6ff", label="spikes ")
plot!([],[], linewidth=2, color="#729fcfff", label="input ")
plot!([],[], linewidth=2, color="#fcaf3eff", linestyle=:dash, label="filtered input ")
plot!([],[], linewidth=2, color="#ad7fa8ff", label="state ")
plot!([],[], linewidth=2, color="#ef2929ff", linestyle=:dash, label="threshold ")
plot!([],[], linewidth=2, color="#729fcfff", linestyle=:dash, label="acausal rec.  ")
plot!([],[], linewidth=2, color="#fcaf3eff", label="causal LTI rec. ")
plot!(legend=(0.0, 0.9))

l = @layout [   a{0.7w} b
c{0.5h}     ]

plot(p2,p3,p1, layout=l)
end

savefig("figures/reconstruction_lif.pdf")
savefig("figures/reconstruction_lif.svg")

max_stars_repo_name: jleugeri/phd

max_issues_repo_name: jleugeri/phd

max_forks_repo_name: jleugeri/phd


text:
# This file was generated, do not modify it. # hide
const dim = 16

db = MatrixDatabase(randn(Float32, dim, 10^5))
Q = MatrixDatabase(randn(Float32, dim, 30))
dist = SqL2Distance()
G = SearchGraph(; dist, db)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
default_printer(v) = string(round(v, digits=3))
@recipe(Slider) do scene
    Theme(
        value = 0,
        sliderlength = 200,
        sliderheight = 50,
        backgroundcolor = (:gray, 0.01),
        strokecolor = (:black, 0.4),
        strokewidth = 1,
        textcolor = :black,
        slidercolor = (:gray, 0.6),
        buttoncolor = :white,
        buttonsize = 15,
        buttonstroke = 1.5,
        textsize = 15,
        buttonstrokecolor = :black,
        valueprinter = default_printer
    )
end

mouseover() = error("not implemented")
export mouseover
convert_arguments(::Type{<: Slider}, x::AbstractRange) = (x,)

function range_label_bb(tplot, printer_func, range)
    bb = boundingbox(tplot, printer_func(first(range)))
    for elem in Iterators.drop(range, 1)
        bb = union(bb, boundingbox(tplot, printer_func(elem)))
    end
    bb
end
function plot!(slider::Slider)
    @extract(slider, (
        backgroundcolor, strokecolor, strokewidth, slidercolor, buttonstroke,
        buttonstrokecolor, buttonsize, buttoncolor, valueprinter,
        sliderlength, sliderheight, textcolor, textsize
    ))
    range = slider[1]
    val = slider[:value]
    push!(val, first(value(range)))
    label = lift((v, f)-> f(v), val, valueprinter)
    lplot = text!(
        slider, label,
        textsize = textsize,
        align = (:left, :center), color = textcolor,
        position = map((w, h)-> Point2f0(w, h/2), sliderlength, sliderheight)
    ).plots[end]
    lbb = lift(range_label_bb, Node(lplot), valueprinter, range)
    bg_rect = lift(sliderlength, sliderheight, lbb) do w, h, bb
        FRect(0, 0, w + 10 + widths(bb)[1], h)
    end
    poly!(
        slider, bg_rect,
        color = backgroundcolor, linecolor = strokecolor,
        linewidth = strokewidth
    )
    line = lift(sliderlength, sliderheight) do w, h
        Point2f0[(10, h / 2), (w - 10, h / 2)]
    end

    linesegments!(slider, line, color = slidercolor)
    button = scatter!(
        slider, map(x-> x[1:1], line),
        markersize = buttonsize, color = buttoncolor, strokewidth = buttonstroke,
        strokecolor = buttonstrokecolor
    ).plots[end]
    dragslider(slider, button)
end

function dragslider(slider, button)
    mpos = events(slider).mouseposition
    drag_started = Ref(false)
    startpos = Base.RefValue((0.0, 0.0))
    range = slider[1]
    @extract slider (value, sliderlength)
    foreach(events(slider).mousedrag) do drag
        if drag == Mouse.down && mouseover(slider, button)
            startpos[] = mpos[]
            drag_started[] = true
        elseif drag == Mouse.pressed && drag_started[]
            diff = startpos[] .- mpos[]
            startpos[] = mpos[]
            spos = translation(button)[][1] - diff[1]
            l = sliderlength[] - 17.5
            if spos >= 0 && spos <= l
                idx = round(Int, ((spos / l) .* (length(range[]) - 1)) + 1)
                value[] = range[][idx]
                translate!(button, spos, 0, 0)
            end
        else
            drag_started[] = false
            startpos[] = (0.0, 0.0)
        end
        return
    end
end

function move!(x::Slider, idx::Integer)
    r = x[1][]
    len = x[:sliderlength][]
    x[:value] = r[idx]
    xpos = ((idx - 1) / (length(r) - 1)) * len
    translate!(x.plots[end], xpos, 0, 0)
    return
end
export move!
@recipe(Button) do scene
    Theme(
        dimensions = (40, 40),
        backgroundcolor = (:white, 0.4),
        strokecolor = (:black, 0.4),
        strokewidth = 1,
        textcolor = :black,
        clicks = 0
    )
end

function button(func::Function, scene::Scene, txt; kw_args...)
    b = button!(scene, txt; kw_args...)[end]
    foreach(b[:clicks]) do clicks
        func(clicks)
        return
    end
    b
end

function plot!(splot::Button)
    @extract(splot, (
        backgroundcolor, strokecolor, strokewidth,
        dimensions, textcolor, clicks
    ))
    txt = splot[1]
    lplot = text!(
        splot, txt,
        align = (:center, :center), color = textcolor,
        textsize = 15,
        position = map((wh)-> Point2f0(wh./2), dimensions)
    ).plots[end]
    lbb = boundingbox(lplot) # on purpose static so we hope text won't become too long?
    # bg_rect = lift(dimensions) do wh
    #     IRect(0, 0, Vec(wh))
    # end
    # p = poly!(
    #     splot, bg_rect,
    #     color = backgroundcolor, linecolor = strokecolor,
    #     linewidth = strokewidth
    # )
    foreach(events(splot).mousebuttons) do mb
        if ispressed(mb, Mouse.left) && mouseover(parent(splot), lplot)
            clicks[] = clicks[] + 1
        end
        return
    end
    splot
end

window_open(scene::Scene) = getscreen(scene) != nothing && isopen(getscreen(scene))

function playbutton(f, scene, range, rate = (1/30))
    b = button(scene, "▶ ")
    isplaying = Ref(false)
    play_idx = Ref(1)
    foreach(b[:clicks]) do x
        if !isplaying[] && x > 0 # check that this isn't before any clicks
            isplaying[] = true
            @async begin
                b.plots[1][1][] = "𝅛𝅛"
                tstart = time()
                while (isplaying[] && window_open(scene))
                    if time() - tstart >= rate
                        f(range[play_idx[]])
                        play_idx[] = mod1(play_idx[] + 1, length(range))
                        tstart = time()
                        force_update!()
                    end
                    yield()
                end
                isplaying[] = false
            end
        else
            b.plots[1][1][] = "▶ "
            isplaying[] = false
        end
        nothing
    end
    b
end

struct Popup
    scene::Scene
    open::Node
    position::Node{Point2f0}
    width::Node{Point2f0}
end

function hbox!(plots::Vector{T}; kw_args...) where T <: AbstractPlot
    N = length(plots)
    h = 0.0
    for idx in 1:N
        p = plots[idx]
        translate!(p, 0.0, h, 0.0)
        swidth = widths(boundingbox(p))
        h += (swidth[2] * 1.2)
    end
end

function textslider(ui, range, label)
    a = slider!(ui, range, raw = true)[end][:value]
    text!(ui, "$label:", raw = true, align = (:left, :center))
    a
end

function sample_color(f, ui, colormesh, v)
    mpos = ui.events.mouseposition
    sub = Scene(ui, transformation = Transformation(), px_area = pixelarea(ui), theme = theme(ui))
    select = scatter!(
        sub, lift((x, r)-> [Point2f0(x) .- minimum(r)], mpos, pixelarea(sub)),
        markersize = 10, color = (:white, 0.2), strokecolor = :white,
        strokewidth = 5, visible = lift(identity, theme(ui, :visible)), raw = true
    )[end]
    translate!(select, 0, 0, 10)

    foreach(mpos, ui.events.mousebuttons) do mp, mb
        bb = FRect2D(modelmatrix(sub)[] * modelmatrix(colormesh)[] * boundingbox(colormesh))
        mp = Point2f0(mp) .- minimum(pixelarea(sub)[])
        if Point2f0(mp) in bb
            select[:visible] = true
            if ispressed(mb, Mouse.left)
                sv = (Point2f0(mp) .- minimum(bb)) ./ widths(bb)
                f(RGBAf0(HSV(v[], sv[1], sv[2])))
            end
        else
            select[:visible] = false
        end
        return
    end
end

function popup(parent, position, width)
    pos_n = Node(Point2f0(position))
    width_n = Node(Point2f0(width))

    harea = lift(pos_n, width_n) do p, wh
        IRect(p[1], p[2] + wh[2] - 20, wh[1], 20)
    end
    parea = lift(pos_n, width_n) do p, wh
        IRect(Point2f0(p), Point2f0(wh))
    end
    popup = Scene(parent, parea)
    theme(popup)[:visible] = Node(false)
    header = Scene(popup, harea)
    theme(popup)[:plot] = Attributes(raw = true)
    campixel!(popup)
    campixel!(header)
    theme(header)[:plot] = Attributes(raw = true)
    theme(header)[:visible] = theme(popup, :visible)
    initialized = Ref(false)
    but = button(header, "x") do click
        if initialized[]
            theme(popup, :visible)[] = !theme(popup, :visible)[]
        else
            initialized[] = true
        end
        return
    end
    foreach(width_n) do wh
        translate!(but, wh[1] - 30, -10, 120)
    end
    poly!(header, lift(r-> FRect(0, 0, widths(r)), harea), color = (:gray, 0.1))
    poly!(popup, lift(wh-> FRect(2, 2, (wh - 4)...), width_n), color = :white, strokecolor = :black, strokewidth = 2)
    Popup(Scene(popup, theme = theme(popup)), theme(popup, :visible), pos_n, width_n)
end

function mouse_selection end
export mouse_selection

function colorswatch(ui)
    pop = popup(ui, (50, 50), (250, 300))
    sub_ui = pop.scene
    hsv_hue = textslider(sub_ui, 1:360, "hue")
    colors = lift(hsv_hue) do V
        [HSV(V, 0, 0), HSV(V, 1, 0), HSV(V, 1, 1), HSV(V, 0, 1)]
    end
    S = 200
    colormesh = mesh!(sub_ui,
        # TODO implement decompose correctly to just have this be IRect(0, 0, S, S)
        [(0, 0), (S, 0), (S, S), (0, S)],
        [1, 2, 3, 3, 4, 1],
        color = colors, raw = true, shading = false
    )[end]
    color = Node(RGBAf0(0,0,0,1))
    sample_color(sub_ui, colormesh, hsv_hue) do c
        color[] = c
    end
    hbox!(sub_ui.plots)
    translate!(sub_ui, 10, 0, 0)
    rect = IRect(0, 0, 50, 50)
    swatch = poly!(ui, rect, color = color, raw = true, visible = true)[end]
    pop.open[] = false
    foreach(ui.events.mousebuttons) do mb
        if ispressed(mb, Mouse.left)
            plot, idx = mouse_selection(ui)
            if plot in swatch.plots
                pop.position[] = Point2f0(events(ui).mouseposition[] .+ 50)
                pop.open[] = true
            end
        end
        return
    end
    hbox!(sub_ui.plots)
    color, pop
end


text:
using Test
using LinearAlgebra
using PolynomialBasis
using IntervalArithmetic
using Revise
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature

function corner_distance_function(x::V,xc) where {V<:AbstractVector}
    v = xc - x
    if all(x .<= xc)
        minimum(v)
    elseif all(x.> xc)
        return -sqrt(v'*v)
    elseif x[2] > xc[2]
        return v[2]
    else
        return v[1]
    end
end

function corner_distance_function(points::M,xc) where {M<:AbstractMatrix}
    return vec(mapslices(x->corner_distance_function(x,xc),points,dims=1))
end

corner = [0.,0.]
polyorder = 2
numqp = 5
poly = InterpolatingPolynomial(1,2,polyorder)
coeffs = corner_distance_function(poly.basis.points,corner)
update!(poly,coeffs)

quad = area_quadrature(poly,+1,[-1.,-1.],[1.,1.],numqp,numsplits=3)
squad = surface_quadrature(poly,[-1.,-1.],[1.,1.],numqp,numsplits=3)

using Plots
xrange = -1:1e-2:1
contour(xrange,xrange,(x,y)->poly(x,y),levels=[0.0],aspect_ratio=:equal)
scatter!(quad.points[1,:],quad.points[2,:])
scatter!(squad.points[1,:],squad.points[2,:])

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
using Test, IonSim, IonSim.analytical

@testset "analytic_functions -- MolmerSorensen" begin
    t = 0:100:10000
    ν = 1
    Ω = 0.025
    δ = 0.95
    η = 0.1
    n̄ = 0

    # test off-resonant limit: ηΩ << (ν - δ)
    # should get rabi oscillations between |ee> and |gg> with Ω̃ = (ηΩ)^2 / (ν - δ)
    ρgg, ρee = two_ion_ms(t, Ω, ν, δ, η, n̄)
    @test ρee ≈ @.(sin(2π * (η * Ω)^2 / 2(ν - δ) * t)^2)

    # test near resonant limit
    # if we set Ω = (ν-δ) / (2η√K); K<:Int, we should get ρee=ρgg=0.5 at time τ = √K / (2ηΩ)
    Ω = (ν - δ) / (2η)
    τ = 1 / (2η * Ω)
    ρgg, ρee = two_ion_ms([τ], Ω, ν, δ, η, n̄)
    @test ρgg[1] + ρee[1] ≈ 1.0

    Ω = (ν - δ) / (2η * √2)
    τ = √2 / (2η * Ω)
    ρgg, ρee = two_ion_ms([τ], Ω, ν, δ, η, n̄)
    @test ρgg[1] + ρee[1] ≈ 1.0

    # test hot gate
    Ω = (ν - δ) / (2η)
    τ = 1 / (2η * Ω)
    n̄ = 50
    ρgg, ρee = two_ion_ms([τ/2, τ], Ω, ν, δ, η, n̄)
    @test ρgg[2] + ρee[2] ≈ 1.0 rtol=1e-2
    @test ρgg[1] + ρee[1] ≈ 0.75 rtol=1e-2
end

@testset "analyticfunctions -- RabiFlop" begin
    # # test carrier with motion in ground state
    t = 0:0.01:1
    Ω = 1
    n̄ = 0
    η = 0
    @test rabi_flop(t, Ω , η, n̄) ≈ @.(sin(2π * t / 2)^2) rtol=1e-5

    # test bluesideband, Ω̃ ≈ ηΩ
    η = 0.1
    @test rabi_flop(t, Ω , η, n̄, s=1) ≈ @.(sin(2π * η * t / 2)^2) rtol=1e-2
end


text:
#!/usr/bin/env julia

## brillouin.jl
#
#    Author: Leandro Seixas <seixasle@gmail.com>
#            Centre for Advanced 2D Materials (CA2DM),
#            National University of Singapore (NUS)
#            Singapore
#

using YAML
#using DataFrames
#using JulieTest
#using Gadfly

include("src/monkhorst-pack.jl")
include("src/reciprocal.jl")
include("src/firstBZ.jl")

lattice_constant = 2.47

bravais = [ sqrt(3)/2  1/2 0.0;
            sqrt(3)/2 -1/2 0.0;
            0.0  0.0 10.0 ]

bravais_lattice = lattice_constant*bravais

reciprocal_lattice = 2π*inv(bravais_lattice)'

G( n::Integer ) = reciprocal_lattice[n,:]

g_neighbors = collect( [ x*G(1)+y*G(2)+z*G(3) for x=-1:1:1, y=-1:1:1, z=-1:1:1 ] )

#for i in 1:length(kvec)
#  @printf("%14.8f%14.8f%14.8f\n", kvec[i]...)
#end


text:
using JuMP
using Cbc
# cria o modelo para a otimização linear, utiliza o cbc como solver
m = Model(with_optimizer(Cbc.Optimizer))

set_time_limit_sec(m , 1)

n = 6

# define as variaveis binárias xij e o u como inteiro
@variable(m , x[1:n , 1:n] , Bin)
@variable(m , u[1:n] , integer=true)

fix(u[1] , 0)

# cria o grafo de pesos para fazer conexão
y = zeros(Int64 , n , n)

# dados do exemplo do slide da prof(trabalho -  material auxiliar)
# é o exemplo das cidades que está logo depois do MTZ
y = [ 0  186 105 208 177 94 ;
     187 0 89.8 223 255 254 ;
     99.9 89 0  203 204 172 ;
     206 220 203 0 377 295 ;
     168 251 201 376 0 156 ;
     86.8 255 173  293 159  0]

#for i in 1:5
#    for j in 1:5
#        y[i,j] = i+j
#    end
#end

# função objetiva
@objective(m , Min , sum(x.*y))

# constraint de que a soma de uma linha vai ser 1
@constraint(m , [i= 1:n] , sum(x[i,1:end .!=i]) == 1)
# constraint de que a soma de uma coluna vai ser 1
@constraint(m , [i= 1:n] , sum(x[1:end.!=i,i]) == 1)

# constraint do u
@constraint(m , [i = 1:n , j =2:n , i != j] , u[j] >= u[i] + x[i,j] - n *(1-x[i,j]))

#printa os dados do modelo
println(m)

# resolve e printa os resultados
optimize!(m)

uAux = value.(u)
xAux = value.(x)

for i in 1:n
    println("u$(i) $(uAux[i])")
end

print("\n\n Xij:\n")

for i in 1:n
    for j in 1:n
        print("$(xAux[i,j]) ")
    end
    println()
end

max_stars_repo_name: jvlsg/julia-tsp

max_issues_repo_name: jvlsg/julia-tsp

max_forks_repo_name: jvlsg/julia-tsp


text:
config = MOI.Test.Config()

@testset "Model" begin
    optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 1.0, 0.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1.0], [-1.0]])
    mock = bridged_mock(optimize!)
    Tests.zero_polynomial_test(mock, config)
end

# The VectorOfVariables-in-ZeroPolynomialSet is bridged by VectorFunctionizeBridge
# since the free variable is bridged. This tests that the MomentsAttribute is
# passed by the VectorFunctionizeBridge.
@testset "NoFreeVariable" begin
    optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 0.0, 1.0, 0.0, 0.0, 0.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1.0], [-1.0]])
    mock = bridged_mock(optimize!, model = NoFreeVariable{Float64}())
    Tests.zero_polynomial_test(mock, config)
end


text:
"""
    build_sysimg(sysimg_path=default_sysimg_path(), cpu_target="native", userimg_path=nothing; force=false)

Rebuild the system image. Store it in `sysimg_path`, which defaults to a file named `sys.ji`
that sits in the same folder as `libjulia.{so,dylib}`, except on Windows where it defaults
to `Sys.BINDIR/../lib/julia/sys.ji`. Use the cpu instruction set given by `cpu_target`.
Valid CPU targets are the same as for the `-C` option to `julia`, or the `-march` option to
`gcc`. Defaults to `native`, which means to use all CPU instructions available on the
current processor. Include the user image file given by `userimg_path`, which should contain
directives such as `using MyPackage` to include that package in the new system image. New
system image will not replace an older image unless `force` is set to true.
"""
function build_sysimg(
        sysimg_path, userimg_path = nothing;
        verbose = false, quiet = false, release = false,
        home = nothing, startup_file = nothing, handle_signals = nothing,
        sysimage_native_code = nothing, compiled_modules = nothing,
        depwarn = nothing, warn_overwrite = nothing,
        compile = nothing, cpu_target = nothing, optimize = nothing, debug = nothing,
        inline = nothing, check_bounds = nothing, math_mode = nothing,
        cc = nothing, cc_flags = nothing
    )
    # build vanilla backup system image
    clean_sysimg = get_backup!(occursin("debug", basename(Base.julia_cmd().exec[1])), cpu_target)
    static_julia(
        userimg_path, verbose = verbose, quiet = quiet, builddir = sysimg_path, outname = "sys",
        autodeps = true, shared = true, release = release,
        sysimage = clean_sysimg, home = home, startup_file = startup_file, handle_signals = handle_signals,
        sysimage_native_code = sysimage_native_code, compiled_modules = compiled_modules,
        depwarn = depwarn, warn_overwrite = warn_overwrite,
        compile = compile, cpu_target = cpu_target, optimize = optimize, debug = debug,
        inline = inline, check_bounds = check_bounds, math_mode = math_mode,
        cc = cc, cc_flags = cc_flags
    )
end

"""
    build_shared_lib(
        julia_program, output_name = nothing;
        snoopfile = nothing, builddir = nothing, verbose = false, quiet = false,
        init_shared = false, copy_julialibs = true, copy_files = nothing, release = false, Release = false,
        sysimage = nothing, home = nothing, startup_file = nothing, handle_signals = nothing,
        sysimage_native_code = nothing, compiled_modules = nothing,
        depwarn = nothing, warn_overwrite = nothing,
        compile = nothing, cpu_target = nothing, optimize = nothing, debug = nothing,
        inline = nothing, check_bounds = nothing, math_mode = nothing,
        cc = nothing, cc_flags = nothing
    )
    `julia_program` needs to be a Julia script containing a `julia_main` function, e.g. like `examples/hello.jl`
    `snoopfile` is optional and can be a Julia script which calls functions that you want to make sure to have precompiled
    `builddir` is where the compiled artifacts will end up
"""
function build_shared_lib(
        julia_program, output_name = nothing;
        snoopfile = nothing, builddir = nothing, verbose = false, quiet = false,
        init_shared = false, copy_julialibs = true, copy_files = nothing, release = false, Release = false,
        sysimage = nothing, home = nothing, startup_file = nothing, handle_signals = nothing,
        sysimage_native_code = nothing, compiled_modules = nothing,
        depwarn = nothing, warn_overwrite = nothing,
        compile = nothing, cpu_target = nothing, optimize = nothing, debug = nothing,
        inline = nothing, check_bounds = nothing, math_mode = nothing,
        cc = nothing, cc_flags = nothing
    )
    static_julia(
        julia_program, verbose = verbose, quiet = quiet,
        builddir = builddir, outname = output_name, snoopfile = snoopfile, autodeps = true, shared = true,
        init_shared = init_shared, copy_julialibs = copy_julialibs, copy_files = copy_files, release = release, Release = release,
        sysimage = sysimage, home = home, startup_file = startup_file, handle_signals = handle_signals,
        sysimage_native_code = sysimage_native_code, compiled_modules = compiled_modules,
        depwarn = depwarn, warn_overwrite = warn_overwrite,
        compile = compile, cpu_target = cpu_target, optimize = optimize, debug = debug,
        inline = inline, check_bounds = check_bounds, math_mode = math_mode,
        cc = cc, cc_flags = cc_flags
    )
end

"""
    build_executable(
        julia_program, output_name = nothing, c_program = nothing;
        snoopfile = nothing, builddir = nothing, verbose = false, quiet = false,
        copy_julialibs = true, copy_files = nothing, release = false, Release = false,
        sysimage = nothing, home = nothing, startup_file = nothing, handle_signals = nothing,
        sysimage_native_code = nothing, compiled_modules = nothing,
        depwarn = nothing, warn_overwrite = nothing,
        compile = nothing, cpu_target = nothing, optimize = nothing, debug = nothing,
        inline = nothing, check_bounds = nothing, math_mode = nothing,
        cc = nothing, cc_flags = nothing
    )
    `julia_program` needs to be a Julia script containing a `julia_main` function, e.g. like `examples/hello.jl`
    `snoopfile` is optional and can be a Julia script which calls functions that you want to make sure to have precompiled
    `builddir` is where the compiled artifacts will end up
"""
function build_executable(
        julia_program, output_name = nothing, c_program = nothing;
        snoopfile = nothing, builddir = nothing, verbose = false, quiet = false,
        copy_julialibs = true, copy_files = nothing, release = false, Release = false,
        sysimage = nothing, home = nothing, startup_file = nothing, handle_signals = nothing,
        sysimage_native_code = nothing, compiled_modules = nothing,
        depwarn = nothing, warn_overwrite = nothing,
        compile = nothing, cpu_target = nothing, optimize = nothing, debug = nothing,
        inline = nothing, check_bounds = nothing, math_mode = nothing,
        cc = nothing, cc_flags = nothing
    )
    static_julia(
        julia_program, cprog = c_program, verbose = verbose, quiet = quiet,
        builddir = builddir, outname = output_name, snoopfile = snoopfile, autodeps = true, executable = true,
        copy_julialibs = copy_julialibs, copy_files = copy_files, release = release, Release = release,
        sysimage = sysimage, home = home, startup_file = startup_file, handle_signals = handle_signals,
        sysimage_native_code = sysimage_native_code, compiled_modules = compiled_modules,
        depwarn = depwarn, warn_overwrite = warn_overwrite,
        compile = compile, cpu_target = cpu_target, optimize = optimize, debug = debug,
        inline = inline, check_bounds = check_bounds, math_mode = math_mode,
        cc = cc, cc_flags = cc_flags
    )
end

"""
    force_native_image!()
Builds a clean system image, similar to a fresh Julia install.
Can also be used to build a native system image for a downloaded cross compiled julia binary.
"""
function force_native_image!(debug = false) # debug is ignored right now
    sysimg = get_backup!(debug, "native")
    copy_system_image(dirname(sysimg), default_sysimg_path(debug))
end


text:
####################################################
"""
    plt = plotstopo(topo::VisClaw.Topo; kwargs...)

    plotstopo!(plt::Plots.Plot, topo::VisClaw.Topo; kwargs...)

plot topography and bathymetry data using Plots
"""
function plotstopo!(plt, topo::VisClaw.Topo; kwargs...)
	# parse keyword args
    kwdict = KWARG(kwargs)
	seriestype, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_seriestype, :heatmap)
	xlims, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_xlims, extrema(topo.x))
	ylims, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_ylims, extrema(topo.y))

	# plot
    plt = Plots.plot!(plt, topo.x, topo.y, topo.elevation; seriestype=seriestype,
                      xlims=xlims, ylims=ylims, axis_ratio=:equal, kwdict...)
    # return
    return plt
end
####################################################
"""
$(@doc plotstopo!)
"""
plotstopo(topo::VisClaw.Topo; kwargs...) = plotstopo!(Plots.plot(), topo; kwargs...)

####################################################
"""
    plt = plotsdtopo(dtopo::VisClaw.DTopo, itime::Integer=0; kwargs...)

    plotsdtopo!(plt::Plots.Plot, dtopo::VisClaw.DTopo, itime::Integer=0; kwargs...)

plot dtopo data using Plots
"""
function plotsdtopo!(plt, dtopo::VisClaw.DTopo, itime::Integer=0; kwargs...)

    ( (itime < 0) || (dtopo.mt < itime) ) && error("Invalid time")

    if dtopo.mt==1;     z = dtopo.deform
    elseif  itime == 0; z = dtopo.deform[:,:,end]
    else;               z = dtopo.deform[:,:,itime]
    end

	# parse keyword args
    kwdict = KWARG(kwargs)
	seriestype, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_seriestype, :heatmap)

    # plot
    plt = Plots.plot!(plt, dtopo.x, dtopo.y, z; seriestype=seriestype,
                      xlims=extrema(dtopo.x), ylims=extrema(dtopo.y),
                      axis_ratio=:equal, kwdict...)
    # return
    return plt
end
####################################################
"""
$(@doc plotsdtopo!)
"""
plotsdtopo(dtopo::VisClaw.DTopo, itime::Integer=0; kwargs...) = plotsdtopo!(Plots.plot(), dtopo, itime; kwargs...)

####################################################
"""
    plt = plotstoporange(geo::VisClaw.AbstractTopo; kwargs...)

    plotstoporange!(plt::Plots.Plot, geo::VisClaw.AbstractTopo; kwargs...)

plot a range of topo/bath using Plots
"""
function plotstoporange!(plt, geo::VisClaw.AbstractTopo; kwargs...)
	# parse keyword args
    kwdict = KWARG(kwargs)
	label, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_label, "")

	xp = [geo.x[1],  geo.x[1]  , geo.x[end], geo.x[end], geo.x[1]]
	yp = [geo.y[1],  geo.y[end], geo.y[end], geo.y[1]  , geo.y[1]]

	# plot
	plt = Plots.plot!(plt, xp, yp; label=label, kwdict...)

	return plt
end
####################################################
"""
$(@doc plotstoporange!)
"""
plotstoporange(geo::VisClaw.AbstractTopo; kwargs...) = plotstoporange!(Plots.plot(), geo; kwargs...)
####################################################

####################################################
"""
    plt = plotscoastline(topo::VisClaw.Topo; kwargs...)
	plt = plotscoastline(topo::VisClaw.Topo, wetdry::VisClaw.Topo; kwargs...)

    plotscoastline!(plt::Plots.Plot, topo::VisClaw.Topo; kwargs...)
    plotscoastline!(plt::Plots.Plot, topo::VisClaw.Topo, wetdry::VisClaw.Topo; kwargs...)

plot coastlines from topography and bathymetry data using Plots
"""
function plotscoastline!(plt, topo::VisClaw.Topo; kwargs...)
	# parse keyword args
    kwdict = KWARG(kwargs)
	label, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_label, "")
	xlims, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_xlims, extrema(topo.x))
	ylims, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_ylims, extrema(topo.y))

	# plot
	plt = Plots.contour!(plt, topo.x, topo.y, topo.elevation; levels=[0], label=label,
	                     xlims=xlims, ylims=ylims, axis_ratio=:equal, kwdict...)
	return plt
end
####################################################
"""
$(@doc plotscoastline!)
"""
plotscoastline(topo::VisClaw.Topo; kwargs...) = plotscoastline!(Plots.plot(), topo; kwargs...)
####################################################

####################################################
function plotscoastline!(plt, topo::VisClaw.Topo, wetdry::VisClaw.Topo; kwargs...)
	## check args
	(size(topo.elevation) == size(wetdry.elevation)) || error("size of wetdry array should correspond to that of topo array.")

    ## make a new VisClaw.Topo for coastline plot
	dry = findall(convert(BitArray, wetdry.elevation))
    topotmp = copy(topo.elevation)
	topotmp[dry] .= 1e3
	topo_coastline = VisClaw.Topo(topo.ncols, topo.nrows, topo.x, topo.y, topo.dx, topo.dy, topotmp)

    ## plot
	plt = plotscoastline!(plt, topo_coastline; kwargs...)
	return plt
end
####################################################
plotscoastline(topo::VisClaw.Topo, wetdry::VisClaw.Topo; kwargs...) = plotscoastline!(Plots.plot(), topo, wetdry; kwargs...)
####################################################


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# test/variable.jl
# Testing for VariableRef
#############################################################################

module TestVariable

using JuMP
import LinearAlgebra: Symmetric
using Test

include(joinpath(@__DIR__, "utilities.jl"))
include(joinpath(@__DIR__, "JuMPExtension.jl"))

function _test_variable_name_util(variable, name)
    @test name == @inferred JuMP.name(variable)
    @test variable == JuMP.variable_by_name(JuMP.owner_model(variable), name)
end

# Slices three-dimensional DenseAxisArray x[I,J,K]
# I,J,K can be singletons, ranges, colons, etc.
function _sliceof_util(VariableRefType, x, I, J, K)
    y = Array{VariableRefType}(undef, length(I), length(J), length(K))

    ii = 1
    jj = 1
    kk = 1
    for i in I
        for j in J
            for k in K
                y[ii, jj, kk] = x[i, j, k]
                kk += 1
            end
            jj += 1
            kk = 1
        end
        ii += 1
        jj = 1
    end
    idx = [length(I) == 1, length(J) == 1, length(K) == 1]
    return dropdims(y, dims = tuple(findall(idx)...))
end

function test_variable_no_bound(ModelType, VariableRefType)
    model = ModelType()
    @variable(model, nobounds)
    @test !JuMP.has_lower_bound(nobounds)
    @test !JuMP.has_upper_bound(nobounds)
    @test !JuMP.is_fixed(nobounds)
    _test_variable_name_util(nobounds, "nobounds")
    @test zero(nobounds) isa JuMP.GenericAffExpr{Float64,VariableRefType}
    @test one(nobounds) isa JuMP.GenericAffExpr{Float64,VariableRefType}
end

function test_variable_lower_bound_rhs(ModelType, ::Any)
    model = ModelType()
    @variable(model, lbonly >= 0, Bin)
    @test JuMP.has_lower_bound(lbonly)
    @test 0.0 == @inferred JuMP.lower_bound(lbonly)
    @test !JuMP.has_upper_bound(lbonly)
    @test !JuMP.is_fixed(lbonly)
    @test JuMP.is_binary(lbonly)
    @test !JuMP.is_integer(lbonly)
    @test isequal(model[:lbonly], lbonly)
    JuMP.delete_lower_bound(lbonly)
    @test !JuMP.has_lower_bound(lbonly)
    # Name already used
    @test_throws ErrorException @variable(model, lbonly)
end

function test_variable_lower_bound_lhs(ModelType, ::Any)
    model = ModelType()
    @variable(model, 0 <= lblhs, Bin)
    @test JuMP.has_lower_bound(lblhs)
    @test 0.0 == @inferred JuMP.lower_bound(lblhs)
    @test !JuMP.has_upper_bound(lblhs)
    @test !JuMP.is_fixed(lblhs)
    @test JuMP.is_binary(lblhs)
    @test !JuMP.is_integer(lblhs)
    @test isequal(model[:lblhs], lblhs)
end

function test_variable_upper_bound_rhs(ModelType, ::Any)
    model = ModelType()
    @variable(model, ubonly <= 1, Int)
    @test !JuMP.has_lower_bound(ubonly)
    @test JuMP.has_upper_bound(ubonly)
    @test 1.0 == @inferred JuMP.upper_bound(ubonly)
    @test !JuMP.is_fixed(ubonly)
    @test !JuMP.is_binary(ubonly)
    @test JuMP.is_integer(ubonly)
    @test isequal(model[:ubonly], ubonly)
    JuMP.delete_upper_bound(ubonly)
    @test !JuMP.has_upper_bound(ubonly)
end

function test_variable_upper_bound_lhs(ModelType, ::Any)
    model = ModelType()
    @variable(model, 1 >= ublhs, Int)
    @test !JuMP.has_lower_bound(ublhs)
    @test JuMP.has_upper_bound(ublhs)
    @test 1.0 == @inferred JuMP.upper_bound(ublhs)
    @test !JuMP.is_fixed(ublhs)
    @test !JuMP.is_binary(ublhs)
    @test JuMP.is_integer(ublhs)
    @test isequal(model[:ublhs], ublhs)
end

function test_variable_interval(ModelType, ::Any)
    function has_bounds(var, lb, ub)
        @test JuMP.has_lower_bound(var)
        @test lb == @inferred JuMP.lower_bound(var)
        @test JuMP.has_upper_bound(var)
        @test ub == @inferred JuMP.upper_bound(var)
        @test !JuMP.is_fixed(var)
    end
    model = ModelType()
    @variable(model, 0 <= bothb1 <= 1)
    has_bounds(bothb1, 0.0, 1.0)
    @variable(model, 0 ≤ bothb2 ≤ 1)
    has_bounds(bothb2, 0.0, 1.0)
    @variable(model, 1 >= bothb3 >= 0)
    has_bounds(bothb3, 0.0, 1.0)
    @variable(model, 1 ≥ bothb4 ≥ 0)
    has_bounds(bothb4, 0.0, 1.0)
    @test_macro_throws ErrorException @variable(model, 1 ≥ bothb5 ≤ 0)
    @test_macro_throws ErrorException @variable(model, 1 ≤ bothb6 ≥ 0)
end

function test_variable_fix(ModelType, ::Any)
    model = ModelType()
    @variable(model, fixed == 1.0)
    @test !JuMP.has_lower_bound(fixed)
    @test !JuMP.has_upper_bound(fixed)
    @test JuMP.is_fixed(fixed)
    @test 1.0 == @inferred JuMP.fix_value(fixed)
    JuMP.unfix(fixed)
    @test !JuMP.is_fixed(fixed)
    JuMP.set_lower_bound(fixed, 0.0)
    @test_throws Exception JuMP.fix(fixed, 1.0)
    JuMP.fix(fixed, 1.0; force = true)
    @test !JuMP.has_lower_bound(fixed)
    @test !JuMP.has_upper_bound(fixed)
    @test JuMP.is_fixed(fixed)
    @test 1.0 == @inferred JuMP.fix_value(fixed)
end

function test_variable_custom_index_sets(ModelType, ::Any)
    model = ModelType()
    @variable(model, onerangeub[-7:1] <= 10, Int)
    @variable(model, manyrangelb[0:1, 10:20, 1:1] >= 2)
    @test JuMP.has_lower_bound(manyrangelb[0, 15, 1])
    @test 2 == @inferred JuMP.lower_bound(manyrangelb[0, 15, 1])
    @test !JuMP.has_upper_bound(manyrangelb[0, 15, 1])

    s = ["Green", "Blue"]
    @variable(model, x[i = -10:10, s] <= 5.5, Int, start = i + 1)
    @test 5.5 == @inferred JuMP.upper_bound(x[-4, "Green"])
    _test_variable_name_util(x[-10, "Green"], "x[-10,Green]")
    @test JuMP.start_value(x[-3, "Blue"]) == -2
    @test isequal(model[:onerangeub][-7], onerangeub[-7])
    @test_throws KeyError model[:foo]
end

function test_variable_anonymous(ModelType, ::Any)
    model = ModelType()
    @test_throws ErrorException @variable(model, [(0, 0)])  # #922
    x = @variable(model, [(0, 2)])
    @test "" == @inferred JuMP.name(x[0])
    @test "" == @inferred JuMP.name(x[2])
end

function test_variable_is_valid_delete(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)
    @test JuMP.is_valid(model, x)
    JuMP.delete(model, x)
    @test !JuMP.is_valid(model, x)
    second_model = ModelType()
    @test_throws Exception JuMP.delete(second_model, x)
end

function test_variable_bounds_set_get(ModelType, ::Any)
    model = ModelType()
    @variable(model, 0 <= x <= 2)
    @test 0 == @inferred JuMP.lower_bound(x)
    @test 2 == @inferred JuMP.upper_bound(x)
    set_lower_bound(x, 1)
    @test 1 == @inferred JuMP.lower_bound(x)
    set_upper_bound(x, 3)
    @test 3 == @inferred JuMP.upper_bound(x)
    @variable(model, q, Bin)
    @test !JuMP.has_lower_bound(q)
    @test !JuMP.has_upper_bound(q)

    @variable(model, 0 <= y <= 1, Bin)
    @test 0 == @inferred JuMP.lower_bound(y)
    @test 1 == @inferred JuMP.upper_bound(y)

    @variable(model, fixedvar == 2)
    @test 2.0 == @inferred JuMP.fix_value(fixedvar)
    JuMP.fix(fixedvar, 5)
    @test 5 == @inferred JuMP.fix_value(fixedvar)
    @test_throws Exception JuMP.lower_bound(fixedvar)
    @test_throws Exception JuMP.upper_bound(fixedvar)
end

function test_variable_starts_set_get(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:3])
    x0 = collect(1:3)
    JuMP.set_start_value.(x, x0)
    @test JuMP.start_value.(x) == x0
    @test JuMP.start_value.([x[1], x[2], x[3]]) == x0

    @variable(model, y[1:3, 1:2])
    @test_throws DimensionMismatch JuMP.set_start_value.(y, collect(1:6))
end

function test_variable_integrality_set_get(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:3])

    JuMP.set_integer(x[2])
    JuMP.set_integer(x[2])  # test duplicated call
    @test JuMP.is_integer(x[2])
    JuMP.unset_integer(x[2])
    @test !JuMP.is_integer(x[2])

    JuMP.set_binary(x[1])
    JuMP.set_binary(x[1])  # test duplicated call
    @test JuMP.is_binary(x[1])
    @test_throws Exception JuMP.set_integer(x[1])
    JuMP.unset_binary(x[1])
    @test !JuMP.is_binary(x[1])

    @variable(model, y, binary = true)
    @test JuMP.is_binary(y)
    @test_throws Exception JuMP.set_integer(y)
    JuMP.unset_binary(y)
    @test !JuMP.is_binary(y)

    @variable(model, z, integer = true)
    @test JuMP.is_integer(z)
    @test_throws Exception JuMP.set_binary(z)
    JuMP.unset_integer(z)
    @test !JuMP.is_integer(z)
end

function test_variable_repeated_elements(ModelType, ::Any)
    # Tests repeated elements in index set throw error (JuMP issue #199).
    model = ModelType()
    index_set = [:x, :x, :y]
    @test_throws ErrorException (@variable(
        model,
        unused_variable[index_set],
        container = DenseAxisArray
    ))
    @test_throws ErrorException (@variable(
        model,
        unused_variable[index_set],
        container = SparseAxisArray
    ))
    @test_throws ErrorException (@variable(
        model,
        unused_variable[index_set, [1]],
        container = DenseAxisArray
    ))
    @test_throws ErrorException (@variable(
        model,
        unused_variable[index_set, [1]],
        container = SparseAxisArray
    ))
end

function test_variable_oneto_index_set(ModelType, VariableRefType)
    # Tests that Base.OneTo can be used in index set (JuMP issue #933).
    model = ModelType()
    auto_var = @variable(model, [Base.OneTo(3), 1:2], container = Auto)
    @test auto_var isa Matrix{VariableRefType}
    @test (3, 2) == @inferred size(auto_var)
    array_var = @variable(model, [Base.OneTo(3), 1:2], container = Array)
    @test array_var isa Matrix{VariableRefType}
    @test (3, 2) == @inferred size(array_var)
    denseaxisarray_var =
        @variable(model, [Base.OneTo(3), 1:2], container = DenseAxisArray)
    @test denseaxisarray_var isa JuMP.Containers.DenseAxisArray{VariableRefType}
    @test length.(axes(denseaxisarray_var)) == (3, 2)
end

function test_variable_base_name_in_macro(ModelType, ::Any)
    model = ModelType()
    @variable(model, normal_var)
    _test_variable_name_util(normal_var, "normal_var")
    no_indices = @variable(model, base_name = "foo")
    _test_variable_name_util(no_indices, "foo")
    # Note that `z` will be ignored in name.
    indices = @variable(model, z[i = 2:3], base_name = "t")
    _test_variable_name_util(indices[2], "t[2]")
    return _test_variable_name_util(indices[3], "t[3]")
end

function test_variable_name(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)
    _test_variable_name_util(x, "x")
    JuMP.set_name(x, "y")
    @test JuMP.variable_by_name(model, "x") isa Nothing
    _test_variable_name_util(x, "y")
    y = @variable(model, base_name = "y")
    err(name) = ErrorException("Multiple variables have the name $name.")
    @test_throws err("y") JuMP.variable_by_name(model, "y")
    JuMP.set_name(y, "x")
    _test_variable_name_util(x, "y")
    _test_variable_name_util(y, "x")
    JuMP.set_name(x, "x")
    @test_throws err("x") JuMP.variable_by_name(model, "x")
    @test JuMP.variable_by_name(model, "y") isa Nothing
    JuMP.set_name(y, "y")
    _test_variable_name_util(x, "x")
    return _test_variable_name_util(y, "y")
end

function test_variable_condition_in_indexing(ModelType, ::Any)
    function test_one_dim(x)
        @test 5 == @inferred length(x)
        for i in 1:10
            if iseven(i)
                @test haskey(x, i)
            else
                @test !haskey(x, i)
            end
        end
    end

    function test_two_dim(y)
        @test 15 == @inferred length(y)
        for j in 1:10, k in 3:2:9
            if isodd(j + k) && k <= 8
                @test haskey(y, (j, k))
            else
                @test !haskey(y, (j, k))
            end
        end
    end

    model = ModelType()
    # Parses as ref on 0.7.
    @variable(model, named_one_dim[i = 1:10; iseven(i)])
    test_one_dim(named_one_dim)
    # Parses as vcat on 0.7.
    anon_one_dim = @variable(model, [i = 1:10; iseven(i)])
    test_one_dim(anon_one_dim)
    # Parses as typed_vcat on 0.7.
    @variable(model, named_two_dim[j = 1:10, k = 3:2:9; isodd(j + k) && k <= 8])
    test_two_dim(named_two_dim)
    # Parses as vect on 0.7.
    anon_two_dim =
        @variable(model, [j = 1:10, k = 3:2:9; isodd(j + k) && k <= 8])
    return test_two_dim(anon_two_dim)
end

function test_variable_macro_return_type(ModelType, VariableRefType)
    model = ModelType()
    @variable(model, x[1:3, 1:4, 1:2], start = 0.0)
    @test typeof(x) == Array{VariableRefType,3}
    @test typeof(JuMP.start_value.(x)) == Array{Float64,3}

    @variable(model, y[1:0], start = 0.0)
    @test typeof(y) == Vector{VariableRefType}
    # No type to infer for an empty collection.
    @test typeof(JuMP.start_value.(y)) == Vector{Union{Nothing,Float64}}

    @variable(model, z[1:4], start = 0.0)
    @test typeof(z) == Vector{VariableRefType}
    @test typeof(JuMP.start_value.(z)) == Vector{Float64}
end

function test_variable_start_value_on_empty(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:4, 1:0, 1:3], start = 0)  # Array{VariableRef}
    @variable(model, y[1:4, 2:1, 1:3], start = 0)  # DenseAxisArray
    @variable(model, z[1:4, Set(), 1:3], start = 0)  # SparseAxisArray

    @test JuMP.start_value.(x) == Array{Float64}(undef, 4, 0, 3)
    # TODO: Decide what to do here. I don't know if we still need to test this
    #       given broadcast syntax.
    # @test typeof(JuMP.start_value(y)) <: JuMP.DenseAxisArray{Float64}
    # @test JuMP.size(JuMP.start_value(y)) == (4,0,3)
    # @test typeof(JuMP.start_value(z)) ==
    #   JuMP.DenseAxisArray{Float64,3,Tuple{UnitRange{Int},Set{Any},UnitRange{Int}}}
    # @test length(JuMP.start_value(z)) == 0
end

function test_variable_denseaxisarray_slices(ModelType, VariableRefType)
    # Test slicing DenseAxisArrays (JuMP issue #684).
    model = ModelType()
    @variable(model, x[1:3, 1:4, 1:2], container = DenseAxisArray)
    @variable(model, y[1:3, -1:2, 3:4])
    @variable(model, z[1:3, -1:2:4, 3:4])
    @variable(model, w[1:3, -1:2, [:red, "blue"]])

    #@test x[:] == vec(_sliceof_util(VariableRefType, x, 1:3, 1:4, 1:2))
    @test x isa JuMP.Containers.DenseAxisArray
    @test x[:, :, :].data == _sliceof_util(VariableRefType, x, 1:3, 1:4, 1:2)
    @test x[1, :, :].data == _sliceof_util(VariableRefType, x, 1, 1:4, 1:2)
    @test x[1, :, 2].data == _sliceof_util(VariableRefType, x, 1, 1:4, 2)
    @test_throws KeyError x[1, :, 3]
    #@test x[1:2,:,:].data == _sliceof_util(VariableRefType, x, 1:2, 1:4, 1:2)
    #@test x[1:2,:,2].data == _sliceof_util(VariableRefType, x, 1:2, 1:4, 2)
    #@test x[1:2,:,1:2].data == _sliceof_util(VariableRefType, x, 1:2, 1:4, 1:2)
    @test_throws KeyError x[1:2, :, 1:3]

    #@test y[:] == vec(_sliceof_util(VariableRefType, y, 1:3, -1:2, 3:4))
    @test y[:, :, :].data == _sliceof_util(VariableRefType, y, 1:3, -1:2, 3:4)
    @test y[1, :, :].data == _sliceof_util(VariableRefType, y, 1, -1:2, 3:4)
    @test y[1, :, 4].data == _sliceof_util(VariableRefType, y, 1, -1:2, 4)
    @test_throws KeyError y[1, :, 5]
    # @test y[1:2,:,:] == _sliceof_util(VariableRefType, y, 1:2, -1:2, 3:4)
    # @test y[1:2,:,4] == _sliceof_util(VariableRefType, y, 1:2, -1:2, 4)
    # @test y[1:2,:,3:4] == _sliceof_util(VariableRefType, y, 1:2, -1:2, 3:4)
    # @test_throws BoundsError y[1:2,:,1:3]

    #@test z[:] == vec(_sliceof_util(VariableRefType, z, 1:3, -1:2:4, 3:4))
    @test z[:, 1, :].data == _sliceof_util(VariableRefType, z, 1:3, 1, 3:4)
    @test z[1, 1, :].data == _sliceof_util(VariableRefType, z, 1, 1, 3:4)
    @test_throws KeyError z[:, 5, 3]
    # @test z[1:2,1,:] == _sliceof_util(VariableRefType, z, 1:2, 1, 3:4)
    # @test z[1:2,1,4] == _sliceof_util(VariableRefType, z, 1:2, 1, 4)
    # @test z[1:2,1,3:4] == _sliceof_util(VariableRefType, z, 1:2, 1, 3:4)
    # @test_throws BoundsError z[1:2,1,1:3]

    #@test w[:] == vec(_sliceof_util(VariableRefType, w, 1:3, -1:2, [:red,"blue"]))
    @test w[:, :, :] == w
    @test w[1, :, "blue"].data ==
          _sliceof_util(VariableRefType, w, 1, -1:2, ["blue"])
    @test w[1, :, :red].data ==
          _sliceof_util(VariableRefType, w, 1, -1:2, [:red])
    @test_throws KeyError w[1, :, "green"]
    # @test w[1:2,:,"blue"] == _sliceof_util(VariableRefType, w, 1:2, -1:2, ["blue"])
    # @test_throws ErrorException w[1:2,:,[:red,"blue"]]
end

function test_variable_end_indexing(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[0:2, 1:4])
    @variable(model, z[0:2])
    @test x[end, 1] == x[2, 1]
    @test x[0, end-1] == x[0, 3]
    @test z[end] == z[2]
    # TODO: It is redirected to x[11] as it is the 11th element but linear
    #       indexing is not supported
    @test_throws BoundsError x[end-1]
end

function test_variable_unsigned_index(ModelType, ::Any)
    # Tests unsigned int can be used to construct index set (JuMP issue #857).
    model = ModelType()
    t = UInt(4)
    @variable(model, x[1:t])
    @test 4 == @inferred num_variables(model)
end

function test_variable_symmetric(ModelType, ::Any)
    model = ModelType()

    @variable model x[1:2, 1:2] Symmetric
    @test x isa Symmetric
    @test x[1, 2] === x[2, 1]
    @test model[:x] === x

    y = @variable model [1:2, 1:2] Symmetric
    @test y isa Symmetric
    @test y[1, 2] === y[2, 1]
end

function test_variable_skewsymmetric(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:2, 1:2] in SkewSymmetricMatrixSpace())
    @test x[1, 2] == -x[2, 1]
    @test iszero(x[1, 1])
    @test iszero(x[2, 2])
    @test model[:x] === x
    @test num_variables(model) == 1
    @variable(model, z[1:3, 1:3] in SkewSymmetricMatrixSpace())
    @test z[1, 2] == -z[2, 1]
    @test z[1, 3] == -z[3, 1]
    @test z[2, 3] == -z[3, 2]
    @test iszero(z[1, 1])
    @test iszero(z[2, 2])
    @test iszero(z[3, 3])
    @test model[:z] === z
    @test num_variables(model) == 4
    y = @variable(model, [1:3, 1:3] in SkewSymmetricMatrixSpace())
    @test y[1, 2] == -y[2, 1]
    @test y[2, 3] == -y[3, 2]
    @test iszero(y[3, 3])
    @test num_variables(model) == 7
end

function test_variables_constrained_on_creation(ModelType, ::Any)
    model = ModelType()

    err = ErrorException(
        "In `@variable(model, x[1:2] in SecondOrderCone(), set = PSDCone())`: Cannot specify set twice, it was already set to `\$(Expr(:escape, :(SecondOrderCone())))` so the `set` keyword argument is not allowed.",
    )
    @test_macro_throws err @variable(
        model,
        x[1:2] in SecondOrderCone(),
        set = PSDCone()
    )
    err = ErrorException(
        "In `@variable(model, x[1:2] in SecondOrderCone(), PSD)`: Cannot specify set twice, it was already set to `\$(Expr(:escape, :(SecondOrderCone())))` so the `PSD` argument is not allowed.",
    )
    @test_macro_throws err @variable(model, x[1:2] in SecondOrderCone(), PSD)
    err = ErrorException(
        "In `@variable(model, x[1:2] in SecondOrderCone(), Symmetric)`: Cannot specify `Symmetric` when the set is already specified, the variable is constrained to belong to `\$(Expr(:escape, :(SecondOrderCone())))`.",
    )
    @test_macro_throws err @variable(
        model,
        x[1:2] in SecondOrderCone(),
        Symmetric
    )
    err = ErrorException(
        "In `@variable(model, x[1:2], set = SecondOrderCone(), set = PSDCone())`: `set` keyword argument was given 2 times.",
    )
    @test_macro_throws err @variable(
        model,
        x[1:2],
        set = SecondOrderCone(),
        set = PSDCone()
    )

    @variable(model, x[1:2] in SecondOrderCone())
    @test num_constraints(model, typeof(x), MOI.SecondOrderCone) == 1
    @test name(x[1]) == "x[1]"
    @test name(x[2]) == "x[2]"

    @variable(model, [1:2] in SecondOrderCone())
    @test num_constraints(model, typeof(x), MOI.SecondOrderCone) == 2

    @variable(model, [1:3] in MOI.SecondOrderCone(3))
    @test num_constraints(model, typeof(x), MOI.SecondOrderCone) == 3

    z = @variable(model, z in MOI.Semiinteger(1.0, 2.0))
    @test num_constraints(model, typeof(z), MOI.Semiinteger{Float64}) == 1

    @variable(model, set = MOI.Semiinteger(1.0, 2.0))
    @test num_constraints(model, typeof(z), MOI.Semiinteger{Float64}) == 2

    @variable(model, [1:3, 1:3] in PSDCone())
    @test num_constraints(
        model,
        typeof(x),
        MOI.PositiveSemidefiniteConeTriangle,
    ) == 1
end

function test_batch_delete_variables(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:3] >= 1)
    @objective(model, Min, sum([1, 2, 3] .* x))
    @test all(is_valid.(model, x))
    delete(model, x[[1, 3]])
    @test all((!is_valid).(model, x[[1, 3]]))
    @test is_valid(model, x[2])
    second_model = ModelType()
    @test_throws Exception JuMP.delete(second_model, x[2])
    @test_throws Exception JuMP.delete(second_model, x[[1, 3]])
end

function test_Model_all_variable(::Any, ::Any)
    model = Model()
    @variable(model, x)
    @variable(model, y)
    @test [x, y] == @inferred JuMP.all_variables(model)
end

function test_Model_macro_variables(::Any, ::Any)
    model = Model()
    @variables model begin
        0 ≤ x[i = 1:2] ≤ i
        y ≥ 2, Int, (start = 0.7)
        z ≤ 3, (start = 10)
        q, (Bin, start = 0.5)
    end

    @test "x[1]" == @inferred JuMP.name(x[1])
    @test 0 == @inferred JuMP.lower_bound(x[1])
    @test 1 == @inferred JuMP.upper_bound(x[1])
    @test !JuMP.is_binary(x[1])
    @test !JuMP.is_integer(x[1])
    @test JuMP.start_value(x[1]) === nothing

    @test "x[2]" == @inferred JuMP.name(x[2])
    @test 0 == @inferred JuMP.lower_bound(x[2])
    @test 2 == @inferred JuMP.upper_bound(x[2])
    @test !JuMP.is_binary(x[2])
    @test !JuMP.is_integer(x[2])
    @test JuMP.start_value(x[2]) === nothing

    @test "y" == @inferred JuMP.name(y)
    @test 2 == @inferred JuMP.lower_bound(y)
    @test !JuMP.has_upper_bound(y)
    @test !JuMP.is_binary(y)
    @test JuMP.is_integer(y)
    @test JuMP.start_value(y) === 0.7

    @test "z" == @inferred JuMP.name(z)
    @test !JuMP.has_lower_bound(z)
    @test 3 == @inferred JuMP.upper_bound(z)
    @test !JuMP.is_binary(z)
    @test !JuMP.is_integer(z)
    @test JuMP.start_value(z) === 10.0

    @test "q" == @inferred JuMP.name(q)
    @test !JuMP.has_lower_bound(q)
    @test !JuMP.has_upper_bound(q)
    @test JuMP.is_binary(q)
    @test !JuMP.is_integer(q)
    @test JuMP.start_value(q) === 0.5
end

# @testset "Variables for JuMPExtension.MyModel" begin
#     variables_test(JuMPExtension.MyModel, JuMPExtension.MyVariableRef)
# end

function test_Model_dual_variable(::Any, ::Any)
    model = Model()
    @variable(model, x == 0)
    exception = ErrorException(
        "To query the dual variables associated with a variable bound, first " *
        "obtain a constraint reference using one of `UpperBoundRef`, `LowerBoundRef`, " *
        "or `FixRef`, and then call `dual` on the returned constraint reference.\nFor " *
        "example, if `x <= 1`, instead of `dual(x)`, call `dual(UpperBoundRef(x))`.",
    )
    @test_throws exception JuMP.dual(x)
end

function test_Model_value_containers(::Any, ::Any)
    model = Model()
    @variable(model, x[1:2])
    exception = ErrorException(
        "`JuMP.value` is not defined for collections of JuMP types. Use " *
        "Julia's broadcast syntax instead: `JuMP.value.(x)`.",
    )
    @test_throws exception JuMP.value(x)
end

function test_Model_get_variable_coefficient(::Any, ::Any)
    m = Model()
    x = @variable(m, x)
    y = @variable(m, y)
    @test coefficient(x, x) == 1.0
    @test coefficient(x, y) == 0.0
    @test coefficient(x, x, x) == 0.0
    @test coefficient(x, y, x) == coefficient(x, x, y) == 0.0
    @test coefficient(x, y, y) == 0.0
end

function _mock_reduced_cost_util(
    obj_sense::MOI.OptimizationSense,
    #obj_value::Float64,
    var_obj_coeff::Float64,
    #var_value,
    var_bound_type::Symbol,
    var_bounds_dual = nothing,
    has_duals::Bool = var_bounds_dual !== nothing,
)
    mockoptimizer =
        MOIU.MockOptimizer(MOIU.Model{Float64}(), eval_objective_value = false)
    m = JuMP.direct_model(mockoptimizer)
    if var_bound_type === :lower
        @variable(m, x >= 0)
        if has_duals
            @assert isa(var_bounds_dual, Float64)
            has_duals && MOI.set(
                mockoptimizer,
                MOI.ConstraintDual(),
                JuMP.optimizer_index(JuMP.LowerBoundRef(x)),
                var_bounds_dual,
            )
        end
    elseif var_bound_type === :upper
        @variable(m, x <= 10)
        if has_duals
            @assert isa(var_bounds_dual, Float64)
            has_duals && MOI.set(
                mockoptimizer,
                MOI.ConstraintDual(),
                JuMP.optimizer_index(JuMP.UpperBoundRef(x)),
                var_bounds_dual,
            )
        end
    elseif var_bound_type === :fixed
        @variable(m, x == 10)
        if has_duals
            @assert isa(var_bounds_dual, Float64)
            MOI.set(
                mockoptimizer,
                MOI.ConstraintDual(),
                JuMP.optimizer_index(JuMP.FixRef(x)),
                var_bounds_dual,
            )
        end
    elseif var_bound_type === :both
        @variable(m, 0 <= x <= 10)
        if has_duals
            @assert length(var_bounds_dual) == 2
            @assert eltype(var_bounds_dual) == Float64
            lb_dual, ub_dual = var_bounds_dual
            MOI.set(
                mockoptimizer,
                MOI.ConstraintDual(),
                JuMP.optimizer_index(JuMP.LowerBoundRef(x)),
                lb_dual,
            )
            MOI.set(
                mockoptimizer,
                MOI.ConstraintDual(),
                JuMP.optimizer_index(JuMP.UpperBoundRef(x)),
                ub_dual,
            )
        end
    elseif var_bound_type === :none
        @variable(m, x)
        @assert var_bounds_dual === nothing
    else
        error("unrecognized bound type")
    end

    @objective(m, obj_sense, var_obj_coeff * x)

    if has_duals
        MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
        MOI.set(mockoptimizer, MOI.ResultCount(), 1)
        MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
        MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    end

    return x
end

function test_Model_reduced_cost(::Any, ::Any)
    Min = MOI.MIN_SENSE
    Max = MOI.MAX_SENSE
    # The method should always fail if duals are not available.
    x = _mock_reduced_cost_util(Min, 1.0, :none)
    @test_throws ErrorException reduced_cost(x)
    x = _mock_reduced_cost_util(Min, 1.0, :fixed)
    @test_throws ErrorException reduced_cost(x)
    x = _mock_reduced_cost_util(Min, 1.0, :lower)
    @test_throws ErrorException reduced_cost(x)
    x = _mock_reduced_cost_util(Min, 1.0, :upper)
    @test_throws ErrorException reduced_cost(x)
    x = _mock_reduced_cost_util(Min, 1.0, :both)
    @test_throws ErrorException reduced_cost(x)
    # My reimplementation of the tests suggested by @odow.
    # Note that the floating point values are compared by equality because
    # there is no risk of the solver messing this up (mocks are being used).
    # First the fixed variable tests.
    x = _mock_reduced_cost_util(Min, 1.0, :none, nothing, true) # free var
    @test reduced_cost(x) == 0.0
    x = _mock_reduced_cost_util(Min, 1.0, :fixed, 1.0) # min x, x == 10
    @test reduced_cost(x) == 1.0
    x = _mock_reduced_cost_util(Max, 1.0, :fixed, -1.0) # max x, x == 10
    @test reduced_cost(x) == 1.0
    x = _mock_reduced_cost_util(Min, -1.0, :fixed, -1.0) # min -x, x == 10
    @test reduced_cost(x) == -1.0
    x = _mock_reduced_cost_util(Max, -1.0, :fixed, 1.0) # max -x, x == 10
    @test reduced_cost(x) == -1.0
    # Then the double bounded variables.
    #x = _mock_reduced_cost_util(Min, 1.0, :both, (0.0, 1.0)) # min x, 0 <= x <= 10
    #@test reduced_cost(x) == 1.0
    #x = _mock_reduced_cost_util(Max, 1.0, :both, (-1.0, 0.0)) # max x, 0 <= x <= 10
    #@test reduced_cost(x) == 1.0
    #x = _mock_reduced_cost_util(Min, -1.0, :both, (-1.0, 0.0)) # min -x, 0 <= x <= 10
    #@test reduced_cost(x) == -1.0
    #x = _mock_reduced_cost_util(Max, -1.0, :both, (0.0, 1.0)) # max -x, 0 <= x <= 10
    #@test reduced_cost(x) == -1.0
    # Test for a single upper bound and a single lower bound.
    x = _mock_reduced_cost_util(Min, 1.0, :lower, 1.0) # min x, 0 <= x
    @test reduced_cost(x) == 1.0
    x = _mock_reduced_cost_util(Max, 1.0, :upper, 1.0) # max x, x <= 10
    @test reduced_cost(x) == 1.0
end

function test_Model_value(::Any, ::Any)
    @test value(1) === 1
    @test value(1.0) === 1.0
    @test value(JuMP._MA.Zero()) === 0.0
end

function test_Model_value_var(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:2])
    vals = Dict(x[1] => 1.0, x[2] => 2.0)
    f = vidx -> vals[vidx]
    @test value(x[1], f) === 1.0
    @test value(x[2], f) === 2.0
end

function test_Model_relax_integrality(::Any, ::Any)
    model = Model()
    @variable(model, x, Bin)
    @variable(model, -1 <= y <= 2, Bin)
    @variable(model, 0.1 <= z <= 0.6, Bin)
    @variable(model, a, Int)
    @variable(model, -1 <= b <= 2, Int)
    unrelax = relax_integrality(model)

    @test !is_binary(x)
    @test lower_bound(x) == 0.0
    @test upper_bound(x) == 1.0

    @test !is_binary(y)
    @test lower_bound(y) == 0.0
    @test upper_bound(y) == 1.0

    @test !is_binary(z)
    @test lower_bound(z) == 0.1
    @test upper_bound(z) == 0.6

    @test !is_integer(a)
    @test !has_lower_bound(a)
    @test !has_upper_bound(a)

    @test !is_integer(b)
    @test lower_bound(b) == -1.0
    @test upper_bound(b) == 2.0

    unrelax()

    @test is_binary(x)
    @test !has_lower_bound(x)
    @test !has_upper_bound(x)

    @test is_binary(y)
    @test lower_bound(y) == -1.0
    @test upper_bound(y) == 2.0

    @test is_binary(z)
    @test lower_bound(z) == 0.1
    @test upper_bound(z) == 0.6

    @test is_integer(a)
    @test !has_lower_bound(a)
    @test !has_upper_bound(a)

    @test is_integer(b)
    @test lower_bound(b) == -1.0
    @test upper_bound(b) == 2.0

    fix(x, 1)
    unrelax = relax_integrality(model)
    @test !is_binary(x)
    @test is_fixed(x)
    unrelax()
    @test is_binary(x)
    @test is_fixed(x)

    fix(x, 0)
    unrelax = relax_integrality(model)
    @test !is_binary(x)
    @test is_fixed(x)
    unrelax()
    @test is_binary(x)
    @test is_fixed(x)

    fix(a, 1)
    unrelax = relax_integrality(model)
    @test !is_integer(a)
    @test is_fixed(a)
    unrelax()
    @test is_integer(a)
    @test is_fixed(a)
end

function test_Model_relax_integrality_error_cases(::Any, ::Any)
    model = Model()
    @variable(model, x)
    @constraint(model, x in MOI.Semicontinuous(1.0, 2.0))
    err = ErrorException(
        "Support for relaxing semicontinuous constraints " *
        "is not yet implemented.",
    )
    @test_throws err relax_integrality(model)

    model = Model()
    @variable(model, x)
    @constraint(model, x in MOI.Semiinteger(1.0, 2.0))
    err = ErrorException(
        "Support for relaxing semi-integer constraints " *
        "is not yet implemented.",
    )
    @test_throws err relax_integrality(model)

    model = Model()
    @variable(model, x, Bin)
    fix(x, 2)
    err = ErrorException(
        "The model has no valid relaxation: binary variable " *
        "fixed out of bounds.",
    )
    @test_throws err relax_integrality(model)

    model = Model()
    @variable(model, x, Bin)
    fix(x, -1)
    err = ErrorException(
        "The model has no valid relaxation: binary variable " *
        "fixed out of bounds.",
    )
    @test_throws err relax_integrality(model)
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        f = getfield(@__MODULE__, name)
        @testset "$(name)" begin
            f(Model, VariableRef)
        end
        if !startswith("$(name)", "test_Model_")
            @testset "$(name)-JuMPExtension" begin
                f(JuMPExtension.MyModel, JuMPExtension.MyVariableRef)
            end
        end
    end
end

end  # module TestVariable

TestVariable.runtests()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# SteinKernel

abstract type SteinKernel end

## METHODS TO IMPLEMENT ##
function k(kernel::SteinKernel, x::Array{Float64, 1}, y::Array{Float64, 1})
    error("Must implement the method k")
end

function gradxk(kernel::SteinKernel, x::Array{Float64, 1}, y::Array{Float64, 1})
    error("Must implement the method gradxk")
end

function gradxyk(kernel::SteinKernel, x::Array{Float64, 1}, y::Array{Float64, 1})
    error("Must implement the method gradxyk")
end

# Computes \grad_y k(x, y) = \grad_x k(y, x)
function gradyk(kernel::SteinKernel, x::Array{Float64, 1}, y::Array{Float64, 1})
    gradxk(kernel, y, x)
end

# This is the Chris Oates Langevin Stein kernel
function k0(kernel::SteinKernel,
            x::Array{Float64, 1},
            y::Array{Float64, 1},
            gradlogpx::Array{Float64, 1},
            gradlogpy::Array{Float64, 1})
    k(kernel, x, y) * dot(gradlogpx, gradlogpy) +
        dot(gradlogpx, gradyk(kernel, x, y)) +
        dot(gradlogpy, gradxk(kernel, x, y)) +
        gradxyk(kernel, x, y)
end


text:
export add_command!

"""
    add_command!(
        c::Client,
        prefix::AbstractString,
        func::Function;
        tag::Symbol=gensym(),
        expiry::Union{Int, Period, Nothing}=nothing,
    )
    add_command!(
        func::Function;
        c::Client,
        prefix::AbstractString,
        tag::Symbol=gensym(),
        expiry::Union{Int, Period, Nothing}=nothing,
    )

Add a text command handler. The handler function should take two arguments: A
[`Client`](@ref) and a [`Message`](@ref). The keyword arguments are identical to
[`add_handler!`](@ref). `do` syntax is also accepted.
"""
function add_command!(
    c::Client,
    prefix::AbstractString,
    func::Function;
    tag::Symbol=gensym(),
    expiry::Union{Int, Period, Nothing}=nothing,
)
    if !hasmethod(func, (Client, Message))
        throw(ArgumentError("Handler function must accept (::Client, ::Message)"))
    end

    function handler(c::Client, e::MessageCreate)
        e.message.author.id == me(c).id && return
        startswith(e.message.content, prefix) || return
        func(c, e.message)
    end

    add_handler!(c, MessageCreate, handler; tag=tag, expiry=expiry)
end

function add_command!(
    func::Function,
    c::Client,
    prefix::AbstractString,
    tag::Symbol=gensym(),
    expiry::Union{Int, Period, Nothing}=nothing,
)
    return add_command!(c, prefix, func; tag=tag, expiry=expiry)
end

# TODO: A much nicer command framework.