text:
################################################################################
#
#  Predicates
#
################################################################################

issimple(::Type{AnticNumberField}) = true

issimple(::AnticNumberField) = true

################################################################################
#
#  Field constructions
#
################################################################################

doc"""
***
  NumberField(f::fmpq_poly)

>  The number field Q[x]/f generated by f.
"""
function NumberField(f::fmpq_poly; cached::Bool = true)
  return NumberField(f, "_a", cached = cached)
end

function NumberField(f::fmpz_poly, s::Symbol; cached::Bool = true)
  Qx, x = PolynomialRing(FlintQQ, string(parent(f).S))
  return NumberField(Qx(f), String(s), cached = cached)
end

function NumberField(f::fmpz_poly, s::AbstractString; cached::Bool = true)
  Qx, x = PolynomialRing(FlintQQ, string(parent(f).S))
  return NumberField(Qx(f), s, cached = cached)
end

function NumberField(f::fmpz_poly; cached::Bool = true)
  Qx, x = PolynomialRing(FlintQQ, string(parent(f).S))
  return NumberField(Qx(f), cached = cached)
end

global const number_field = NumberField

doc"""
    pure_extension(n::Int, gen::Integer) -> AnticNumberField, nf_elem
    pure_extension(n::Int, gen::fmpz) -> AnticNumberField, nf_elem

The number field with defining polynomial $x^n-gen$.
"""
function pure_extension(n::Int, gen::Integer; cached::Bool = true)
  return pure_extension(n, fmpz(gen), cached = cached)
end

function pure_extension(n::Int, gen::fmpz; cached::Bool = true)
  kx, x = FlintQQ["x"]
  return number_field(x^n - gen, cached = cached)
end

doc"""
    cyclotomic_field(n::Int) -> AnticNumberField, nf_elem

The $n$-th cyclotomic field defined by the $n$-the cyclotomic polynomial.
"""
function cyclotomic_field(n::Int; cached::Bool = true)
  Zx, x = FlintZZ["x"]
  return number_field(cyclotomic(n, x), "z_$n", cached = cached)
end

#fields with usually large class groups...
# TODO: Some sort of reference?
doc"""
    wildanger_field(n::Int, B::fmpz) -> AnticNumberField, nf_elem
    wildanger_field(n::Int, B::Integer) -> AnticNumberField, nf_elem

Returns the field with defining polynomial $x^n + \sum_{i=0}^n (-1)^iBx^i$.
These fields tend to have non-trivial class groups.
"""
function wildanger_field(n::Int, B::fmpz; cached::Bool = true)
  Qx, x = PolynomialRing(FlintQQ)
  f = x^n
  for i=0:n-1
    f += (-1)^(n-i)*B*x^i
  end
  return NumberField(f, "_\$", cached = cached)
end

function wildanger_field(n::Int, B::Integer; cached::Bool = true)
  return wildanger_field(n, fmpz(B), cached = cached)
end

################################################################################
#
#  Characteristic
#
################################################################################

characteristic(::AnticNumberField) = 0

################################################################################
#
#  Predicates
#
################################################################################

doc"""
    ispure_extension(K::AnticNumberField) -> Bool

Tests if $K$ is pure, that is, if the defining polynomial is $x^n - g$ for some
rational number $g$.
"""
function ispure_extension(K::AnticNumberField)
  if !ismonic(K.pol)
    return false
  end
  return all(i->iszero(coeff(K.pol, i)), 1:degree(K)-1)
end

doc"""
    iskummer_extension(K::AnticNumberField) -> Bool

Tests if $K$ is a Kummer extension of $\mathbb Q$, that is, if the defining
polynomial is of the form $x^2 - g$.
"""
function iskummer_extension(K::AnticNumberField)
  if degree(K) != 2
    return false
  end
  return ispure_extension(K)
end

################################################################################
#
#  Class group
#
################################################################################

doc"""
    class_group(K::AnticNumberField) -> GrpAbFinGen, Map

Shortcut for {{{class_group(maximal_order(K))}}}: returns the class
group as an abelian group and a map from this group to the set
of ideals of the maximal order.
"""
function class_group(K::AnticNumberField)
  return class_group(maximal_order(K))
end

################################################################################
#
#  Basis
#
################################################################################

doc"""
    basis(K::AnticNumberField) -> Vector{nf_elem}

Returns the power basis of $K$, that is, the elements $1,a,\dotsc,a^{d - 1}$,
where $d$ is the degree of $K$.
"""
function basis(K::AnticNumberField)
  n = degree(K)
  g = gen(K);
  d = Array{typeof(g)}(n)
  b = K(1)
  for i = 1:n-1
    d[i] = b
    b *= g
  end
  d[n] = b
  return d
end

################################################################################
#
#  Torsion units and related functions
#
################################################################################

