text:
"""
Solves the equations of fluid flow using a "Stable Fluids"-like algorithm
based on the FFT in a Kolmogorov Flow scenario. This setup results in forced
isotropic turbulence.

Momentum:           ∂u/∂t + (u ⋅ ∇) u = − ∇p + ν ∇²u + f

Incompressibility:  ∇ ⋅ u = 0

u:  Velocity (2d vector)
p:  Pressure
f:  Forcing (here Kolmogorov Forcing)
ν:  Kinematic Viscosity
t:  Time
∇:  Nabla operator (defining nonlinear convection, gradient and divergence)
∇²: Laplace Operator

----

Setup:

A rectangular domain (16x9) with periodic boundary conditions in both axes.

            +-------------------------------------------------+
            |                                                 |
            |             *                      *            |
            |          *           *    *    *                |
            |                                                 |
            |                                 *               |
            |     *       *                                   |
            |                      *     *                    |
            |                                            *    |
            |                                                 |
            |  *          *      *                 *       *  |
            |                            *                    |
            +-------------------------------------------------+

-> The fluid is initialized at rest

-> Then, a forcing is applied according to

        ⎡   c * sin(k * y)  ⎤
    f = ⎢                   ⎥
        ⎣         0         ⎦

        with scaling c and wavenumber k

-----

Solution Strategy:

-> Start with zero velocity everywhere: u = [0, 0]

1. Add forces

    w₁ = u + Δt f

2. Convect by self-advection (set the value at the current
   location to be the value at the position backtraced
   on the streamline.) -> unconditionally stable

    w₂ = w₁(p(x, −Δt))

3. Subtract mean velocity for stabilization

4. Diffuse and Project in Fourier Domain

    4.1 Forward Transformation into Fourier Domain

        w₂ → w₃

    4.2 Diffuse by "low-pass filtering" (convolution
        is multiplication in the Fourier Domain)

        w₄ = exp(− k² ν Δt) w₃

    4.3 Compute the (pseudo-) pressure in the Fourier Domain
        by evaluating the divergence in the Fourier Domain

        q = w₄ ⋅ k / ||k||₂

    4.4 Correct the velocities such that they are incompressible

        w₅ = w₄ − q k / ||k||₂

    4.5 Inverse Transformation back into spatial domain

        w₆ ← w₅

5. Subtract mean velocity for stabilization

6. Repeat

k = [ k_x, k_y ] are the spatial frequencies (= wavenumbers)

The Fourier Transformation implicitly prescribes the periodic
Boundary Conditions

-----

To create a video/gif out of the saved images run

    ffmpeg -f image2 -framerate 24 -i kolmogorov_%05d.png kolmogorov_animation.gif

"""

using FFTW
using Plots
using ProgressMeter
using Interpolations
using LinearAlgebra
using Printf
using Statistics

N_POINTS_Y = 360
ASPECT_RATIO = 16/9
KINEMATIC_VISCOSITY = 1.0 / 1000.0
TIME_STEP_LENGTH = 0.001
N_TIME_STEPS = 750

FORCING_WAVENUMBER = 8
FORCING_SCALE = 100.0

function main()
    n_points_x = Int(ASPECT_RATIO * N_POINTS_Y)
    x_extent = n_points_x / N_POINTS_Y - 1.0e-7
    x_interval = range(0.0, x_extent, length=n_points_x)
    y_interval = range(0.0, 1.0, length=N_POINTS_Y)

    coordinates_x = [x for x in x_interval, y in y_interval]
    coordinates_y = [y for x in x_interval, y in y_interval]

    wavenumbers_1d_x = rfftfreq(n_points_x) .* n_points_x
    n_fft_points_x = length(wavenumbers_1d_x)
    wavenumbers_1d_y = fftfreq(N_POINTS_Y) .* N_POINTS_Y
    n_fft_points_y = length(wavenumbers_1d_y)

    wavenumbers_x = [k_x for k_x in wavenumbers_1d_x, k_y in wavenumbers_1d_y]
    wavenumbers_y = [k_y for k_x in wavenumbers_1d_x, k_y in wavenumbers_1d_y]
    wavenumbers_norm = [norm([k_x, k_y]) for k_x in wavenumbers_1d_x, k_y in wavenumbers_1d_y]

    decay = exp.(- TIME_STEP_LENGTH .* KINEMATIC_VISCOSITY .* wavenumbers_norm.^2)
    wavenumbers_norm[iszero.(wavenumbers_norm)] .= 1.0
    normalized_wavenumbers_x = wavenumbers_x ./ wavenumbers_norm
    normalized_wavenumbers_y = wavenumbers_y ./ wavenumbers_norm

    force_x = FORCING_SCALE * sin.(FORCING_WAVENUMBER * pi * coordinates_y)

    # Preallocate arrays
    backtraced_coordinates_x = zeros(Float32, n_points_x, N_POINTS_Y)
    backtraced_coordinates_y = zeros(Float32, n_points_x, N_POINTS_Y)

    velocity_x = zeros(Float32, n_points_x, N_POINTS_Y)
    velocity_y = zeros(Float32, n_points_x, N_POINTS_Y)

    velocity_x_prev = zeros(Float32, n_points_x, N_POINTS_Y)
    velocity_y_prev = zeros(Float32, n_points_x, N_POINTS_Y)

    velocity_x_fft = zeros(Complex{Float32}, n_fft_points_x, n_fft_points_y)
    velocity_y_fft = zeros(Complex{Float32}, n_fft_points_x, n_fft_points_y)
    pressure_fft = zeros(Complex{Float32}, n_fft_points_x, n_fft_points_y)

    d_u_d_y_fft = zeros(Complex{Float32}, n_fft_points_x, n_fft_points_y)
    d_v_d_x_fft = zeros(Complex{Float32}, n_fft_points_x, n_fft_points_y)
    curl_fft = zeros(Complex{Float32}, n_fft_points_x, n_fft_points_y)
    curl = zeros(Float32, n_points_x, N_POINTS_Y)

    interpolator = interpolate(
        (x_interval, y_interval),
        velocity_x,
        Gridded(Linear()),
    )

    theme(:dark)

    @showprogress "Timestepping ..." for iter in 1:N_TIME_STEPS
        # (1) Apply the forces
        velocity_x_prev .+= force_x * TIME_STEP_LENGTH

        # (2) Self-Advection by backtracing and interpolation
        backtraced_coordinates_x .= mod1.(
            coordinates_x - TIME_STEP_LENGTH * velocity_x_prev,
            x_extent,
        )
        backtraced_coordinates_y .= mod1.(
            coordinates_y - TIME_STEP_LENGTH * velocity_y_prev,
            1.0,
        )
        interpolator.coefs .= velocity_x_prev
        velocity_x .= interpolator.(backtraced_coordinates_x, backtraced_coordinates_y)
        interpolator.coefs .= velocity_y_prev
        velocity_y .= interpolator.(backtraced_coordinates_x, backtraced_coordinates_y)

        # (3) Stabilize by subtracting the mean velocities
        velocity_x .-= mean(vec(velocity_x))
        velocity_y .-= mean(vec(velocity_y))

        # (4.1) Transform into Fourier domain
        velocity_x_fft = rfft(velocity_x)
        velocity_y_fft = rfft(velocity_y)

        # (4.2) Diffuse by low-pass filtering
        velocity_x_fft .*= decay
        velocity_y_fft .*= decay

        # (4.3) Compute the Pseudo-Pressure by Divergence in Fourier Domain
        pressure_fft = (
            velocity_x_fft .* normalized_wavenumbers_x
            +
            velocity_y_fft .* normalized_wavenumbers_y
        )

        # (4.4) Project the velocities to be incompressible
        velocity_x_fft -= pressure_fft .* normalized_wavenumbers_x
        velocity_y_fft -= pressure_fft .* normalized_wavenumbers_y

        # (4.5) Transform back into spatial domain
        velocity_x = irfft(velocity_x_fft, n_points_x)
        velocity_y = irfft(velocity_y_fft, n_points_x)

        # (5) Stabilize by subtracting the mean velocities
        velocity_x .-= mean(vec(velocity_x))
        velocity_y .-= mean(vec(velocity_y))

        # (6) Advance in time
        velocity_x_prev = velocity_x
        velocity_y_prev = velocity_y

        # Visualize
        d_u_d_y_fft = im .* wavenumbers_y .* velocity_x_fft
        d_v_d_x_fft = im .* wavenumbers_x .* velocity_y_fft
        curl_fft = d_v_d_x_fft - d_u_d_y_fft
        curl = irfft(curl_fft, n_points_x)

        curl = sign.(curl) .* sqrt.(abs.(curl) ./ quantile(vec(curl), 0.8))

        heatmap(
            x_interval,
            y_interval,
            curl',
            c = :seaborn_icefire_gradient,
            size = (1920, 1080),
            clim = (-5.0, 5.0),
            axis = false,
            showaxis = false,
            legend = :none,
            ticks = false,
            margin = 0.0Plots.mm,
            annotations = (
                0.06,
                1.0 - 0.06,
                Plots.text(
                    "$(@sprintf("iter: %05d", iter))",
                    pointsize = 40,
                    color = :white,
                    halign = :left,
                )
            )
        )
        savefig("kolmogorov_$(@sprintf("%05d", iter)).png")

    end

end

main()


text:
# Part 1
layout = Array{Int8}(hcat(map(line -> 'L' .== collect(line), readlines("day11_input.txt"))...)')

function singlestep(current)
    updated = copy(current)
    region = CartesianIndices(current)
    ifirst, ilast = first(region), last(region)
    iunit = oneunit(ifirst)
    for i in region
        neighbours = [current[j] for j in max(ifirst, i - iunit):min(ilast, i + iunit) if j != i]
        if (current[i] == 1 && count(neighbours .== 2) == 0) updated[i] = 2 end
        if (current[i] == 2 && count(neighbours .== 2) >= 4) updated[i] = 1 end
    end
    updated
end

function steadystate(func, current)
    previous = zeros(size(current))
    while current != previous
        previous, current = current, func(current)
    end
    current
end

println("Part 1: occupied count = $(count(steadystate(singlestep, layout) .== 2))")

# Part 2
function singlestep2(current)
    updated = copy(current)
    region = CartesianIndices(current)
    iunit = oneunit(first(region))
    for i in region
        neighbours = []
        for direction in -iunit:iunit
            if !iszero(direction)
                j = i + direction
                while j in region
                    if current[j] != 0
                        push!(neighbours, current[j])
                        break
                    end
                    j += direction
                end
            end
        end
        if (current[i] == 1 && count(neighbours .== 2) == 0) updated[i] = 2 end
        if (current[i] == 2 && count(neighbours .== 2) >= 5) updated[i] = 1 end
    end
    updated
end

println("Part 2: occupied count = $(count(steadystate(singlestep2, layout) .== 2)) ")

ext: jl

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


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

using ArchGDAL
using GDAL

function set_band_names(filename::String, band_names)
    GC.gc()
    local_ods = GDAL.gdalopen(filename, GDAL.GA_Update)
    for _b in 1:(length(band_names))
        GDAL.gdalsetdescription(GDAL.gdalgetrasterband(local_ods, _b), band_names[_b])
    end
    local_ods = nothing
end

function initiate_output_datasets(output_files, x_len, y_len, output_bands, reference_dataset)
    GC.gc()
    for _o in 1:length(output_files)
        outDataset = ArchGDAL.create(output_files[_o], driver=ArchGDAL.getdriver("ENVI"), width=x_len,
                                     height=y_len, nbands=output_bands[_o], dtype=Float64)
        ArchGDAL.setproj!(outDataset, ArchGDAL.getproj(reference_dataset))
        @info "Output Image Size (x,y,b): $x_len, $y_len, $output_bands.  Creating output fractional cover dataset."
        try
            ArchGDAL.setgeotransform!(outDataset, ArchGDAL.getgeotransform(reference_dataset))
        catch e
            println("No geotransorm avaialble, proceeding without")
        end

        for _b in 1:length(output_bands)
            ArchGDAL.setnodatavalue!(ArchGDAL.getband(outDataset,_b), -9999)
        end

        output = zeros(x_len, y_len, output_bands[_o]) .- 9999
        ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])

        outDataset = nothing

    end
end

function write_results(output_files, output_bands, x_len, y_len, results, args)
    GC.gc()

    @info "Write primary output"
    output = zeros(y_len, x_len, output_bands[1]) .- 9999
    for res in results
        if isnothing(res[2]) == false
            output[res[1],res[3], :] = res[2]
        end
    end
    output = permutedims( output, (2,1,3))
    outDataset = ArchGDAL.read(output_files[1], flags=1)
    ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
    outDataset = nothing

    ods_idx = 2
    if args.n_mc > 1
        @info "Write uncertainty output"
        output = zeros(y_len, x_len, output_bands[ods_idx]) .- 9999
        for res in results
            if isnothing(res[4]) == false
                output[res[1],res[3], :] = res[4]
            end
        end

        output = permutedims( output, (2,1,3))
        outDataset = ArchGDAL.read(output_files[ods_idx], flags=1,alloweddrivers =["ENVI"])
        ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
        outDataset = nothing
        ods_idx += 1
    end

    if args.write_complete_fractions == 1
        @info "Write complete fractions output"
        output = zeros(y_len, x_len, output_bands[ods_idx]) .- 9999
        for res in results
            if isnothing(res[5]) == false
                output[res[1],res[3], :] = res[5]
            end
        end

        output = permutedims( output, (2,1,3))
        outDataset = ArchGDAL.read(output_files[ods_idx], flags=1,alloweddrivers =["ENVI"])
        ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
        outDataset = nothing
        ods_idx += 1
    end

