text:
# Pairwise Alignment
# ==================
#
# Pairwise alignment type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioAlignments.jl/blob/master/LICENSE.md

"""
Pairwise alignment
"""
mutable struct PairwiseAlignment{S1,S2}
    a::AlignedSequence{S1}
    b::S2
end

function Base.iterate(aln::PairwiseAlignment, ij=(2,1))
    i, j = ij
    if i > lastindex(aln.a.aln.anchors)
        return nothing
    end

    anchors = aln.a.aln.anchors
    anchor = anchors[i]
    seq = aln.a.seq
    ref = aln.b
    seqpos = anchors[i-1].seqpos
    refpos = anchors[i-1].refpos

    if ismatchop(anchor.op)
        x = seq[seqpos + j]
        y = ref[refpos + j]
    elseif isinsertop(anchor.op)
        x = seq[seqpos + j]
        y = BioSymbols.gap(eltype(ref))
    elseif isdeleteop(anchor.op)
        x = BioSymbols.gap(eltype(seq))
        y = ref[refpos + j]
    else
        @assert false
    end

    if ismatchop(anchor.op) || isinsertop(anchor.op)
        if j < anchor.seqpos - seqpos
            j += 1
        else
            i += 1
            j = 1
        end
    else
        if j < anchor.refpos - refpos
            j += 1
        else
            i += 1
            j = 1
        end
    end

    return (x, y), (i, j)
end

Base.length(aln::PairwiseAlignment) = count_aligned(aln)
Base.eltype(::Type{PairwiseAlignment{S1,S2}}) where {S1,S2} = Tuple{eltype(S1),eltype(S2)}

"""
    count(aln::PairwiseAlignment, target::Operation)

Count the number of positions where the `target` operation is applied.
"""
function Base.count(aln::PairwiseAlignment, target::Operation)
    anchors = aln.a.aln.anchors
    n = 0
    for i in 2:lastindex(anchors)
        op = anchors[i].op
        if op == target
            if ismatchop(op) || isinsertop(op)
                n += anchors[i].seqpos - anchors[i-1].seqpos
            elseif isdeleteop(op)
                n += anchors[i].refpos - anchors[i-1].refpos
            end
        end
    end
    return n
end

"""
    count_matches(aln)

Count the number of matching positions.
"""
count_matches(aln::PairwiseAlignment) = count(aln, OP_SEQ_MATCH)

"""
    count_mismatches(aln)

Count the number of mismatching positions.
"""
count_mismatches(aln::PairwiseAlignment) = count(aln, OP_SEQ_MISMATCH)

"""
    count_insertions(aln)

Count the number of inserting positions.
"""
count_insertions(aln::PairwiseAlignment) = count(aln, OP_INSERT)

"""
    count_deletions(aln)

Count the number of deleting positions.
"""
count_deletions(aln::PairwiseAlignment) = count(aln, OP_DELETE)

"""
    count_aligned(aln)

Count the number of aligned positions.
"""
function count_aligned(aln::PairwiseAlignment)
    anchors = aln.a.aln.anchors
    n = 0
    for i in 2:lastindex(anchors)
        op = anchors[i].op
        if ismatchop(op) || isinsertop(op)
            n += anchors[i].seqpos - anchors[i-1].seqpos
        elseif isdeleteop(op)
            n += anchors[i].refpos - anchors[i-1].refpos
        end
    end
    return n
end

"""
    seq2ref(aln::PairwiseAlignment, i::Integer)::Tuple{Int,Operation}

Map a position `i` from the first sequence to the second.
"""
function seq2ref(aln::PairwiseAlignment, i::Integer)::Tuple{Int,Operation}
    return seq2ref(aln.a, i)
end

"""
    ref2seq(aln::PairwiseAlignment, i::Integer)::Tuple{Int,Operation}

Map a position `i` from the second sequence to the first.
"""
function ref2seq(aln::PairwiseAlignment, i::Integer)::Tuple{Int,Operation}
    return ref2seq(aln.a, i)
end

# Printers
# --------

function Base.show(io::IO, aln::PairwiseAlignment)
    println(io, summary(aln), ':')
    print(io, aln)
end

function Base.print(io::IO, aln::PairwiseAlignment)
    print_pairwise_alignment(io, aln)
end

function print_pairwise_alignment(io::IO, aln::PairwiseAlignment; width::Integer=60)
    seq = aln.a.seq
    ref = aln.b
    anchors = aln.a.aln.anchors
    # width of position numbers
    posw = ndigits(max(anchors[end].seqpos, anchors[end].refpos)) + 1

    i = 0
    seqpos = anchors[1].seqpos
    refpos = anchors[1].refpos
    seqbuf = IOBuffer()
    refbuf = IOBuffer()
    matbuf = IOBuffer()
    next_xy = iterate(aln)
    while next_xy !== nothing
        (x, y), s = next_xy
        next_xy = iterate(aln ,s)

        i += 1
        if x != BioSymbols.gap(eltype(seq))
            seqpos += 1
        end
        if y != BioSymbols.gap(eltype(ref))
            refpos += 1
        end

        if i % width == 1
            print(seqbuf, "  seq:", lpad(seqpos, posw), ' ')
            print(refbuf, "  ref:", lpad(refpos, posw), ' ')
            print(matbuf, " "^(posw + 7))
        end

        print(seqbuf, x)
        print(refbuf, y)
        print(matbuf, x == y ? '|' : ' ')

        if i % width == 0
            print(seqbuf, lpad(seqpos, posw))
            print(refbuf, lpad(refpos, posw))
            print(matbuf)

            println(io, String(take!(seqbuf)))
            println(io, String(take!(matbuf)))
            println(io, String(take!(refbuf)))

            if next_xy !== nothing
                println(io)
                seek(seqbuf, 0)
                seek(matbuf, 0)
                seek(refbuf, 0)
            end
        end
    end

    if i % width != 0
        print(seqbuf, lpad(seqpos, posw))
        print(refbuf, lpad(refpos, posw))
        print(matbuf)

        println(io, String(take!(seqbuf)))
        println(io, String(take!(matbuf)))
        println(io, String(take!(refbuf)))
    end
end


text:
"""
`findindices(data::Vector)`

Find indices of peaks in `data`. For flat peaks, return only the lowest index.
"""
function findindices(data::Vector)

    # Print warning if the data vector contains NaN values
    if any(isnan.(data))
        @warn "Data vector contains NaN values!"
    end

    n = length(data)  # length of data
    indices = Int[]  # index vector for found peaks
    index_candidate = 0  # index of a potential peak

    # iterate data points
    for i = 2:n-1

        # if previous data point was smaller,
        # the current index becomes a candidate
        if data[i-1] < data[i]
            index_candidate = i
        end

        # if a candidate exists and future data point is smaller,
        # then the point is a peak and its index saved
        if index_candidate > 0 && data[i] > data[i+1]
            push!(indices, index_candidate)
            index_candidate = 0
        end

    end

    return indices

end

"""
`findprominences(data::Vector, indices::Vector)`

Find prominence for peaks at `indices` of `data` vector.
"""
function findprominences(data::Vector, indices::Vector)

    # number of peaks and vector of peak prominences
    n_indices = length(indices)
    prominences = typeof(data)(undef, n_indices)

    # iterate peak indices
    for i_index = 1:n_indices

        # current peak index and peak value
        current_index = indices[i_index]
        peak_value = data[current_index]

        # run left until a value higher than the current peak is reached
        left_index = 1
        for i = i_index-1:-1:1
            if peak_value < data[indices[i]]
                left_index = indices[i]
                break
            end
        end

        # find minimum of the left interval
        left_minimum = minimum(data[left_index:indices[i_index]])

        # run right until a value higher than the current peak is reached
        right_index = length(data)
        for i = i_index+1:1:n_indices
            if peak_value < data[indices[i]]
                right_index = indices[i]
                break
            end
        end

        # find minimum of the right interval
        right_minimum = minimum(data[indices[i_index]:right_index])

        # the higest of the two minima is a reference point
        reference_point = max(left_minimum, right_minimum)

        # prominence
        prominences[i_index] = peak_value - reference_point

    end

    return prominences

end

"""
`findwidthbounds(data, x, indices, prominences)`

Find width bounds for peaks in `data` vector. `x` is location vector for `data`.
"""
function findwidthbounds(data, x, indices, prominences)

    n_indices = length(indices)
    width_bounds = Matrix{Float64}(undef, n_indices, 2)

    for i_index = 1:n_indices

        current_index = indices[i_index]
        peak = data[current_index]  # peak

        # reference height at which to measure the peak width
        height_reference = peak - prominences[i_index] * 0.5

        # find point to the left from the peak where w_ref intersect data
        left_index = 1
        for i = current_index:-1:2
            y1 = data[i]
            y2 = data[i-1]
            if y2 <= height_reference <= y1
                # interpolate the precise point
                x1 = x[i]
                xdiff = x1 - x[i-1]
                left_index = x1 - xdiff / (y2 - y1) * (height_reference - y1)
                break
            end
        end

        # find point to the right from the peak where w_ref intersect data
        right_index = length(data)
        for i = current_index:1:length(data)-1
            y1 = data[i]
            y2 = data[i+1]
            if y2 <= height_reference <= y1
                # interpolate the precise point
                x1 = x[i]
                xdiff = x[i+1] - x1
                right_index = x1 + xdiff / (y2 - y1) * (height_reference - y1)
                break
            end
        end

        width_bounds[i_index, :] = [left_index, right_index]

    end

    return width_bounds

end