doc"""
***
    istorsion_unit(x::nf_elem, checkisunit::Bool = false) -> Bool

Returns whether $x$ is a torsion unit, that is, whether there exists $n$ such
that $x^n = 1$.

If `checkisunit` is `true`, it is first checked whether $x$ is a unit of the
maximal order of the number field $x$ is lying in.
"""
function istorsion_unit(x::nf_elem, checkisunit::Bool = false)
  if checkisunit
    _isunit(x) ? nothing : return false
  end

  K = parent(x)
  d = degree(K)
  c = conjugate_data_arb(K)
  r, s = signature(K)

  while true
    @vprint :UnitGroup 2 "Precision is now $(c.prec) \n"
    l = 0
    @vprint :UnitGroup 2 "Computing conjugates ... \n"
    cx = conjugates_arb(x, c.prec)
    A = ArbField(c.prec, false)
    for i in 1:r
      k = abs(cx[i])
      if k > A(1)
        return false
      elseif isnonnegative(A(1) + A(1)//A(6) * log(A(d))//A(d^2) - k)
        l = l + 1
      end
    end
    for i in 1:s
      k = abs(cx[r + i])
      if k > A(1)
        return false
      elseif isnonnegative(A(1) + A(1)//A(6) * log(A(d))//A(d^2) - k)
        l = l + 1
      end
    end

    if l == r + s
      return true
    end
    refine(c)
  end
end

doc"""
***
    torsion_unit_order(x::nf_elem, n::Int)

Given a torsion unit $x$ together with a multiple $n$ of its order, compute
the order of $x$, that is, the smallest $k \in \mathbb Z_{\geq 1}$ such
that $x^`k` = 1$.

It is not checked whether $x$ is a torsion unit.
"""
function torsion_unit_order(x::nf_elem, n::Int)
  # TODO:
  # This is lazy
  # Someone please change this
  y = deepcopy(x)
  for i in 1:n
    if isone(y)
      return i
    end
    mul!(y, y, x)
  end
  error("Something odd in the torsion unit order computation")
end

#################################################################################################
#
#  Normal Basis
#
#################################################################################################

doc"""
***
    normal_basis(K::Nemo.AnticNumberField) -> nf_elem

Given a number field K which is normal over Q, return
an element generating a normal basis of K over Q.
"""
function normal_basis(K::Nemo.AnticNumberField)
  n = degree(K)
  Aut = automorphisms(K)

  length(Aut) != degree(K) && error("The field is not normal over the rationals!")

  A = zero_matrix(FlintQQ, n, n)
  r = one(K)
  while true
    r = rand(basis(K), -n:n)
    for i = 1:n
      y = Aut[i](r)
      for j = 1:n
        A[i,j] = coeff(y, j - 1)
      end
    end
    if !iszero(det(A))
      break
    end
  end
  return r
end

################################################################################
#
#  Subfield check
#
################################################################################

function _issubfield(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  g = L.pol
  Lx, x = PolynomialRing(L, "x", cached = false)
  fL = Lx()
  rts = roots(f, L)
  for r in rts
  #for i = 0:degree(f)
  #  setcoeff!(fL, i, L(coeff(f, i)))
  #end
  #fac = factor(fL)
  #for (a, b) in fac
  #  if degree(a) == 1
  #    c1 = coeff(a, 0)
  #    c2 = coeff(a, 1)
  #    h = parent(L.pol)(-c1*inv(c2))
  #    return true, NfToNfMor(K, L, h(gen(L)))
    h = parent(L.pol)(r)
    return true, NfToNfMor(K, L, h(gen(L)))
  end
  return false, NfToNfMor(K, L, L())
end

doc"""
***
      issubfield(K::AnticNumberField, L::AnticNumberField) -> Bool, NfToNfMor

> Returns "true" and an injection from $K$ to $L$ if $K$ is a subfield of $L$.
> Otherwise the function returns "false" and a morphism mapping everything to 0.
"""
function issubfield(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  g = L.pol
  if mod(degree(g), degree(f)) != 0
    return false, NfToNfMor(K, L, L())
  end
  t = divexact(degree(g), degree(f))
  try
    OK = _get_maximal_order_of_nf(K)
    OL = _get_maximal_order_of_nf(L)
    if mod(discriminant(OL), discriminant(OK)^t) != 0
      return false, NfToNfMor(K, L, L())
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    # We could factorize the discriminant of f, but we only test small primes.
    p = 3
    df = discriminant(f)
    dg = discriminant(g)
    while p < 10000
      if p > df || p > dg
        break
      end
      if mod(valuation(df, p), 2) == 0
        p = next_prime(p)
        continue
      end
      if mod(dg, p^t) != 0
        return false, NfToNfMor(K, L, L())
      end
      p = next_prime(p)
    end
  end
  return _issubfield(K, L)
end

################################################################################
#
#  Isomorphism
#
################################################################################

doc"""
***
    isisomorphic(K::AnticNumberField, L::AnticNumberField) -> Bool, NfToNfMor

Returns "true" and an isomorphism from $K$ to $L$ if $K$ and $L$ are isomorphic.
Otherwise the function returns "false" and a morphism mapping everything to 0.
"""
function isisomorphic(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  g = L.pol
  if degree(f) != degree(g)
    return false, NfToNfMor(K, L, L())
  end
  if signature(K) != signature(L)
    return false, NfToNfMor(K, L, L())
  end
  try
    OK = _get_maximal_order_of_nf(K)
    OL = _get_maximal_order_of_nf(L)
    if discriminant(OK) != discriminant(OL)
      return false, NfToNfMor(K, L, L())
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    t = discriminant(f)//discriminant(g)
    if !issquare(numerator(t)) || !issquare(denominator(t))
      return false, NfToNfMor(K, L, L())
    end
  end
  return _issubfield(K, L)
end

################################################################################
#
#  Compositum
#
################################################################################

doc"""
   compositum(K::AnticNumberField, L::AnticNumberField) -> AnticNumberField, Map, Map
> Assuming $L$ is normal (which is not checked), compute the compositum $C$ of the
> 2 fields together with the embedding of $K \to C$ and $L \to C$.
"""
function compositum(K::AnticNumberField, L::AnticNumberField)
  lf = factor(K.pol, L)
  d = degree(first(lf.fac)[1])
  if any(x->degree(x) != d, keys(lf.fac))
    error("2nd field cannot be normal")
  end
  KK = NumberField(first(lf.fac)[1])[1]
  Ka, m1, m2 = absolute_field(KK)
  return Ka, hom(K, Ka, m1(gen(KK))), m2
end

################################################################################
#
#  Serialization
#
################################################################################

# This function can be improved by directly accessing the numerator
# of the fmpq_poly representing the nf_elem
doc"""
***
    write(io::IO, A::Array{nf_elem, 1}) -> Void

> Writes the elements of `A` to `io`. The first line are the coefficients of
> the defining polynomial of the ambient number field. The following lines
> contain the coefficients of the elements of `A` with respect to the power
> basis of the ambient number field.
"""
function write(io::IO, A::Array{nf_elem, 1})
  if length(A) == 0
    return
  else
    # print some useful(?) information
    print(io, "# File created by Hecke $VERSION_NUMBER, $(Base.Dates.now()), by function 'write'\n")
    K = parent(A[1])
    polring = parent(K.pol)

    # print the defining polynomial
    g = K.pol
    d = denominator(g)

    for j in 0:degree(g)
      print(io, coeff(g, j)*d)
      print(io, " ")
    end
    print(io, d)
    print(io, "\n")

    # print the elements
    for i in 1:length(A)

      f = polring(A[i])
      d = denominator(f)

      for j in 0:degree(K)-1
        print(io, coeff(f, j)*d)
        print(io, " ")
      end

      print(io, d)

      print(io, "\n")
    end
  end
end

doc"""
***
    write(file::String, A::Array{nf_elem, 1}, flag::ASCIString = "w") -> Void

> Writes the elements of `A` to the file `file`. The first line are the coefficients of
> the defining polynomial of the ambient number field. The following lines
> contain the coefficients of the elements of `A` with respect to the power
> basis of the ambient number field.
>
> Unless otherwise specified by the parameter `flag`, the content of `file` will be
> overwritten.
"""
function write(file::String, A::Array{nf_elem, 1}, flag::String = "w")
  f = open(file, flag)
  write(f, A)
  close(f)
end

# This function has a bad memory footprint
doc"""
***
    read(io::IO, K::AnticNumberField, ::Type{nf_elem}) -> Array{nf_elem, 1}

> Given a file with content adhering the format of the `write` procedure,
> this functions returns the corresponding object of type `Array{nf_elem, 1}` such that
> all elements have parent $K$.

**Example**

    julia> Qx, x = FlintQQ["x"]
    julia> K, a = NumberField(x^3 + 2, "a")
    julia> write("interesting_elements", [1, a, a^2])
    julia> A = read("interesting_elements", K, Hecke.nf_elem)
"""
function read(io::IO, K::AnticNumberField, ::Type{Hecke.nf_elem})
  Qx = parent(K.pol)

  A = Array{nf_elem, 1}()

  i = 1

  for ln in eachline(io)
    if ln[1] == '#'
      continue
    elseif i == 1
      # the first line read should contain the number field and will be ignored
      i = i + 1
    else
      coe = map(Hecke.fmpz, split(ln, " "))
      t = fmpz_poly(Array(slice(coe, 1:(length(coe) - 1))))
      t = Qx(t)
      t = divexact(t, coe[end])
      push!(A, K(t))
      i = i + 1
    end
  end

  return A
end

doc"""
***
    read(file::String, K::AnticNumberField, ::Type{nf_elem}) -> Array{nf_elem, 1}

> Given a file with content adhering the format of the `write` procedure,
> this functions returns the corresponding object of type `Array{nf_elem, 1}` such that
> all elements have parent $K$.

**Example**

    julia> Qx, x = FlintQQ["x"]
    julia> K, a = NumberField(x^3 + 2, "a")
    julia> write("interesting_elements", [1, a, a^2])
    julia> A = read("interesting_elements", K, Hecke.nf_elem)
"""
function read(file::String, K::AnticNumberField, ::Type{Hecke.nf_elem})
  f = open(file, "r")
  A = read(f, K, Hecke.nf_elem)
  close(f)
  return A
end


text:
export StoppingCriterion
include("problem.jl")
include("options.jl")

include("debugOptions.jl")
include("recordOptions.jl")

include("gradientPlan.jl")
include("proximalPlan.jl")
include("subGradientPlan.jl")

include("stoppingCriterion.jl")
include("stepsize.jl")

ext: jl


text:
using KernelAbstractions
using KernelAbstractions.Extras
using CUDA
using CUDAKernels
using Random123
using WormlikeChain
using StaticArrays

const N= 256

function getrow(A,row)
    ntuple(i -> A[row,i], Val(3))
end

@kernel function foo(A,B)
    I = @index(Global)
    myA= getrow(A,I)
    myB= getrow(B,I)
    myA= myA .+ myB
    @inbounds A[I,1] +=  myA[1]

end

fookh = foo(CPU(), 8)
fookd = foo(CUDADevice(), N)

A_h = ones(Float32,N,3)
A_d = CUDA.ones(Float32,N,3)
B_h = ones(Float32,N,3)
B_d = CUDA.ones(Float32,N,3)

#@time wait(fookd(A_d, B_d, ndrange=N))
#@time wait(fookh(A_h, B_h, ndrange=N))
show(A_d)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Spot_julia_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Spot_julia")
JLLWrappers.@generate_main_file("Spot_julia", UUID("b9ff107b-b78e-5f57-8772-57a339630ba9"))
end  # module Spot_julia_jll


text:
"Adapted From: https://machinelearningmastery.com/differential-evolution-from-scratch-in-python/"
#Written By Johnathan Bizzano

export diffevorandomsinglebin

#Mutate the genes
function mutate(a, b, c, scalefactor)
    return a + scalefactor * (b - c)
end

#Perform gene binomial crossover & clip the mutated genes
function clipNCrossover(mutated, lower_bound, upper_bound, target, cr)
    if rand() >= cr
        return target
    end
    if mutated < lower_bound
        return lower_bound
    elseif mutated > upper_bound
        return upper_bound
    end
    return mutated
end

function choose_individual(population, excluded_individuals...)
    choice = rand(1:length(population))
    while choice in excluded_individuals
        choice = rand(1:length(population))
    end
    choice
end

init_rand(lower_bound, upper_bound, rnginterval) = rand(lower_bound:rnginterval:upper_bound)

"Differential Evolution of DeRand1Bin algorithm. ScaleFactor ∈ [0, 2]. CrossoverRate ∈ [0, 1]"
function diffevorandomsinglebin(device::Device, error_func::Problem, initial_values::AbstractArray{T},
            lower_bound::T, upper_bound::T, population_size::Int, iterations::Int;
            scalefactor = .8, crossoverrate = .7, rnginterval = 1E-8, IsVerbose=true) where T

    (scalefactor > 2 || scalefactor < 0) && error("Scale Factor ∈ [0, 2]")
    (crossoverrate > 1 || crossoverrate < 0) && error("Cross Over Rate ∈ [0, 1]")
    (population_size < 4) && error("Population Size Must Be > 3")

    dims = length(initial_values)
    mutated = alloc(device, T, dims)

    #Create Population with Initial Values as first entry
    population = alloc(device, population_size, dims)
    copy!(view(population, 1, 1:dims), initial_values)

    #Initialize Random
    broadcast!(init_rand, view(population, 2:population_size, dims), lower_bound, upper_bound, rnginterval)

    #Init Eval
    errors = [error_func(ind) for ind in population]

    #Find the Best Init Performer
    best_idx = argmin(errors)
    best_error = errors[best_idx]
    prev_error = best_error

    for i in 1:iterations
        #Iterate all current Solutions
        for j in 1:population_size
            selected = view(population, j, 1:dims)
            #Choose Three Individuals in Population (Not Current)
            idx1 = choose_individual(population, j)
            idx2 = choose_individual(population, j, idx1)
            idx3 = choose_individual(population, j, idx1, idx2)

            broadcast!(mutate, mutated, view(population, idx1, 1:dims), view(population, idx2, 1:dims), view(population, idx3, 1:dims), scalefactor)
            broadcast!(clipNCrossover, mutated, lower_bound, upper_bound, selected, crossoverrate)

            #Check to see if crossover was better then current individual
            crossover_error = error_func(mutated)
            if crossover_error < error_func(selected)
                copy!(selected, mutated)
                errors[j] = crossover_error
            end
        end

        #Check Best Performer in this iteratation
        idx = argmin(errors)
        if errors[idx] < prev_error
            prev_error = errors[idx]
            best_error = prev_error
            best_idx = idx
        end

        (IsVerbose) && println("Iter:$i. Best Error:$best_error")
    end

    return (population[best_idx], best_error)
end


text:
ENV["GKS_ENCODING"] = "utf-8" # Allows the use of unicode characters in Plots.jl
using Plots
using StatsPlots
using Turing
using Bijectors
using Random
using DynamicPPL: getlogp, settrans!, getval, reconstruct, vectorize, setval!

# Set a seed.
Random.seed!(0)

# Define a strange model.
@model function gdemo(x)
    s² ~ InverseGamma(2, 3)
    m ~ Normal(0, sqrt(s²))
    bumps = sin(m) + cos(m)
    m = m + 5*bumps
    for i in eachindex(x)
      x[i] ~ Normal(m, sqrt(s²))
    end
    return s², m
end

# Define our data points.
x = [1.5, 2.0, 13.0, 2.1, 0.0]

# Set up the model call, sample from the prior.
model = gdemo(x)
vi = Turing.VarInfo(model)

# Convert the variance parameter to the real line before sampling.
# Note: We only have to do this here because we are being very hands-on.
# Turing will handle all of this for you during normal sampling.
dist = InverseGamma(2,3)
svn = vi.metadata.s.vns[1]
mvn = vi.metadata.m.vns[1]
setval!(vi, vectorize(dist, Bijectors.link(dist, reconstruct(dist, getval(vi, svn)))), svn)
settrans!(vi, true, svn)

# Evaluate surface at coordinates.
function evaluate(m1, m2)
    spl = Turing.SampleFromPrior()
    vi[svn] = [m1]
    vi[mvn] = [m2]
    model(vi, spl)
    getlogp(vi)
end

function plot_sampler(chain; label="")
    # Extract values from chain.
    val = get(chain, [:s, :m, :lp])
    ss = link.(Ref(InverseGamma(2, 3)), val.s)
    ms = val.m
    lps = val.lp

    # How many surface points to sample.
    granularity = 100

    # Range start/stop points.
    spread = 0.5
    σ_start = minimum(ss) - spread * std(ss); σ_stop = maximum(ss) + spread * std(ss);
    μ_start = minimum(ms) - spread * std(ms); μ_stop = maximum(ms) + spread * std(ms);
    σ_rng = collect(range(σ_start, stop=σ_stop, length=granularity))
    μ_rng = collect(range(μ_start, stop=μ_stop, length=granularity))

    # Make surface plot.
    p = surface(σ_rng, μ_rng, evaluate,
          camera=(30, 65),
        #   ticks=nothing,
          colorbar=false,
          color=:inferno,
          title=label)

    line_range = 1:length(ms)

    scatter3d!(ss[line_range], ms[line_range], lps[line_range],
        mc =:viridis, marker_z=collect(line_range), msw=0,
        legend=false, colorbar=false, alpha=0.5,
        xlabel="σ", ylabel="μ", zlabel="Log probability",
        title=label)

    return p
end;

samplers = [
    (Gibbs(HMC(0.01, 5, :s), PG(20, :m)), "Gibbs{HMC, PG}"),
    (HMC(0.01, 10), "HMC"),
    (HMCDA(200, 0.65, 0.3), "HMCDA"),
    (MH(), "MH()"),
    (NUTS(0.65), "NUTS(0.65)"),
    (NUTS(0.95), "NUTS(0.95)"),
    (NUTS(0.2), "NUTS(0.2)"),
    (PG(20), "PG(20)"),
    (PG(50), "PG(50)")]

for (i, (spl, spl_name)) in enumerate(samplers)
    c = sample(model, spl, 1000)
    p = plot_sampler(c, label="$spl_name")
    savefig(joinpath(@__DIR__, "sampler-figs", "samplers-$i.svg"))
end


text:
using Compat

@windows_only begin
    using WinRPM
    WinRPM.install("libfreetype6", yes=true)
end

const freetype = Libdl.find_library(
    ["libfreetype", "libfreetype-6", "libfreetype.6", "freetype"],
    [
        "/usr/lib/i386-linux-gnu", "/usr/lib/x86_64-linux-gnu",
        Pkg.dir("WinRPM", "deps", "usr", "$(Sys.ARCH)-w64-mingw32", "sys-root", "mingw", "bin"),
        "/Developer/SDKs/MacOSX10.7.sdk/usr/X11/include/freetype2", "/usr/X11/lib/"
    ]
)

if isempty(freetype)
    error("FreeType library could not be found. Please file an issue here: https://github.com/jhasse/FreeType.jl/issues/new")
end

# Save lib and version into seperate file to statically include the library with precompile
f = open("deps.jl", "w")
write(f, """
const freetype = \"$(escape_string(freetype))\"
""")
close(f)


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

@with_kw mutable struct ParametricGomoryData{T <: AbstractFloat}
    relaxation::Union{Nothing, SecondStageLP{T}} = nothing
    cglp::Union{Nothing, JuMP.Model} = nothing
    unrelax::Function = () -> nothing
end

@with_kw mutable struct ParametricGomoryParameters
    maximum_iterations::Int = 1
    optimizer = nothing
    strategy::ParametricGomoryStrategy = Gomory()
end
"""
    ParametricGomoryMaster

Master functor object for using weak optimality cuts in an integer L-shaped algorithm. Requires all first-stage decisions to be binary. Create by supplying a [`UseParametricGomory`](@ref) object through `integer_strategy` in `LShaped.Optimizer` or set the [`IntegerStrategy`](@ref) attribute.

"""
struct ParametricGomoryMaster <: AbstractIntegerAlgorithm
    integer_variables::Vector{MOI.VariableIndex}
end

function initialize_integer_algorithm!(::ParametricGomoryMaster, first_stage::JuMP.Model)
    # Sanity check
    if all(is_binary, all_decision_variables(first_stage, 1))
       @warn "All first-stage decisions are binary. Convexification or CombinatorialCuts are more approriate integer strategies"
    end
    # Gather integer variables
    append!(master.integer_variables, gather_integer_variables(subproblem))
    # Sanity check
    if isempty(master.integer_variables)
        @warn "No integer variables in subproblem $(subproblem.id). Integer strategy is superfluous."
    end
    return nothing
end

function handle_integrality!(lshaped::AbstractLShaped, ::ParametricGomoryMaster)
    @unpack τ = lshaped.parameters
    if gap(lshaped) <= τ
        set_metadata!(lshaped.execution.metadata, :converged, true)
    end
    return nothing
end
"""
    ParametricGomoryWorker

Worker functor object for using weak optimality cuts in an integer L-shaped algorithm. Create by supplying a [`UseParametricGomory`](@ref) object through `integer_strategy` in `LShaped.Optimizer` or set the [`IntegerStrategy`](@ref) attribute.

"""
struct ParametricGomoryWorker{T <: AbstractFloat} <: AbstractIntegerAlgorithm
    data::ParametricGomoryData{T}
    parameters::ParametricGomoryParameters
    integer_variables::Vector{MOI.VariableIndex}
    cutting_planes::Vector{MOI.ConstraintIndex}

    function ParametricGomoryWorker(::Type{T}; kw...) where T <: AbstractFloat
        T_ = promote_type(T, Float32)
        worker = new{T_}(ParametricGomoryData{T_}(),
                         ParametricGomoryParameters(; kw...),
                         Vector{MOI.VariableIndex}(),
                         Vector{MOI.ConstraintIndex}())
        return worker
    end
end

function initialize_integer_algorithm!(worker::ParametricGomoryWorker, subproblem::SubProblem)
    # Sanity check
    if any(is_integer, all_decision_variables(subproblem.model, StochasticPrograms.stage(subproblem.model))) ||
       any(is_integer, all_auxiliary_variables(subproblem.model))
        if !(worker.parameters.strategy isa Gomory)
            if subproblem.id == 1
                @warn "Switching convexification strategy to Gomory to handle integer variables."
            end
            worker.parameters.strategy = Gomory()
        end
    end
    # Gather integer variables
    append!(worker.integer_variables, gather_integer_variables(subproblem))
    # Sanity check
    if isempty(worker.integer_variables)
        @warn "No integer variables in subproblem $(subproblem.id). Integer strategy is superfluous."
    end
    # Relax integer restrictions
    worker.data.unrelax = relax_decision_integrality(subproblem.model)
    # Extract model information
    standard_form = worker.parameters.strategy isa Gomory
    worker.data.relaxation = SecondStageLP(subproblem; standard_form = standard_form)
    return nothing
end

function restore!(model::MOI.ModelLike, worker::ParametricGomoryWorker)
    # Delete any added cutting planes
    if !isempty(worker.cutting_planes)
        MOI.delete(model, worker.cutting_planes)
    end
    # Restore integral restrictions
    worker.data.unrelax()
    return nothing
end

function integer_variables(worker::ParametricGomoryWorker)
    return worker.integer_variables
end

function solve_subproblem(subproblem::SubProblem{T},
                          metadata,
                          ::NoFeasibilityAlgorithm,
                          worker::ParametricGomoryWorker,
                          x::AbstractVector) where T <: AbstractFloat
    # Initial solve
    cut = solve_subproblem(subproblem, x)
    # Check if integer restrictions are satisfied
    # and collect any fractional variables
    fractional = MOI.VariableIndex[]
    set_metadata!(metadata,
                  subproblem.id,
                  :integral_solution,
                  check_integrality_restrictions(subproblem, fractional))
    # Check if solution satisfies integer requirements
    if isempty(fractional)
        # No cutting planes needed, just return cut
        return cut
    end
    # Check if L-shaped has converged
    if has_metadata(metadata, :converged)
        # Solve subproblem with integer restrictions
        return solve_unrelaxed(subproblem, metadata, worker, x)
    end
    # Otherwise, add cutting planes until solution is integer
    # or user-limit is reached
    num_iterations = 1
    while true
        inds = MOI.get(subproblem.optimizer, MOI.ListOfVariableIndices())
        ȳ = MOI.get.(subproblem.optimizer, MOI.VariablePrimal(), inds[length(x)+1:end])
        try
            # Generate and add cutting plane
            πx, πy, g = cutting_plane(worker,
                                      worker.parameters.strategy,
                                      subproblem,
                                      x,
                                      ȳ,
                                      fractional)
            f, set, coeffs = moi_constraint(subproblem, πx, πy, g)
            # Add cutting plane to subproblem
            ci = MOI.add_constraint(subproblem.optimizer, f, set)
            push!(worker.cutting_planes, ci)
            # Update subproblem
            if !isempty(coeffs)
                push!(subproblem.linking_constraints, ci)
                push!(subproblem.masterterms, coeffs)
            end
        catch
            @warn "Cutting plane could not be generated in subproblem $(subproblem.id)"
            # If cutting plane could not be generated,
            # return latest working cut and try again next iteration
            return cut
        end
        # Re-solve subproblem
        cut = solve_subproblem(subproblem, x)
        fractional = MOI.VariableIndex[]
        set_metadata!(metadata,
                      subproblem.id,
                      :integral_solution,
                      check_integrality_restrictions(subproblem, fractional))
        if isempty(fractional) || num_iterations >= worker.parameters.maximum_iterations
            # Return cut
            return cut
        end
        # Increment iteration count
        num_iterations += 1
    end
end

function solve_subproblem(subproblem::SubProblem,
                          metadata,
                          feasibility_algorithm::FeasibilityCutsWorker,
                          worker::ParametricGomoryWorker,
                          x::AbstractVector)
    # Prepare auxiliary problem
    model = subproblem.optimizer
    if !prepared(feasibility_algorithm)
        prepare!(model, feasibility_algorithm)
    end
    # Relax integrality
    unrelax = relax_decision_integrality(subproblem.model)
    # Optimize auxiliary problem
    MOI.optimize!(model)
    # Sanity check that aux problem could be solved
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if !(status ∈ AcceptableTermination)
        error("Subproblem $(subproblem.id) was not solved properly during feasibility check, returned status code: $status")
    end
    # check objective
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    w = correction * MOI.get(model, MOI.ObjectiveValue())
    # Ensure correction is available in master
    set_metadata!(metadata, subproblem.id, :correction, correction)
    # Check feasibility
    if w > sqrt(eps())
        # Subproblem is infeasible, create feasibility cut
        cut = FeasibilityCut(subproblem, x)
        # Unrelax integrality restrictions again
        unrelax()
        return cut
    end
    # Unrelax integrality restrictions again
    unrelax()
    # Restore subproblem and solve as usual
    restore_subproblem!(subproblem)
    return solve_subproblem(subproblem, metadata, NoFeasibilityAlgorithm(), worker, x)
end

function solve_unrelaxed(subproblem::SubProblem, metadata, worker::ParametricGomoryWorker, x::AbstractVector)
    # Delete any added cutting planes
    if !isempty(worker.cutting_planes)
        MOI.delete(subproblem.optimizer, worker.cutting_planes)
    end
    empty!(worker.cutting_planes)
    # Restore integral restrictions
    worker.data.unrelax()
    worker.data.unrelax = () -> nothing
    # Solve subproblem with integer restrictions
    MOI.optimize!(subproblem.optimizer)
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if status ∈ AcceptableTermination
        # Integer restrictions are satisfied if optimal
        set_metadata!(metadata,
                      subproblem.id,
                      :integral_solution,
                      true)
        # Get sense
        sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
        correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
        # Create sense-corrected optimality cut
        Q = correction * MOI.get(subproblem.optimizer, MOI.ObjectiveValue())
        return OptimalityCut(spzeros(length(x)), Q, subproblem.id)
    elseif status == MOI.INFEASIBLE
        return Infeasible(subproblem)
    elseif status == MOI.DUAL_INFEASIBLE
        return Unbounded(subproblem)
    else
        error("Subproblem $(subproblem.id) was not solved properly, returned status code: $status")
    end
end

function moi_constraint(subproblem::SubProblem{T}, πx::AbstractVector, πy::AbstractVector, g::AbstractFloat) where T <: AbstractFloat
    # Create cutting plane
    f = zero(AffineDecisionFunction{T})
    coeffs = Vector{Tuple{Int, Int, T}}()
    N = StochasticPrograms.stage(subproblem.model)
    all_known = mapreduce(vcat, 1:N-1) do s
        index.(all_known_decision_variables(subproblem.model, s))
    end
    all_decisions = index.(all_decision_variables(subproblem.model, N))
    ix = 0
    iy = 0
    for vi in MOI.get(subproblem.optimizer, MOI.ListOfVariableIndices())
        if vi in all_known
            ix += 1
            coeff = T(πx[ix])
            if abs(coeff) >= sqrt(eps())
                push!(f.known_part.terms, MOI.ScalarAffineTerm(coeff, vi))
                push!(coeffs, (1, vi.value, coeff))
            else
                πx[ix] = zero(T)
            end
        elseif vi in all_decisions
            iy += 1
            coeff = T(πy[iy])
            if abs(coeff) >= sqrt(eps())
                push!(f.decision_part.terms, MOI.ScalarAffineTerm(coeff, vi))
            else
                πy[iy] = zero(T)
            end
        else
            iy += 1
            coeff = T(πy[iy])
            if abs(coeff) >= sqrt(eps())
                push!(f.variable_part.terms, MOI.ScalarAffineTerm(coeff, vi))
            else
                πy[iy] = zero(T)
            end
        end
    end
    set = MOI.GreaterThan{Float64}(g)
    return f, set, coeffs
end

function cutting_plane(worker::ParametricGomoryWorker,
                       ::Gomory,
                       subproblem::SubProblem,
                       x̄::AbstractVector,
                       ȳ::AbstractVector,
                       fractional::Vector{MOI.VariableIndex})
    # Pick smallest fractional index
    fractional_index = fractional[1].value - length(x̄)
    # Get matrices
    T = worker.data.relaxation.T
    W = worker.data.relaxation.W
    h = worker.data.relaxation.h
    ny = length(ȳ)
    # Get current optimal basis
    basis = optimal_basis(subproblem)
    nonbasis = setdiff(1:size(W,2),basis)
    # Compute row of basis inverse corresponding
    # to fractional index
    i = something(findfirst(j -> j == fractional_index, basis), 0)
    if i == 0
        error("Fractional variable not in basis.")
    end
    ei = zeros(size(W, 1))
    ei[i] = 1.0
    Bᵢ⁻¹ = W[:,basis]' \ ei
    # Get source rows
    nonbasic = setdiff(1:ny, basis)
    nonbasic_slack = setdiff(ny+1:size(W,2), basis)
    w = (Bᵢ⁻¹'*W[:,nonbasis])
    γ = (Bᵢ⁻¹'*T)[:]
    ρ = Bᵢ⁻¹'*h
    # Translate x to the origin
    for (j,x) in enumerate(x̄)
        if x == 1.
            ρ -= γ[j]
            γ[j] *= -1
        end
    end
    # Calculate Gomory coefficients
    ϕ(x) = x - floor(x)
    #πx = ϕ.(γ)
    #w = ϕ.(w)
    #ρ = ϕ(ρ)
    ϕγ = ϕ.(γ)
    ϕw = ϕ.(w)
    ρ  = ϕ(ρ)
    πx = min.(ϕγ, (ρ .* (1 .- ϕγ) ./ (1 .- ρ)))
    w  = min.(ϕw, (ρ .* (1 .- ϕw) ./ (1 .- ρ)))
    # Translate back
    for (j,x) in enumerate(x̄)
        if x == 1.
            ρ -= πx[j]
            πx[j] *= -1
        end
    end
    πy = zero(ȳ)
    πy[nonbasic] .= w[nonbasic]
    # Map slacks to rows
    for (coeff,slack) in zip(w[length(nonbasic)+1:end], nonbasic_slack)
        row = slack - ny
        πx .-= W[row,slack] * coeff * T[row,:]
        πy .-= W[row,slack] * coeff * W[row,1:ny]
        ρ -= W[row,slack] * coeff * h[row]
    end
    # Extend standard form LP with new cuts
    add_row!(worker.data.relaxation, πx, πy, ρ; standard_form = true)
    # Return cutting plane
    return πx, πy, ρ
end

function cutting_plane(worker::ParametricGomoryWorker,
                       ::LiftAndProject,
                       subproblem::SubProblem{T},
                       x̄::AbstractVector,
                       ȳ::AbstractVector,
                       fractional::Vector{MOI.VariableIndex}) where T <: AbstractFloat
    # Check optimizer
    if worker.parameters.optimizer === nothing
        error("Cannot generate CGLP without optimizer.")
    end
    # Largest fractional index
    fractional_index = fractional[end].value - length(x̄)
    model = subproblem.model
    nx = length(x̄)
    ny = length(ȳ)
    nw = nx + ny - 1
    # Prepare data
    Gᵢ = Vector{Int}()
    Gⱼ = Vector{Int}()
    Gᵥ = Vector{T}()
    Hᵢ = Vector{Int}()
    Hⱼ = Vector{Int}()
    Hᵥ = Vector{T}()
    Fᵢ = Vector{Int}()
    Fⱼ = Vector{Int}()
    Fᵥ = Vector{T}()
    lhs = Vector{T}()
    rhs = Vector{T}()
    i = 0
    # Loop over constraints
    for (F,S) in list_of_constraint_types(model)
        for cref in all_constraints(model, F, S)
            if S <: StochasticPrograms.SingleDecisionSet ||
                S <: StochasticPrograms.MultipleDecisionSet ||
                S <: StochasticPrograms.SingleKnownSet ||
                S <: StochasticPrograms.MultipleKnownSet
                continue
            elseif F <: VariableRef || F <: DecisionRef
                continue
            elseif F <: AffExpr
                i += 1
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                for term in f.terms
                    idx = term.variable_index.value - nx
                    if idx == fractional_index
                        # yⱼ
                        push!(Hᵢ, i)
                        push!(Hᵥ, term.coefficient)
                        push!(Hⱼ, idx)
                    else
                        # yⱼ
                        push!(Fᵢ, i)
                        push!(Fᵥ, term.coefficient)
                        push!(Fⱼ, nx + ny + nx + idx)
                        # 1-yⱼ
                        push!(Hᵢ, i+1)
                        push!(Hᵥ, term.coefficient)
                        push!(Hⱼ, idx)
                        push!(Fᵢ, i+1)
                        push!(Fᵥ, -term.coefficient)
                        push!(Fⱼ, nx + ny + nx + idx)
                    end
                end
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                if set isa MOI.LessThan
                    # yⱼ
                    push!(lhs, -Inf)
                    push!(rhs, zero(T))
                    push!(Hᵢ, i)
                    push!(Hᵥ, -set.upper)
                    push!(Hⱼ, fractional_index)
                    # 1-yⱼ
                    push!(lhs, -Inf)
                    push!(rhs, set.upper)
                    push!(Hᵢ, i+1)
                    push!(Hᵥ, set.upper)
                    push!(Hⱼ, fractional_index)
                elseif set isa MOI.GreaterThan
                    # yⱼ
                    push!(lhs, zero(T))
                    push!(rhs, Inf)
                    push!(Hᵢ, i)
                    push!(Hᵥ, -set.lower)
                    push!(Hⱼ, fractional_index)
                    # 1-yⱼ
                    push!(lhs, set.lower)
                    push!(rhs, Inf)
                    push!(Hᵢ, i+1)
                    push!(Hᵥ, set.lower)
                    push!(Hⱼ, fractional_index)
                elseif set isa MOI.EqualTo
                    # yⱼ
                    push!(lhs, zero(T))
                    push!(rhs, zero(T))
                    push!(Hᵢ, i)
                    push!(Hᵥ, -set.value)
                    push!(Hⱼ, fractional_index)
                    # 1-yⱼ
                    push!(lhs, set.value)
                    push!(rhs, set.value)
                    push!(Hᵢ, i+1)
                    push!(Hᵥ, set.value)
                    push!(Hⱼ, fractional_index)
                end
                i += 1
            elseif F <: DecisionAffExpr
                i += 1
                f = MOI.get(model, MOI.ConstraintFunction(), cref)
                for term in f.variable_part.terms
                    idx = term.variable_index.value - nx
                    if idx == fractional_index
                        # yⱼ
                        push!(Hᵢ, i)
                        push!(Hᵥ, term.coefficient)
                        push!(Hⱼ, idx)
                    else
                        # yⱼ
                        push!(Fᵢ, i)
                        push!(Fᵥ, term.coefficient)
                        push!(Fⱼ, nx + idx - 1)
                        # 1-yⱼ
                        push!(Hᵢ, i+1)
                        push!(Hᵥ, term.coefficient)
                        push!(Hⱼ, idx)
                        push!(Fᵢ, i+1)
                        push!(Fᵥ, -term.coefficient)
                        push!(Fⱼ, nx + idx - 1)
                    end
                end
                for term in f.decision_part.terms
                    idx = term.variable_index.value - nx
                    if idx == fractional_index
                        # yⱼ
                        push!(Hᵢ, i)
                        push!(Hᵥ, term.coefficient)
                        push!(Hⱼ, idx)
                    else
                        # yⱼ
                        push!(Fᵢ, i)
                        push!(Fᵥ, term.coefficient)
                        push!(Fⱼ, nx + idx - 1)
                        # 1-yⱼ
                        push!(Hᵢ, i+1)
                        push!(Hᵥ, term.coefficient)
                        push!(Hⱼ, idx)
                        push!(Fᵢ, i+1)
                        push!(Fᵥ, -term.coefficient)
                        push!(Fⱼ, nx + idx - 1)
                    end
                end
                for term in f.known_part.terms
                    idx = term.variable_index.value
                    # yⱼ
                    push!(Fᵢ, i)
                    push!(Fᵥ, term.coefficient)
                    push!(Fⱼ, idx)
                    # 1-yⱼ
                    push!(Gᵢ, i+1)
                    push!(Gᵥ, term.coefficient)
                    push!(Gⱼ, idx)
                    push!(Fᵢ, i+1)
                    push!(Fᵥ, -term.coefficient)
                    push!(Fⱼ, idx)
                end
                set = MOI.get(model, MOI.ConstraintSet(), cref)
                if set isa MOI.LessThan
                    # yⱼ
                    push!(lhs, -Inf)
                    push!(rhs, zero(T))
                    push!(Hᵢ, i)
                    push!(Hᵥ, -set.upper)
                    push!(Hⱼ, fractional_index)
                    # 1-yⱼ
                    push!(lhs, -Inf)
                    push!(rhs, set.upper)
                    push!(Hᵢ, i+1)
                    push!(Hᵥ, set.upper)
                    push!(Hⱼ, fractional_index)
                elseif set isa MOI.GreaterThan
                    # yⱼ
                    push!(lhs, zero(T))
                    push!(rhs, Inf)
                    push!(Hᵢ, i)
                    push!(Hᵥ, -set.lower)
                    push!(Hⱼ, fractional_index)
                    # 1-yⱼ
                    push!(lhs, set.lower)
                    push!(rhs, Inf)
                    push!(Hᵢ, i+1)
                    push!(Hᵥ, set.lower)
                    push!(Hⱼ, fractional_index)
                elseif set isa MOI.EqualTo
                    # yⱼ
                    push!(lhs, zero(T))
                    push!(rhs, zero(T))
                    push!(Hᵢ, i)
                    push!(Hᵥ, -set.value)
                    push!(Hⱼ, fractional_index)
                    # 1-yⱼ
                    push!(lhs, set.value)
                    push!(rhs, set.value)
                    push!(Hᵢ, i+1)
                    push!(Hᵥ, set.value)
                    push!(Hⱼ, fractional_index)
                end
                i += 1
            else
                error("Cannot extract linear representation from $F-$S constraint.")
            end
        end
    end
    # x ≥ 0
    for j in 1:length(x̄)
        i += 1
        # yⱼ
        push!(Fᵢ, i)
        push!(Fᵥ, 1.)
        push!(Fⱼ, j)
        push!(lhs, zero(T))
        push!(rhs, Inf)
        # 1-yⱼ
        push!(Gᵢ, i+1)
        push!(Gᵥ, 1.)
        push!(Gⱼ, j)
        push!(Fᵢ, i+1)
        push!(Fᵥ, -1.)
        push!(Fⱼ, j)
        push!(lhs, zero(T))
        push!(rhs, Inf)
        i += 1
    end
    # x ≤ 1
    for j in 1:length(x̄)
        i += 1
        # yⱼ
        push!(Fᵢ, i)
        push!(Fᵥ, 1.)
        push!(Fⱼ, j)
        push!(Hᵢ, i)
        push!(Hᵥ, -1.)
        push!(Hⱼ, fractional_index)
        push!(lhs, -Inf)
        push!(rhs, zero(T))
        # 1-yⱼ
        push!(Gᵢ, i+1)
        push!(Gᵥ, 1.)
        push!(Gⱼ, j)
        push!(Fᵢ, i+1)
        push!(Fᵥ, -1.)
        push!(Fⱼ, j)
        push!(Hᵢ, i+1)
        push!(Hᵥ, 1.)
        push!(Hⱼ, fractional_index)
        push!(lhs, -Inf)
        push!(rhs, 1.)
        i += 1
    end
    # y ≥ 0
    for j in 1:length(ȳ)
        i += 1
        if j == fractional_index
            # yⱼ
            push!(Hᵢ, i)
            push!(Hᵥ, 1.)
            push!(Hⱼ, j)
            push!(lhs, zero(T))
            push!(rhs, Inf)
        else
            # yⱼ
            push!(Fᵢ, i)
            push!(Fᵥ, 1.)
            push!(Fⱼ, nx + j - 1)
            push!(lhs, zero(T))
            push!(rhs, Inf)
            # 1-yⱼ
            push!(Hᵢ, i+1)
            push!(Hᵥ, 1.)
            push!(Hⱼ, j)
            push!(Fᵢ, i+1)
            push!(Fᵥ, -1.)
            push!(Fⱼ, nx + j - 1)
            push!(lhs, zero(T))
            push!(rhs, Inf)
            i += 1
        end
    end
    # y ≤ 1
    for j in 1:length(ȳ)
        i += 1
        if j == fractional_index
            # 1-yⱼ
            push!(Hᵢ, i)
            push!(Hᵥ, 1.)
            push!(Hⱼ, j)
            push!(lhs, -Inf)
            push!(rhs, 1.)
        else
            # yⱼ
            push!(Fᵢ, i)
            push!(Fᵥ, 1.)
            push!(Fⱼ, nx + j - 1)
            push!(Hᵢ, i)
            push!(Hᵥ, -1.)
            push!(Hⱼ, fractional_index)
            push!(lhs, -Inf)
            push!(rhs, zero(T))
            # 1-yⱼ
            push!(Hᵢ, i+1)
            push!(Hᵥ, 1.)
            push!(Hⱼ, j)
            push!(Fᵢ, i+1)
            push!(Fᵥ, -1.)
            push!(Fⱼ, nx + j - 1)
            push!(Hᵢ, i+1)
            push!(Hᵥ, 1.)
            push!(Hⱼ, fractional_index)
            push!(lhs, -Inf)
            push!(rhs, 1.)
            i += 1
        end
    end
    G̃ = sparse(Gᵢ, Gⱼ, Gᵥ, i, nx)
    H̃ = sparse(Hᵢ, Hⱼ, Hᵥ, i, ny)
    F̃ = sparse(Fᵢ, Fⱼ, Fᵥ, i, nw)
    G, f = StochasticPrograms.canonical(G̃, lhs, rhs; direction = :geq)
    H, f = StochasticPrograms.canonical(H̃, lhs, rhs; direction = :geq)
    F, f = StochasticPrograms.canonical(F̃, lhs, rhs; direction = :geq)
    # Generate CGLP
    CGLP = Model(worker.parameters.optimizer)
    n = size(F, 1)
    @variable(CGLP, π[1:n] >= 0)
    @objective(CGLP, Min, π⋅(H*ȳ) - π⋅(f - G*x̄))
    @constraint(CGLP, F'*π .== 0)
    @constraint(CGLP, sum(π) == 1)
    # Solve CGLP
    optimize!(CGLP)
    if termination_status(CGLP) == MOI.OPTIMAL
        π = value.(π)
        πx = G'*π
        πy = H'*π
        g = π⋅f
        # Extend standard form LP with new cuts
        add_row!(worker.data.relaxation, πx, πy, g)
        # Return cutting plane
        return πx, πy, g
    else
        error("CGLP could not be solved.")
    end
end

function cutting_plane(worker::ParametricGomoryWorker,
                       ::CuttingPlaneTree,
                       subproblem::SubProblem,
                       x̄::AbstractVector,
                       ȳ::AbstractVector,
                       fractional::Vector{MOI.VariableIndex})
    # Check optimizer
    if worker.parameters.optimizer === nothing
        error("Cannot generate CGLP without optimizer.")
    end
    # Largest fractional index
    fractional_index = fractional[end].value - length(x̄)
    # Dimensions
    lp = worker.data.relaxation
    m  = size(lp.T, 1)
    nx = size(lp.T, 2)
    ny = size(lp.W, 2)
    length(x̄) == nx || error("Given first-stage decision $x̄ does not match second-stage problem dimensions.")
    length(ȳ) == ny || error("Given second-stage decision $ȳ does not match second-stage problem dimensions.")
    ȳ[abs.(ȳ) .<= sqrt(eps())] .= 0.0
    x_lb = fill(0., nx)
    x_ub = fill(1., nx)
    lb_finite = isfinite.(lp.lb)
    ub_finite = isfinite.(lp.ub)
    # Generate CGLP
    CGLP = Model(worker.parameters.optimizer)
    @variable(CGLP, πx[1:nx])
    @variable(CGLP, πy[1:ny])
    @variable(CGLP, tx >= 0)
    @variable(CGLP, ty >= 0)
    @variable(CGLP, λ₂₁[1:m] >= 0)
    @variable(CGLP, λ₂₂[1:m] >= 0)
    @variable(CGLP, ν₁₁[1:nx] >= 0)
    @variable(CGLP, ν₁₂[1:nx] >= 0)
    @variable(CGLP, ν₂₁[1:ny] >= 0)
    @variable(CGLP, ν₂₂[1:ny] >= 0)
    @variable(CGLP, μ₁₁[1:nx] >= 0)
    @variable(CGLP, μ₁₂[1:nx] >= 0)
    @variable(CGLP, μ₂₁[1:ny] >= 0)
    @variable(CGLP, μ₂₂[1:ny] >= 0)
    @objective(CGLP, Min, tx + ty)
    @constraint(CGLP, vcat(tx, πx) in MOI.NormOneCone(nx+1))
    @constraint(CGLP, vcat(ty, πy) in MOI.NormOneCone(ny+1))
    @constraint(CGLP, πx .== lp.T'*λ₂₁ + μ₁₁ - ν₁₁)
    @constraint(CGLP, πx .== lp.T'*λ₂₂ + μ₁₂ - ν₁₂)
    @constraint(CGLP, πy .== lp.W'*λ₂₁ + μ₂₁ - ν₂₁)
    @constraint(CGLP, πy .== lp.W'*λ₂₂ + μ₂₂ - ν₂₂)
    lp.ub[fractional_index] = 0.
    @constraint(CGLP, lp.h⋅λ₂₁ + x_lb⋅μ₁₁ + lp.lb[lb_finite]⋅μ₂₁[lb_finite] - x_ub⋅ν₁₁ - lp.ub[ub_finite]⋅ν₂₁[ub_finite] >= πx⋅x̄ + πy⋅ȳ + 1)
    lp.ub[fractional_index] = 1.
    lp.lb[fractional_index] = 1.
    @constraint(CGLP, lp.h⋅λ₂₂ + x_lb⋅μ₁₂ + lp.lb[lb_finite]⋅μ₂₂[lb_finite] - x_ub⋅ν₁₂ - lp.ub[ub_finite]⋅ν₂₂[ub_finite] >= πx⋅x̄ + πy⋅ȳ + 1)
    lp.lb[fractional_index] = 0.
    # Solve CGLP
    optimize!(CGLP)
    if termination_status(CGLP) == MOI.OPTIMAL
        # Extract cutting plane parameters
        πx = value.(CGLP[:πx])
        πy = value.(CGLP[:πy])
        g  = πx⋅x̄ + πy⋅ȳ + 1.
        # Extend standard form LP with new cuts
        add_row!(worker.data.relaxation, πx, πy, g)
        # Return cutting plane
        return πx, πy, g
    else
        error("CGLP could not be solved to optimality.")
    end
end

# API
# ------------------------------------------------------------
"""
    ParametricGomory

Factory object for using convexification to handle integer recourse. Pass to `integer_strategy` in `LShaped.Optimizer` or set the [`IntegerStrategy`](@ref) attribute.

"""
struct ParametricGomory <: AbstractIntegerStrategy
    parameters::ParametricGomoryParameters
end
ParametricGomory(; kw...) = ParametricGomory(ParametricGomoryParameters(; kw...))

function master(wc::ParametricGomory, ::Type{T}) where T <: AbstractFloat
    return ParametricGomoryMaster()
end

function worker(wc::ParametricGomory, ::Type{T}) where T <: AbstractFloat
    return ParametricGomoryWorker(T; type2dict(wc.parameters)...)
end
function worker_type(::ParametricGomory)
    return ParametricGomoryWorker{Float64}
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

include("graph.jl")

"""
    LazyBridgeOptimizer(model::MOI.ModelLike)

The `LazyBridgeOptimizer` is a bridge optimizer that supports multiple bridges,
and only bridges things which are not supported by the internal model.

Internally, the `LazyBridgeOptimizer` solves a shortest hyper-path problem to
determine which bridges to use.

In general, you should use [`full_bridge_optimizer`](@ref) instead of this
constructor because [`full_bridge_optimizer`](@ref) automatically adds a large
number of supported bridges.

See also: [`add_bridge`](@ref), [`remove_bridge`](@ref), [`has_bridge`](@ref)
and [`full_bridge_optimizer`](@ref).

## Example

```jldoctest; setup=:(using MathOptInterface; const MOI = MathOptInterface)
julia> model = MOI.Bridges.LazyBridgeOptimizer(MOI.Utilities.Model{Float64}())
MOIB.LazyBridgeOptimizer{MOIU.Model{Float64}}
with 0 variable bridges
with 0 constraint bridges
with 0 objective bridges
with inner model MOIU.Model{Float64}

julia> MOI.Bridges.add_bridge(model, MOI.Bridges.Variable.FreeBridge{Float64})

julia> MOI.Bridges.has_bridge(model, MOI.Bridges.Variable.FreeBridge{Float64})
true
```
"""
mutable struct LazyBridgeOptimizer{OT<:MOI.ModelLike} <: AbstractBridgeOptimizer
    # Internal model
    model::OT
    # Bridged variables
    variable_map::Variable.Map
    var_to_name::Dict{MOI.VariableIndex,String}
    name_to_var::Union{Dict{String,MOI.VariableIndex},Nothing}
    # Bridged constraints
    constraint_map::Constraint.Map
    con_to_name::Dict{MOI.ConstraintIndex,String}
    name_to_con::Union{Dict{String,MOI.ConstraintIndex},Nothing}
    # Bridged objective
    objective_map::Objective.Map
    # Bellman-Ford graph
    graph::Graph
    # List of types of available bridges
    variable_bridge_types::Vector{Any}
    variable_node::OrderedDict{Tuple{Type},VariableNode}
    variable_types::Vector{Tuple{Type}}
    # List of types of available bridges
    constraint_bridge_types::Vector{Any}
    constraint_node::OrderedDict{Tuple{Type,Type},ConstraintNode}
    constraint_types::Vector{Tuple{Type,Type}}
    # List of types of available bridges
    objective_bridge_types::Vector{Any}
    objective_node::OrderedDict{Tuple{Type},ObjectiveNode}
    objective_types::Vector{Tuple{Type}}
    # Cache for (F, S) -> BridgeType. Avoids having to look up
    # `concrete_bridge_type` at runtime, which is slow.
    cached_bridge_type::Dict{Any,Type}

    function LazyBridgeOptimizer(model::MOI.ModelLike)
        return new{typeof(model)}(
            model,
            Variable.Map(),
            Dict{MOI.VariableIndex,String}(),
            nothing,
            Constraint.Map(),
            Dict{MOI.ConstraintIndex,String}(),
            nothing,
            Objective.Map(),
            Graph(),
            Any[],
            OrderedDict{Tuple{Type},VariableNode}(),
            Tuple{Type}[],
            Any[],
            OrderedDict{Tuple{Type,Type},ConstraintNode}(),
            Tuple{Type,Type}[],
            Any[],
            OrderedDict{Tuple{Type},ObjectiveNode}(),
            Tuple{Type}[],
            Dict{Any,Type}(),
        )
    end
end

Variable.bridges(bridge::LazyBridgeOptimizer) = bridge.variable_map

Constraint.bridges(bridge::LazyBridgeOptimizer) = bridge.constraint_map

Objective.bridges(b::LazyBridgeOptimizer) = b.objective_map

"""
    _reset_bridge_graph(b::LazyBridgeOptimizer)

Reset the hyper-graph maintained by the [LazyBridgeOptimizer](@ref), but do not
remove the bridges that have been added.

We need to reset the graph whenever we add or delete bridges from `b`, because
it might change the optimal hyper-paths. Recomputing from scratch is a lot of
work, but it will only happen if we interleave calls to [`add_bridge`](@ref),
[`remove_bridge`](@ref), and model builders like [`MOI.add_constraint`](@ref).
"""
function _reset_bridge_graph(b::LazyBridgeOptimizer)
    empty!(b.variable_node)
    empty!(b.variable_types)
    empty!(b.constraint_node)
    empty!(b.constraint_types)
    empty!(b.objective_node)
    empty!(b.objective_types)
    empty!(b.graph)
    empty!(b.cached_bridge_type)
    return
end

"""
    _variable_nodes(b::LazyBridgeOptimizer, ::Type{BT}) where {BT}

Return the list of `VariableNode` that would be added if `BT` is used in `b`.
"""
function _variable_nodes(b::LazyBridgeOptimizer, ::Type{BT}) where {BT}
    return VariableNode[
        node(b, S) for (S,) in added_constrained_variable_types(BT)
    ]
end

"""
    _constraint_nodes(b::LazyBridgeOptimizer, ::Type{BT}) where {BT}

Return the list of `ConstraintNode` that would be added if `BT` is used in `b`.
"""
function _constraint_nodes(b::LazyBridgeOptimizer, ::Type{BT}) where {BT}
    return ConstraintNode[
        node(b, F, S) for (F, S) in added_constraint_types(BT)
    ]
end

"""
    _edge(b::LazyBridgeOptimizer, index::Int, BT::Type{<:AbstractBridge})

Return the `Edge` or `ObjectiveEdge` in the hyper-graph associated with the
bridge `BT`, where `index` is the index of `BT` in the list of bridges.
"""
function _edge(b::LazyBridgeOptimizer, index::Int, BT::Type{<:AbstractBridge})
    return Edge(index, _variable_nodes(b, BT), _constraint_nodes(b, BT))
end

# Method for objective bridges because they produce ObjectiveEdge.
function _edge(
    b::LazyBridgeOptimizer,
    index::Int,
    BT::Type{<:Objective.AbstractBridge},
)
    return ObjectiveEdge(
        index,
        _variable_nodes(b, BT),
        _constraint_nodes(b, BT),
        node(b, set_objective_function_type(BT)),
    )
end

_functionized_type(::Nothing) = nothing

function _functionized_type(
    ::Type{<:Constraint.ScalarFunctionizeBridge{T}},
) where {T}
    return MOI.ScalarAffineFunction{T}
end

function _functionized_type(
    ::Type{<:Constraint.VectorFunctionizeBridge{T}},
) where {T}
    return MOI.VectorAffineFunction{T}
end

function _functionized_type(b::LazyBridgeOptimizer, ::Type{MOI.VariableIndex})
    return _functionized_type(
        _first_functionize_bridge(
            b.constraint_bridge_types,
            Constraint.ScalarFunctionizeBridge,
        ),
    )
end

function _functionized_type(
    b::LazyBridgeOptimizer,
    ::Type{MOI.VectorOfVariables},
)
    return _functionized_type(
        _first_functionize_bridge(
            b.constraint_bridge_types,
            Constraint.VectorFunctionizeBridge,
        ),
    )
end

"""
    node(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractSet})

Return the `VariableNode` associated with set `S` in `b`.
"""
function node(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractSet})
    # If we support the set, the node is 0.
    if (
        S <: MOI.AbstractScalarSet &&
        MOI.supports_add_constrained_variable(b.model, S)
    ) || (
        S <: MOI.AbstractVectorSet &&
        MOI.supports_add_constrained_variables(b.model, S)
    )
        return VariableNode(0)
    end
    # If (S,) is stored in .variable_node, we've already added the node
    # previously.
    variable_node = get(b.variable_node, (S,), nothing)
    if variable_node !== nothing
        return variable_node
    end
    # This is a new (S,). We need to add it to the graph.
    variable_node = add_variable_node(b.graph)
    b.variable_node[(S,)] = variable_node
    push!(b.variable_types, (S,))
    F = MOI.Utilities.variable_function_type(S)
    if is_bridged(b, MOI.Reals)
        # The solver doesn't support adding free variables.
        FF = _functionized_type(b, F)
        if FF !== nothing
            # The cost is `3` because we assume the distance of the variable
            # node `MOI.Reals` is `1` (that is, it is bridged by
            # `Variable.FreeBridge`. We also use the functionize bridge which
            # has cost 1. And we add `+1` as we treat constrained variables as
            # constraints.
            set_variable_constraint_node(
                b.graph,
                variable_node,
                node(b, FF, S),
                3,
            )
        end
    else
        # We add `+1` as we treat constrained variables as constraints.
        set_variable_constraint_node(b.graph, variable_node, node(b, F, S), 1)
    end
    for (i, BT) in enumerate(b.variable_bridge_types)
        if Variable.supports_constrained_variable(BT, S)
            edge = _edge(b, i, Variable.concrete_bridge_type(BT, S))::Edge
            add_edge(b.graph, variable_node, edge)
        end
    end
    return variable_node
end

"""
    node(b::LazyBridgeOptimizer, ::Type{F}, ::Type{S})

Return the `ConstraintNode` associated with constraint `F`-in-`S` in `b`.
"""
function node(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    # If we support the constraint type, the node is 0.
    if MOI.supports_constraint(b.model, F, S)
        return ConstraintNode(0)
    end
    # If (F, S) is stored in .constraint_node, we've already added the node
    # previously.
    constraint_node = get(b.constraint_node, (F, S), nothing)
    if constraint_node !== nothing
        return constraint_node
    end
    # This is a new (F, S). We need to add it to the graph.
    constraint_node = add_constraint_node(b.graph)
    b.constraint_node[(F, S)] = constraint_node
    push!(b.constraint_types, (F, S))
    for (i, BT) in enumerate(b.constraint_bridge_types)
        if MOI.supports_constraint(BT, F, S)
            edge = _edge(b, i, Constraint.concrete_bridge_type(BT, F, S))::Edge
            add_edge(b.graph, constraint_node, edge)
        end
    end
    return constraint_node
end

"""
    node(b::LazyBridgeOptimizer, ::Type{F})

Return the `ObjectiveNode` associated with constraint `F` in `b`.
"""
function node(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractScalarFunction})
    # If we support the objective function, the node is 0.
    if MOI.supports(b.model, MOI.ObjectiveFunction{F}())
        return ObjectiveNode(0)
    end
    # If (f,) is stored in .objective_node, we've already added the node
    # previously.
    objective_node = get(b.objective_node, (F,), nothing)
    if objective_node !== nothing
        return objective_node
    end
    # This is a new (F,). We need to add it to the graph.
    objective_node = add_objective_node(b.graph)
    b.objective_node[(F,)] = objective_node
    push!(b.objective_types, (F,))
    for (i, BT) in enumerate(b.objective_bridge_types)
        if Objective.supports_objective_function(BT, F)
            bridge_type = Objective.concrete_bridge_type(BT, F)
            edge = _edge(b, i, bridge_type)::ObjectiveEdge
            add_edge(b.graph, objective_node, edge)
        end
    end
    return objective_node
end

"""
    _bridge_types(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge}),

A fallback helper to return the correct list of bridges given the type of `BT`.
"""
function _bridge_types(
    b::LazyBridgeOptimizer,
    ::Type{<:Variable.AbstractBridge},
)
    return b.variable_bridge_types
end

function _bridge_types(
    b::LazyBridgeOptimizer,
    ::Type{<:Constraint.AbstractBridge},
)
    return b.constraint_bridge_types
end

function _bridge_types(
    b::LazyBridgeOptimizer,
    ::Type{<:Objective.AbstractBridge},
)
    return b.objective_bridge_types
end

"""
    add_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})

Enable the use of the bridges of type `BT` by `b`.
"""
function add_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})
    if !has_bridge(b, BT)
        push!(_bridge_types(b, BT), BT)
        _reset_bridge_graph(b)
    end
    return
end

"""
    remove_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})

Disable the use of the bridges of type `BT` by `b`.
"""
function remove_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})
    bridge_types = _bridge_types(b, BT)
    i = findfirst(isequal(BT), bridge_types)
    if i === nothing
        error(
            "Cannot remove bridge `$BT` as it was never added or was already",
            " removed.",
        )
    end
    deleteat!(bridge_types, i)
    _reset_bridge_graph(b)
    return
end

"""
    has_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})

Return a `Bool` indicating whether the bridges of type `BT` are used by `b`.
"""
function has_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})
    return findfirst(isequal(BT), _bridge_types(b, BT)) !== nothing
end

# It only bridges when the constraint is not supporting, hence the name "Lazy"
function is_bridged(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractScalarSet})
    return !MOI.supports_add_constrained_variable(b.model, S)
end

function is_bridged(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractVectorSet})
    return !MOI.supports_add_constrained_variables(b.model, S)
end

function is_bridged(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return !MOI.supports_constraint(b.model, F, S)
end

function is_bridged(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractScalarFunction},
)
    return !MOI.supports(b.model, MOI.ObjectiveFunction{F}())
end

function bridge_index(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractSet})
    return bridge_index(b.graph, node(b, S))
end

function bridge_index(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return bridge_index(b.graph, node(b, F, S))
end

function bridge_index(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction})
    return bridge_index(b.graph, node(b, F))
end

supports_constraint_bridges(::LazyBridgeOptimizer) = true

function supports_bridging_constrained_variable(
    b::LazyBridgeOptimizer,
    S::Type{<:MOI.AbstractSet},
)
    variable_node = node(b, S)
    if !iszero(bridge_index(b.graph, variable_node))
        return true
    end
    constraint_node = b.graph.variable_constraint_node[variable_node.index]
    return constraint_node.index != INVALID_NODE_INDEX &&
           !iszero(bridge_index(b.graph, constraint_node))
end

function supports_bridging_constraint(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return !iszero(bridge_index(b, F, S))
end

function supports_bridging_objective_function(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractScalarFunction},
)
    return !iszero(bridge_index(b, F))
end

function is_variable_bridged(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractSet})
    return is_bridged(b, S) && is_variable_edge_best(b.graph, node(b, S))
end

function bridge_type(b::LazyBridgeOptimizer, S::Type{<:MOI.AbstractSet})
    bt = get(b.cached_bridge_type, (S,), nothing)
    if bt !== nothing
        return bt::Type
    end
    index = bridge_index(b, S)
    if iszero(index)
        throw(MOI.UnsupportedConstraint{MOIU.variable_function_type(S),S}())
    end
    new_bt = Variable.concrete_bridge_type(b.variable_bridge_types[index], S)
    b.cached_bridge_type[(S,)] = new_bt
    return new_bt::Type
end

function bridge_type(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    bt = get(b.cached_bridge_type, (F, S), nothing)
    if bt !== nothing
        return bt::Type
    end
    index = bridge_index(b, F, S)
    if iszero(index)
        throw(MOI.UnsupportedConstraint{F,S}())
    end
    new_bt =
        Constraint.concrete_bridge_type(b.constraint_bridge_types[index], F, S)
    b.cached_bridge_type[(F, S)] = new_bt
    return new_bt::Type
end

function bridge_type(
    b::LazyBridgeOptimizer,
    F::Type{<:MOI.AbstractScalarFunction},
)
    bt = get(b.cached_bridge_type, (F,), nothing)
    if bt !== nothing
        return bt::Type
    end
    index = bridge_index(b, F)
    if iszero(index)
        throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction{F}()))
    end
    new_bt = Objective.concrete_bridge_type(b.objective_bridge_types[index], F)
    b.cached_bridge_type[(F,)] = new_bt
    return new_bt::Type