end

function load_line(reflectance_file::String, reflectance_uncertainty_file::String, line::Int64,
                       good_bands::Array{Bool}, refl_nodata::Float64)

        reflectance_dataset = ArchGDAL.read(reflectance_file,alloweddrivers =["ENVI"])
        img_dat = convert(Array{Float64},ArchGDAL.readraster(reflectance_file,alloweddrivers =["ENVI"])[:,line,:])
        img_dat = img_dat[:, good_bands]
        good_data = .!all(img_dat .== refl_nodata, dims=2)[:,1]
        img_dat = img_dat[good_data,:]

        if sum(good_data) > 1
            if reflectance_uncertainty_file != ""
                unc_dat = convert(Array{Float64},ArchGDAL.readraster(reflectance_uncertainty_file,alloweddrivers =["ENVI"])[:,line,:])
                unc_dat = unc_dat[:, good_bands]
                unc_dat = unc_dat[good_data,:]
            else
                unc_dat = nothing
            end
        else
            return nothing, nothing, good_data
        end

        return img_dat, unc_dat, good_data
end


text:
#==============================================================================

          Solving the HJB with a state constraint using
                the explicit method

               Translated from matlab code on Ben Moll's website:
                              http://www.princeton.edu/~moll/HACTproject.htm

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

using LinearAlgebra, SparseArrays, Plots

γ = 2 # parameter from CRRA utility
ρ = 0.05 # the discount rate

z1=.1
z2 =.2
z= [z1 z2]

λ1 = 1.5
λ2 = 1.0

λ = [λ1 λ2]

H = 500 # number of points in the grid space
amin= -0.15
amax= 10

a = LinRange(amin,amax,H)
da = (amax-amin)/(H-1)

aa = [a a]
zz = ones(H,1)*z

maxit = 100 # the maximum number of iterations we allow the finite differencing algorithm
crit = 10^(-6) # our critical value
Δ = 40

dVf=zeros(H,2)
dVb = zeros(H,2)
c= zeros(H,2)

# Create the matrix that captures the systems dynamic from λ process
Aswitch = [-sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[1]; sparse(I,H,H)*λ[2] -sparse(I,H,H)*λ[2]]

# Set up parameters for outer loop
Ir = 20 # number of iterations, values to check
rmin= -0.05
rmax=0.04
r_grid =LinRange(rmin,rmax,Ir)

# Inital Guess
r = r_grid[1]

v0=zeros(H,2)
v0[:,1] = (z[1] .+ max(r,0.01).*a).^(1-γ)/(1-γ)/ρ
v0[:,2] = (z[2] .+ max(r,0.01).*a).^(1-γ)/(1-γ)/ρ

global v0=v0

dist=[]
V_n=[]
g_r=[]
dV_r=[]
V_r=[]
c_r=[]
adot=[]

S =zeros(Ir,1)

# Loop that finds the equilibrium interest rate
for ir in 1:Ir
    global r = r_grid[ir]

    if ir>1
        global v0 = V_r[ir-1]
    end

    global v = v0

    # The finite differeing loop
    for n in 1:maxit
        global V=v
        push!(V_n,V)
        # Forward differencing
        dVf[1:H-1,:]=(V[2:H,:]-V[1:H-1,:])/da
        dVf[H,:] = (z .+ r.*amax).^(-γ) # impose state constraint at the max, just in case
        # Backward differencing
        dVb[2:H,:]=(V[2:H,:]-V[1:H-1,:])/da
        dVb[1,:]= (z .+r.*amin).^(-γ)

        # Find consumption and savings with the forward difference
        cf = dVf.^(-1/γ)
        sf = zz +r.*aa -cf

        # Find consumption and savings with the backward difference
        cb= dVb.^(-1/γ)
        sb = zz +r.*aa -cb

        # Find consumption and savings at steady state
        c0 = zz + r.*aa
        dV0 = c0.^(-γ)

        # Now implement the upwind scheme in order to select the best differencing method
        If = sf .> 0
        Ib = sb .< 0
        I0 = (1 .- If - Ib)

        # State constraint at amin is automatically implemented

        global dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0
        global c = dV_Upwind.^(-1/γ)

        u = c.^(1-γ)/(1-γ)
        # Construct matrix for the evolution of the system
        X = -min.(sb,0)/da
        Y= -max.(sf,0)/da + min.(sb,0)/da
        Z = max.(sf,0)/da

        A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
        A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
        AA = [A1 spzeros(H,H); spzeros(H,H) A2]

        global A = AA + Aswitch

        B = (ρ +1/Δ)*sparse(I,2*H,2*H) -A
        u_stacked = [u[:,1]; u[:,2]]
        V_stacked = [V[:,1]; V[:,2]]

        b = u_stacked + V_stacked/Δ
        V_stacked = B\b # Solves the system of equations

        V = [V_stacked[1:H] V_stacked[H+1:2*H]]

        V_change = V - v

        global v = V

        push!(dist, findmax(abs.(V_change))[1])

        if dist[n] < crit
            println("Value Function converged, Iteration=")
            println(n)
            break
        end
    end

    #====================================================
            Fokker-Planck Equation
    ====================================================#

    AT = A' # we need the transpose of our A matrix
    b =zeros(2*H,1)

    # We need to fix one value to prevent the matrix from being singular
    i_fix =1
    b[i_fix] = .1
    row = [zeros(1,i_fix-1) 1 zeros(1,2*H-i_fix)]
    AT[i_fix,:] = row

    # Using this solve the linear system
    gg = AT\b
    g_sum = gg'*ones(2*H,1)*da
    gg = gg./g_sum # we need gg to sum to one

    g = [gg[1:H] gg[H+1:2*H]]

    # Graphs
    push!(g_r, g)
    push!(adot, zz+r.*aa-c)
    push!(V_r, V)
    push!(c_r, c)
    push!(dV_r, dV_Upwind)

    S[ir] = g[:,1]'*a*da + g[:,2]'*a*da

#end of outer loop
end

# Graphs
ir =1
plot(a,adot[ir][:,1], label="\$ s_1(a,r_1)\$", xlabel="a", ylabel="\$s_i(a)\$", xlims=(amin,amax))
plot!(a,adot[ir][:,2], label="\$ s_2(a,r_1)\$")
plot!(a,adot[ir+1][:,1], label="\$ s_1(a,r_2)\$")
plot!(a,adot[ir+2][:,2], label="\$ s_2(a,r_2)\$")

ir =9
plot(a,c_r[ir][:,1], label="\$ c_1(a,r_1)\$", xlabel="a", ylabel="\$c_i(a)\$", xlims=(amin,amax))
plot!(a,c_r[ir][:,2], label="\$ c_2(a,r_1)\$")
plot!(a,c_r[ir+1][:,1], label="\$ c_1(a,r_2)\$")
plot!(a,c_r[ir+2][:,2], label="\$ c_2(a,r_2)\$")

amax1=.8
plot(a,g_r[ir][:,1], label="\$ g_1(a)\$", xlabel="a", ylabel="\$g_i(a)\$", xlims=(amin,amax1))
plot!(a,adot[ir][:,2], label="\$ g_2(a)\$")

Smax=findmax(S)[1]
amin1=1.1*amin
aaa =LinRange(amin1,Smax,Ir)
rrr = LinRange(rmin, 0.06, Ir)

plot(S,r_grid, ylims=(rmin,0.06), xlims=(amin1,Smax), ylabel="\$r\$", xlabel="S(r)", label="\$S(r)\$")
plot!(zeros(Ir,1), rrr, legend=:bottomright, color=:red, label="\$a=0\$")
plot!(zeros(Ir,1).+amin, rrr, label="\$a=\\underline{a}\$", color=:yellow, line=:dash)
plot!(aaa, ones(Ir,1).*ρ, label="\$r=\\rho \$", line=:dash)


text:
@testset "Graph coloring" begin

@testset "49 US states + DC" begin
    com = CS.init()

    num_colors = 8
    washington = add_var!(com, 1, num_colors)
    montana = add_var!(com, 1, num_colors)
    maine = add_var!(com, 1, num_colors)
    north_dakota = add_var!(com, 1, num_colors)
    south_dakota = add_var!(com, 1, num_colors)
    wyoming = add_var!(com, 1, num_colors)
    wisconsin = add_var!(com, 1, num_colors)
    idaho = add_var!(com, 1, num_colors)
    vermont = add_var!(com, 1, num_colors)
    minnesota = add_var!(com, 1, num_colors)
    oregon = add_var!(com, 1, num_colors)
    new_hampshire = add_var!(com, 1, num_colors)
    iowa = add_var!(com, 1, num_colors)
    massachusetts = add_var!(com, 1, num_colors)
    nebraska = add_var!(com, 1, num_colors)
    new_york = add_var!(com, 1, num_colors)
    pennsylvania = add_var!(com, 1, num_colors)
    connecticut = add_var!(com, 1, num_colors)
    rhode_island = add_var!(com, 1, num_colors)
    new_jersey = add_var!(com, 1, num_colors)
    indiana = add_var!(com, 1, num_colors)
    nevada = add_var!(com, 1, num_colors)
    utah = add_var!(com, 1, num_colors)
    california = add_var!(com, 1, num_colors)
    ohio = add_var!(com, 1, num_colors)
    illinois = add_var!(com, 1, num_colors)
    washington_dc = add_var!(com, 1, num_colors)
    delaware = add_var!(com, 1, num_colors)
    west_virginia = add_var!(com, 1, num_colors)
    maryland = add_var!(com, 1, num_colors)
    colorado = add_var!(com, 1, num_colors)
    kentucky = add_var!(com, 1, num_colors)
    kansas = add_var!(com, 1, num_colors)
    virginia = add_var!(com, 1, num_colors)
    missouri = add_var!(com, 1, num_colors)
    arizona = add_var!(com, 1, num_colors)
    oklahoma = add_var!(com, 1, num_colors)
    north_carolina = add_var!(com, 1, num_colors)
    tennessee = add_var!(com, 1, num_colors)
    texas = add_var!(com, 1, num_colors)
    new_mexico = add_var!(com, 1, num_colors)
    alabama = add_var!(com, 1, num_colors)
    mississippi = add_var!(com, 1, num_colors)
    georgia = add_var!(com, 1, num_colors)
    south_carolina = add_var!(com, 1, num_colors)
    arkansas = add_var!(com, 1, num_colors)
    louisiana = add_var!(com, 1, num_colors)
    florida = add_var!(com, 1, num_colors)
    michigan = add_var!(com, 1, num_colors)

    states = [washington,montana,maine,north_dakota,south_dakota,wyoming,wisconsin,idaho,vermont,minnesota,oregon,new_hampshire,
    iowa,massachusetts,nebraska,new_york,pennsylvania,connecticut,rhode_island,new_jersey,indiana,nevada,utah,california,ohio,
    illinois,washington_dc,delaware,west_virginia,maryland,colorado,kentucky,kansas,virginia,missouri,arizona,oklahoma,north_carolina,
    tennessee,texas,new_mexico,alabama,mississippi,georgia,south_carolina,arkansas,louisiana,florida,michigan]

    add_constraint!(com, washington != oregon)
    add_constraint!(com, washington != idaho)
    add_constraint!(com, oregon != idaho)
    add_constraint!(com, oregon != nevada)
    add_constraint!(com, oregon != california)
    add_constraint!(com, california != nevada)
    add_constraint!(com, california != arizona)
    add_constraint!(com, nevada != idaho)
    add_constraint!(com, nevada != utah)
    add_constraint!(com, nevada != arizona)
    add_constraint!(com, idaho != montana)
    add_constraint!(com, idaho != wyoming)
    add_constraint!(com, idaho != utah)
    add_constraint!(com, utah != wyoming)
    add_constraint!(com, utah != colorado)
    add_constraint!(com, utah != new_mexico)
    add_constraint!(com, utah != arizona)
    add_constraint!(com, arizona != colorado)
    add_constraint!(com, arizona != new_mexico)
    add_constraint!(com, montana != north_dakota)
    add_constraint!(com, montana != south_dakota)
    add_constraint!(com, montana != wyoming)
    add_constraint!(com, wyoming != south_dakota)
    add_constraint!(com, wyoming != nebraska)
    add_constraint!(com, wyoming != colorado)
    add_constraint!(com, colorado != nebraska)
    add_constraint!(com, colorado != kansas)
    add_constraint!(com, colorado != oklahoma)
    add_constraint!(com, colorado != new_mexico)
    add_constraint!(com, new_mexico != oklahoma)
    add_constraint!(com, new_mexico != texas)
    add_constraint!(com, north_dakota != minnesota)
    add_constraint!(com, north_dakota != south_dakota)
    add_constraint!(com, south_dakota != minnesota)
    add_constraint!(com, south_dakota != iowa)
    add_constraint!(com, south_dakota != nebraska)
    add_constraint!(com, nebraska != iowa)
    add_constraint!(com, nebraska != missouri)
    add_constraint!(com, nebraska != kansas)
    add_constraint!(com, kansas != missouri)
    add_constraint!(com, kansas != oklahoma)
    add_constraint!(com, oklahoma != arkansas)
    add_constraint!(com, oklahoma != texas)
    add_constraint!(com, texas != arkansas)
    add_constraint!(com, texas != louisiana)
    add_constraint!(com, minnesota != wisconsin)
    add_constraint!(com, minnesota != iowa)
    add_constraint!(com, iowa != wisconsin)
    add_constraint!(com, iowa != illinois)
    add_constraint!(com, iowa != missouri)
    add_constraint!(com, missouri != illinois)
    add_constraint!(com, missouri != kentucky)
    add_constraint!(com, missouri != tennessee)
    add_constraint!(com, missouri != arkansas)
    add_constraint!(com, arkansas != tennessee)
    add_constraint!(com, arkansas != mississippi)
    add_constraint!(com, arkansas != louisiana)
    add_constraint!(com, louisiana != mississippi)
    add_constraint!(com, wisconsin != illinois)
    add_constraint!(com, wisconsin != michigan)
    add_constraint!(com, illinois != indiana)
    add_constraint!(com, illinois != kentucky)
    add_constraint!(com, kentucky != indiana)
    add_constraint!(com, kentucky != ohio)
    add_constraint!(com, kentucky != west_virginia)
    add_constraint!(com, kentucky != virginia)
    add_constraint!(com, kentucky != tennessee)
    add_constraint!(com, tennessee != virginia)
    add_constraint!(com, tennessee != north_carolina)
    add_constraint!(com, tennessee != georgia)
    add_constraint!(com, tennessee != alabama)
    add_constraint!(com, tennessee != mississippi)
    add_constraint!(com, mississippi != alabama)
    add_constraint!(com, michigan != indiana)
    add_constraint!(com, michigan != ohio)
    add_constraint!(com, indiana != ohio)
    add_constraint!(com, alabama != georgia)
    add_constraint!(com, alabama != florida)
    add_constraint!(com, ohio != pennsylvania)
    add_constraint!(com, ohio != west_virginia)
    add_constraint!(com, maine != new_hampshire)
    add_constraint!(com, new_hampshire != vermont)
    add_constraint!(com, new_hampshire != massachusetts)
    add_constraint!(com, vermont != massachusetts)
    add_constraint!(com, vermont != new_york)
    add_constraint!(com, massachusetts != new_york)
    add_constraint!(com, massachusetts != rhode_island)
    add_constraint!(com, massachusetts != connecticut)
    add_constraint!(com, connecticut != rhode_island)
    add_constraint!(com, new_york != connecticut)
    add_constraint!(com, new_york != pennsylvania)
    add_constraint!(com, new_york != new_jersey)
    add_constraint!(com, pennsylvania != new_jersey)
    add_constraint!(com, pennsylvania != delaware)
    add_constraint!(com, pennsylvania != maryland)
    add_constraint!(com, pennsylvania != west_virginia)
    add_constraint!(com, new_jersey != delaware)
    add_constraint!(com, maryland != washington_dc)
    add_constraint!(com, maryland != west_virginia)
    add_constraint!(com, maryland != virginia)
    add_constraint!(com, washington_dc != virginia)
    add_constraint!(com, west_virginia != virginia)
    add_constraint!(com, virginia != north_carolina)
    add_constraint!(com, north_carolina != south_carolina)
    add_constraint!(com, south_carolina != georgia)
    add_constraint!(com, georgia != florida)

    # :should be :Min or :Max
    @test_throws ErrorException set_objective!(com, :Minimize, CS.vars_max(states))

    set_objective!(com, :Min, CS.vars_max(states))

    status = solve!(com; keep_logs=true)
    CS.save_logs(com, "graph_color_optimize.json")
    rm("graph_color_optimize.json")

    @test status == :Solved
    @test com.best_sol == 4
    @test all([CS.isfixed(var) for var in states])
    @test maximum([value(var) for var in states]) == 4