"""
`findwidthbounds(data, x, indices)`

Find width bounds using a half height of a peak as a reference.
"""
function findwidthbounds(data, x, indices)

    n_indices = length(indices)
    width_bounds = Matrix{Float64}(undef, n_indices, 2)

    for i_index = 1:n_indices

        current_index = indices[i_index]
        peak = data[current_index]  # peak

        # for negative peaks return NaNs
        if peak < 0
            width_bounds[i_index, :] = [NaN, NaN]
            continue
        end

        # reference height at which to measure the peak width
        height_reference = 0.5 * peak

        # find point to the left from the peak where w_ref intersect data
        if i_index == 1
            left_index = 1
        else
            i_minimum = argmin(data[indices[i_index-1]:current_index])
            left_index = indices[i_index-1] + i_minimum -1
        end
        left_x = x[left_index]
        for i = current_index:-1:(left_index+1)
            y1 = data[i]
            y2 = data[i-1]
            if y2 <= height_reference <= y1
                # interpolate the precise point
                x1 = x[i]
                xdiff = x1 - x[i-1]
                left_x = x1 - xdiff / (y2 - y1) * (height_reference - y1)
                break
            end
        end

        # find point to the right from the peak where w_ref intersect data
        if i_index == n_indices
            right_index = length(data)
        else
            i_minimum = argmin(data[current_index:indices[i_index+1]])
            right_index = current_index + i_minimum -1
        end
        right_x = x[right_index]
        for i = current_index:1:(right_index-1)
            y1 = data[i]
            y2 = data[i+1]
            if y2 <= height_reference <= y1
                # interpolate the precise point
                x1 = x[i]
                xdiff = x[i+1] - x1
                right_x = x1 + xdiff / (y2 - y1) * (height_reference - y1)
                break
            end
        end

        width_bounds[i_index, :] = [left_x, right_x]

    end

    return width_bounds

end

"""
`findpeaks(data::Vector, x=1:length(data); kwargs...)`

**Argument**
- `data`: Input signal vector.

**Optional Argument**
- `x`: Locations for the input signal vector.

**Keyword Arguments**
- `npeaks`: Maximum number of peaks to return.
- `sortstr`: Sort peaks. Possible options are "ascend" (the smallest
    peak first) or "descend" (the largest peak first).
- `sortref`: Reference property used to sort the peaks. Possible options are
    "height" (default) and "prominence". This option is ignored if `sortstr` is
    not specified.
- `minheight`: Minimum peak height.
- `minprominence`: Minimum peak prominence.
- `threshold`: Minimum height difference between a peak and the neighboring
    points.
- `mindistance`: Minimum distance between neighboring peaks.
- `widthreference`: "halfheight" will use half heights, instead of half
    prominences, as reference heights to calculate the peak widths.
- `minwidth`: Minimum peak width.
- `maxwidth`: Maximum peak width.

**Returns**
- Struct `PeakResults` holding an array of all peaks and their properties.
"""
function findpeaks(data::Vector, x=1:length(data); kwargs...)

    # find peak indices
    indices = findindices(data)

    # apply minimum height
    if :minheight in keys(kwargs)
        minheight = kwargs[:minheight]
        peaks = data[indices]
        inds = peaks .>= minheight
        indices = indices[inds]
    end

    # apply threshold
    if :threshold in keys(kwargs)
        threshold = kwargs[:threshold]
        inds = []
        n = length(indices)
        for i = 1:n
            index = indices[i]
            height_difference = min(data[index] - data[index-1],
                data[index] - data[index+1])
            if height_difference >= threshold
                push!(inds, i)
            end
        end
        indices = indices[inds]
    end

    # find peak prominences
    prominences = findprominences(data, indices)

    # apply minimal prominence
    if :minprominence in keys(kwargs)
        minprominence = kwargs[:minprominence]
        inds = prominences .>= minprominence
        indices = indices[inds]
        prominences = prominences[inds]
    end

    # find peak width bounds
    if :widthreference in keys(kwargs)
        widthreference = kwargs[:widthreference]
        if widthreference == "halfheight"
            width_bounds = findwidthbounds(data, x, indices)
        else
            width_bounds = findwidthbounds(data, x, indices, prominences)
        end
    else
        width_bounds = findwidthbounds(data, x, indices, prominences)
    end

    # create Peaks struct
    peaks = data[indices]  # peak heights
    locations = x[indices]  # peak locations
    pr = PeakResults(indices, peaks, locations, prominences, width_bounds)

    # apply minwidth
    if :minwidth in keys(kwargs)
        minwidth = kwargs[:minwidth]
        widths = peakwidths(pr)
        inds = []
        n = length(widths)
        for i = 1:n
            width = widths[i]
            if width >= minwidth
                push!(inds, i)
            end
        end
        pr = pr[inds]
    end

    # apply maxwidth
    if :maxwidth in keys(kwargs)
        maxwidth = kwargs[:maxwidth]
        widths = peakwidths(pr)
        inds = []
        n = length(widths)
        for i = 1:n
            width = widths[i]
            if width <= maxwidth
                push!(inds, i)
            end
        end
        pr = pr[inds]
    end

    # apply minimal distance criterium
    if :mindistance in keys(kwargs)

        mindistance = kwargs[:mindistance]

        pr_remaining = pr
        selected_indices = []  # indices of selected peaks
        while !isempty(pr_remaining)

            # extract needed peak properties
            peaks = pr_remaining.peaks
            indices = pr_remaining.indices
            locations = pr_remaining.locations

            # select the highest peak
            index_maxpeak = argmax(peaks)
            push!(selected_indices, indices[index_maxpeak])

            # remove all peaks that are closer than mindistance to the selected peak
            location = locations[index_maxpeak]
            inds = abs.(locations .- location) .>= mindistance
            pr_remaining = pr_remaining[inds]

        end

        # select a subset of pr based on the selected indices
        sort!(selected_indices)
        inds = [x in selected_indices for x in pr.indices]
        pr = pr[inds]
    end

    # sort peaks
    if (:sortref in keys(kwargs)) && !(:sortstr in keys(kwargs))
        @warn "sortstr is not specified, sortref will be ignored"
    end

    if :sortstr in keys(kwargs)
        sortstr = kwargs[:sortstr]
        if :sortref in keys(kwargs)
            sortref = kwargs[:sortref]
        else
            sortref = "height"
        end
        if sortstr == "ascend"
            pr = sort(pr; ref=sortref)
        elseif sortstr == "descend"
            pr = sort(pr; rev=true, ref=sortref)
        end
    end

    # select demanded number of peaks
    if :npeaks in keys(kwargs)
        npeaks = kwargs[:npeaks]
        n = length(indices)
        maxind = min(n, npeaks)
        pr = pr[1:maxind]
    end

    return pr

end

ext: jl

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
"""
Computes linear polynomial basis
OUTPUT: 2D array
"""
function linear_polynomial_basis(x)#linear polynomial basis
    @assert typeof(x) <:Array{T, 1} where T  ||  typeof(x) <:Array{T, 2} where T
    if typeof(x)<:Array{T, 1} where T
        return hcat(1, reshape(x, 1, length(x)))
    else
        return hcat(ones(size(x, 1), 1), x)
    end
end

"""
Computes constant polynomial basis
OUTPUT: 2D Array
"""
function constant_basis(x)
    @assert typeof(x) <:Array{T, 1} where T  ||  typeof(x) <:Array{T, 2} where T
    if typeof(x) <:Array{T, 1} where T
        return reshape([1], 1, 1)
    else
        return reshape(ones(size(x, 1), 1), size(x, 1), 1)
    end
end

if false #test
a = constant_basis([1, 2, 3])
b = constant_basis(reshape([1; 2; 3], 3, 1))
c = linear_polynomial_basis([1 2; 3 4; 5 6])
d = linear_polynomial_basis([1;2;3])
end

"""
Choose covariate type based on p
"""
function covariate_fun(x, p)
    n = size(x, 1)
    d = size(x, 2)
    if p == 1
        return ones(n, 1)
    elseif p == 1 + d
        return hcat(ones(n), x)
    else
        throw(ArgumentError("Only support constant or linear convariate."))
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
## Element type
struct IIDMGenerator <: AbstractNodeLabel
  id::String
  minP::Float64 # MW
  maxP::Float64 # MW
  targetV::Float64 # p.u.
  volt_regulator_on::Bool
  S_bounds::Array{Tuple{Complex, Complex}} # MW+im*MVar
end

## 1. Variables creation
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String) where T <: IIDMGenerator
  bus_vars[elemid] = Variable(variable_name("Sgen", bus, elemid, scenario), Complex)
  return
end

# function create_vars!(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T <: MyType
#   return
# end

## 2. Power balance
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: IIDMGenerator
  cstrname = "UNIT"
  return [cstrname, Polynomial(bus_vars[elemid]), 0, 0]
end

# function Sorig(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T<:MyType
#   return Polynomial()
# end

# function Sdest(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T<:MyType
#   return Polynomial()
# end

## 3. Constraints creation
function constraint(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: IIDMGenerator
  bounds = sort(element.S_bounds, by=x->real(x[1]))
  length(bounds) == 2 || warn("$(length(bounds))x2 ≠ 4 points for power domain of $elemid, $busid")

  Sgen = bus_vars[elemid]
  constraints = Dict{String, Constraint}()

  constraints["Cstr_$(busid)_$(elemid)_real"] = real(bounds[1][1]) << Sgen << real(bounds[2][1])

  a, θ = bounds[1][2], angle(bounds[2][2] - bounds[1][2])
  p =  real((Sgen-a)*exp(im*(π/2-θ)))
  constraints["Cstr_$(busid)_$(elemid)_upper"] = 0 << p << Inf

  a, θ = bounds[1][1], angle(bounds[2][1] - bounds[1][1])
  p = real((Sgen-a)*exp(im*(-π/2-θ)))
  constraints["Cstr_$(busid)_$(elemid)_lower"] = 0 << p << Inf

  return constraints
end

# function constraint(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T <: MyType
#   return Dict{String, Constraint}()
# end

## 4. Generator cost
function cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_var::Dict{String, Variable}, Snode::Polynomial, lb, ub) where T <: IIDMGenerator
  Sgen = gencost = Polynomial()
  if elem_formulation == :NbMinVar
    _, S, lb_gen, ub_gen = Snodal(element, bus, elemid, elem_formulation, bus_elems_var)
    abs((lb-lb_gen) - (ub-ub_gen)) ≤ 5e-5 || warn("Power missing at $bus... $(abs((lb-lb_gen) - (ub-ub_gen)))")
    Sgen = add(Snode, -S - (ub-ub_gen))
  elseif elem_formulation == :NewVar
    Sgen = Polynomial(bus_elems_var[elemid])
  else
    warn("No cost contribution for gen $elemid, $bus")
  end

  add!(Sgen, conj(Sgen))
  gencost = (Sgen)*0.5
  return gencost
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Dierckx