end

function _func_name(::Type{Constraint.ScalarFunctionizeBridge})
    return "VariableIndex", "constraint"
end

function _func_name(::Type{Constraint.VectorFunctionizeBridge})
    return "VectorOfVariables", "constraint"
end

_func_name(::Type{Objective.FunctionizeBridge}) = "VariableIndex", "objective"

function _first_functionize_bridge(bridge_types, target_type)
    index = findfirst(bridge_type -> bridge_type <: target_type, bridge_types)
    if index === nothing
        return nothing
    end
    return bridge_types[index]
end

function _functionize_bridge(bridge_types, target_type)
    bridge_type = _first_functionize_bridge(bridge_types, target_type)
    if bridge_type === nothing
        func, name = _func_name(target_type)
        error(
            "Need to apply a `$target_type` to a `$func` $name because the " *
            "variable is bridged but no such $name bridge type was added. " *
            "Add one with `add_bridge`.",
        )
    end
    return bridge_type
end

function constraint_vector_functionize_bridge(b::LazyBridgeOptimizer)
    return _functionize_bridge(
        b.constraint_bridge_types,
        Constraint.VectorFunctionizeBridge,
    )
end

function constraint_scalar_functionize_bridge(b::LazyBridgeOptimizer)
    return _functionize_bridge(
        b.constraint_bridge_types,
        Constraint.ScalarFunctionizeBridge,
    )