end

@testset "49 US states + DC only 3 colors" begin
    com = CS.init()

    num_colors = 3
    washington = add_var!(com, 1, num_colors)
    montana = add_var!(com, 1, num_colors)
    maine = add_var!(com, 1, num_colors)
    north_dakota = add_var!(com, 1, num_colors)
    south_dakota = add_var!(com, 1, num_colors)
    wyoming = add_var!(com, 1, num_colors)
    wisconsin = add_var!(com, 1, num_colors)
    idaho = add_var!(com, 1, num_colors)
    vermont = add_var!(com, 1, num_colors)
    minnesota = add_var!(com, 1, num_colors)
    oregon = add_var!(com, 1, num_colors)
    new_hampshire = add_var!(com, 1, num_colors)
    iowa = add_var!(com, 1, num_colors)
    massachusetts = add_var!(com, 1, num_colors)
    nebraska = add_var!(com, 1, num_colors)
    new_york = add_var!(com, 1, num_colors)
    pennsylvania = add_var!(com, 1, num_colors)
    connecticut = add_var!(com, 1, num_colors)
    rhode_island = add_var!(com, 1, num_colors)
    new_jersey = add_var!(com, 1, num_colors)
    indiana = add_var!(com, 1, num_colors)
    nevada = add_var!(com, 1, num_colors)
    utah = add_var!(com, 1, num_colors)
    california = add_var!(com, 1, num_colors)
    ohio = add_var!(com, 1, num_colors)
    illinois = add_var!(com, 1, num_colors)
    washington_dc = add_var!(com, 1, num_colors)
    delaware = add_var!(com, 1, num_colors)
    west_virginia = add_var!(com, 1, num_colors)
    maryland = add_var!(com, 1, num_colors)
    colorado = add_var!(com, 1, num_colors)
    kentucky = add_var!(com, 1, num_colors)
    kansas = add_var!(com, 1, num_colors)
    virginia = add_var!(com, 1, num_colors)
    missouri = add_var!(com, 1, num_colors)
    arizona = add_var!(com, 1, num_colors)
    oklahoma = add_var!(com, 1, num_colors)
    north_carolina = add_var!(com, 1, num_colors)
    tennessee = add_var!(com, 1, num_colors)
    texas = add_var!(com, 1, num_colors)
    new_mexico = add_var!(com, 1, num_colors)
    alabama = add_var!(com, 1, num_colors)
    mississippi = add_var!(com, 1, num_colors)
    georgia = add_var!(com, 1, num_colors)
    south_carolina = add_var!(com, 1, num_colors)
    arkansas = add_var!(com, 1, num_colors)
    louisiana = add_var!(com, 1, num_colors)
    florida = add_var!(com, 1, num_colors)
    michigan = add_var!(com, 1, num_colors)

    states = [washington,montana,maine,north_dakota,south_dakota,wyoming,wisconsin,idaho,vermont,minnesota,oregon,new_hampshire,
    iowa,massachusetts,nebraska,new_york,pennsylvania,connecticut,rhode_island,new_jersey,indiana,nevada,utah,california,ohio,
    illinois,washington_dc,delaware,west_virginia,maryland,colorado,kentucky,kansas,virginia,missouri,arizona,oklahoma,north_carolina,
    tennessee,texas,new_mexico,alabama,mississippi,georgia,south_carolina,arkansas,louisiana,florida,michigan]

    add_constraint!(com, washington != oregon)
    add_constraint!(com, washington != idaho)
    add_constraint!(com, oregon != idaho)
    add_constraint!(com, oregon != nevada)
    add_constraint!(com, oregon != california)
    add_constraint!(com, california != nevada)
    add_constraint!(com, california != arizona)
    add_constraint!(com, nevada != idaho)
    add_constraint!(com, nevada != utah)
    add_constraint!(com, nevada != arizona)
    add_constraint!(com, idaho != montana)
    add_constraint!(com, idaho != wyoming)
    add_constraint!(com, idaho != utah)
    add_constraint!(com, utah != wyoming)
    add_constraint!(com, utah != colorado)
    add_constraint!(com, utah != new_mexico)
    add_constraint!(com, utah != arizona)
    add_constraint!(com, arizona != colorado)
    add_constraint!(com, arizona != new_mexico)
    add_constraint!(com, montana != north_dakota)
    add_constraint!(com, montana != south_dakota)
    add_constraint!(com, montana != wyoming)
    add_constraint!(com, wyoming != south_dakota)
    add_constraint!(com, wyoming != nebraska)
    add_constraint!(com, wyoming != colorado)
    add_constraint!(com, colorado != nebraska)
    add_constraint!(com, colorado != kansas)
    add_constraint!(com, colorado != oklahoma)
    add_constraint!(com, colorado != new_mexico)
    add_constraint!(com, new_mexico != oklahoma)
    add_constraint!(com, new_mexico != texas)
    add_constraint!(com, north_dakota != minnesota)
    add_constraint!(com, north_dakota != south_dakota)
    add_constraint!(com, south_dakota != minnesota)
    add_constraint!(com, south_dakota != iowa)
    add_constraint!(com, south_dakota != nebraska)
    add_constraint!(com, nebraska != iowa)
    add_constraint!(com, nebraska != missouri)
    add_constraint!(com, nebraska != kansas)
    add_constraint!(com, kansas != missouri)
    add_constraint!(com, kansas != oklahoma)
    add_constraint!(com, oklahoma != arkansas)
    add_constraint!(com, oklahoma != texas)
    add_constraint!(com, texas != arkansas)
    add_constraint!(com, texas != louisiana)
    add_constraint!(com, minnesota != wisconsin)
    add_constraint!(com, minnesota != iowa)
    add_constraint!(com, iowa != wisconsin)
    add_constraint!(com, iowa != illinois)
    add_constraint!(com, iowa != missouri)
    add_constraint!(com, missouri != illinois)
    add_constraint!(com, missouri != kentucky)
    add_constraint!(com, missouri != tennessee)
    add_constraint!(com, missouri != arkansas)
    add_constraint!(com, arkansas != tennessee)
    add_constraint!(com, arkansas != mississippi)
    add_constraint!(com, arkansas != louisiana)
    add_constraint!(com, louisiana != mississippi)
    add_constraint!(com, wisconsin != illinois)
    add_constraint!(com, wisconsin != michigan)
    add_constraint!(com, illinois != indiana)
    add_constraint!(com, illinois != kentucky)
    add_constraint!(com, kentucky != indiana)
    add_constraint!(com, kentucky != ohio)
    add_constraint!(com, kentucky != west_virginia)
    add_constraint!(com, kentucky != virginia)
    add_constraint!(com, kentucky != tennessee)
    add_constraint!(com, tennessee != virginia)
    add_constraint!(com, tennessee != north_carolina)
    add_constraint!(com, tennessee != georgia)
    add_constraint!(com, tennessee != alabama)
    add_constraint!(com, tennessee != mississippi)
    add_constraint!(com, mississippi != alabama)
    add_constraint!(com, michigan != indiana)
    add_constraint!(com, michigan != ohio)
    add_constraint!(com, indiana != ohio)
    add_constraint!(com, alabama != georgia)
    add_constraint!(com, alabama != florida)
    add_constraint!(com, ohio != pennsylvania)
    add_constraint!(com, ohio != west_virginia)
    add_constraint!(com, maine != new_hampshire)
    add_constraint!(com, new_hampshire != vermont)
    add_constraint!(com, new_hampshire != massachusetts)
    add_constraint!(com, vermont != massachusetts)
    add_constraint!(com, vermont != new_york)
    add_constraint!(com, massachusetts != new_york)
    add_constraint!(com, massachusetts != rhode_island)
    add_constraint!(com, massachusetts != connecticut)
    add_constraint!(com, connecticut != rhode_island)
    add_constraint!(com, new_york != connecticut)
    add_constraint!(com, new_york != pennsylvania)
    add_constraint!(com, new_york != new_jersey)
    add_constraint!(com, pennsylvania != new_jersey)
    add_constraint!(com, pennsylvania != delaware)
    add_constraint!(com, pennsylvania != maryland)
    add_constraint!(com, pennsylvania != west_virginia)
    add_constraint!(com, new_jersey != delaware)
    add_constraint!(com, maryland != washington_dc)
    add_constraint!(com, maryland != west_virginia)
    add_constraint!(com, maryland != virginia)
    add_constraint!(com, washington_dc != virginia)
    add_constraint!(com, west_virginia != virginia)
    add_constraint!(com, virginia != north_carolina)
    add_constraint!(com, north_carolina != south_carolina)
    add_constraint!(com, south_carolina != georgia)
    add_constraint!(com, georgia != florida)

    set_objective!(com, :Min, CS.vars_max(states))

    status = solve!(com)

    @test status == :Infeasible
end

end


text:
export Mixture

"""
    mutable struct Mixture{I,O} <: SFunc{I,O}

`Mixture` defines an *sfunc* representing mixtures of other *sfuncs*.  It contains a
vector of *sfuncs* and a vector of probabilities that those *sfuncs* are selected,
which indices are keys associating the two.  The output type of a `Mixture` is defined by
the output type of its internal components.  The parameters of a `Mixture` are
its probabilities followed by parameters of all its internal components, in order.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `cpdf`
- `expectation`
- `compute_pi`
- `send_lambda`

# Type parameters
- `I`: the input type(s) of the `Mixture`
- `O`: the shared output type(s) of its internal components and the output type(s) of the `Mixture`
"""
mutable struct Mixture{I,O} <: SFunc{I,O}
    components::Vector{<:SFunc{I,O}}
    probabilities::Vector{Float64}