# Convenience function for wavelength conversion
@inline aa_to_invum(wave::Real) = 10000 / wave

# --------------------------------------------------------------------------------

# Optical coefficients
const ccm89_ca = [1.0, 0.17699, -0.50447, -0.02427, 0.72085, 0.01979, -0.7753, 0.32999, 0.0]
const ccm89_cb = [0.0, 1.41338, 2.28305, 1.07233, -5.38434, -0.62251, 5.3026, -2.09002, 0.0]

const od94_ca = [1.0, 0.104, -0.609, 0.701, 1.137, -1.718, -0.827, 1.647, -0.505]
const od94_cb = [0.0, 1.952, 2.908, -3.989, -7.985, 11.102, 5.491, -10.805, 3.347]

"""
    CCM89(;Rv=3.1)

Clayton, Cardelli and Mathis (1989) dust law.

Returns E(B-V) in magnitudes at the given wavelength relative to the extinction
at 5494.5 Å. The default support is [1000, 33333]. Outside of that range this will return 0. `Rv` is the selective extinction and is valid over [2, 6]. A typical value for the Milky Way is 3.1.

# References
[Clayton,Cardelli and Mathis (1989)](https://ui.adsabs.harvard.edu/abs/1989ApJ...345..245C)
"""
@with_kw struct CCM89 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::CCM89)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return ccm89_invum(x, law.Rv, ccm89_ca, ccm89_cb)
end

bounds(::Type{CCM89}) = (1000, 33333)

"""
    OD94(;Rv=3.1)

O'Donnell (1994) dust law.

This is identical to the Clayton, Cardelli and Mathis (1989) dust law, except
for different coefficients used in the optical (3030.3 Å to 9090.9 Å).

# References
[O'Donnell (1994)](https://ui.adsabs.harvard.edu/abs/1994ApJ...422..158O)

# See Also
[`CCM89`](@ref)
"""
@with_kw struct OD94 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::OD94)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return ccm89_invum(x, law.Rv, od94_ca, od94_cb)
end

bounds(::Type{OD94}) = (1000, 33333)

"""
    DustExtinction.ccm89_invum(x, Rv, c_a, c_b)

The algorithm used for the [`CCM89`](@ref) extinction law, given inverse microns, Rv, and a set of coefficients for use in the optical (only difference between ccm89 and od94). For more information, seek the original paper.
"""
function ccm89_invum(x::Real, Rv::Real, c_a::Vector{<:Real}, c_b::Vector{<:Real})
    if x < 0.3
        error("out of bounds of CCM89, support is over $(bounds(CCM89)) angstrom")
    elseif x < 1.1  # Near IR
        y = x^1.61
        a = 0.574y
        b = -0.527y
    elseif x < 3.3  # Optical
        y = x - 1.82
        yn = 1.0
        a = @evalpoly y c_a[1] c_a[2] c_a[3] c_a[4] c_a[5] c_a[6] c_a[7] c_a[8] c_a[9]
        b = @evalpoly y c_b[1] c_b[2] c_b[3] c_b[4] c_b[5] c_b[6] c_b[7] c_b[8] c_b[9]
    elseif x < 8.0  # NUV
        a =  1.752 - 0.316x - (0.104 / ((x - 4.67)^2 + 0.341))
        b = -3.090 + 1.825x + (1.206 / ((x - 4.62)^2 + 0.263))
        if x > 5.9 # Far NUV
            y = x - 5.9
            a += @evalpoly y 0.0 0.0 -0.04473 -0.009779
            b += @evalpoly y 0.0 0.0 0.213 0.1207
        end
    elseif x ≤ 10.0 # FUV
        y = x - 8.0
        a = @evalpoly y -1.073 -0.628 0.137 -0.07
        b = @evalpoly y 13.67 4.257 -0.42 0.374
    else
        error("out of bounds of CCM89, support is over $(bounds(CCM89)) angstrom")
    end
    return a + b / Rv
end

# --------------------------------------------------------------------------------

"""
    CAL00(;Rv=4.05)

Calzetti et al. (2000) Dust Law.

Returns E(B-V) in magnitudes at the given wavelength. `λ` is the wavelength in Å
 and has support over [1200, 22000]. Outside of that range this will return 0.

Calzetti et al. (2000) developed a recipe for dereddening the spectra of
galaxies where massive stars dominate the radiation output. They found the best
 fit value for such galaxies was 4.05±0.80.

# References
[Calzetti et al. (2000)](https://ui.adsabs.harvard.edu/abs/2000ApJ...533..682C)
"""
@with_kw struct CAL00 <: ExtinctionLaw
    Rv::Float64 = 4.05
end
function (law::CAL00)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return cal00_invum(x, law.Rv)
end

bounds(::Type{CAL00}) = (1200, 22000)

"""
    DustExtinction.cal00_invum(x, Rv)

The algorithm used for the [`CAL00`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function cal00_invum(x::Real, Rv::Real)
    if x > 1 / 0.12
        error("out of bounds of CAL00, support is over $(bounds(CAL00)) angstrom")
    elseif x > 1 / 0.63
        k = @evalpoly x -2.156 1.509 -0.198 0.011
    elseif x > 1 / 2.2
        k = @evalpoly x -1.857 1.040
    else
        error("out of bounds of CAL00, support is over $(bounds(CAL00)) angstrom")
    end
    return 1.0 + 2.659 * k / Rv
end

"""
    VCG04(;Rv=3.1)

Valencic, Clayton, & Gordon (2004) dust law.

This model applies to the UV spectral region all the way to 912 Å.
This model was not derived for the optical or NIR.