end

function objective_functionize_bridge(b::LazyBridgeOptimizer)
    return _functionize_bridge(
        b.objective_bridge_types,
        Objective.FunctionizeBridge,
    )
end

function bridging_cost(b::LazyBridgeOptimizer, args...)
    return bridging_cost(b.graph, node(b, args...))
end

recursive_model(b::LazyBridgeOptimizer) = b

function MOI.compute_conflict!(model::LazyBridgeOptimizer)
    return MOI.compute_conflict!(model.model)
end


text:
using FractionalCalculus, Plots

tspan=collect(0:0.01:15)

examplederivative = []

for i in range(0.1, 0.9, step=0.1)
    push!(examplederivative, fracdiff(sin, i, tspan, 0.001, CaputoDiethelm()))
end

plot(tspan, examplederivative, title="Different order of derivative", linewidth=3, label=["α=0.1" "α=0.2" "α=0.3" "α=0.4" "α=0.5" "α=0.6" "α=0.7" "α=0.8" "α=0.9" "α=1"], legend=:bottomright)


text:
module RevisedAffineArithmetic

using FastRounding, IntervalArithmetic, StaticArrays

import Base: +, -, *, /, ==, convert, promote_rule
import IntervalArithmetic: ±

export RevisedAffineForm, ±, @rafvars, @affinize