end

@impl begin
    struct MixtureSupport end
    function support(sf::Mixture{I,O},
                    parranges::NTuple{N,Vector},
                    size::Integer,
                    curr::Vector{<:O}) where {I,O,N}

        subsize = Int(ceil(size / length(sf.components)))
        ranges = [support(comp, parranges, subsize, curr) for comp in sf.components]
        result = []
        for range in ranges
            append!(result, range)
        end
        result = unique(result)
        tresult = convert(Vector{output_type(sf)}, result)
        # tresult = Vector{output_type(sf)}(undef, length(result))
        # copyto!(tresult, result)
        sort!(tresult)
        return tresult
    end
end

#=
@impl begin
    struct MixtureInitialStats end
    initial_stats(sf::Mixture) = [initial_stats(c) for c in sf.components]
end

@impl begin
    struct MixtureAccumulateStats end
    function accumulate_stats(sf::Mixture, existing_stats, new_stats)
        [accumulate_stats(sf.components[i], existing_stats[i], new_stats[i]) for i in 1:length(sf.components)]
    end
end

@impl begin
    struct MixtureExpectedStats end
    function expected_stats(sf::Mixture{I,O},
                          range::__OptVec{<:O},
                          parranges::NTuple{N,Vector},
                          pis::NTuple{M,Dist},
                          child_lambda::Score{<:O}) where {I,O,N,M}
        # The statistics organize mutually exclusive cases.
        # Each case consists of a component i, a parent value p, and a child value c,
        # and represents the #(i selected, c | p, evidence).
        # This is equal to P(i selected | p, evidence) #(c | i selected, p, evidence).
        # P(i selected | p, evidence) is proportional to m.probabilities[i] * \sum_c #(c | i selected, p, evidence).
        # P(C | i selected, p, evidence) is equal to expected_stats(m.components[i], same arguments...)
        compstats = [expected_stats(comp, range, parranges, pis, child_lambda) for comp in sf.components]
        # This is bad. It assumes stats is a Dict, which is okay for tables
        # summed = [Dict((k,sum(v)) for (k,v) in stats) for stats in compstats]
        summed = [sum(sum(v) for v in values(stats)) for stats in compstats]
        pselected = sf.probabilities .* summed
        # return [mult_through(compstats[i], pselected[i]) for i in 1:length(m.components)]
        return [mult_through(compstats[i], sf.probabilities[i]) for i in 1:length(sf.components)]
    end
end

@impl begin
    struct MixtureMaximizeStats end
    function maximize_stats(sf::Mixture, stats)
        probparams = normalize([sum(sum(values(st))) for st in stats])
        compparams = [maximize_stats(sf.components[i], stats[i]) for i in 1:length(sf.components)]
        return (probparams, compparams...)
    end
end
=#

@impl begin
    struct MixtureSupportQuality end

    function support_quality(sf::Mixture{I,O}, parranges) where {I,O}
        q = support_quality_rank(:CompleteSupport)
        for comp in sf.components
            imp = get_imp(MultiInterface.get_policy(), Support, sf, parranges, 0, O[])
            q = min(q, support_quality_rank(support_quality(imp, comp, parranges)))
        end
        return support_quality_from_rank(q)
    end
end

@impl begin
    mutable struct MixtureMakeFactors
        numpartitions::Dict{SFunc, Int64} = Dict{SFunc, Int64}()
    end

    function make_factors(sf::Mixture{I,O},
                        range::__OptVec{<:O},
                        parranges::NTuple{N,Vector},
                        id,
                        parids::Tuple)::Tuple{Vector{<:Scruff.Utils.Factor}, Vector{<:Scruff.Utils.Factor}} where {I,O,N}
        lfactors = Vector{Scruff.Utils.Factor}()
        ufactors = Vector{Scruff.Utils.Factor}()
        numcomps = length(sf.components)
        mixkey = nextkey()
        for (i,comp) in enumerate(sf.components)
            (lcompfactors, ucompfactors) = make_factors(comp, range, parranges, id, parids)
            function process(factors, target)
                for fact in factors
                    dims = [d for d in fact.dims]
                    push!(dims, numcomps)
                    dims = Tuple(dims)
                    keys = [k for k in fact.keys]
                    push!(keys, mixkey)
                    keys = Tuple(keys)
                    entries = Float64[]
                    for e in fact.entries
                        for j = 1:numcomps
                            push!(entries, i == j ? e : 1.0) # 1.0 means irrelevant
                        end
                    end
                    relevantfact = Factor(dims, keys, entries)
                    push!(target, relevantfact)
                end
            end
            process(lcompfactors, lfactors)
            process(ucompfactors, ufactors)
        end
        mixdims = (numcomps,)
        mixkeys = (mixkey,)
        mixentries = sf.probabilities
        mixfact = Factor(mixdims, mixkeys, mixentries)
        push!(lfactors, mixfact)
        push!(ufactors, mixfact)
        return (lfactors, ufactors)
    end
end

@impl begin
    struct MixtureComputePi end
    function compute_pi(sf::Mixture{I,O},
                     range::__OptVec{<:O},
                     parranges::NTuple{N,Vector},
                     incoming_pis::Tuple)::Dist{<:O} where {N,I,O}
        function f(i)
            cp = compute_pi(sf.components[i], range, parranges, incoming_pis)
            sf.probabilities[i] .* [cpdf(cp, (), x) for x in range]
        end
        scaled = [f(i) for i in 1:length(sf.components)]
        result = sum(scaled)
        return Cat(range, normalize(result))
    end
end

@impl begin
    struct MixtureSendLambda end
    function send_lambda(sf::Mixture{I,O},
                       lambda::Score{<:O},
                       range::__OptVec{<:O},
                       parranges::NTuple{N,Vector},
                       incoming_pis::Tuple,
                       parent_ix::Integer)::Score where {N,I,O}

        # Need to make sure the target parent range is a Vector{T} rather than a Vector{Any}
        T = typeof(parranges[parent_ix][1])
        target_parrange :: Vector{T} = parranges[parent_ix]
        lams = [send_lambda(comp, lambda, range, parranges, incoming_pis, parent_ix) for comp in sf.components]
        scores = [[get_score(lams[j], target_parrange[i]) for i in 1:length(target_parrange)] for j in 1:length(sf.components)]
        scaled = sf.probabilities .* scores
        result = zeros(Float64, length(target_parrange))
        for sc in scaled
            result .+= sc
        end
        return SoftScore(target_parrange, result)
    end
end

# This does not seem to fit
@impl begin
    struct MixtureSample end
    function sample(sf::Mixture{I,O}, x::I)::O where {I,O}
        probs = sf.probabilities/sum(sf.probabilities)
        cat = Categorical(probs)
        which_component = rand(cat)
        component = sf.components[which_component]
        return sample(component, x)
    end
end

@impl begin
    struct MixtureCpdf end
    function cpdf(sf::Mixture{I,O}, i::I, o::O)::AbstractFloat where {I,O}
        complpdf = [cpdf(comp, i, o) for comp in sf.components]
        probs = sf.probabilities/sum(sf.probabilities)
        return sum(probs .* complpdf)
    end
end

@impl begin
    struct MixtureExpectation end

    function expectation(sf::Mixture{I,O}, x::I)::O where {I,O}
        probs = sf.probabilities/sum(sf.probabilities)
        cat = Categorical(probs)
        which_component = rand(cat)
        component = sf.components[which_component]
        return expectation(component, x)
    end
end


text:
#######
# sum #
#######

# basic sum #
#-----------#

function Base.sum(x::TrackedArray{V,D}) where {V,D}
    tp = tape(x)
    out = track(sum(value(x)), D, tp)
    record!(tp, SpecialInstruction, sum, x, out)
    return out
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(sum)})
    input = instruction.input
    output = instruction.output
    istracked(input) && increment_deriv!(input, deriv(output))
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(sum)})
    input = instruction.input
    value!(instruction.output, sum(value(input)))
    return nothing
end

# sum over dimensions #
#---------------------#

function record_sum!(y::TrackedArray, x)
    sum!(value(y), value(x))
    record!(tape(x, y), SpecialInstruction, sum!, x, y, index_bound(y, x))
    return y
end

Base.sum!(y::TrackedArray, x::TrackedArray) = record_sum!(y, x)
Base.sum!(y::TrackedArray, x::AbstractArray) = record_sum!(y, x)

function Base.sum(x::TrackedArray{V,D}, dims) where {V,D}
    tp = tape(x)
    out = track(sum(value(x), dims), D, tp)
    record!(tp, SpecialInstruction, sum!, x, out, index_bound(out, x))
    return out
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(sum!)})
    input = instruction.input
    output = instruction.output
    bound = instruction.cache
    istracked(input) && reduction_increment_deriv!(input, deriv(output), bound)
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(sum!)})
    sum!(value(instruction.output), value(instruction.input))
    return nothing
end

########
# mean #
########

function Statistics.mean(x::TrackedArray{V,D}) where {V,D}
    tp = tape(x)
    out = track(mean(value(x)), D, tp)
    record!(tp, SpecialInstruction, mean, x, out)
    return out
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(mean)})
    input = instruction.input
    output = instruction.output
    istracked(input) && increment_deriv!(input, inv(length(input)) * deriv(output))
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(mean)})
    input = instruction.input
    value!(instruction.output, mean(value(input)))
    return nothing
end

#######
# dot #
#######

function record_dot(x, y, ::Type{D}) where D
    tp = tape(x, y)
    out = track(dot(value(x), value(y)), D, tp)
    cache = (similar(x, D), similar(y, D))
    record!(tp, SpecialInstruction, dot, (x, y), out, cache)
    return out
end

LinearAlgebra.dot(x::TrackedArray{X,D}, y::TrackedArray{Y,D}) where {X,Y,D} = record_dot(x, y, D)

for A in ARRAY_TYPES
    @eval LinearAlgebra.dot(x::TrackedArray{X,D}, y::$A) where {X,D} = record_dot(x, y, D)
    @eval LinearAlgebra.dot(x::$A, y::TrackedArray{Y,D}) where {Y,D} = record_dot(x, y, D)
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(dot)})
    a, b = instruction.input
    pull_value!(a)
    pull_value!(b)
    value!(instruction.output, dot(value(a), value(b)))
    return nothing
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(dot)})
    a, b = instruction.input
    a_tmp, b_tmp = instruction.cache
    output = instruction.output
    output_deriv = deriv(output)
    if istracked(a)
        copyto!(b_tmp, value(b))
        lmul!(output_deriv, b_tmp)
        increment_deriv!(a, b_tmp)
    end
    if istracked(b)
        copyto!(a_tmp, value(a))
        lmul!(output_deriv, a_tmp)
        increment_deriv!(b, a_tmp)
    end
    unseed!(output)
    return nothing
end


text:
export invariant_set, invariant_sets, invariant_sets!, algebraiclift

using FillArrays

#function r(A::Matrix{T}, c::Vector{T}=zeros(T, size(A, 1))) where T
#    [one(T) zeros(T, 1, size(A, 2))
#     c A]
#end

# x+ = Ax + Bu, x ∈ X, u ∈ U
# y+ = [A B; 0 0]y [0; I]w, y ∈ X × U
# [I 0]y+ = [A B]y, y ∈ X × U
function algebraiclift(s::ConstrainedLinearControlDiscreteSystem)
    A = [s.A s.B]
    E = [Matrix(one(eltype(s.A)) * I, size(s.A)...) zeros(eltype(s.A), size(s.B)...)]
    return ConstrainedLinearAlgebraicDiscreteSystem(A, E, stateset(s) * inputset(s))
end
function algebraiclift(s::ConstrainedLinearControlDiscreteSystem{T, MTA, MTB, ST, FullSpace}) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST}
    as = algebraiclift(LinearControlDiscreteSystem(s.A, s.B))
    return ConstrainedLinearAlgebraicDiscreteSystem(as.A, as.E, stateset(s))
end
function algebraiclift(s::LinearControlDiscreteSystem)
    n = statedim(s)
    # We threat zero rows separately to avoid having approximate identity in `E`
    # as it can lead to numerical difficuly for the SDP solver down the road.
    z = findall(i -> all(j -> iszero(s.B[i, j]), 1:size(s.B, 2)), 1:n)
    nz = setdiff(1:n, z)
    Ez = Matrix(one(eltype(s.A)) * I, n, n)[z, :]
    null = nullspace(s.B[nz, :]')'
    Enz = zeros(eltype(null), size(null, 1), n)
    Enz[:, nz] = null
    return LinearAlgebraicDiscreteSystem([s.A[z, :]; Enz * s.A], [Ez; Enz])
end
algebraiclift(s::ConstrainedDiscreteIdentitySystem) = s
algebraiclift(S::AbstractVector) = algebraiclift.(S)
algebraiclift(S::Fill) = Fill(algebraiclift(first(S)), length(S))
function algebraiclift(h::HybridSystem)
    HybridSystem(h.automaton, algebraiclift(h.modes), algebraiclift(h.resetmaps), h.switchings)
end

function constrain_invariance(model, ::ContinuousIdentitySystem, set) end
function constrain_invariance(model, s::Union{ConstrainedContinuousIdentitySystem, ConstrainedDiscreteIdentitySystem}, set)
    return @constraint(model, set ⊆ stateset(s))
end

function constrain_invariance(model, ::IdentityMap, source_set, target_set, λ) end
function constrain_invariance(model, s::LinearMap, source_set, target_set, λ)
    return @constraint(model, s.A * source_set ⊆ target_set)
end
function constrain_invariance(model, s::LinearAlgebraicDiscreteSystem,
                              source_set, target_set, λ)
    return @constraint(model, s.A * source_set ⊆ s.E * target_set,
                       S_procedure_scaling = λ)
end

function Polyhedra.translate(system::ConstrainedLinearAlgebraicDiscreteSystem, c)
    return ConstrainedAffineAlgebraicDiscreteSystem(
        system.A, system.E, system.E * c - system.A * c,
        Polyhedra.translate(stateset(system), c))
end

struct ConstrainedAffineAlgebraicDiscreteSystem{T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}, VTB <: AbstractVector{T}, ST} <: AbstractDiscreteSystem
    A::MTA
    E::MTE
    b::VTB
    X::ST