# References
[Valencic, Clayton, & Gordon (2004)](https://ui.adsabs.harvard.edu/abs/2004ApJ...616..912V/)
"""
@with_kw struct VCG04 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::VCG04)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return vcg04_invum(x, law.Rv)
end

bounds(::Type{VCG04}) = (1250, 3030.3)

"""
    DustExtinction.vcg04_invum(x, Rv)

The algorithm used for the [`VCG04`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function vcg04_invum(x::Real, Rv::Real)
    if 3.3 ≤ x ≤ 8.0  # NUV
        a = 1.808 - 0.215 * x - 0.134 / ((x - 4.558)^2 + 0.566)
        b = -2.350 + 1.403 * x + 1.103 / ((x - 4.587)^2 + 0.263)
    else
        error("out of bounds of VCG04, support is over $(bounds(VCG04)) angstrom")
    end
    if 5.9 ≤ x ≤ 8.0  # far-NUV
        y = x - 5.9
        a += @evalpoly y 0.0 0.0 -0.0077 -0.0030
        b += @evalpoly y 0.0 0.0 0.2060 0.0550
    end

    return a + b / Rv
end

"""
    GCC09(;Rv=3.1)

Gordon, Cartledge, & Clayton (2009) dust law.

This model applies to the UV spectral region all the way to 909.09 Å.
This model was not derived for the optical or NIR.

# References
[Gordon, Cartledge, & Clayton (2009)](https://ui.adsabs.harvard.edu/abs/2009ApJ...705.1320G/)
"""
@with_kw struct GCC09 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::GCC09)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return gcc09_invum(x, law.Rv)
end

bounds(::Type{GCC09}) = (909.09, 3030.3)

"""
    DustExtinction.gcc09_invum(x, Rv)

The algorithm used for the [`GCC09`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function gcc09_invum(x::Real, Rv::Real)
    if 3.3 ≤ x ≤ 11.0  # NUV
        a = 1.894 - 0.373 * x - 0.0101 / ((x - 4.57)^2 + 0.0384)
        b = -3.490 + 2.057 * x + 0.706 / ((x - 4.59)^2 + 0.169)
    else # out of bounds
        error("out of bounds of GCC09, support is over $(bounds(GCC09)) angstrom")
    end
    if 5.9 ≤ x ≤ 11.0  # far-NUV
        y = x - 5.9
        a += @evalpoly y 0.0 0.0 -0.110 -0.0100
        b += @evalpoly y 0.0 0.0 0.531 0.0544
    end

    return a + b / Rv
end

# x value above which FM90 parametrization used
const f99_x_cutval_uv = aa_to_invum(2700)
# required UV points for spline interpolation
const f99_x_splineval_uv = aa_to_invum.((2700, 2600))

# Shape models used by F99 and F04
function _curve_F99_method(
    x,
    Rv,
    c1,
    c2,
    c3,
    c4,
    x0,
    gamma,
    optnir_axav_x,
    optnir_axav_y,
    )

    # add in required spline points, otherwise just spline points
    if x >= f99_x_cutval_uv
        xuv = (f99_x_splineval_uv..., x)
    else
        xuv = f99_x_splineval_uv
    end

    # FM90 model and values
    fm90_model = FM90(c1=c1, c2=c2, c3=c3, c4=c4, x0=x0, gamma=gamma)
    # evaluate model and get results in A(x)/A(V)
    axav_fm90 = @. fm90_model(aa_to_invum(xuv)) / Rv + 1

    # ignore the spline points
    if x >= f99_x_cutval_uv
        axav = last(axav_fm90)
    else
        # **Optical Portion**

        # save spline points
        y_splineval_uv = axav_fm90

        # spline points
        x_splineval_optir = (0.0, optnir_axav_x...)

        # determine optical/IR values at spline points
        y_splineval_optir = (0.0, optnir_axav_y...)
        spline_x = (x_splineval_optir..., f99_x_splineval_uv...)
        spline_y = (y_splineval_optir..., y_splineval_uv...)
        spl = Spline1D(collect(spline_x), collect(spline_y), k=3)
        axav = spl(x)
    end

    # return A(x)/A(V)
    return axav
end

# spline points
const f99_optnir_axav_x = aa_to_invum.((26500, 12200, 6000, 5470, 4670, 4110))
const f99_nir_axebv_y_params = @. (0.265, 0.829) / 3.1

# c1-c2 correlation terms
const f99_c3 = 3.23
const f99_c4 = 0.41
const f99_x0 = 4.596
const f99_gamma = 0.99

"""
    F99(;Rv=3.1)
Fitzpatrick (1999) dust law.
Returns E(B-V) in magnitudes at the given wavelength relative to the
extinction. This model applies to the UV and optical to NIR spectral range.
The default support is [1000, 33333] Å. Outside of that range this will return
0. Rv is the selective extinction and is valid over [2, 6]. A typical value for
the Milky Way is 3.1.
# References
[Fitzpatrick (1999)](https://ui.adsabs.harvard.edu/abs/1999PASP..111...63F/)
"""
@with_kw struct F99 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::F99)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return f99_invum(x, law.Rv)
end

bounds(::Type{F99}) = (1000.0, 33333.3)

"""
    DustExtinction.f99_invum(x, Rv)
The algorithm used for the [`F99`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function f99_invum(x::Real, Rv::Real)
    if !(0.3 <= x <= 10.0)
        error("out of bounds of F99, support is over $(bounds(F99)) angstrom")
    end

    # terms depending on Rv
    c2 = @evalpoly (1. / Rv) -0.824 4.717
    # original F99 c1-c2 correlation
    c1 = @evalpoly c2 2.030 -3.007

    # determine optical/IR values at spline points
    #    Final optical spline point has a leading "-1.208" in Table 4
    #    of F99, but that does not reproduce Table 3.
    #    Additional indication that this is not correct is from
    #    fm_unred.pro
    #    which is based on FMRCURVE.pro distributed by Fitzpatrick.
    #    --> confirmation needed?
    #
    #    Also, fm_unred.pro has different coeff and # of terms,
    #    but later work does not include these terms
    #    --> check with Fitzpatrick?
    opt_axebv_y =
        (@evalpoly Rv -0.426 1.0044),
        (@evalpoly Rv -0.050 1.0016),
        (@evalpoly Rv  0.701 1.0016),
        (@evalpoly Rv  1.208 1.0032 -0.00033)

    nir_axebv_y = @. f99_nir_axebv_y_params * Rv
    optnir_axebv_y = @. (nir_axebv_y..., opt_axebv_y...) / Rv

    return _curve_F99_method(
            x,
            Rv,
            c1,
            c2,
            f99_c3,
            f99_c4,
            f99_x0,
            f99_gamma,
            f99_optnir_axav_x,
            optnir_axebv_y,
        )
end

# spline points
const f04_opt_axav_x = aa_to_invum.((6000, 5470, 4670, 4110))
const f04_nir_axav_x = (0.50, 0.75, 1.0)

# **Use NIR spline x values in FM07, clipped to K band for now
const f04_optnir_axav_x = (f04_nir_axav_x..., f04_opt_axav_x...)

# c1-c2 correlation terms
const f04_c3 = 2.991
const f04_c4 = 0.319
const f04_x0 = 4.592
const f04_gamma = 0.922

"""
    F04(;Rv=3.1)
Fitzpatrick (2004) dust law.
Returns E(B-V) in magnitudes at the given wavelength relative to the
extinction. This model applies to the UV and optical to NIR spectral range.
The default support is [1000, 33333] Å. Outside of that range this will return
0. Rv is the selective extinction and is valid over [2, 6]. A typical value for
the Milky Way is 3.1.
Equivalent to the F99 model with an updated NIR Rv dependence
See also Fitzpatrick & Massa (2007, ApJ, 663, 320)
# References
[Fitzpatrick (2004)](https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F/)
"""
@with_kw struct F04 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::F04)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return f04_invum(x, law.Rv)
end

bounds(::Type{F04}) = (1000.0, 33333.3)

"""
    DustExtinction.f04_invum(x, Rv)
The algorithm used for the [`F04`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function f04_invum(x::Real, Rv::Real)
    if !(0.3 <= x <= 10.0)
        error("out of bounds of F04, support is over $(bounds(F04)) angstrom")
    end

    # original F99 Rv dependence
    c2 = @evalpoly (1. / Rv) -0.824 4.717
    # updated F04 C1-C2 correlation
    c1 = @evalpoly c2 2.18 -2.91

    # **Keep optical spline points from F99:
    #    Final optical spline point has a leading "-1.208" in Table 4
    #    of F99, but that does not reproduce Table 3.
    #    Additional indication that this is not correct is from
    #    fm_unred.pro
    #    which is based on FMRCURVE.pro distributed by Fitzpatrick.
    #    --> confirmation needed?
    opt_axebv_y =
        (@evalpoly Rv -0.426 1.0044),
        (@evalpoly Rv -0.050 1.0016),
        (@evalpoly Rv  0.701 1.0016),
        (@evalpoly Rv  1.208 1.0032 -0.00033)

    # updated NIR curve from F04, note R dependence
    # Julia v1.0 workaround: https://github.com/JuliaAstro/DustExtinction.jl/pull/31#commitcomment-40605778
    nir_axebv_y_coeff = (0.63 * Rv - 0.84)
    nir_axebv_y = @. nir_axebv_y_coeff * f04_nir_axav_x^1.84

    optnir_axebv_y = @. (nir_axebv_y..., opt_axebv_y...) / Rv

    return _curve_F99_method(
            x,
            Rv,
            c1,
            c2,
            f04_c3,
            f04_c4,
            f04_x0,
            f04_gamma,
            f04_optnir_axav_x,
            optnir_axebv_y,
        )
end


text:
if isdir(Pkg.dir("OhMyREPL"))
    @eval using OhMyREPL
else
    warn("OhMyREPL not installed")
end

max_stars_repo_name: sanga/dotfiles

max_issues_repo_name: sanga/dotfiles

max_forks_repo_name: sanga/dotfiles


text:
#!/usr/bin/env julia

using Luxor

width, height = 4000, 4000
margin = 500

fname = "/tmp/test-image.pdf"
Drawing(width, height, fname)
origin()
background("grey25")

setline(5)
sethue("green")

image = readpng(dirname(@__FILE__) * "/../docs/figures/julia-logo-mask.png")

w = image.width
h = image.height

pagetiles = Tiler(width, height, 7, 9)
tw = pagetiles.tilewidth/2
for (pos, n) in pagetiles
    circle(pos, tw, :stroke)
    circle(pos, tw, :clip)
    gsave()
    translate(pos)
    scale(.95, .95)
    rotate(rand(0.0:pi/8:2pi))
    placeimage(image, -w/2, -h/2)
    grestore()
    clipreset()
end

finish()
println("finished test: output in $(fname)")


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.

"""
    test_modification_set_function_single_variable(
        model::MOI.ModelLike,
        config::Config,
    )

Test that modifying the function of a `VariableIndex`-in-`LessThan` constraint
throws a [`SettingVariableIndexNotAllowed`](@ref) error.
"""
function test_modification_set_function_single_variable(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    x, y = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x + T(1) * y
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    err = MOI.SettingVariableIndexNotAllowed()
    @test_throws err MOI.set(model, MOI.ConstraintFunction(), c, y)
    return
end

"""
    test_modification_set_singlevariable_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test set modification VariableIndex-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_set_singlevariable_lessthan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1))],
    )
    MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(T(2)))
    @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.LessThan(T(2))
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        variable_primal = [(x, T(2))],
        constraint_primal = [(c, T(2))],
        constraint_dual = [(c, T(-1))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_set_singlevariable_lessthan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[2]),
            MOI.FEASIBLE_POINT,
        ),
    )
    return
end

"""
    test_modification_transform_singlevariable_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test set transformation of a VariableIndex-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_transform_singlevariable_lessthan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1))],
    )
    c2 = MOI.transform(model, c, MOI.GreaterThan(T(2)))
    @test !MOI.is_valid(model, c)
    @test MOI.is_valid(model, c2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        variable_primal = [(x, T(2))],
        constraint_primal = [(c2, T(2))],
        constraint_dual = [(c2, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_transform_singlevariable_lessthan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[2]),
            MOI.FEASIBLE_POINT,
        ),
    )
    return
end

"""
    test_modification_set_scalaraffine_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying set of ScalarAffineFunction-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_set_scalaraffine_lessthan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(model, f, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1))],
    )
    MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(T(2)))
    @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.LessThan(T(2))
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        variable_primal = [(x, T(2))],
        constraint_primal = [(c, T(2))],
        constraint_dual = [(c, T(-1))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_set_scalaraffine_lessthan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1],
        ),
    )
    return
end

"""
    test_modification_coef_scalaraffine_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying a variable coefficient in a ScalarAffineFunction-in-LessThan
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_coef_scalaraffine_lessthan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(model, f, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1))],
    )
    MOI.modify(model, c, MOI.ScalarCoefficientChange(x, T(2)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_coef_scalaraffine_lessthan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1//2],
        ),
    )
    return
end

"""
    test_modification_func_scalaraffine_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test setting the function in a ScalarAffineFunction-in-LessThan
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_func_scalaraffine_lessthan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(model, f, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1))],
    )
    MOI.set(
        model,
        MOI.ConstraintFunction(),
        c,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0)),
    )
    foo = MOI.get(model, MOI.ConstraintFunction(), c)
    @test foo ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0))
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        # constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_func_scalaraffine_lessthan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1//2],
        ),
    )
    return
end