include("affineform.jl")
include("arithmetic.jl")
include("conversions.jl")
include("vars.jl")
end


text:
module SingleColumnModels

include(joinpath("Grids", "FiniteDifferenceGrids.jl"))
include(joinpath("DomainDecompositions", "DomainDecompositions.jl"))
include(joinpath("StateVecs", "StateVecs.jl"))
include(joinpath("EDMF", "EDMF.jl"))

end # module


text:
module HackMoji

export runme

include("interpreter.jl")

end


text:
"""
    adjustedFire(t, proj, weapon, aero,tank)

Adjsut the elevation and azimuth angle to hit the target within a precision of 1 mm. The function returns the adjusted
elevation and azimuth angles in degrees.
"""
function adjustedFire!(t::AbstractTarget, proj::AbstractPenetrator, weapon::Gun, aero::DataFrame,tank::Tank;w=Wind(0.0, 0.0),atmosphere=nothing)
    precision=0.001
    missDist = 1e3
    QE = 0.0
    AZ = 0.0
    windIn = wind(tank, w)
    while missDist >precision
       # println("Find QE AZ")
    QE,AZ=QEfinderMPMM(t, proj, weapon,aero,w_bar=windIn,atm=atmosphere)
    #println("QE", " ", QE)
    #println("AZ", " ", AZ)
    if isnan(QE)
        println("NaN")
        break
    end
    weapon.QE = QE#rad2deg(QE)
    weapon.AZ = AZ#rad2deg(AZ)
    tank.canon.θ = QE#rad2deg(QE)
    tank.turret.ξ = AZ#rad2deg(AZ)
    windIn = wind(tank, w)
    proj.position=muzzlePos(tank)
    proj.velocity=muzzleVel(tank)
    #println("compute impact point")
    impactP= trajectoryMPMM(proj, t, weapon,aero,w_bar=windIn,atm=atmosphere)[1]
    #println("computed", " ", impactP)
    missDist = euclidean(impactP,t.position)
    #println("miss distance", " ", missDist)
    end
    return QE,AZ