end
MathematicalSystems.statedim(s::ConstrainedAffineAlgebraicDiscreteSystem) = size(s.A, 1)
MathematicalSystems.stateset(s::ConstrainedAffineAlgebraicDiscreteSystem) = s.X
MathematicalSystems.inputdim(::ConstrainedAffineAlgebraicDiscreteSystem) = 0
MathematicalSystems.islinear(::ConstrainedAffineAlgebraicDiscreteSystem) = false
MathematicalSystems.isaffine(::ConstrainedAffineAlgebraicDiscreteSystem) = true

struct AffineAlgebraicDiscreteSystem{T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}, VTB <: AbstractVector{T}} <: AbstractDiscreteSystem
    A::MTA
    E::MTE
    b::VTB
end
MathematicalSystems.statedim(s::AffineAlgebraicDiscreteSystem) = size(s.A, 1)
MathematicalSystems.inputdim(::AffineAlgebraicDiscreteSystem) = 0
MathematicalSystems.islinear(::AffineAlgebraicDiscreteSystem) = false
MathematicalSystems.isaffine(::AffineAlgebraicDiscreteSystem) = true

function constrain_invariance(model, s::AffineAlgebraicDiscreteSystem,
                              source_set, target_set, λ)
    return @constraint(model, s.A * source_set + s.b ⊆ s.E * target_set,
                       S_procedure_scaling = λ)
end

function isinfeasible(status::Tuple{MOI.TerminationStatusCode, MOI.ResultStatusCode, MOI.ResultStatusCode, String}, no_objective::Bool)
    # Mosek may find an INFEASIBILITY_CERTIFICATE with a STALL/MOI.SLOW_PROGRESS status.
    # SDPA may prove that the problem is INFEASIBILITY_OR_UNBOUNDED which means INFEASIBLE if there is no objective but it has not INFEASIBILITY_CERTIFICATE.
    status[3] == MOI.INFEASIBILITY_CERTIFICATE || status[1] == MOI.INFEASIBLE || (no_objective && status[1] == MOI.INFEASIBLE_OR_UNBOUNDED)
end
function isfeasible(status::Tuple{MOI.TerminationStatusCode, MOI.ResultStatusCode, MOI.ResultStatusCode, String}, no_objective::Bool)
    status[2] == MOI.FEASIBLE_POINT
end
function isdecided(status::Tuple{MOI.TerminationStatusCode, MOI.ResultStatusCode, MOI.ResultStatusCode, String}, no_objective::Bool)
    return isinfeasible(status, no_objective) || isfeasible(status, no_objective)
end

"""
    invariant_sets!(sets, modes_to_compute, system::AbstractHybridSystem,
                    args...; kwargs...)

Similar to `invariant_sets(system, args...; kwargs...)` but stores the result
in `sets` and only compute the modes in `modes_to_compute`. The other sets in
`enabled` that are not in `modes_to_compute` are assumed to have the value given
in `sets`.
"""
function invariant_sets! end

"""
    invariant_sets(system::AbstractHybridSystem, factory::JuMP.OptimizerFactory,
                   set_variables::AbstractVector{<:SetProg.AbstractVariable};
                   volume_heuristic = nth_root,
                   infeasibility_certificates = nothing,
                   verbose=1,
                   λ=Dict{transitiontype(system), Float64}(),
                   enabled = states(system))

Compute maximal invariant sets of the family `set_variables` for the modes of
`system` using the solver provided by `factory`. The volume of the sets is
estimated using `volume_heuristic`. If the program is infeasible, the
certificates for each transition are stored in `infeasibility_certificates`.
For the containment of non-homogeneous, the S-procedure might be a Bilinear
Matrix Inequality (BMI) which is NP-hard to solve. To avoid that, provides the
value of `λ` to use in the dictionary `λ`. To ignore some state and the
transitions involving these states in the computation, give an `enabled` vector
without them.
"""
function invariant_sets end

function default_variable(system::AbstractSystem, factory::JuMP.OptimizerFactory)
    return default_variable(stateset(system), factory)
end
function default_variable(p::Polyhedra.Rep, factory::JuMP.OptimizerFactory)
    solver = something(Polyhedra.default_solver(p), factory)
    center, radius = Polyhedra.chebyshevcenter(p, solver)
    return SetProg.Ellipsoid(point = SetProg.InteriorPoint(center))
end

function invariant_sets!(
    sets, modes_to_compute, s::AbstractHybridSystem, factory::JuMP.OptimizerFactory,
    set_variables::AbstractVector{<:SetProg.AbstractVariable} = [default_variable(mode, factory) for mode in s.modes];
    λ=Dict{transitiontype(s), Float64}(),
    enabled = states(s),
    volume_heuristic = nth_root,
    infeasibility_certificates = nothing,
    verbose=1
)
    model = SOSModel(factory)
    #set_vrefs = @variable(model, [q in modes_to_compute], Ellipsoid(point=h[q]))
    # FIXME calls JuMP.variable_type(model, Ellipsoid(point=h[q])) then
    # complains that q is not defined
    set_vrefs = Dict(state => @variable(model, variable_type=set_variables[state]) for state in modes_to_compute)

    if volume_heuristic !== nothing
        @objective(model, Max,
                   sum(set -> volume_heuristic(volume(set)), values(set_vrefs)))
    end

    λouts = Dict{transitiontype(s), JuMP.AffExpr}()
    transition_constraints = HybridSystems.transition_property(
        s, ConstraintRef{Model, SetProg.ConstraintIndex, SetProg.SetShape})

    for q in modes_to_compute
        source_set = set_vrefs[q]
        # Invariance constraint for mode
        constrain_invariance(model, HybridSystems.mode(s, q), source_set)
        # Invariance constraint for transitions
        for t in out_transitions(s, q)
            λin = get(λ, t, nothing)
            target_state = target(s, t)
            if target_state in enabled
                target_set = target_state in modes_to_compute ? set_vrefs[target_state] : sets[target_state]
                if λin !== nothing
                    λouts[t] = λin
                end
                con_ref = constrain_invariance(
                    model, resetmap(s, t), source_set, target_set, λin)
                if con_ref !== nothing
                    transition_constraints[t] = con_ref
                end
            end
        end
    end

    if verbose >= 3
        print(model)
    end

    JuMP.optimize!(model)

    if verbose >= 4
        print(model)
    end

    if verbose >= 1
        try
            @show MOI.get(model, MOI.SolveTime())
        catch err
            if !(err isa ArgumentError)
                rethrow(err)
            end
        end
        @show JuMP.raw_status(model)
        @show JuMP.termination_status(model)
        @show JuMP.primal_status(model)
        @show JuMP.dual_status(model)
        @show JuMP.objective_value(model)
    end

    if verbose >= 2
        for (t, λout) in λouts
            println("λ for $t is $(JuMP.value(λout))")
        end
    end

    status = (JuMP.termination_status(model),
              JuMP.primal_status(model),
              JuMP.dual_status(model),
              JuMP.raw_status(model))

    if isfeasible(status, volume_heuristic === nothing)
        for q in modes_to_compute
            sets[q] = JuMP.value(set_vrefs[q])
        end
    elseif isinfeasible(status, volume_heuristic === nothing) &&
            infeasibility_certificates !== nothing &&
            status[3] != MOI.NO_SOLUTION # if `status[1]` is `MOI.INFEASIBLE_OR_UNBOUNDED`, there might be no certificate.
        for q in modes_to_compute
            for t in out_transitions(s, q)
                if target(s, t) in enabled
                    infeasibility_certificates[t] = SetProg.SumOfSquares.moment_matrix(transition_constraints[t])
                end
            end
        end
    end

    return status
end
function invariant_sets(s::AbstractHybridSystem, args...; kws...)
    sets = HybridSystems.state_property(s, SetProg.Sets.AbstractSet{Float64})
    invariant_sets!(sets, states(s), s, args...; kws...)
    return sets
end

const HybridSystemWithControlResetMap = HybridSystem{<:AbstractAutomaton, <:AbstractSystem, <:LinearControlDiscreteSystem}
function invariant_sets!(sets, modes_to_compute, s::HybridSystemWithControlResetMap, args...; kws...)
    return invariant_sets!(sets, modes_to_compute, algebraiclift(s), args...;
                           kws...)
end
function invariant_sets(s::HybridSystemWithControlResetMap, args...; kws...)
    return invariant_sets(algebraiclift(s), args...; kws...)
end

_mode(s::AbstractContinuousSystem) = s
_mode(s::AbstractDiscreteSystem) = ConstrainedContinuousIdentitySystem(statedim(s), stateset(s))

_resetmap(s::AbstractContinuousSystem) = IdentityMap(statedim(s))
_resetmap(s::ConstrainedDiscreteIdentitySystem) = IdentityMap(statedim(s))
_resetmap(s::ConstrainedLinearDiscreteSystem) = LinearMap(s.A)
# TODO need to add LinearAlgebraicMap to MathematicalSystems
_resetmap(s::ConstrainedLinearAlgebraicDiscreteSystem) = LinearAlgebraicDiscreteSystem(s.A, s.E)
_resetmap(s::ConstrainedAffineAlgebraicDiscreteSystem) = AffineAlgebraicDiscreteSystem(s.A, s.E, s.b)

function invariant_set(
    system::Union{
        ConstrainedContinuousIdentitySystem,
        ConstrainedDiscreteIdentitySystem,
        ConstrainedLinearDiscreteSystem,
        ConstrainedLinearAlgebraicDiscreteSystem,
        ConstrainedAffineAlgebraicDiscreteSystem},
    factory::JuMP.OptimizerFactory,
    set_variable::SetProg.AbstractVariable=default_variable(system, factory);
    λ=nothing,
    kws...)
    hs = HybridSystem(
        HybridSystems.OneStateAutomaton(1),
        [_mode(system)], [_resetmap(system)], [AutonomousSwitching()])
    λs = Dict{transitiontype(hs), Float64}()
    if λ != nothing
        λs[HybridSystems.OneStateTransition(1)] = λ
    end
    sets = invariant_sets(hs, factory, [set_variable]; λ = λs, kws...)
    return sets[1]
end
function invariant_set(s::ConstrainedLinearControlDiscreteSystem{T, MTA, MTB, ST, FullSpace}, args...; kws...) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST}
    return invariant_set(algebraiclift(s), args...; kws...)
end
function invariant_set(s::ConstrainedLinearControlDiscreteSystem, args...; kws...)
    set = invariant_set(algebraiclift(s), args...; kws...)
    return project(set, 1:statedim(s))
end


text:
o=Dict(:h5=>false,
       :mfile=>"train4_p2.jld2",
       :epochs=>30,
       :lr=>0.0001,
       :p=>12,
       :ema=>0.999f0,
       :batchsize=>64,
       :selfattn=>true,
       :gating=>true,
       :shuffle=>true,
       :sorted=>false,
       :prefix=>"train4_p2",
       :vocab_size=>90,
       :embed_size=>300,
       :dhome=>"data/",
       :loadresnet=>false,
       :d=>512)

max_stars_repo_name: feyzaakyurek/MAC-Network

max_issues_repo_name: feyzaakyurek/MAC-Network

max_forks_repo_name: feyzaakyurek/MAC-Network


text:
function countAndFilter!(val::T, list::AbstractArray{T}) where {T <: Integer}
    nocc = 0
    function doFilter(v)
        if v == val
            nocc += 1
            false
        else
            true
        end
    end
    filter!(doFilter, list)
    return nocc
end
function test_countAndFilter()
    a = vcat(-1:4, 1:5, 2:10)
    countAndFilter!(1, a)
    a = vcat(-1:4, 1:5, 2:10)
    countAndFilter!(6, a)
    a
end

function countAndFilter2!(val::T, list::AbstractArray{T}) where {T <: Integer}
    n1 = count(v->v==val, list)
    filter!(v->v!=val, list)
    n1
end
function test_countAndFilter2()
    a = vcat(-1:4, 1:5, 2:10)
    countAndFilter2!(1, a)
    a = vcat(-1:4, 1:5, 2:10)
    countAndFilter2!(6, a)
    a
end

function countAndFilter3!(val::T, list::AbstractArray{T}) where {T <: Integer}
    length(list) - length(filter!(v->v != val, list))