"""
    test_modification_func_vectoraffine_nonneg(
        model::MOI.ModelLike,
        config::Config,
    )

Test setting the function in a VectorAffineFunction-in-Nonnegatives
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_func_vectoraffine_nonneg(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x, y = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    f = T(1) * x + T(2) * y
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    g = MOI.Utilities.operate(vcat, T, T(1) * x, T(2) * y)
    c = MOI.add_constraint(model, g, MOI.Nonnegatives(2))
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        variable_primal = [(x, T(0)), (y, T(0))],
        constraint_primal = [(c, T[0, 0])],
    )
    MOI.set(
        model,
        MOI.ConstraintFunction(),
        c,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), y)),
            ],
            T[-1, -3//2],
        ),
    )
    foo = MOI.get(model, MOI.ConstraintFunction(), c)
    @test foo ≈ MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x)),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), y)),
        ],
        T[-1, -3//2],
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(5 // 2),
        variable_primal = [(x, T(1)), (y, T(3 // 4))],
        constraint_primal = [(c, T[0, 0])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_func_vectoraffine_nonneg),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[0, 0]),
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1, 3//4]),
        ),
    )
    return
end

"""
    test_modification_const_vectoraffine_zeros(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying the constant term in a VectorAffineFunction-in-Zeros constraint.
If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_const_vectoraffine_zeros(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x, y = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x + T(2) * y
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), y)),
            ],
            T[0, 0],
        ),
        MOI.Zeros(2),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        variable_primal = [(x, T(0)), (y, T(0))],
        constraint_primal = [(c, T[0, 0])],
    )
    MOI.modify(model, c, MOI.VectorConstantChange(T[-1, -3//2]))
    _test_model_solution(
        model,
        config;
        objective_value = T(5 // 2),
        variable_primal = [(x, T(1)), (y, T(3 // 4))],
        constraint_primal = [(c, T[0, 0])],
    )
    return
end

version_added(::typeof(test_modification_const_vectoraffine_zeros)) = v"0.10.7"

function setup_test(
    ::typeof(test_modification_const_vectoraffine_zeros),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[0, 0]),
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1, 3//4]),
        ),
    )
    return
end

"""
    test_modification_const_vectoraffine_nonpos(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying the constant term in a VectorAffineFunction-in-Nonpositives
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_const_vectoraffine_nonpos(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x, y = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(1) * x + T(2) * y
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), y)),
            ],
            T[0, 0],
        ),
        MOI.Nonpositives(2),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        variable_primal = [(x, T(0)), (y, T(0))],
        constraint_primal = [(c, T[0, 0])],
    )
    MOI.modify(model, c, MOI.VectorConstantChange(T[-1, -3//2]))
    _test_model_solution(
        model,
        config;
        objective_value = T(5 // 2),
        variable_primal = [(x, T(1)), (y, T(3 // 4))],
        constraint_primal = [(c, T[0, 0])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_const_vectoraffine_nonpos),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[0, 0]),
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1, 3//4]),
        ),
    )
    return
end

"""
    test_modification_multirow_vectoraffine_nonpos(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying the variable coefficients in a
VectorAffineFunction-in-Nonpositives constraint. If `config.solve=true` confirm
that it solves correctly.
"""
function test_modification_multirow_vectoraffine_nonpos(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), x)),
            ],
            T[-1, -1],
        ),
        MOI.Nonpositives(2),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T[-1//2, 0])],
    )
    MOI.modify(model, c, MOI.MultirowChange(x, [(1, T(4)), (2, T(3))]))
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 4),
        variable_primal = [(x, T(1 // 4))],
        constraint_primal = [(c, T[0, -1//4])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_multirow_vectoraffine_nonpos),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//4]),
        ),
    )
    return
end

"""
    test_modification_const_scalar_objective(
        model::MOI.ModelLike,
        config::Config,
    )

Test the constant of a scalaraffine objective. If `config.solve=true` confirm
that it solves correctly.
"""
function test_modification_const_scalar_objective(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(2))
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(3),
        variable_primal = [(x, T(1))],
    )
    MOI.modify(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarConstantChange(T(3)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(4),
        variable_primal = [(x, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_const_scalar_objective),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
    )
    return
end

"""
    test_modification_coef_scalar_objective(
        model::MOI.ModelLike,
        config::Config
    )

Test modifying a variable coefficient in a scalaraffine objective. If
`config.solve=true` confirm that it solves correctly.
"""
function test_modification_coef_scalar_objective(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
    )
    MOI.modify(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarCoefficientChange(x, T(3)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(3),
        variable_primal = [(x, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_coef_scalar_objective),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
    )
    return
end

"""
    test_modification_delete_variable_with_single_variable_obj(
        model::MOI.ModelLike,
        config::Config,
    )

Test deleting a variable that is the objective function.
"""
function test_modification_delete_variable_with_single_variable_obj(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    c = MOI.add_constraint(model, x, MOI.GreaterThan(T(1)))
    MOI.delete(model, y)
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_modification_delete_variable_with_single_variable_obj),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.VariableIndex, MOI.GreaterThan{T}) => T[1],
        ),
    )
    return
end

"""
    test_modification_delete_variables_in_a_batch(
        model::MOI.ModelLike,
        config::Config,
    )

Test deleting many variables in a batch (i.e. using the delete method which
takes a vector of variable references). If `config.solve=true` confirm that it
solves correctly.
"""
function test_modification_delete_variables_in_a_batch(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x, y, z = MOI.add_variables(model, 3)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    f = T(1) * x + T(2) * y + T(3) * z
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.add_constraint.(model, [x, y, z], MOI.GreaterThan(T(1)))
    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test MOI.is_valid(model, z)
    if _supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), T(6), config)
    end
    MOI.delete(model, [x, z])
    @test !MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test !MOI.is_valid(model, z)
    if _supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), T(2), config)
    end
    return
end

function setup_test(
    ::typeof(test_modification_delete_variables_in_a_batch),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1, 1, 1]),
            MOI.FEASIBLE_POINT,
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
        ),
    )
    return
end

"""
    test_modification_affine_deletion_edge_cases(
        model::MOI.ModelLike,
        config::Config,
    )

Test various edge cases relating to deleting affine constraints. This requires
    + ScalarAffineFunction-in-LessThan; and
    + VectorAffineFunction-in-Nonpositives.
"""
function test_modification_affine_deletion_edge_cases(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    # helpers. The function 1.0x + 0.0
    saf = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    vaf = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x))],
        T[0],
    )
    vaf2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x))],
        T[-2],
    )
    # max x
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), saf)
    # test adding a VectorAffineFunction -in- LessThan
    c1 = MOI.add_constraint(model, vaf, MOI.Nonpositives(1))
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        constraint_primal = [(c1, T[0])],
    )
    # test adding a ScalarAffineFunction -in- LessThan
    c2 = MOI.add_constraint(model, saf, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        constraint_primal = [(c1, T[0]), (c2, T(0))],
    )
    # now delete the VectorAffineFunction
    MOI.delete(model, c1)
    @test_throws MOI.InvalidIndex{typeof(c1)} MOI.delete(model, c1)
    try
        MOI.delete(model, c1)
    catch err
        @test err.index == c1
    end
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        constraint_primal = [(c2, T(1))],
    )
    # add a different VectorAffineFunction constraint
    c3 = MOI.add_constraint(model, vaf2, MOI.Nonpositives(1))
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        constraint_primal = [(c2, T(1)), (c3, T[-1])],
    )
    # delete the ScalarAffineFunction
    MOI.delete(model, c2)
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        constraint_primal = [(c3, T[0])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_affine_deletion_edge_cases),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[2])),
    )
    return
end

"""
    test_modification_incorrect(model::MOI.ModelLike, config::Config{T}) where {T}

Test that constraint sets and functions cannot be set for the wrong type.
"""
function test_modification_incorrect(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0)),
        MOI.EqualTo(T(1)),
    )
    @test_throws(
        MOI.SetTypeMismatch{MOI.EqualTo{T},MOI.LessThan{T}},
        MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(T(1))),
    )
    @test_throws(
        MOI.FunctionTypeMismatch{MOI.ScalarAffineFunction{T},MOI.VariableIndex},
        MOI.set(model, MOI.ConstraintFunction(), c, x),
    )
    return
end

"""
    test_modification_incorrect_VariableIndex(
        model::MOI.ModelLike,
        config::Config,
    )

Test that constraint sets cannot be set for the wrong set type, and that
VariableIndex functions cannot be modified.
"""
function test_modification_incorrect_VariableIndex(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.LessThan{T})
    @requires(
        MOI.supports_constraint(model, MOI.VariableIndex, MOI.GreaterThan{T}),
    )
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    @test_throws(
        MOI.SetTypeMismatch{MOI.GreaterThan{T},MOI.LessThan{T}},
        MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(T(1))),
    )
    @test_throws(
        MOI.FunctionTypeMismatch{MOI.VariableIndex,MOI.ScalarAffineFunction{T}},
        MOI.set(model, MOI.ConstraintFunction(), c, T(1) * x),
    )
    y = MOI.add_variable(model)
    @test_throws(
        MOI.SettingVariableIndexNotAllowed(),
        MOI.set(model, MOI.ConstraintFunction(), c, y),
    )
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
function test_conditioning_actions_on_state()
    # Condition available actions on next state....
    # GridWorld(sx=2,sy=3) w reward at (2,3):
    #
    # Here's our grid, with some actions missing:
    #
    # |state (x,y)____available actions__|
    # |5 (1,3)__l,r__|6 (2,3)__u,d,l,r+REWARD__|
    # |3 (1,2)__l,u__|4 (2,2)_______l,r________|
    # |1 (1,1)__l,r__|2 (2,1)_______u,r________|
    # 7 (0,0) is absorbing state
    mdp = SpecialGridWorld(LegacyGridWorld(sx=2, sy=3, rs = [GridWorldState(2,3)], rv = [10.0]))

    solver = ValueIterationSolver(verbose=true)
    policy = solve(solver, mdp)

    println(policy.policy)

    # up: 1, down: 2, left: 3, right: 4
    correct_policy = [4,1,1,3,4,1,1] # alternative policies
    # for state 6 are possible, but since they are ordered
    # such that 1 comes first, 1 will always be the policy
    return policy.policy == correct_policy
end

@test test_conditioning_actions_on_state() == true


text:
struct TeXParseError <: Exception
    msg::String
    stack::Stack
    position::Int
    data
end

function Base.showerror(io::IO, e::TeXParseError)
    println(io, "TeXParseError: ",  e.msg)
    show_state(io, e.stack, e.position, e.data)
end

function show_state(io::IO, stack, position, data)
    # Everything is shifted by one
    position -= 1
    if position > lastindex(data)
        println(io, "after the end of the data")
        println(io, data)
    else
        # Compute the index of the char from the byte index
        position_to_id = zeros(Int, lastindex(data))
        id = 0
        for p in 1:lastindex(data)
            if isvalid(data, p)
                id += 1
            end
            position_to_id[p] = id
        end
        p = position_to_id[position]
        println(io, "at position ", position, " (string index ", p, ")")
        println(io, data)
        println(io, " " ^ (p - 1), "^")
    end

    println(io, "Stack before")
    show_stack(io, stack)
    println(io)
end

show_state(stack, position, data) = show_state(stdout, stack, position, data)

function show_stack(io, stack)
    for (k, level) in enumerate(stack)
        K = length(stack) - k + 1
        print(io, "[$K] ", level)
    end
end

show_stack(stack) = show_stack(stdout, stack)

function show_debug_info(stack, position, data, action_name)
    @info "Action " * String(action_name)
    show_state(stack, position, data)
end

# Super and subscript
super = re"\^"
super.actions[:exit] = [:end_command_builder, :setup_decorated, :begin_super]

sub = re"_"
sub.actions[:exit] = [:end_command_builder, :setup_decorated, :begin_sub]

# Groups
lbrace = re"{"
lbrace.actions[:exit] = [:end_command_builder, :begin_group]

rbrace = re"}"
rbrace.actions[:exit] = [:end_command_builder, :end_group, :end_token]

# Commands
bslash = re"\\"
bslash.actions[:exit] = [:end_command_builder, :begin_command_builder]

command_char = re"[A-Za-z]"
command_char.actions[:exit] = [:push_char, :end_token]

# Characters
space = re" "
space.actions[:exit] = [:end_command_builder]
special_char = lbrace | rbrace | bslash | super | sub | command_char | space
other_char = re"." \ special_char
other_char.actions[:exit] = [:end_command_builder, :push_char, :end_token]

mathexpr = re.rep(special_char | other_char)
mathexpr.actions[:exit] = [:end_command_builder]

machine = Automa.compile(mathexpr)

current(stack) = first(stack)
current_head(stack) = head(current(stack))
push_to_current!(stack, arg) = push!(current(stack).args, arg)

const require_token = [
    :subscript,
    :superscript,
    :left_delimiter,
    :right_delimiter
]

function requirement(stack)
    current_head(stack) == :argument_gatherer && return :argument
    current_head(stack) in require_token && return :token
    return :none
end

function end_token!(stack)
    requirement(stack) != :token && return
    token = pop!(stack)

    if current_head(stack) in [:decorated, :underover, :integral]
        decorated = pop!(stack)
        id = head(token) == :subscript ? 2 : 3

        !isnothing(decorated.args[id]) && throw(
            TeXParseError("multiple subscripts or superscripts", stack, p, data))
        decorated.args[id] = first(token.args)
        push_to_current!(stack, decorated)
    elseif head(token) == :left_delimiter
        push_to_current!(stack, token)
    elseif head(token) == :right_delimiter
        push_to_current!(stack, token)
        delimited = pop!(stack)

        # Simplify by constructing a single :delimited expr
        left = delimited.args[1]
        right = delimited.args[end]
        args = delimited.args[2:end-1]
        content = length(args) == 1 ? first(args) : TeXExpr(:group, args)
        simplified = TeXExpr(:delimited, [left.args[1], content, right.args[1]])
        push_to_current!(stack, simplified)
    end
end

# Set all command as function to play more nicely with Revise.jl

_begin_group!(stack, p, data) = push!(stack, TeXExpr(:group))

function _end_group!(stack, p, data)
    current_head(stack) != :group && throw(
        TeXParseError("Unexpected '}' at position $(p-1)", stack, p, data))
    group = pop!(stack)

    # Remplace empty groups by a zero-width space
    if isempty(group.args)
        group = TeXExpr(:space, 0.0)
    # Remove nestedness for group with a single element
    elseif length(group.args) == 1
        group = first(group.args)
    end

    if requirement(stack) == :argument
        push_to_current!(stack, group)

        command_builder = current(stack)
        head, required_n_args = command_builder.args[1:2]
        args = command_builder.args[3:end]

        if required_n_args == length(args)
            pop!(stack)
            command = TeXExpr(head, args)
            push_to_current!(stack, command)
        end
    else
        push_to_current!(stack, group)
    end
end

function _push_char!(stack, p, data)
    if current_head(stack) == :skip_char
        pop!(stack)
    elseif isvalid(data, p-1)
        char = data[prevind(data, p)]
        push_to_current!(stack, canonical_expr(char))
    end
end

function _end_token!(stack, p, data)
    if isvalid(data, p-1)
        end_token!(stack)
    end
end

_begin_sub!(stack, p, data) = push!(stack, TeXExpr(:subscript))
_begin_super!(stack, p, data) = push!(stack, TeXExpr(:superscript))

function _setup_decorated!(stack, p, data)
    core = pop!(current(stack).args)
    if head(core) ∉ [:decorated, :integral, :underover]
        push!(stack, TeXExpr(:decorated, Any[core, nothing, nothing]))
    else
        push!(stack, core)
    end
end

_begin_command_builder!(stack, p, data) = push!(stack, TeXExpr(:command_builder))

function _end_command_builder!(stack, p, data)
    if current_head(stack) == :command_builder
        command_builder = pop!(stack)

        args = command_builder.args
        skip_char = false

        # One character command are always tested even if the char is not
        # a letter
        if length(args) == 0
            current_char = data[prevind(data, p)]
            push!(args, current_char)
            skip_char = true
        end

        command_name = String(Char.(args))
        command = "\\" * command_name

        if command_name == "left"
            push!(stack, TeXExpr(:delimited))
            push!(stack, TeXExpr(:left_delimiter))
        elseif command_name == "right"
            current_head(stack) != :delimited && throw(
                TeXParseError("unexpected '\\right' at position $(p-1)",
                stack, p, data))
            push!(stack, TeXExpr(:right_delimiter))
        elseif haskey(command_to_canonical, command)
            expr = command_to_canonical[command]

            if head(expr) == :argument_gatherer
                push!(stack, expr)
            else
                push_to_current!(stack, expr)
                end_token!(stack)
            end
        else
            throw(
                TeXParseError("unsupported command $command",
                stack, p, data))
        end

        if skip_char
            push!(stack, TeXExpr(:skip_char))
        end
    end
end

action_names = [
    :begin_group, :end_group, :push_char, :end_token, :begin_sub, :begin_super,
    :setup_decorated, :begin_command_builder, :end_command_builder]

actions = map(action_names) do action_name
    function_name = Symbol("_" * String(action_name) * "!")

    return action_name => quote
        if showdebug
            show_debug_info(stack, p, data, $(QuoteNode(action_name)))
        end

        $function_name(stack, p, data)

        if showdebug
            println("Stack after")
            show_stack(stack)
            println()
        end
    end
end

actions = Dict(actions...)

context = Automa.CodeGenContext()
@eval function texparse(data ; showdebug=false)
    # Allows string to start with _ or ^
    if !isempty(data) && (data[1] == '_' || data[1] == '^')
        data = "{}" * data
    end

    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)

    # Needed to avoid problem with multi bytes unicode chars
    stack = Stack{Any}()
    push!(stack, TeXExpr(:expr))

    try
        $(Automa.generate_exec_code(context, machine, actions))
    catch
        throw(TeXParseError("unexpected error while parsing", stack, p, data))
    end

    while current_head(stack) == :skip_char
        pop!(stack)
    end

    if length(stack) > 1
        err = TeXParseError(
            "end of string reached with unfinished $(current(stack).head)",
            stack, p_eof, data)
        throw(err)
    end

    final_expr = current(stack)
    # Never return an empty expression
    if isempty(final_expr.args)
        push!(final_expr.args, TeXExpr(:space, 0.0))
    end
    return final_expr
end

@doc """
    texparse(data::String ; showdebug=false)

Parse a string representing a single LaTeX expression into nested TeXExpr.

See the documentation for the possible combinations of expression head and
arguments.

Setting `showdebug` to `true` show a very verbose break down of the parsing.
""" texparse

"""
    texparse(data::LaTeXString ; showdebug=false)

Parse a LaTeXString composed of a single LaTeX math expression into nested
TeXExpr.
"""
texparse(data::LaTeXString ; showdebug=false) = texparse(data[2:end-1] ; showdebug=showdebug)


text:
using Test
using GRASP.Symmetries

@testset "Symmetries" begin

#
# Parity
# ------------------------------------------------------------------------------------------
@test isdefined(Symmetries, :even)
@test isdefined(Symmetries, :odd)

@test Parity('+') == Symmetries.even
@test Parity('-') == Symmetries.odd
@test_throws ArgumentError Parity(' ')

@test Parity(1) == Symmetries.even
@test Parity(-1) == Symmetries.odd
@test_throws ArgumentError Parity(0)
@test_throws ArgumentError Parity(-2)
@test_throws ArgumentError Parity(42)

@test string(Symmetries.even) == "+"
@test string(Symmetries.odd)  == "-"

@test parse(Parity, "+")   == Symmetries.even
@test parse(Parity, "-")   == Symmetries.odd
@test parse(Parity, " + ") == Symmetries.even
@test parse(Parity, "\t+") == Symmetries.even

@test_throws Meta.ParseError parse(Parity, "")
@test_throws Meta.ParseError parse(Parity, " ")
@test_throws Meta.ParseError parse(Parity, "++")
@test_throws Meta.ParseError parse(Parity, "α")

#
# Angular momentum
# ------------------------------------------------------------------------------------------
@test_throws ArgumentError AngularMomentum(-2)
@test_throws ArgumentError AngularMomentum(-1//2)

@test string(AngularMomentum(0))     == "0"
@test string(AngularMomentum(1//2))  == "1/2"
@test string(AngularMomentum(2//2))  == "1"
@test string(AngularMomentum(3//2))  == "3/2"
@test string(AngularMomentum(5))     == "5"
@test string(AngularMomentum(21//2)) == "21/2"
@test string(AngularMomentum(40//4)) == "10"

@test parse(AngularMomentum, "")     == AngularMomentum(0)
@test parse(AngularMomentum, "  ")   == AngularMomentum(0)
@test parse(AngularMomentum, "0")    == AngularMomentum(0)
@test parse(AngularMomentum, "1")    == AngularMomentum(1)
@test parse(AngularMomentum, "22")   == AngularMomentum(22)
@test parse(AngularMomentum, "1/2")  == AngularMomentum(1//2)
@test parse(AngularMomentum, "11/2") == AngularMomentum(11//2)
@test parse(AngularMomentum, "10/2") == AngularMomentum(10//2)

@test convert(Rational, AngularMomentum(0)) == 0//1
@test convert(Rational, AngularMomentum(0)) == 0
@test convert(Rational, AngularMomentum(22)) == 22
@test convert(Rational, AngularMomentum(22)) == 22//1
@test convert(Rational, AngularMomentum(1//2)) == 1//2

@test AngularMomentum(0) + AngularMomentum(1//2) == AngularMomentum(1//2)
@test AngularMomentum(1//2) + AngularMomentum(0) == AngularMomentum(1//2)
@test AngularMomentum(1//2) + AngularMomentum(1//2) == AngularMomentum(1)
@test AngularMomentum(1) + AngularMomentum(1//2) == AngularMomentum(3//2)
@test AngularMomentum(9//2) + AngularMomentum(10) == AngularMomentum(29//2)

@test absdiff(AngularMomentum(0), AngularMomentum(0)) == AngularMomentum(0)
@test absdiff(AngularMomentum(0), AngularMomentum(5)) == AngularMomentum(5)
@test absdiff(AngularMomentum(5), AngularMomentum(0)) == AngularMomentum(5)
@test absdiff(AngularMomentum(1), AngularMomentum(1)) == AngularMomentum(0)
@test absdiff(AngularMomentum(3//2), AngularMomentum(3//2)) == AngularMomentum(0)
@test absdiff(AngularMomentum(1//2), AngularMomentum(5)) == AngularMomentum(9//2)

@test AngularMomentum(0) < AngularMomentum(1//2)
@test AngularMomentum(0) <= AngularMomentum(0)
@test AngularMomentum(3//2) <= AngularMomentum(3//2)
@test AngularMomentum(3//2) < AngularMomentum(2)
@test AngularMomentum(3//2) <= AngularMomentum(2)
@test !(AngularMomentum(3//2) < AngularMomentum(3//2))
@test !(AngularMomentum(4//2) <= AngularMomentum(3//2))
@test AngularMomentum(4//2) > AngularMomentum(3//2)

#
# Angular symmetries
# ------------------------------------------------------------------------------------------
@test string(AngularSymmetry(AngularMomentum(0), Symmetries.even))   == "0+"
@test string(AngularSymmetry(AngularMomentum(1//2), Symmetries.odd)) == "1/2-"

@test string(AngularSymmetry(0, Symmetries.even))   == "0+"
@test string(AngularSymmetry(1//2, Symmetries.odd)) == "1/2-"

@test string(AngularSymmetry(AngularMomentum(0), '-'))    == "0-"
@test string(AngularSymmetry(AngularMomentum(1//2), '+')) == "1/2+"

@test string(AngularSymmetry(0, true))     == "0+"
@test string(AngularSymmetry(1//2, false)) == "1/2-"
@test string(AngularSymmetry(1, '-'))      == "1-"
@test string(AngularSymmetry(3//2, '+'))   == "3/2+"

end # @testset "Symmetries"


text:
# ------------------------------------------------------------------------------------------
# # Paquetes
#
# Julia tiene màs de 1686 paquetes registrados, conformando una gran parte del ecosistema de
# Julia.
#
# Para ver todos los paquetes, visita
#
# https://pkg.julialang.org/
#
# o bien
#
# https://juliaobserver.com/
#
# Ahora vamos a aprender a usarlos
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# La primera vez que quieres usar un paquete en Julia, hay que agregarlo
# ------------------------------------------------------------------------------------------

#Pkg.add("Example")

# ------------------------------------------------------------------------------------------
# Cada vez que usas Julia (empezar una nueva sesión en el REPL, abrir un notebook por
# primera vez, por ejemplo), tienes que cargar el paquete usando la palabra reservada
# `using`
# ------------------------------------------------------------------------------------------

using Example

# ------------------------------------------------------------------------------------------
# En el código fuente de  `Example.jl` en
#
# https://github.com/JuliaLang/Example.jl/blob/master/src/Example.jl
#
# Vemos una función declarada como
#
# ```
# hello(who::String) = "Hello, $who"
# ```
#
# Si cargamos `Example`, debemos poder llamar `hello`
# ------------------------------------------------------------------------------------------

hello("it's me. I was wondering if after all these years you'd like to meet.")

# ------------------------------------------------------------------------------------------
# Ahora vamos a jugar con el paquete de Colors
# ------------------------------------------------------------------------------------------

#Pkg.add("Colors")

using Colors

# ------------------------------------------------------------------------------------------
# Creemos una bandeja de 100 colores
# ------------------------------------------------------------------------------------------

bandeja = distinguishable_colors(100)

# ------------------------------------------------------------------------------------------
# y podemos crear una matriz colorida aleatoriamente con rand
# ------------------------------------------------------------------------------------------

rand(bandeja, 3, 3)

# ------------------------------------------------------------------------------------------
# En el próximo notebook, vamos a usar un nuevo paquete para graficar datos
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
# 7.1 Usa el paquete de  (código fuente en https://github.com/JuliaMath/Primes.jl) para
# encontrar el número primer más grande menor a 1,000,000
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module WiringDiagrams

using Reexport

include("Directed.jl")
include("CPortGraphs.jl")
include("Undirected.jl")
include("MonoidalDirected.jl")
include("MonoidalUndirected.jl")
include("Algebras.jl")
include("Algorithms.jl")
include("Expressions.jl")
include("ScheduleUndirected.jl")

@reexport using .DirectedWiringDiagrams
@reexport using .CPortGraphs
@reexport using .UndirectedWiringDiagrams
@reexport using .MonoidalDirectedWiringDiagrams
@reexport using .MonoidalUndirectedWiringDiagrams
@reexport using .WiringDiagramAlgebras
@reexport using .WiringDiagramAlgorithms
@reexport using .WiringDiagramExpressions
@reexport using .ScheduleUndirectedWiringDiagrams

include("Serialization.jl")
include("GraphML.jl")
include("JSON.jl")

using .WiringDiagramSerialization
export convert_from_graph_data, convert_to_graph_data
@reexport using .GraphMLWiringDiagrams
@reexport using .JSONWiringDiagrams

end


text:
"""
# module kernels

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

functions to create LAMINART kernels for lgn -> L6, L4 and W kernel
"""

module LaminartKernels

export kern_d, kern_b

# lgn to l6/l4
using NNlib, ImageFiltering, Images, OffsetArrays

function kern_d_ph(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l >> 1
    g = σ == 0 ? [1f0] : [exp(-x * cos(θ) / (2f0 * σ)) for x = -w:w]
    centered(g / sum(g))
end

function kern_d_pv(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l >> 1
    g = σ == 0 ? [1f0] : [exp(-x * sin(θ) / (2f0 * σ)) for x = -w:w]
    centered(g / sum(g))
end

function kern_d_mh(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l >> 1
    g = σ == 0 ? [1f0] : [exp(x * cos(θ) / (2f0 * σ)) for x = -w:w]
    centered(g / sum(g))
end

function kern_d_mv(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l >> 1
    g = σ == 0 ? [1f0] : [exp(x * sin(θ) / (2f0 * σ)) for x = -w:w]
    centered(g / sum(g))
end

function kern_d_p(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    kern_d_pv(σ, θ, l) .* transpose(kern_d_ph(σ, θ, l))
end

function kern_d_m(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    kern_d_mv(σ, θ, l) .* transpose(kern_d_mh(σ, θ, l))
end

function kern_d(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    kern_d_p(σ, θ, l) .- kern_d_m(σ, θ, l)
end

# function kern_d(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
#     relu.(kern_d(σ, θ, l)) .- relu.(-1 .*(kern_d(σ, θ, l)))
# end

function kern_b(σ::Real, θ::Real, l = 4 * ceil(Int, σ) + 1)
    relu.(kern_d(σ, θ, l)) .+ relu.(-1f0 .* (kern_d(σ, θ, l)))
end

"""
used to make W kernels
returns 2D kernel with elongated Gaussian rotated
"""
function gaussian_rot(
    σ_x::Real,
    σ_y::Real,
    θ::Real,
    l = 4 * ceil(Int, max(σ_a, σ_b)) + 1,
)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l >> 1
    g = OffsetArray(fill(0.0f0, l, l), -w:w, -w:w)
    #todo add when σ_x or/and σ_y == 0
    for x ∈ -w:w, y ∈ -w:w
        g[x, y] = exp(
            -1f0 / 2f0 * (
                (((x * cos(θ) - y * sin(θ)) / σ_x)^2f0) +
                (((x * sin(θ) + y * cos(θ)) / σ_y)^2f0)
            ),
        )
    end
    centered(g / sum(g))
end

# Rotating function
function fun_R(x, y, θ)
    x * cos(θ) - y * sin(θ), x * sin(θ) + y * cos(θ)
end

end


text:
"""
    PersistenceDiagram <: AbstractVector{PersistenceInterval}

Type for representing persistence diagrams. Behaves exactly like a vector of
`PersistenceInterval`s, but can have additional metadata attached to it. It supports pretty
printing and plotting.

Can be used as a table with any function that uses the
[`Tables.jl`](https://github.com/JuliaData/Tables.jl) interface. Note that using it as a
table will only keep interval endpoints and the `dim` and `threshold` attributes.

# Example

```jldoctest
julia> diagram = PersistenceDiagram([(1, 3), (3, 4), (1, Inf)]; dim=1, custom_metadata=:a)
3-element 1-dimensional PersistenceDiagram:
 [1.0, 3.0)
 [3.0, 4.0)
 [1.0, ∞)

julia> diagram[1]
[1.0, 3.0)

julia> sort(diagram; by=persistence, rev=true)
3-element 1-dimensional PersistenceDiagram:
 [1.0, ∞)
 [1.0, 3.0)
 [3.0, 4.0)

julia> propertynames(diagram)
(:dim, :custom_metadata)

julia> dim(diagram)
1

julia> diagram.custom_metadata
:a
```
"""
struct PersistenceDiagram <: AbstractVector{PersistenceInterval}
    intervals::Vector{PersistenceInterval}
    meta::NamedTuple
end

function PersistenceDiagram(intervals::Vector{PersistenceInterval}; kwargs...)
    meta = (; kwargs...)
    return PersistenceDiagram(intervals, meta)
end
function PersistenceDiagram(intervals::AbstractVector{PersistenceInterval}; kwargs...)
    return PersistenceDiagram(collect(intervals); kwargs...)
end
function PersistenceDiagram(pairs::AbstractVector{<:Tuple}; kwargs...)
    return PersistenceDiagram(PersistenceInterval.(pairs); kwargs...)
end
function PersistenceDiagram(table)
    rows = Tables.rows(table)
    if isempty(rows)
        return PersistenceDiagram(PersistenceInterval[])
    else
        firstrow = first(rows)
        dim = hasproperty(firstrow, :dim) ? firstrow.dim : missing
        threshold = hasproperty(firstrow, :threshold) ? firstrow.threshold : missing
        intervals = map(rows) do row
            d = hasproperty(row, :dim) ? row.dim : missing
            t = hasproperty(row, :threshold) ? row.threshold : missing
            if !isequal(d, dim)
                error("different `dim`s detected. Try splitting the table first.")
            end
            if !isequal(t, threshold)
                error("different `threshold`s detected. Try splitting the table first.")
            end
            PersistenceInterval(row.birth, row.death)
        end
        return PersistenceDiagram(intervals; dim=dim, threshold=threshold)
    end
end

function Base.show(io::IO, diag::PersistenceDiagram)
    return summary(io, diag)
end
function Base.summary(io::IO, diag::PersistenceDiagram)
    if haskey(diag.meta, :dim)
        print(io, length(diag), "-element ", dim(diag), "-dimensional PersistenceDiagram")
    else
        print(io, length(diag), "-element PersistenceDiagram")
    end
end

###
### Array interface
###
Base.size(diag::PersistenceDiagram) = size(diag.intervals)
Base.getindex(diag::PersistenceDiagram, i::Integer) = diag.intervals[i]
Base.setindex!(diag::PersistenceDiagram, x, i::Integer) = diag.intervals[i] = x
Base.firstindex(diag::PersistenceDiagram) = 1
Base.lastindex(diag::PersistenceDiagram) = length(diag.intervals)

function Base.similar(diag::PersistenceDiagram)
    return PersistenceDiagram(similar(diag.intervals); diag.meta...)
end
function Base.similar(diag::PersistenceDiagram, dims::Tuple)
    return PersistenceDiagram(similar(diag.intervals, dims); diag.meta...)
end

###
### Meta
###
function Base.getproperty(diag::PersistenceDiagram, key::Symbol)
    if hasfield(typeof(diag), key)
        return getfield(diag, key)
    elseif haskey(diag.meta, key)
        return diag.meta[key]
    else
        error("$diag has no $key")
    end
end
function Base.propertynames(diag::PersistenceDiagram, private::Bool=false)
    if private
        return tuple(:intervals, :meta, propertynames(diag.meta)...)
    else
        return tuple(:intervals, propertynames(diag.meta)...)
    end
end

"""
    threshold(diagram::PersistenceDiagram)

Get the threshold of persistence diagram. Equivalent to `diagram.threshold`.
"""
threshold(diag::PersistenceDiagram) = diag.threshold

"""
    dim(diagram::PersistenceDiagram)

Get the dimension of persistence diagram. Equivalent to `diagram.dim`.
"""
dim(diag::PersistenceDiagram) = diag.dim


text:
plotvector(x, str="k-",label="", alpha=1.0) = plot3D([0;x[1]],[0;x[2]],[0;x[3]], str,alpha=alpha,label=label)

proj(v, x) = vec( dot(v, x) / dot(v, v) * v )

unitvector(x) = x / norm(x)


text:
module Old_Impl

export MultiheadAttention

const Abstract3DTensor{T} = AbstractArray{T, 3}

include("batched_tril.jl")
include("batchedmul.jl")

# an multi-head attention used in Transformers.jl v0.1.13
include("./mh_atten.jl")
include("./getmask.jl")

end


text:
module ScientificTypesBase

# Type exports
export Convention

# re-export-able types and methods
export Scientific, Found, Unknown, Known, Finite, Infinite,
    OrderedFactor, Multiclass, Count, Continuous,
    ScientificTimeType, ScientificDate, ScientificDateTime,
    ScientificTime,
    Textual, Binary,
    ColorImage, GrayImage, Image, Table,
    Density, Sampleable,
    ManifoldPoint,
    Annotated, AnnotationOf, Multiset, Iterator
export scitype, scitype_union, elscitype, nonmissing, trait

# utils (should not be re-exported)
export TRAIT_FUNCTION_GIVEN_NAME, set_convention

# -------------------------------------------------------------------
# Scientific Types

abstract type Found          end
abstract type Known <: Found end
abstract type Unknown <: Found end

abstract type Annotated{S} <: Known end
abstract type AnnotationOf{S} <: Known end
abstract type Multiset{S} <: Known end

# for iterators over objects with scitype Ω that do not have some
# AbstractVector scitype:
abstract type Iterator{Ω} end

abstract type           Infinite <: Known end
abstract type          Finite{N} <: Known end
abstract type         Image{W,H} <: Known end
abstract type ScientificTimeType <: Known end
abstract type            Textual <: Known end
abstract type           Table{K} <: Known end

abstract type Continuous <: Infinite end
abstract type  Count <: Infinite end

abstract type    Multiclass{N} <: Finite{N} end
abstract type OrderedFactor{N} <: Finite{N} end

abstract type ScientificDate <: ScientificTimeType end
abstract type ScientificTime <: ScientificTimeType end
abstract type ScientificDateTime <: ScientificTimeType end

abstract type  GrayImage{W,H} <: Image{W,H} end
abstract type ColorImage{W,H} <: Image{W,H} end

# when sampled, objects with these scitypes return objects of scitype Ω:
abstract type Sampleable{Ω} end
abstract type Density{Ω} <: Sampleable{Ω} end

abstract type ManifoldPoint{M} <: Known end

# aliases:
const Binary     = Finite{2}
const Scientific = Union{Missing,Found} # deprecated (no longer publicized)

# for internal use
const Arr = AbstractArray

# -------------------------------------------------------------------
# Convention

abstract type Convention end
struct NoConvention <: Convention end

const CONVENTION = Ref{Convention}(NoConvention())

"""
    set_convention(C)

Set the current convention to `C`.
"""
set_convention(C::Convention) = (CONVENTION[] = C; nothing)

"""
    convention()

Return the current convention.
"""
function convention()::Convention
    C = CONVENTION[]
    if C isa NoConvention
        @warn "No convention specified. Did you forget to use the " *
              "`set_convention` function?"
    end
    return C
end

# -------------------------------------------------------------------
# trait

const TRAIT_FUNCTION_GIVEN_NAME = Dict{Symbol,Function}()

"""
    trait(X)

Check `X` against traits specified in `TRAIT_FUNCTION_GIVEN_NAME` and
return a symbol corresponding to the first matching trait, or `:other`
if `X` didn't match any of the trait functions. If more than one trait
matches, throw an error.

"""
function trait(X)::Symbol
    ret = :other
    found = false
    for (name, f) in TRAIT_FUNCTION_GIVEN_NAME
        if f(X)
            found && error("Bad dictionary "*
                           "`ScientificTypesBase.TRAIT_FUNCTION_GIVEN_NAME`. ")
            ret = name
            found = true
        end
    end
    return ret
end

# -----------------------------------------------------------------
# nonmissing

if VERSION < v"1.3"
    # see also discourse.julialang.org/t/get-non-missing-type-in-the-case-of-parametric-type/29109
    """
        nonmissingtype(TT)

    Return the type `T` if `TT = Union{Missing,T}` for some `T` and return `TT`
    otherwise.
    """
    function nonmissingtype(::Type{T}) where T
        return T isa Union ? ifelse(T.a == Missing, T.b, T.a) : T
    end
end
nonmissing = nonmissingtype

# -----------------------------------------------------------------
# Constructor for table scientific type

"""
Table(...)

Constructor for the `Table` scientific type with:

```
Table(S1, S2, ..., Sn) <: Table
```

where  `S1, ..., Sn` are the scientific type of the table's columns which
are expected to be represented by abstract vectors.
"""
Table(Ts::Type...) = Table{<:Union{(Arr{<:T,1} for T in Ts)...}}

# -----------------------------------------------------------------
# scitype

include("scitype.jl")

end # module


text:
function train_model(data, pools, model_type, model_C, model_gamma)
    model = instantiate(eval(Symbol(model_type)), data, pools, Dict{Symbol, Any}())
    initialize!(model, FixedParameterInitialization(MLKernels.GaussianKernel(model_gamma), model_C))
    set_adjust_K!(model, true)
    status = fit!(model, DEFAULT_SOLVER)
    if status ∉ [JuMP.MathOptInterface.OPTIMAL, JuMP.MathOptInterface.LOCALLY_SOLVED, JuMP.MathOptInterface.ALMOST_LOCALLY_SOLVED]
        throw(ArgumentError("Cannot fit model due to solver error: '$(status)')."))
    end
    return model
end

predict_data(model, data) = SVDD.classify.(SVDD.predict(model, data))
score_data(model, data) = SVDD.predict(model, data)

function run_query_strategy(data, pools, model, qs_type, query_history)
    qs = initialize_qs(eval(Symbol(qs_type)), model, data, Dict{Symbol, Any}())
    pool_map = labelmap(pools)
    haskey(pool_map, :U) || throw(ArgumentError("No more points that are unlabeled."))
    scores = qs_score(qs, data, pool_map)
    @assert length(scores) == size(data, 2)
    candidates = [i for i in pool_map[:U] if i ∉ query_history]
    return [argmax(scores[candidates])]
end

function run_ocal(r)
    data, pools, params = r["data"], r["labels"], r["params"]

    # train global model
    model_global = train_model(data, pools, params["classifier"], params["C"], params["gamma"])
    prediction_global = predict_data(model_global, data)

    # execute query strategy
    query_ids = run_query_strategy(data, pools, model_global, params["query_strategy"], vcat(r["query_history"]...))

    # train subspace models
    prediction_subspaces = []
    score_subspace_grids = []
    for (k, s) in enumerate(r["subspaces"])
        model_subspace = train_model(data[s, :], pools, params["classifier"], params["C"], params["gamma"])
        push!(prediction_subspaces, predict_data(model_subspace, data[s, :]))
        push!(score_subspace_grids, score_data(model_subspace, r["subspace_grids"][k]))
    end

    # rank subspaces
    ranking_subspaces = collect(1:length(r["subspaces"]))

    return Dict(
        :prediction_global => prediction_global,
        :prediction_subspaces => prediction_subspaces,
        :score_subspace_grids => score_subspace_grids,
        :query_ids => query_ids,
        :ranking_subspaces => ranking_subspaces)
end