end


text:
module TestCore

using FLoops
using Test

function test_two_states()
    xs = 1:10
    @floop begin
        s = 0
        p = 1
        for x in xs
            s += x
            p *= x
        end
    end
    @test s == sum(xs)
    @test p == prod(xs)
end

function test_complex_assignments()
    xs = 1:10
    @floop begin
        s = 0
        p = s + 1
        for x in xs
            s += x
            p *= x
        end
    end
    @test s == sum(xs)
    @test p == prod(xs)
end

function test_no_states()
    xs = 1:10
    ys = similar(xs)
    @floop for (i, x) in enumerate(xs)
        ys[i] = 2x
    end
    @test ys == 2 .* xs
end

function test_at_inbounds()
    xs = 1:10
    @floop begin
        s = 0
        @inbounds for x in xs
            s += x
        end
    end
    @test s == sum(xs)
end

function test_function()
    xs = 1:10
    @floop begin
        s = 0
        @inbounds for x in xs
            f() = return s + x
            s = f()
        end
    end
    @test s == sum(xs)
end

function test_let_and_break()
    function demo(xs)
        @floop begin
            s = 0
            p = 1
            for x in xs
                s += x
                p *= x
                let p = x
                    if p == 5
                        s = 0
                        break
                    end
                end
            end
        end
        return (s, p)
    end
    @test demo(1:10) == (0, prod(1:5))
    @test demo(1:4) == (sum(1:4), prod(1:4))