end
function test_countAndFilter3()
    a = vcat(-1:4, 1:5, 2:10)
    countAndFilter3!(1, a)
    a = vcat(-1:4, 1:5, 2:10)
    countAndFilter3!(6, a)
    a
end
using BenchmarkTools

@show @benchmark test_countAndFilter()
@show @benchmark test_countAndFilter2()
@show @benchmark test_countAndFilter3()

t1 = test_countAndFilter()
t2 = test_countAndFilter2()
t3 = test_countAndFilter3()
@assert t1 == t2 == t3


text:
using Revise
using Plots
using Basins
using DynamicalSystems
using RegionTrees

function newton_map(dz,z, p, n)
    f(x) = x^p[1]-1
    df(x)= p[1]*x^(p[1]-1)
    z1 = z[1] + im*z[2]
    dz1 = f(z1)/df(z1)
    z1 = z1 - dz1
    dz[1]=real(z1)
    dz[2]=imag(z1)
    return
end

# dummy function to keep the initializator happy
function newton_map_J(J,z0, p, n)
   return
end

ds = DiscreteDynamicalSystem(newton_map,[0.1, 0.2], [3] , newton_map_J)
integ  = integrator(ds)
res=10
xg=range(-1.,1.,length=res)
yg=range(-1.,1.,length=res)

@time bsn=Basins.basins_map2D_tree(xg, yg, integ; r_init=0.125, r_max=0.125/4)

function evaluate(basin,v)
    return findleaf(basin,v).data
end

function compute_frac_dim(bsn)

    root_n(x,n) = foldl((y,_) -> parent(y), 1:n; init=x)

    N = zeros(Int, 20)
    wd=1e3
    k_mx = 0
    for lf in allcells(bsn.basin)

        k = 1
        while !isnothing(root_n(lf,k))
            k += 1
        end
        N[k] += 1
        wd = min(wd,lf.boundary.widths[1])
        k_mx = max(k_mx,k)
    end
    @show N
    epsi = reverse(wd * 2 .^ range(0,k_mx-1,step=1))
    return epsi,N[1:k_mx]
end

epsi,N = compute_frac_dim(bsn)

elog = log10.(1 ./ epsi)
Nlog = log10.(N)
D = linear_region(elog[4:end],Nlog[4:end])
@show D

plt = plot(xlim=(xg[1], xg[end]), ylim=(yg[1], yg[end]), legend=nothing)

x = range(xg[1], xg[end], length=1000)
y = range(yg[1], yg[end], length=1000)
heatmap!(plt, x, y, (x, y) -> evaluate(bsn.basin, SVector(x, y)), fill=true)

# plt = plot(xlim=(-1.5, 1.5), ylim=(-1.5, 1.5), legend=nothing)
 for leaf in allleaves(bsn.basin)
     v = hcat(collect(vertices(leaf.boundary))...)
     #plot!(plt, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]])
     #@show leaf.data
     plot!(plt, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]], linecolor=:black, lw=0.4)

 end
 plt


text:
using JuMP, JuMPeR, Gurobi, Random, Distributions, LinearAlgebra, DataFrames, Plots

n = 10 # Number of facilities
m = 50 # Number of customers

# Generating random data (please don't change the seeds.)
facilities = 0.6.*rand(MersenneTwister(5), n,2) .+ 0.2;
customers = rand(MersenneTwister(2), m, 2);
c = [LinearAlgebra.norm(customers[i, :] .- facilities[j, :])[1] for j=1:n, i=1:m];
f = rand(MersenneTwister(3), n)*1 .+ 5;
s = rand(MersenneTwister(4), n)*2 .+ 15;
d = rand(MersenneTwister(5), m)*0.5 .+ 0.75

# P matrix
R_D = 0.25
P = [0.2*exp(-1/R_D .*LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
P = (P .>= 0.2*exp(-1/R_D .* R_D)) .* P

""" Cutting plane facility location model.
    Note this is just the nominal problem, but with linear policies y(z) = u + Vz when z = 0. """
function CP_facility_model(c::Matrix, f::Vector)
    n, m = size(c)
    @assert length(f) == n
    model = Model(solver = GurobiSolver(OutputFlag = 0))

    # VARIABLES
    @variable(model, x[1:n], Bin)     # Facility locations
    @variable(model, V[1:n, 1:m, 1:m])  # Linear decision rule w.r.t. z
    @variable(model, u[1:n, 1:m])  # Fixed policy

    # CONSTRAINTS
    for j=1:m # Demand constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for i=1:n];
        @constraint(model, sum(dt .* u) >= d[j])
    end
    for i=1:n # Capacity constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for j=1:m];
        @constraint(model, sum(dt .* u) <= s[i] * x[i])
    end
    # OBJECTIVE
    @variable(model, F)
    @constraint(model, F >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * u[i, j] for i=1:n, j=1:m))
    @objective(model, Min, F)
    return model, x, u, V
end

function apply_heuristic(model, x, u, V)
    @constraint(model, u .>= 0)
    @constraint(model, V .>= 0)
    return
end

""" Finds and adds worst case cuts for the facility location problem. """
function find_wc_cuts(model, x, u, V, xvals, uvals, Vvals, rho, Gamm)
    wc_model = Model(solver = GurobiSolver(OutputFlag = 0)) # suppressing printouts.
    @variable(wc_model, -rho <= z[1:m] <= rho)
    @variable(wc_model, normdummy[1:m] >= 0)
    @constraint(wc_model, [i=1:m], normdummy[i] >= z[i])
    @constraint(wc_model, [i=1:m], normdummy[i] >= -z[i])
    @constraint(wc_model, sum(normdummy) <= Gamm)

    count = length(model.linconstr)
    for i = 1:n # nonnegativity constraints
        for j=1:m
            @objective(wc_model, Min, uvals[i,j] + sum(Vvals[i,j,:] .* z))
            solve(wc_model)
            if getobjectivevalue(wc_model) < -1e-5 # if constraint is violated
                new_z = getvalue(z)
                @constraint(model, u[i,j] + sum(V[i,j,:] .* new_z) >= 0)
            end
        end
    end
    @info("$(length(model.linconstr) - count) nonnegativity cuts added. ")
    count = length(model.linconstr)
    for j = 1:m # demand constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for i=1:n];
        @objective(wc_model, Min, sum(dt .* uvals) + sum(Vvals[:,j,:] * z) - d[j] - (P*z)[j])
        solve(wc_model)
        if getobjectivevalue(wc_model) < -1e-5 # if constraint is violated
            new_z = getvalue(z)
            @constraint(model, sum(dt .* u) + sum(V[:,j,:] * new_z) >= d[j] + (P*new_z)[j])
        end
    end
    @info("$(length(model.linconstr) - count) demand constraint cuts added.")
    count = length(model.linconstr)
    for i = 1:n # capacity constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for j=1:m];
        @objective(wc_model, Min, - sum(dt .* uvals) - sum(Vvals[i,:,:] * z) + s[i] * xvals[i])
        solve(wc_model)
        if getobjectivevalue(wc_model) < -1e-5 # if constraint is violated
            new_z = getvalue(z)
            @constraint(model, sum(dt .* u) + sum(V[i,:,:] * new_z) <= s[i] * x[i])
        end
    end
    @info("$(length(model.linconstr) - count) capacity constraint cuts added.")
    # objective
    @objective(wc_model, Max, sum(f[j] * xvals[j] for j = 1:n) + sum(c[i, j] * (uvals[i, j] + sum(Vvals[i,j,:] .* z)) for i=1:n, j=1:m))
    solve(wc_model)
    if getobjectivevalue(wc_model) > getvalue(model.obj) + 1e-5 # if constraint is violated
        new_z = getvalue(z)
        @constraint(model, model.obj >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * (u[i, j] + sum(V[i,j,:] .* new_z)) for i=1:n, j=1:m))
        @info("Objective cut added.")
    end
end

""" Plots the solution of the facility location model.
Blue circles are active facilities with different capacities.
Orange plus signs are other potential facility locations.
Rays describe connections between facilities and demand nodes.
"""
function plot_solution(model, x, y, cost = nothing)
    plt = scatter(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s .* getvalue(x))
    scatter!(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s, markershape = :+)
    for i=1:n
        for j=1:m
            if getvalue(y[i,j]) >= 1e-10
                plot!([customers[j, 1], facilities[i,1]], [customers[j,2], facilities[i,2]], linewidth = getvalue(y[i,j]), legend=false)

            end
        end
    end
    if cost == nothing
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
                title = "Total cost: $(round(getobjectivevalue(model), sigdigits=5))")
    else
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
        title = "Total cost: $(round(cost,sigdigits=5))")
    end
    println("Facility cost: $(getvalue(sum(f[j] * x[j] for j = 1:n)))")
    println("Transportation cost: $(getvalue(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)))")
    return plt
end

rho = 1
Gamm = 5
model, x, u, V = CP_facility_model(c, f)
apply_heuristic(model, x, u, V)

# Do 15 iterations of cuts
for i=1:15
    @info("Iteration $(i).")
    modelsize = length(model.linconstr)
    solve(model)
    xvals, uvals, Vvals = getvalue(x), getvalue(u), getvalue(V)
    find_wc_cuts(model, x, u, V, xvals, uvals, Vvals , rho, Gamm)
    if length(model.linconstr) == modelsize
        @info("Optimum reached.")
        @info("Optimal cost:$(getobjectivevalue(model)).")
        break
    end
end

# Problem should converge in ~100 adversarial iterations.

plot_solution(model, x, u)


text:
using CUDA

function hello_world()
    @cuprintln("Greetings from block $(blockIdx().x), thread $(threadIdx().x)!")
    return
end

@cuda blocks=2 threads=2 hello_world()
synchronize()


text:
module Taxonomy
    using NCBITaxonomy
    using Phylo
    using AbstractTrees

    include(joinpath(".", "types.jl"))
    export AbstractTaxonomy, AbstractRank, AbstractTaxon, AbstractClassificationTree, TaxonomicRank, NCBITaxonWrapper

    include(joinpath(".", "nodes.jl"))
    export TaxonNode, FixedRankClassificationTree, FlexibleRankClassificationTree
    export taxon

    include(joinpath(".", "taxonomies.jl"))
    export FixedRankTaxonomy, FlexibleRankTaxonomy
    export ranks, tree, root, defaultranks

    include(joinpath(".", "builders.jl"))
    export taxonomy, _addchild!, _addnode!

    include(joinpath(".", "setoperations.jl"))
    export ∈,⊊,⊂,⊆,∪,∩
end


text:
# Contents:
# Sprite type
# Board type

type Sprite
    y::Int
    x::Int
    is_alive::Bool

    Sprite(ypos,xpos) = new(ypos,xpos,true)
end

function move!(s::Sprite, direction::ASCIIString)
    contains(direction, "n") && (s.y -= 1)
    contains(direction, "s") && (s.y += 1)
    contains(direction, "e") && (s.x += 1)
    contains(direction, "w") && (s.x -= 1)
end

function teleport!(s::Sprite, h, w)
    s.y = rand(1:h)
    s.x = rand(1:w)
end

skip(s::Sprite) = nothing

type Board
    height::Int
    width::Int
    init_num_robots::Int
    live_robots::Int
    wait_mode::Bool
    robots_when_waiting::Int

    #=
    The matrix representations is a pair of 3D matrices (1 4D matrix)
    in which each "slice" shows
    the existence of a sprite, robot, or scrap pile at each coordinate
    e.g. matrix_rep[1,2,5,3] will be 1 iff there is a scrap pile at
    coordinate 2,5 in the first 3D matrix.
    Likewise, matrix_rep[2,2,5,1] and [2,2,5,2] denote the existence
    of a sprite or robot at [2,5] in the second 3D matrix.
    After mutations are finished, switch_active_board is called.
    =#

    matrix_rep::Array{Int,4}
    active::Int
    inactive::Int

    sprite::Sprite

    # Constructors
    Board(h=24,w=60,nrobots=10) = begin
        this = new( h , w , nrobots , nrobots , false , 0 , zeros(Int,2,h,w,3) , 1 , 2 )

        #random initialization of sprite and robots
        rand_coords = sample(1:h*w ,nrobots+1, replace = false)
        this.sprite = Sprite((rand_coords[1]-1) ÷ w + 1, (rand_coords[1]-1) ÷ h + 1)
        this.matrix_rep[this.active,this.sprite.y,this.sprite.x,1] = 1
        for i = 2:length(rand_coords)
            slice(this.matrix_rep,this.active,:,:,2)[rand_coords[i]] = 1
        end
        return this
    end

    Board(nrobots::Int) = Board(24,60,nrobots)
end

has_robots(b::Board) = b.live_robots != 0

function robots_chase_sprite!(s::Sprite,old_robot_field::AbstractArray{Int,2}, new_robot_field::AbstractArray{Int,2})
    height, width = size(old_robot_field)
    for y = 1:height, x = 1:width
        if old_robot_field[y,x] == 1
            old_robot_field[y,x] = 0
            new_robot_field[towards(y,s.y) , towards(x,s.x)] += 1
        end
    end
end

function scrap_robots!(robot_field::AbstractArray{Int,2},scrap_field::AbstractArray{Int,2})
    for i = 1:length(robot_field)
        if scrap_field[i] == 1
            robot_field[i] = 0
        elseif robot_field[i] > 1
            robot_field[i] = 0
            scrap_field[i] = 1
        end
    end
end

function copy_scrap_field!(b::Board)
    b.matrix_rep[b.active,:,:,3] = b.matrix_rep[b.inactive,:,:,3]
end

function process_robot_turn!(b::Board)
    old_robot_field = slice(b.matrix_rep,b.active,:,:,2)
    new_robot_field = slice(b.matrix_rep,b.inactive,:,:,2)
    old_scrap_field = slice(b.matrix_rep,b.active,:,:,3)
    new_scrap_field = slice(b.matrix_rep,b.inactive,:,:,3)

    new_scrap_field[:] = old_scrap_field[:]

    robots_chase_sprite!(b.sprite, old_robot_field, new_robot_field)

    scrap_robots!(new_robot_field,new_scrap_field)

    b.live_robots = sum(b.matrix_rep[b.inactive,:,:,2])
end

function switch_active_board!(b::Board)
    b.active, b.inactive = b.inactive, b.active
end

unset_sprite_pos!(b::Board) = b.matrix_rep[b.active,b.sprite.y,b.sprite.x,1] = 0
set_sprite_pos!(b::Board) = b.matrix_rep[b.inactive,b.sprite.y,b.sprite.x,1] = 1

function move_sprite!(b::Board, direction::ASCIIString)
    unset_sprite_pos!(b)
    move!(b.sprite, direction)
    set_sprite_pos!(b)
end

function teleport_sprite!(b::Board)
    unset_sprite_pos!(b)
    teleport!(b.sprite,b.height,b.width)
    set_sprite_pos!(b)
end

function skip_sprite(b::Board) # equivalent to nothing for now
    unset_sprite_pos!(b)
    skip(b.sprite)
    set_sprite_pos!(b)
end

robot_on_square(b::Board,y,x) = b.matrix_rep[b.inactive,y,x,2] == 1
scrap_on_square(b::Board,y,x) = b.matrix_rep[b.inactive,y,x,3] == 1
is_sprite_on_robot(b::Board) = robot_on_square(b,b.sprite.y,b.sprite.x)
is_sprite_on_scrap(b::Board) = scrap_on_square(b,b.sprite.y,b.sprite.x)

function scrap_sprite!(b::Board)
    (is_sprite_on_scrap(b) || is_sprite_on_robot(b)) && (b.sprite.is_alive = false)
end

function enter_wait_mode!(b::Board)
    b.wait_mode = true
    b.robots_when_waiting = b.live_robots
end

function is_inbounds(b::Board, y ,x)
    (y < 1 || y > b.height) && return false
    (x < 1 || x > b.width) && return false
    return true
end

function is_valid(m::Char, b::Board)
    y = b.sprite.y
    x = b.sprite.x
    if m == 'h'
        checky, checkx = y,x-1
    elseif m == 'j'
        checky, checkx = y+1,x
    elseif m == 'k'
        checky, checkx = y-1,x
    elseif m == 'l'
        checky, checkx = y,x+1
    elseif m == 'y'
        checky, checkx = y-1,x-1
    elseif m == 'u'
        checky, checkx = y-1,x+1
    elseif m == 'b'
        checky, checkx = y+1,x-1
    elseif m == 'n'
        checky, checkx = y+1,x+1
    elseif m == ' '
        return !robot_in_dist_one(b,y,x)
    elseif m == 't'
        return true
    elseif m == 'w'
        return true
    else
        return false
    end

    #=if is_inbounds(b,checky, checkx) && !robot_in_dist_one(b,checky, checkx) && !scrap_on_square(b,checky, checkx)=#
    #=else=#
        #=@show is_inbounds(b,checky, checkx)=#
        #=@show !robot_in_dist_one(b,checky, checkx)=#
        #=@show !scrap_on_square(b,checky, checkx)=#
    #=end=#
    return is_inbounds(b,checky, checkx) &&
    !robot_in_dist_one(b,checky, checkx) &&
    !scrap_on_square(b,checky, checkx)
end

function robot_in_dist_one(b::Board, y::Int,x::Int)
    for i = -1:1, j = -1:1
        if is_inbounds(b,y+i,x+j)
            if b.matrix_rep[b.active,y+i,x+j,2] == 1
                return true
            end
        end
    end
    return false
end

function get_score(b::Board)
    retVal = 10 * (b.init_num_robots - b.live_robots)
    if b.live_robots == 0 && b.sprite.is_alive && b.wait_mode
        retVal += b.robots_when_waiting - b.live_robots
    end
    retVal
end

max_stars_repo_name: RatanRSur/robots-from-scratch

max_issues_repo_name: RatanRSur/robots-from-scratch

max_forks_repo_name: RatanRSur/robots-from-scratch


text:
"""
    skipmissingpairwise(fn::Function)::Function
Given a function `fn` with method `fn(x::RealVector,y::RealVector)::Float64` returns a
function `fn_out(X,Y)` that calls `fn` iteratively on the columns of `X` and `Y` and where
missing elements in those columns are pairwise skipped.

`fn_out` has five methods implementing the single-matrix, two-matrix, matrix-vector,
vector-matrix and vector-vector cases.

For example, if `X` is of type `RealOrMissingMatrix` then `fn_out(X)` returns a matrix
`C` where `C[i,j] == fn(skipmissingpairwise(X[:,i],X[:,j])...)`

# Example
```
julia> x = [missing;1;2;5;4;missing];

julia> y = [1;3;missing;2;3;6];

julia> hcat(x,y)
6×2 Matrix{Union{Missing, Int64}}:
  missing  1
 1         3
 2          missing
 5         2
 4         3
  missing  6

julia> KendallTau.skipmissingpairwise(StatsBase.corkendall)(x,y)
-0.8164965809277261

julia> StatsBase.corkendall([1;5;4],[3;2;3])
-0.8164965809277261
```
"""
function skipmissingpairwise(fn::Function)::Function

    function fn_out(X::RealOrMissingMatrix)
        n = size(X, 2)
        #= TODO this hard-wires on diagonal elements to 1.0,
        OK for cor but wrong in general. =#
        C = Matrix{Float64}(I, n, n)
        for j = 2:n
            for i = 1:j - 1
                C[i,j] = C[j,i] = fn(skipmissingpairwise(X[:,j], X[:,i])...)
            end
        end
        return C
    end

    function fn_out(X::RealOrMissingMatrix,Y::RealOrMissingMatrix)
        nr = size(X, 2)
        nc = size(Y, 2)
        C = Matrix{Float64}(undef, nr, nc)
        for j = 1:nr
            for i = 1:nc
                C[j,i] = fn(skipmissingpairwise(X[:,j], Y[:,i])...)
            end
        end
        return C
    end

    #= Return matrix. Consistent with Statistics.cor and corspearman but not with corkendall,
    but maybe that doesn't matter? =#
    function fn_out(X::RealOrMissingMatrix,y::RealOrMissingVector)
        size(X, 1) == length(y) ||
                throw(DimensionMismatch("X and y have inconsistent dimensions"))
            n = size(X, 2)
            return(reshape([fn(skipmissingpairwise(X[:,i], y)...) for i in 1:n], n, 1))
        end

    function fn_out(x::RealOrMissingVector,Y::RealOrMissingMatrix)
        size(Y, 1) == length(x) ||
            throw(DimensionMismatch("x and Y have inconsistent dimensions"))
        n = size(Y, 2)
        return(reshape([fn(skipmissingpairwise(x, Y[:,i])...) for i in 1:n], 1, n))
    end

    function fn_out(x::RealOrMissingVector,y::RealOrMissingVector)
        length(x) == length(y) || throw(DimensionMismatch("Vectors must have same length"))
        fn(skipmissingpairwise(x, y)...)
    end

    return(fn_out)
end

function skipmissingpairwise(x::RealVector, y::RealVector)
    x, y
end

"""
    skipmissingpairwise(x::RealOrMissingVector, y::RealOrMissingVector)
Returns a pair `(a,b)`, filtered copies of `x` and `y`, in which elements `x[i]` and `y[i]`
are "skipped" (filtered out) if either `ismissing(x[i])` or `ismissing(y[i])`.
"""
function skipmissingpairwise(x::RealOrMissingVector{T}, y::RealOrMissingVector{U}) where T where U

    length(x) == length(y) || error("Vectors must have same length")

    T2 = x isa Vector{Missing} ? Missing : T
    U2 = y isa Vector{Missing} ? Missing : U

    nout::Int = 0
    @inbounds for i = 1:length(x)
        if !(ismissing(x[i]) || ismissing(y[i]))
            nout += 1
        end
    end

    res1 = Vector{T2}(undef, nout)
    res2 = Vector{U2}(undef, nout)
    j::Int = 0

    @inbounds for i = 1:length(x)
        if !(ismissing(x[i]) || ismissing(y[i]))
            j += 1
            res1[j] = x[i]
            res2[j] = y[i]
        end
    end

    return(res1, res2)
end


text:
using LinearAlgebra: normalize

"""
    compute_changes_per_element(req_data, res_data)

Computes the number of changes occuring for each element by comparing all the
elements from the request with the received ones.

[Location based metric] Where are the changes occuring in the model?
"""
function compute_changes_per_element(req_data, res_data)::Vector{Float32}
    return sum([
        weights .!= req_data
        for weights in res_data
    ]) / length(res_data)
end

"""
    compute_round_changes(req_data, res_data)

Computes the percentage of changes occuring for the given round by comparing
all the element from the request with the received ones.
"""
function compute_round_changes(req_data, res_data)::Float32
    result = sum(compute_changes_per_element(req_data, res_data))
    result /= length(req_data)
    return result
end

"""
    compute_elements_difference(req_data, res_data, values_range)

Computes the importance of changes occuring for each element by computing the
average absolute difference between the element from the request with the
received ones.

[Location based metric] Where are the changes occuring in the model?
"""
function compute_elements_difference(req_data, res_data)::Vector{Float32}
    return normalize(sum([
        abs.(Float32.(weights) - Float32.(req_data))
        for weights in res_data
    ]) / length(res_data))
end


text:
module PolygonalFem

using LinearAlgebra, StatsBase, SparseArrays, GroupSlices
using NearestNeighbors
using ReverseDiff, CoordinateTransformations, Interpolations
using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile
using FileIO, JLD2, GeometryBasics, Makie, GLMakie, Colors, PyCall

include("assembly.jl")
include("assembly_elas.jl")
include("adjoint.jl")
include("utils.jl")
include("plot.jl")

#-------------------------------------------------------------------------------
if :GraphicTools ∈ names(Main, all = true, imported = true)
    include("solveP1.jl")
    include("_genmeshes.jl")
end
eye(n) = Matrix(I, n, n)
#-------------------------------------------------------------------------------
"""
    u, p, t, meshboundary = vem(filename::String = "squarepolmesh_coarse",
                                nc::Int64 = 1_00;
                                resolution::Int64 = 400)

N.B.
* nc = 100, 1_000 or 10_000
* adapted code from the article / matlab's code  "The virtual element method in 50 lines of matlab". See https://arxiv.org/pdf/1604.06021.pdf
"""
function vem(filename::String = "squarepolmesh_coarse", nc::Int64 = 1_00;
             resolution::Int64 = 400, visible::Bool = false)
    # computes the virtual element solution of the Poisson problem on
    # the specified mesh
    # load the mesh  + pv       : vertices of the cells
    #                + cells    : polygonal cells (⚠ orientation is crucial)
    #                + elem     : triangulated cells
    #                + (pb, tb) : restricted delaunay mesh
    mesh_filename = "$(@__DIR__)/../test/data/$(filename)_$(nc).jld2"
    println("read file $(mesh_filename)")
    JLD2.@load(mesh_filename, ps, Iv, pv, bbelem, elem, pb, tb)
    cells, bcells = bbelem, bbelem
    # all triangles
    t = vcat(elem...)
    # boundary points / rhs
    rhs, boundary_condition = if occursin("square", mesh_filename)
        rhs_sqr, boundary_condition_sqr
    else
        rhs_L, boundary_condition_L
    end
    n_dofs, n_polys = size(pv, 1), 3 # method has 1 degree of freedom per vertex
    # identify boundary points
    meshboundary = unique(btri(t)[:])
    u = zeros(n_dofs) # degrees of freedom of the virtual element solution
    # call assemble function
    IK, JK, SK, IF, SF = assembKM_vemKsource(pv, cells, rhs,
                                             boundary_condition, meshboundary)
    K = sparse(IK, JK, SK, n_dofs, n_dofs)
    F = Vector(sparsevec(IF, SF, n_dofs))
    boundary_vals = boundary_condition(pv[meshboundary, :])
    # vertices which aren’t on the boundary
    internal_dofs =  setdiff(1:n_dofs, meshboundary)
    # solve
    u[internal_dofs] = K[internal_dofs, internal_dofs] \ F[internal_dofs]
    u[meshboundary] .= boundary_vals # set the boundary values
    # plot
    if resolution > 0
        Iv = visible ? Iv : 1:length(cells)
        plotsolution(u, pv, cells, bcells, Iv, ps,
                     mesh_filename = mesh_filename, resolution = resolution)
    end
    return u, pv, t, meshboundary
end
#-------------------------------------------------------------------------------
function rhs_sqr(points)
    x, y = points[:, 1], points[:, 2]
    return 15 * sin.(π * x) .* sin.(π * y)
end
#-------------------------------------------------------------------------------
function boundary_condition_sqr(points)
    x, y = points[:, 1], points[:, 2]
    return (1 .- x) .* y .* sin.(π * x)
    #return zeros(length(x))
end
#-------------------------------------------------------------------------------
function rhs_L(points)
    return zeros(size(points, 1))
end
#-------------------------------------------------------------------------------
function boundary_condition_L(points)
    x, y = points[:, 1], points[:, 2]
    pV = [[x[k], y[k]] for k = 1:length(x)]
    rθ = PolarFromCartesian().(pV)
    r, θ = [rθ[k].r for k = 1:length(x)],  [rθ[k].θ for k = 1:length(x)]
    return  r .^ (2 / 3.) .* sin.( 6 * (θ .- π / 2) / 3)
end
end # module


text:
using DiscreteFiltering
using OrdinaryDiffEq: ODEProblem, solve, Tsit5
using LinearAlgebra
using Plots
using SparseArrays

## Domain
a = 0.0
b = 2π
domain = PeriodicIntervalDomain(a, b)
# domain = ClosedIntervalDomain(a, b)

## Discretization
N = 100
M = N
x = discretize(domain, M)
ξ = discretize(domain, N)
Δx = (b - a) / N

## Filter
h₀ = 2.1Δx
h(x) = h₀ * (1 - 1 / 2 * cos(x))
dh(x) = h₀ / 2 * sin(x)
α(x) = 1 / 3 * dh(x) * h(x)
f = TopHatFilter(h)

# h₀ = 5.1Δx
# h(x) = h₀ # * (1 - 1 / 2 * cos(x))
# dh(x) = 0.0 # h₀ / 2 * sin(x)
# α(x) = 1 / 3 * dh(x) * h(x)
# f = GaussianFilter(h, Δx / 2)

## Time
T = 1.5
t = T

## Plot filter
plot(x, h)
ylims!((0, ylims()[2]))

## Get matrices
C = advection_matrix(domain, N)
D = diffusion_matrix(domain, N)
W = filter_matrix(f, domain, M, N)
R = reconstruction_matrix(f, domain, M, N)
A = spdiagm(α.(ξ))

## Inspect matrices
spy(W)
spy(R)

## Exact solutions
u(ξ, t) = sin(ξ - t) + 3 / 5 * cos(5(ξ - t)) + 1 / 25 * sin(20(ξ - 1 - t))
u_int(ξ, t) = -cos(ξ - t) + 3 / 25 * sin(5(ξ - t)) - 1 / 25 / 20 * cos(20(ξ - 1 - t))
ū(x, t) = 1 / 2h(x) * (u_int(x + h(x), t) - u_int(x - h(x), t))

## Discrete initial conditions
uₕ = u.(ξ, 0.0)
ūₕ = ū.(x, 0.0)
uₕ_allbar = W * uₕ

plot(x, uₕ, label = "Discretized")
plot!(x, ūₕ, label = "Filtered-then-discretized")
plot!(x, uₕ_allbar, label = "Discretized-then-filtered")

## Solve discretized problem
duₕ(uₕ, p, t) = -C * uₕ
prob = ODEProblem(duₕ, uₕ, (0, T))
sol = solve(prob, Tsit5(), abstol = 1e-6, reltol = 1e-4)

plot(x, uₕ, label = "Initial conditions")
plot!(x, sol(t), label = "Discretized")
plot!(x, u.(x, t), label = "Exact")

## Solve filtered-and-then-discretized problem
dūₕ(ūₕ, p, t) = (-C + A * D) * ūₕ
prob_bar = ODEProblem(dūₕ, ūₕ, (0, T))
sol_bar = solve(prob_bar, Tsit5(), abstol = 1e-6, reltol = 1e-4)

plot(x, uₕ, label = "Initial")
plot!(x, ūₕ, label = "Initial filtered")
plot!(x, sol_bar(t), label = "Filtered-then-discretized")
plot!(x, ū.(x, t), label = "Filtered exact")
# plot!(x, 500*α.(ξ))

## Solve discretized-and-then-filtered problem
duₕ_allbar(uₕ_allbar, p, t) = -W * (C * (R * uₕ_allbar))
# duₕ_allbar(uₕ_allbar, p, t) = -W * (C * (W \ uₕ_allbar))
prob_allbar = ODEProblem(duₕ_allbar, W * uₕ, (0, T))
sol_allbar = solve(prob_allbar, Tsit5(), abstol = 1e-6, reltol = 1e-4)

plot(x, uₕ, label = "Initial")
plot!(x, uₕ_allbar, label = "Initial discretized-then-filtered")
plot!(x, sol_allbar(t), label = "Discretized-then-filtered")
plot!(x, W * u.(ξ, t), label = "Exact")

## PGFPlotsX
pgfplotsx()

## Comparison
p = plot(
    size = (400, 300),
    xlabel = "\$x\$",
    legend = :topright,
)
plot!(p, x, uₕ, label = "\$u(x, t = 0.0)\$")
# plot!(p, x, ūₕ, label = "Initial filtered")
# plot!(p, x, sol(t), label = "Discretized")
# plot!(p, x, sol_bar(t), label = "Filtered-then-discretized")
# plot!(p, x, sol_allbar(t), label = "Discretized-then-filtered")
# plot!(p, x, [u.(x, t), ū.(x, t)], label = "Exact")
plot!(p, x, u.(x, t), label = "\$u(x, t = $t)\$")
plot!(p, x, ū.(x, t), label = "\$u(x, t = $t)\$")
# ylims!(p, minimum(uₕ), maximum(uₕ))
display(p)
savefig(p, "output/advection/solution.tikz")

## Relative error
u_exact = u.(ξ, t)
ū_exact = ū.(x, t)
err = abs.(sol(t) - u_exact) ./ maximum(abs.(u_exact))
err_bar = abs.(sol_bar(t) - u_exact) ./ maximum(abs.(u_exact))
err_allbar = abs.(sol_allbar(t) - u_exact) ./ maximum(abs.(u_exact))

plot()
plot!(ξ, err, label = "Unfiltered discretized")
plot!(x, err_bar, label = "Filtered-then-discretized")
plot!(x, err_allbar, label = "Discretized-then-filtered")


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

const T = OrdinaryDiffEq.ODECompositeSolution

@testset "integrator" begin
  prob = DDEProblemLibrary.prob_dde_constant_1delay_ip

  integrator = init(prob, MethodOfSteps(AutoTsit5(Rosenbrock23())))
  @test integrator.sol isa T
  @test integrator.integrator.sol isa T

  sol1 = solve!(integrator)
  @test sol1 isa T

  # compare integration grid
  sol2 = solve(prob, MethodOfSteps(Tsit5()))
  @test sol1.t == sol2.t
  @test sol1.u == sol2.u

  # compare interpolation
  @test sol1(0:0.1:1) == sol2(0:0.1:1)
end

@testset "issue #148" begin
  alg = MethodOfSteps(Tsit5())
  compositealg = MethodOfSteps(CompositeAlgorithm((Tsit5(), RK4()), integrator -> 1))

  prob = DDEProblem((du, u, h, p, t) -> (du[1] = -h(p, t - 1)[1]; nothing),
                    [1.0], (p, t) -> [1.0], (0.0, 5.0))
  sol = solve(prob, alg)
  compositesol = solve(prob, compositealg)
  @test compositesol isa T
  @test sol.t == compositesol.t
  @test sol.u == compositesol.u

  prob_oop = DDEProblem((u, h, p, t) -> -h(p, t - 1), [1.0], (p, t) -> [1.0],
                        (0.0, 5.0))
  sol_oop = solve(prob_oop, alg)
  compositesol_oop = solve(prob_oop, compositealg)
  @test compositesol_oop isa T
  @test sol_oop.t == compositesol_oop.t
  @test sol_oop.u == compositesol_oop.u
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
"""
Handles the functionality for applying various gramatical rules.
"""
module Inflector

import Unicode
import SearchLight

const vowels = ["a", "e", "i", "o", "u"]

"""
    tosingular(word::String; is_irregular::Bool = false) :: Union{Nothing,String}

Returns the singural form of `word`.
"""
function tosingular(word::String; is_irregular::Bool = Inflector.is_irregular(word)) :: Union{Nothing,String}
  ( is_irregular || ! endswith(word, "s") ) && return to_singular_irregular(word)
  endswith(word, "ies") && ! in(word[end-3], vowels) && return (word[1:end-3] * "y")
  endswith(word, "s") && return word[1:end-1]

  nothing
end

"""
    to_singular_irregular(word::String) :: Union{Nothing,String}

Returns the singular form of the irregular word `word`.
"""
function to_singular_irregular(word::String) :: Union{Nothing,String}
  irr = irregular(word)
  irr !== nothing ? irr[1] : nothing
end

"""
    to_plural(word::String; is_irregular::Bool = false) :: Union{Nothing,String}

Returns the plural form of `word`.
"""
function to_plural(word::String; is_irregular::Bool = Inflector.is_irregular(word)) :: Union{Nothing,String}
  is_irregular && return to_plural_irregular(word)
  endswith(word, "y") && ! in(word[end-1], vowels) && return (word[1:end-1] * "ies") # category -> categories // story -> stories
  is_singular(word) ? (word * "s") : (word)
end

"""
    to_plural_irregular(word::String) :: Union{Nothing,String}

Returns the plural form of the irregular word `word`.
"""
function to_plural_irregular(word::String) :: Union{Nothing,String}
  irr = irregular(word)
  irr !== nothing ? (irr)[2] : nothing
end

"""
    from_underscores(word::String) :: String

Generates `SnakeCase` form of `word` from `underscore_case`.
"""
function from_underscores(word::String) :: String
  mapreduce(x -> uppercasefirst(x), *, split(word, "_"))
end

"""
    is_singular(word::String) :: Bool

Returns wether or not `word` is a singular.
"""
function is_singular(word::String) :: Bool
  ! is_plural(word)
end

"""
    is_plural(word::String) :: Bool

Returns wether or not `word` is a plural.
"""
function is_plural(word::String) :: Bool
  word = Unicode.normalize(word, casefold = true)
  irr_word = irregular(word)
  (irr_word !== nothing && word != (irr_word)[1]) ||
    (irr_word !== nothing && word == (irr_word)[2]) ||
      endswith(word, "s")
end

"""
    irregulars() :: Vector{Tuple{String,String}}

Returns a `vector` of words with irregular singular or plural forms.
"""
function irregulars() :: Vector{Tuple{String,String}}
  vcat(IRREGULAR_NOUNS, SearchLight.config.inflector_irregulars)
end

"""
    irregular(word::String) :: Union{Nothing,Tuple{String,String}}

Wether or not `word` has an irregular singular or plural form.
"""
function irregular(word::String) :: Union{Nothing,Tuple{String,String}}
  word = Unicode.normalize(word, casefold = true)

  for (k, v) in IRREGULAR_NOUNS
    (word == k || word == v) && return (k::String, v::String)
  end

  nothing
end

"""
    function is_irregular(word::String) :: Bool

Whether or not `word` has a singular or plural irregular form.
"""
function is_irregular(word::String) :: Bool
  irregular(word) === nothing ? false : true
end

"""
Vector of nouns with irregular forms of singular and/or plural.
Each tuple contains pairs of singular and plural forms.
"""
const IRREGULAR_NOUNS = Vector{Tuple{String,String}}([
  ("alumnus", "alumni"),
  ("cactus", "cacti"),
  ("focus", "foci"),
  ("fungus", "fungi"),
  ("nucleus", "nuclei"),
  ("radius", "radii"),
  ("stimulus", "stimuli"),
  ("axis", "axes"),
  ("analysis", "analyses"),
  ("basis", "bases"),
  ("crisis", "crises"),
  ("diagnosis", "diagnoses"),
  ("ellipsis", "ellipses"),
  ("hypothesis", "hypotheses"),
  ("oasis", "oases"),
  ("paralysis", "paralyses"),
  ("parenthesis", "parentheses"),
  ("synthesis", "syntheses"),
  ("synopsis", "synopses"),
  ("thesis", "theses"),
  ("appendix", "appendices"),
  ("index", "indeces"),
  ("matrix", "matrices"),
  ("beau", "beaux"),
  ("bureau", "bureaus"),
  ("tableau", "tableaux"),
  ("child", "children"),
  ("man", "men"),
  ("ox", "oxen"),
  ("woman", "women"),
  ("bacterium", "bacteria"),
  ("corpus", "corpora"),
  ("criterion", "criteria"),
  ("curriculum", "curricula"),
  ("datum", "data"),
  ("genus", "genera"),
  ("medium", "media"),
  ("memorandum", "memoranda"),
  ("phenomenon", "phenomena"),
  ("stratum", "strata"),
  ("deer", "deer"),
  ("fish", "fish"),
  ("means", "means"),
  ("offspring", "offspring"),
  ("series", "series"),
  ("sheep", "sheep"),
  ("species", "species"),
  ("foot", "feet"),
  ("goose", "geese"),
  ("tooth", "teeth"),
  ("antenna", "antennae"),
  ("formula", "formulae"),
  ("nebula", "nebulae"),
  ("vertebra", "vertebrae"),
  ("vita", "vitae"),
  ("louse", "lice"),
  ("mouse", "mice"),
  ("quiz", "quizzes"),
  ("search", "searches")
])
end