end

function test_let_and_continue()
    function demo(xs)
        @floop begin
            s = 0
            p = 1
            for x in xs
                let p = x
                    if p == 5
                        s = 0
                        continue
                    end
                end
                s += x
                p *= x
            end
        end
        return (s, p)
    end
    @test demo(1:10) == (sum(6:10), prod(1:10) ÷ 5)
    @test demo(1:4) == (sum(1:4), prod(1:4))
end

function test_continue()
    xs = 1:3
    ys = []
    @floop for x in xs
        x == 1 && continue
        push!(ys, x)
    end
    @test ys == [2, 3]
end

function test_break()
    xs = 1:3
    ys = []
    @floop for x in xs
        push!(ys, x)
        x == 2 && break
    end
    @test ys == [1, 2]
end

function test_return()
    function demo()
        ys = []
        @floop for x in 1:3
            push!(ys, x)
            if x == 2
                return ys
            end
        end
    end
    @test demo() == [1, 2]
end

function test_product()
    xs = 1:10
    ys = 2:3:15
    actual = []
    @floop for x in xs, y in ys
        push!(actual, (y, x))
    end
    desired = []
    for x in xs, y in ys
        push!(desired, (y, x))
    end
    @test actual == desired
end

function test_triangular_2()
    xs = 1:10
    actual = []
    @floop for x in xs, y in x:2:10
        push!(actual, (y, x))
    end
    desired = []
    for x in xs, y in x:2:10
        push!(desired, (y, x))
    end
    @test actual == desired
end

function test_triangular_3()
    xs = 1:10
    actual = []
    @floop for x in xs, y in x:2:10, z in x:y:20
        push!(actual, (y, x, z))
    end
    desired = []
    for x in xs, y in x:2:10, z in x:y:20
        push!(desired, (y, x, z))
    end
    @test !isempty(actual)  # make sure it tests _something_
    @test actual == desired
end

function test_internal_goto()
    xs = 1:10
    @floop begin
        s = 0
        for x in xs
            @goto a
            s += x
            @label a
        end
    end
    @test s == 0
end

function test_external_goto()
    xs = 1:10
    @floop begin
        s = 0
        for x in xs
            s += x
            s > 5 && @goto a
        end
    end
    @label a
    @test s == 6
end

end  # module


text:
s = "stressed"
println(reverse(s))


text:
module BinaryRecordReader

using MacroTools

export @construct_reader, construct_reader_deep, construct_reader_shallow

global const def_reads = Set{DataType}()

macro construct_reader(Struct, dims)
    array_specs = eval(dims)
    s = eval(:($__module__.$Struct))
    esc(construct_reader_exp(s, array_specs))
end

"""
    @construct_reader S array_specs

Generates a `Base.read(io::IOStream, ::Type{S})` method that can used read type `S`, which is not an `isbitstype`.

# Example
```julia-repl
julia> struct CompoundStruct
           long::Float32
           COG::UInt16
           data::Matrix{Float64}
       end
julia> @construct_reader CompoundStruct (data=(10,10),)
```
The above will generate a suitable `read` method assuming that the `data` field is a 10x10 matrix.

See also: `construct_reader_deep`, `@construct_reader`
"""
macro construct_reader(Struct)
    esc(:(@construct_reader $Struct NamedTuple()))
end

function construct_reader_exp(s, array_specs)
    isbitstype(s) && (@warn("Attempt to create a reader for plain data type $s"); return nothing) # short ciruite for isbitstype
    expr = quote
        function Base.read(io::IOStream, ::Type{$s})
        end end |> prettify
    for (f, t) in  zip(fieldnames(s), fieldtypes(s))
        if haskey(array_specs, f)
            t <: Array || throw(TypeError(:construct_reader_exp, "construction of intenral reader", Array, t))
            et = eltype(t)
            isbitstype(et) || et ∈ def_reads || throw(ArgumentError("Array element $et has no properly defined read function"))
            push!(expr.args[2].args, :($f = read!(io, $t(undef, $(array_specs[f])))))
        else
            isbitstype(t) || t ∈ def_reads || throw(ArgumentError("$t has no properly defined read function"))
            push!(expr.args[2].args, :($f = read(io, $t)))
        end
    end
    fields = [f for f in fieldnames(s)]
    push!(expr.args[2].args, Expr(:call, s, fields...))
    push!(def_reads, s) # marking it as defined
    expr
end

function strip_array_specs(s, array_specs)
    a = Dict{Symbol,Tuple}()
    b = Vector{Symbol}()
    for k in keys(array_specs)
        parts = split(string(k), "__")
        if length(parts) == 1
            a[k] = array_specs[k]
        elseif parts[1] == string(s)
            a[Symbol(parts[2])] = array_specs[k]
        else
            push!(b, Symbol(parts[1]))
        end
    end
    a, b
end

"""
    construct_reader_shallow(S,array_specs)

Generates a `Base.read(io::IOStream, ::Type{S})` method that can used read type `S`, which is not an `isbitstype`.
The named tuble, `array_spec`, is used to define the dimensionality of the field in question.

# Example
```julia-repl
julia> struct CompoundStruct
           long::Float32
           COG::UInt16
           data::Matrix{Float64}
       end
julia> construct_reader_shallow(CompoundStruct, (data=(10,10),))
```
The above will generate a suitable `read` method assuming that the `data` field is a 10x10 matrix.

See also: `construct_reader_deep`, `@construct_reader`
"""
function construct_reader_shallow(s, array_specs)
    eval(construct_reader_exp(s, array_specs))
end

"""
    construct_reader_deep(S,array_specs)

Generates a `Base.read(io::IOStream, ::Type{S})` method that can used read type `S`, which is not an `isbitstype`.
If the class in question has other fields that made of up of other non-`isbitstype`s a full specification put in place
in the `array_spec`. The syntax for   `types` is `SubTypeName__fieldname`.

# Example
```julia-repl
julia> struct C
            lat::Float32
            time::Matrix{Int32}
        end
julia> struct B
            fun::UInt8
            funMat::Matrix{UInt16}
        end
julia> struct A
            long::Float32
            COG::UInt16
            funky::AnotherSimpleRecMat
            data::Matrix{SimpleRecMat}
        end

julia> construct_reader_deep(A, (data=(5,5), B__funMat=(2,2), C__time =(3,3)))

```
The above will generate three `read` methods for all the three types defined.

See also: `construct_reader_shallow`, `@construct_reader`
"""
function construct_reader_deep(s, array_specs)
    isbitstype(s) && (@warn("Attempt to create a reader for plain data type $s"); return nothing) # short ciruite for isbitstype
    a, b = strip_array_specs(s, array_specs)
    for (f, t) in  zip(fieldnames(s), fieldtypes(s))
        if haskey(a, f)
            t <: Array || throw(TypeError(:construct_reader_deep, "construction of intenral reader", Array, t))
            et = eltype(t)
            if !isbitstype(et) && (Symbol(et) ∈ b || et ∉ def_reads )
                construct_reader_deep(et, array_specs)
            end
        else
            isbitstype(t) || t ∈ def_reads || t <: Array && throw(ArgumentError("$t has no properly defined read function"))
            if !isbitstype(t) && (Symbol(t) ∈ b || t ∉ def_reads )
                construct_reader_deep(t, array_specs)
            end
        end
    end
    eval(construct_reader_exp(s, a))
end

end


text:
module TSMLextra

using Reexport

@reexport using TSML
export fit!, transform!

include("system.jl")
using .System

include("datareader.jl")
using .DataReaders
export DataReader

include("datawriter.jl")
using .DataWriters
export DataWriter

if LIB_CRT_AVAILABLE # from System module
    include("caret.jl")
    using .CaretLearners
    export CaretLearner
end

if LIB_SKL_AVAILABLE # from System module
    include("scikitlearn.jl")
    using .SKLearners
    export SKLearner
    include("skpreprocessor.jl")
    using .SKPreprocessors
    export SKPreprocessor
    include("skcrossvalidator.jl")
    using .SKCrossValidators
    export crossvalidate
end

end # module

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra