text:
using BandedMatrices, LinearAlgebra, Test, Random
import BandedMatrices: banded_qr!
Random.seed!(0)

@testset "QR tests" begin
    for T in (Float64,ComplexF64,Float32,ComplexF32)
        A=brand(T,10,10,3,2)
        Q,R=qr(A)
        @test Matrix(Q)*Matrix(R) ≈ A
        b=rand(T,10)
        @test mul!(similar(b),Q,mul!(similar(b),Q',b)) ≈ b
        for j=1:size(A,2)
            @test Q' * A[:,j] ≈ R[:,j]
        end
        A=brand(T,14,10,3,2)

        Q,R=qr(A)
        @test Matrix(Q)*Matrix(R) ≈ A

        for k=1:size(A,1),j=1:size(A,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        A=brand(T,10,14,3,2)
        Q,R=qr(A)
        @test Matrix(Q)*Matrix(R) ≈ A

        for k=1:size(Q,1),j=1:size(Q,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        A=brand(T,100,100,3,4)
        @test qr(A).factors ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).factors
        @test qr(A).τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
        b=rand(T,100)
        @test qr(A)\b ≈ Matrix(A)\b
        b=rand(T,100,2)
        @test qr(A)\b ≈ Matrix(A)\b
        @test_throws DimensionMismatch qr(A) \ randn(3)
        @test_throws DimensionMismatch qr(A).Q'randn(3)

        A=brand(T,102,100,3,4)
        @test qr(A).factors ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).factors
        @test qr(A).τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
        b=rand(T,102)
        @test qr(A)\b ≈ Matrix(A)\b
        b=rand(T,102,2)
        @test qr(A)\b ≈ Matrix(A)\b
        @test_throws DimensionMismatch qr(A) \ randn(3)
        @test_throws DimensionMismatch qr(A).Q'randn(3)

        A=brand(T,100,102,3,4)
        @test qr(A).factors ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).factors
        @test qr(A).τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
        b=rand(T,100)
        @test_broken qr(A)\b ≈ Matrix(A)\b

        A = Tridiagonal(randn(T,99), randn(T,100), randn(T,99))
        @test qr(A).factors ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).factors
        @test qr(A).τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
        b=rand(T,100)
        @test qr(A)\b ≈ Matrix(A)\b
        b=rand(T,100,2)
        @test qr(A)\b ≈ Matrix(A)\b
        @test_throws DimensionMismatch qr(A) \ randn(3)
        @test_throws DimensionMismatch qr(A).Q'randn(3)
    end

    @testset "Mixed types" begin
        A=brand(10,10,3,2)
        b=rand(ComplexF64,10)
        Q,R=qr(A)
        @test R\(Q'*b) ≈ qr(A)\b ≈ Matrix(A)\b

        A=brand(ComplexF64,10,10,3,2)
        b=rand(10)
        Q,R=qr(A)
        @test R\(Q'*b) ≈ qr(A)\b ≈ Matrix(A)\b

        A = BandedMatrix{Int}(undef, (2,1), (4,4))
        A.data .= 1:length(A.data)
        Q, R = qr(A)
        @test Q*R ≈ A
    end

    @testset "banded Views" begin
        # this test adaptive QR essentially as implemented in InfiniteLinearAlgebra.jl
        A = brand(100,100,1,1)
        V = view(copy(A),1:5,1:6)
        @test qr(V) isa QR{Float64,<:BandedMatrix{Float64}}
        @test qr(V).R ≈ qr(Matrix(V)).R
        @test qr(V).τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(V)).τ
        @test qr(V).Q ≈ qr(Matrix(V)).Q
        @test Matrix(qr(V)) ≈ V
        B = BandedMatrix(A,(1,2)) # pad
        V = view(copy(B),1:5,1:6)
        @test qr!(V) isa QR{Float64,<:SubArray{Float64,2,<:BandedMatrix{Float64}}}
        τ = Array{Float64}(undef,100)
        B = BandedMatrix(A,(1,2)) # pad
        V = view(B,1:6,1:5)
        F1 = qr(V)
        F2 = banded_qr!(V, view(τ,1:size(V,2)))
        @test F1.factors ≈ F2.factors ≈ B[1:6,1:5]
        @test F1.τ ≈ F2.τ
        F = qr(A)
        @test τ[1:size(V,2)]  ≈ F.τ[1:5]

        lmul!(F1.Q', view(B,1:6,6:7))
        @test B[1:5,6:7] ≈ F.factors[1:5,6:7]
        banded_qr!(view(B,6:100,6:100), view(τ,6:100))
        τ[end] = 0
        @test B ≈ F.factors
        @test τ ≈ F.τ
    end
end


text:
using Dates
#using ProfileView

include("../btg.jl")
include("../datasets/load_abalone.jl")

#ind = 350:370
ind = 3300:3359
#posx = 1:3 #
posx = [1;4]
posc = 1:1
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean
Fx = data[ind, posc]
y = float(target[ind])
pind = 10:10 #prediction index
#pind = i:i
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
#prediction data
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))
#rangeθ = [100.0 200]
rangeθ = [100.0 200]
rangeλ = [-1.0 1] #we will always used 1 range scale for lambda
btg1 = btg(trainingData1, rangeθ, rangeλ; quadtype = ["MonteCarlo", "MonteCarlo"])
#θ1 = btg1.nodesWeightsθ.nodes[1, 6] #pick some theta value which doubles as quadrature node
#λ1 = btg1.nodesWeightsλ.nodes[3]
##################################################################################################
PyPlot.close("all") #close existing windows

#(pdf, cdf, dpdf) = solve(btg1); #initialize training_buffer_dicts, solve once so can use fast techiques to extrapolate submatrix determinants, etc.
#pdf(x0, Fx0, .4);
#(pdf, cdf, dpdf) = solve(btg1; validate = i) #do LOOCV this time around
#a1 = y-> pdf(x0, Fx0, y)

"""
leave-one-out-training-data
"""
function lootd(td::AbstractTrainingData, i::Int64)
    x = getPosition(td)
    Fx = getCovariates(td)
    z = getLabel(td)
    x_minus_i = x[[1:i-1;i+1:end], :]
    Fx_minus_i = Fx[[1:i-1;i+1:end], :]
    z_minus_i = z[[1:i-1;i+1:end]]
    x_i = x[i:i, :]
    Fx_i = Fx[i:i, :]
    z_i = z[i:i, :]
    return trainingData(x_minus_i, Fx_minus_i, z_minus_i), x_i, Fx_i, z_i
end
"""
Run
"""
function test_LOOCV(btg1::btg, m::Int64, n::Int64; fast)
    fast = fast
    m = m; n=n

    #plt, axs = PyPlot.subplots(m, n,figsize=[20,11])
    plt, axs = PyPlot.subplots(m, n)

    #figure(1)
    before = Dates.now()
    println("LOOCV type is fast?: $fast")
    for j = 8:8
        ind1 = Int64(ceil(j/n))
        ind2 = Int64(j - n*(floor((j-.1)/n)))
        println("iteration $j")
        if j <m*n
            if fast
                (pdf, cdf, dpdf) = solve(btg1, validate = j);
                b1 = y0 -> pdf(x0, Fx0, y0);
                c1 = y0 -> cdf(x0, Fx0, y0);
                (x, y) = plt_data(b1, .01, 1.2, 100)
                (xc, yc) = plt_data(c1, .01, 1.2, 100)
            else
                (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, j)
                btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = ["MonteCarlo", "MonteCarlo"])
                (pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg2)
                b1 = y -> pdf_minus_i(x_i, Fx_i, y)
                c1 = y -> cdf_minus_i(x_i, Fx_i, y)
                (x, y) = plt_data(b1, .01, 1.2, 100)
                (xc, yc) = plt_data(c1, .01, 1.2, 100)
            end
            axs[ind1, ind2].plot(x, y, color = "red", linewidth = 2.0, linestyle = ":")
            axs[ind1, ind2].plot(xc, yc, color = "orange", linewidth = 2.0, linestyle = ":")
            #axs[ind1, ind2].plot(x1, y1, color = "blue", linewidth = 1.0, linestyle = "--")
            axs[ind1, ind2].axvline(x =  getLabel(btg1.trainingData)[j])
        else #do annotation in last box
            #println("annotating...")
            after = Dates.now()
            elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)
            qtype =  btg1.quadType
            axs[ind1, ind2].annotate(
            L"$\theta: $" * "$rangeθ" * "\n" *
            L"$\lambda: $" * "$rangeλ" * "\n" *
            "$qtype \n" *
            "ind: " * "$ind \n" *
            "x: " * "$posx " *", cov: " * "$posc\n" *
            "time: " * "$elapsedmin",
            xy=[.03;.03],
            fontsize = 5.0)
            # *
            #L"$\lambda: $" * "$rangeλ" * ".\n" *
            #"quad: " * "$btg1.quadType",
        end
    end
    #for ax in axs
    #    ax.set(xlabel="x-label", ylabel="y-label")
    #end
    # Hide x labels and tick labels for top plots and y ticks for right plots.
    for ax in axs
        ax.label_outer()
    end
end

test_LOOCV(btg1, 6, 10; fast = true)

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using DocumenterCitations
using Test

@testset "tex2unicode" begin
    @test DocumenterCitations.tex2unicode("-- ---") == "– —"
    @test DocumenterCitations.tex2unicode(
        raw"\`{o}\'{o}\^{o}\~{o}\={o}\u{o}\.{o}\\\"{o}\r{a}\H{o}\v{s}\d{u}\c{c}\k{a}\b{b}\~{a}") == "òóôõōŏȯöåőšụçąḇã"
    @test DocumenterCitations.tex2unicode(
        raw"\i{}\o{}\O{}\l{}\L{}\i\o\O\l\L") == "ıøØłŁıøØłŁ"
    @test DocumenterCitations.tex2unicode(
        raw"\t{oo}{testText}\t{az}") == "o͡otestTexta͡z"
    @test DocumenterCitations.tex2unicode(
        raw"{\o}verline") == "øverline"
    @test DocumenterCitations.tex2unicode(
        raw"\overline") == "\\overline"
end


text:
"""
   Delegates loading functions & traits for a generic theoretical spectrograph

Author: Eric Ford and collaborators
Created: August 2020
"""

"""
Module providing types and traits and customized functions for a generic theoretical spectrograph.
"""
module TheoreticalInstrument
using ..RvSpectMLBase
using ..RvSpectMLBase.InstrumentsCommon
#import ..RvSpectML: AbstractInstrument, AbstractInstrument1D, AbstractInstrument2D
using DataFrames # , FITSIO

default_theoretical_instrument_resolution = 137500
default_theoretical_instrument_λ_min = 4700.0
default_theoretical_instrument_λ_max = 6000.0

""" Trait for 1D spectra from a theoretical instrument """
struct TheoreticalInstrument1D <: AbstractInstrument1D
    λ_min::Float64
    λ_max::Float64
    resolution::Float64

    function TheoreticalInstrument1D(λl::Real, λh::Real, res::Real)
        @assert 3000 <= λl < λh
        @assert λl < λh <= 10000
        @assert 10000 <= res <= 1000000
        new(λl,λh,res)
    end
end

function TheoreticalInstrument1D(;λ_min::Real = default_theoretical_instrument_λ_min, λ_max::Real = default_theoretical_instrument_λ_max, resolution::Real = default_theoretical_instrument_resolution )
    TheoreticalInstrument1D(λ_min,λ_max,resolution)
end

λ_min(x::TheoreticalInstrument1D) = x.λ_min
λ_max(x::TheoreticalInstrument1D) = x.λ_max

""" Trait for a 2D spectra from a theoretical instrument """
struct TheoreticalInstrument2D <: AbstractInstrument2D
    λ_min::AbstractArray{Float64,1}
    λ_max::AbstractArray{Float64,1}
    pixels_per_order::Int64

    function TheoreticalInstrument2D(λ_in::AbstractArray{Float64,2}; pixels_per_order::Integer = size(λ_in,1))
        @assert size(λ_in,1) >= 6
        @assert size(λ_in,2) >= 1
        @assert 1 <= pixels_per_order <= 16386

        lambda_min = vec(minimum(λ_in,dims=1))
        lambda_max = vec(maximum(λ_in,dims=1))
        new(lambda_min, lambda_max, pixels_per_order)
    end
end

# Trait for any spectra from TheoreticalInstrument (could improve by using SimpleTraits)
const AnyTheoreticalInstrument = Union{TheoreticalInstrument1D,TheoreticalInstrument2D}
export TheoreticalInstrument, TheoreticalInstrument1D, TheoreticalInstrument2D, AnyTheoreticalInstrument

λ_min(x::TheoreticalInstrument2D) = minimum(x.λ_min)
λ_max(x::TheoreticalInstrument2D) = maximum(x.λ_max)

include("traits.jl")
export min_pixel, max_pixel, min_order, max_order, min_pixel_in_order, max_pixel_in_order
export orders_to_use_default, min_col_default, max_col_default
export metadata_symbols_default, metadata_strings_default
export default_ccf_mask_v_width

# Need to think about what to do for IO, if anything.
#include("io.jl")
#export make_manifest, read_metadata, read_data

include("util.jl")

include("gen_data.jl")
export generate_spectrum, generate_spectra_timeseries

end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
using LightGraphs
using MetaGraphs
using JuMP
using Cbc

#=
optimization.jl:

models(using JuMP) some variations around the vanilla LPBasic optimization problem.

=#

export ConstructModel

"""
    function ConstructModel(instance::ProblemInstance; <keyword arguments>)::Model

Construct a JuMP model from the LPBasic problem instance and return it.

# Arguments
- `useSlack::Bool=false`: whether to use slackened constraints. If used, slack_penalty must be defined
- `slack_penalty::Float64`: the penalty taken to the objective-value for each unit of slack taken
- `useBinaryVariant::Bool=false`: whether to use binary decision variables
- `minFrequencyMultiplier::Number=0.5`: the multiplier used to estimate an edge's minimum frequency requirement from the edge's base frequency
- `maxFrequencyMultiplier::Number=2.0`: the multiplier used to estimate an edge's maximum frequency requirement from the edge's base frequency
- `useNonIntegerDecisionVariables=false`: whether to use relaxed, non-integer decision variables

"""
function ConstructModel(instance::ProblemInstance; useSlack::Bool = false, slack_penalty::Union{Nothing, Float64} = nothing, useBinaryVariant::Bool = false, minFrequencyMultiplier::Number = 0.5, maxFrequencyMultiplier::Number = 2.0, useNonIntegerDecisionVariables = false)::Model
    if slack_penalty == 0.0
        @warn "Setting slack to zero should cause the optimization to not choose any line for the solution."
    end
    if slack_penalty != nothing && ~useSlack
        @warn "Slack penalty is set, but useSlack is not set to true. If the use of slackened constraints is intended, please set the useSlack argument to true"
    end
    if minFrequencyMultiplier > maxFrequencyMultiplier
        @warn "min/max base-frequency multiplier choice may result in an unsolvable model: $(minFrequencyMultiplier), $(maxFrequencyMultiplier)"
    end
    if useNonIntegerDecisionVariables && useBinaryVariant
        error("Unable to use both binary and non-integer decision variables simultaneously")
    end

    #creates JuMP model
    model = Model()

    #add decision variables
    if useBinaryVariant
        #we use binary decision variables: each line can either be fully taken(with a preset frequency), or not taken at all
        @variable(model, y[1:length(instance.lines)], Bin, container = Array)
    elseif useNonIntegerDecisionVariables
        #variables are not constrained as integers
        @variable(model, y[1:length(instance.lines)] >= 0, container = Array)
    else
        #if flexible frequencies are used, the default behaviour, we use integer decision variables:
        #thus, each line can be taken with an arbitrary frequency (including 0 for not taken)
        @variable(model, y[1:length(instance.lines)] >= 0, Int, container = Array)

    end
    #rethinking this restriction, I thought: why cant I just use float values here? The representation with frequencies over a T time span is already an aproximation of real operations anyway...

    if useSlack && !useNonIntegerDecisionVariables
        #slack variables for the altered formulation
        @variable(model, l_slack[1:ne(instance.network)] >= 0, Int, container = Array) #lower slack
        @variable(model, u_slack[1:ne(instance.network)] >= 0, Int, container = Array) #upper slack
        #the cost penalty should stop the optimizer from using both simultaneously for any single edge
    elseif useSlack && useNonIntegerDecisionVariables
        #slack variables not constrained to integer values
        @variable(model, l_slack[1:ne(instance.network)] >= 0, container = Array) #lower slack
        @variable(model, u_slack[1:ne(instance.network)] >= 0, container = Array) #upper slack

    end

    #shortcut auxiliary function for clearer reading:
    IsEdgeOnLine_A(edge::LightGraphs.SimpleGraphs.SimpleEdge, lineIndex::Int64) = IsEdgeOnLine(instance, edge, lineIndex, numerical_returns = true)

    #collect edges
    edges_array = [e for e in edges(instance.network)]

    #add frequency restrictions
    for i in 1:length(edges_array)
        e = edges_array[i]

        #f_min and f_max estimation: turns out to be pretty important in practical experimentation
        f_min = Int64(ceil(get_prop(instance.network, e, :base_freq) * minFrequencyMultiplier)) #ceil is important so we pretty much never get a 0 value, which could cause the line set to not be connected
        f_max = Int64(round(get_prop(instance.network, e, :base_freq) * maxFrequencyMultiplier))

        #=
        #this has sometimes raised a warning about to many unperformant additions beeing called. Maybe restructure it? to-do
        @constraint(model, f_min <=

        if useBinaryVariant
            sum(instance.lines[j].frequency * y[j] * IsEdgeOnLine_A(e, j) for j in 1:length(instance.lines))
        else
            #y[j] already represent a variable line frequency. thus we do not multiply it by instance.lines[j].frequency
            sum(y[j] * IsEdgeOnLine_A(e, j) for j in 1:length(instance.lines))
        end
        +
        if ~useSlack
            0
        else
            #plus l_slack, so having some l_slack != 0 helps accomodate the lower limit
            #minus u_slack, helping with the upper limit
            l_slack[i] - u_slack[i]
        end

        <= f_max)
        =#

        ex = @expression(model,
            if useBinaryVariant
                sum(instance.lines[j].frequency * y[j] * IsEdgeOnLine_A(e, j) for j in 1:length(instance.lines))
            else
                #y[j] already represent a variable line frequency. thus we do not multiply it by instance.lines[j].frequency
                sum(y[j] * IsEdgeOnLine_A(e, j) for j in 1:length(instance.lines))
            end

            +
            if ~useSlack
                0
            else
                #plus l_slack, so having some l_slack != 0 helps accomodate the lower limit
                #minus u_slack, helping with the upper limit
                l_slack[i] - u_slack[i]
            end

        )

        @constraint(model, f_min <= ex <= f_max)
    end

    @objective(model, Min,
        if useBinaryVariant
            sum(y[i] * instance.lines[i].cost * instance.lines[i].frequency for i=1:length(instance.lines))
        else
            #y[i] already represents the taken frequency
            sum(y[i] * instance.lines[i].cost for i=1:length(instance.lines))
        end
        +
        if ~useSlack
            0 #if not using slack, slack penalty is 0
        else
            #sum of taken slack * slack_penalty
            sum( (l_slack[i] + u_slack[i]) * slack_penalty for i=1:ne(instance.network))
        end
    )

    return model
end


text:
using EngEconomics, Roots

# Given
MARR = 0.08

aInit = -45000
aAnnual = -10000 + 18000
aScrap = 3000
aLife = 5

bInit = -25000
bAnnual = -4000 + 13000
bScrap = 6000
bLife = 5

cInit = -20000
cAnnual = -1900 + 9000
cScrap = 4600
cLife = 5

# AW
aAW = aInit * capitalRecoveryFactor(MARR, aLife) + aAnnual + aScrap * sinkingFundFactor(MARR, aLife)
bAW = bInit * capitalRecoveryFactor(MARR, bLife) + bAnnual + bScrap * sinkingFundFactor(MARR, bLife)
cAW = cInit * capitalRecoveryFactor(MARR, cLife) + cAnnual + cScrap * sinkingFundFactor(MARR, cLife)


text:
#=
Compute the running median of a sequence of numbers. That is, given a stream of numbers, print out the median of the list so far on each new element.

Recall that the median of an even-numbered list is the average of the two middle numbers.

For example, given the sequence [2, 1, 5, 7, 2, 0, 5], your algorithm should print out:

2
1.5
2
3.5
2
2
2
=#

function print_running_median(arr::Array; store = false)
    if store
        median_arr = []
    end
    for index in 1:length(arr)
        current_arr = sort(arr[1:index])
        if length(current_arr) % 2 == 0
            median = (current_arr[Int8(length(current_arr)/2)] + current_arr[Int8(length(current_arr)/2 + 1)]) / 2
            println(median)
            if store
                push!(median_arr, median)
            end
        else
            median = current_arr[ceil(Int8, length(current_arr)/2)]
            println(median)
            if store
                push!(median_arr, median)
            end
        end
    end
    if store
        return median_arr
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using Revise
using Soss
using ResumableFunctions

hmmStep = @model (s0, step, noise) begin
    s1 ~ EqualMix(step.(s0))
    y ~ noise(s1.x)
end;

function step(s)
    m = @model s begin
        ν ~ HalfCauchy()
        x ~ StudentT(ν, s,1)
    end
    m(s=s)
end;

function noise(s)
    m = @model s begin
        y ~ Normal(s,1)
    end
    m(s=s)
end;

s0 = rand(Normal(0,10), 100);
particles(s0)

rand(hmmStep(s0=s0, step=step, noise=noise)) |> pairs

dynamicHMC(hmmStep(s0=s0, step=step, noise=noise), (y=(y=1.0,),)) |> particles

step(s) = StudentT(3,s,1);
noise(s) = Normal(s,1);

struct Chain
    init
    step
end

@resumable function Base.rand(c::Chain)
    x = rand(c.init)
    @yield x
    while true
        x = rand(c.step(x))
        @yield x
    end
end

# function Distributions.logpdf(c::Chain, xs)

c = Chain(s0, step)

for(n,x) in enumerate(rand(c))
    n > 10 && break
    println(x)
end

m = @model ch begin
    c ~ ch(Normal(), x -> Normal(x,1))
end


text:
if basename(pwd()) == "aoc"
    cd("2019/6")
end

struct Planet
    name
    parent
    children
end
Planet(name) = Planet(name, nothing, Planet[])

function Base.getindex(planet::Planet, name::AbstractString)
    if planet.name == name
        planet
    else
        for child in planet.children
            result = child[name]
            if !isnothing(result)
                return result
            end
        end
    end
end

Base.haskey(planet::Planet, name::AbstractString) = !isnothing(planet[name])

function loadsystem(filename::AbstractString)
    system = Planet("COM")
    lines = readlines(filename)
    while length(lines) > 0
        line = popfirst!(lines)
        parent, child = split(line, ')')
        if haskey(system, parent)
            push!(system[parent].children, Planet(child, system[parent], Planet[]))
        else
            push!(lines, line)
        end
    end
    system
end

function orbits(planet::Planet)
    depth(planet) * length(planet.children) + sum(orbits.(planet.children), init = 0)
end

depth(::Nothing) = 0
depth(planet::Planet) = 1 + depth(planet.parent)

function parents(planet::Planet)
    if isnothing(planet.parent)
        [planet]
    else
        [planet, parents(planet.parent)...]
    end
end

part1(filename::AbstractString) = orbits(loadsystem(filename))
@assert part1("example.txt") == 42

function part2(filename::AbstractString)
    system = loadsystem(filename)
    you = system["YOU"].parent
    santa = system["SAN"].parent
    common = first(intersect(parents(you), parents(santa)))
    depth(you) + depth(santa) - 2 * depth(common)
end

@assert part2("example2.txt") == 4
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module TopicModels

#Imports
import Base.length

using Random, Distributions
using SpecialFunctions: loggamma
using Clustering: randindex

#Exports
export Corpus,
       LdaCorpus,
       Model,
       State,
       readDocs,
       readLexicon,
       termToWordSequence,
       topTopicWords,
       trainModel,
       CorpusTopics,
       CorpusARI,
       DocsARI

#Data that we make or find in real life:
include("Data.jl")

#Bayesian learning and inference:
include("Computation.jl")

#Stuff like perplexity and ARI:
include("Validation.jl")
end #module


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Simulate Schwarzschild spacetime
#--------------------------------------------------------------------
using NLsolve, Printf, Roots

function find_t_of_UV(U::T, V::T, M::Number)::T where {T<:Number}
    @assert V > 0   # ensure you're in region I or II
    @assert U*V < T(1) # ensure you don't hit the singularity
    if U*V == 0     # r = 2M
        t = V       # enforce uniqueness
    elseif U > 0    # r < 2M
        t = -2M*log(U/V)
    elseif U < 0    # r > 2M
        t = -2M*log(-U/V)
    else
        error("Domain error")
    end
    return t
end

function find_r_of_UV(U::T, V::T, M::Number)::T where {T<:Number}
    @assert V > T(0)       # ensure you're in region I or II
    @assert U*V < T(1)     # ensure you don't hit the singularity
    if U*V == T(0)         # r = 2M
        r = 2M
    else                # r < 2M or r > 2M
        f(r) = (r/2M - 1)*exp(r/2M) + U*V
        try
            # r    = find_zero(f, 2M)
            # NOTE: Introduce bracketing
            if U*V > T(0)
                r = find_zero(f, (eps(T), 2*M))
            else
                # WARNING: Would fail if r > 1000*M
                r = find_zero(f, (2M, 1000*M))
            end
        catch
            @show U, V
            exit()
        end
    end
    @assert r > 0
    return r
end

function find_U_of_tr(t::T, r::T, M::Number)::T where {T<:Number}
    @assert r > 0
    if r == 2M
        U = 0
    else
        rstar = r + 2M*log(abs((r/2M)-1))
        u = t - rstar
        r > 2M ?  U = -exp(-u/4M) : U = +exp(-u/4M)
    end
    return U
end

function find_V_of_tr(t::T, r::T, M::Number)::T where {T<:Number}
    @assert r > 0
    if r == 2M
        V = t           # make the inverse mapping unique.
    else
        rstar = r + 2M*log(abs((r/2M)-1))
        v = t + rstar
        V = exp(v/4M)
    end
    return V
end

function schwarzschild(PS::ProductSpace{S1, S2})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
    η0 = Field(PS, (u,v)->find_r_of_UV(u,v,M))
    ϕ0 = Field(PS, (u,v)->0)
    f0 = ((16*M^3)/η0)*exp(-η0/2M)
    a0 = -sqrt(2*f0)
    return (a0, η0)
end

function singularity(PS::ProductSpace{S1, S2})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
    a = Field(PS, (u,v)->NaN)
    η = Field(PS, (u,v)->NaN)
    return (a, η)
end

function excision(PS::ProductSpace{S1, S2})::Bool where {S1, S2}
    uv = Field(PS, (u,v)->u*v)
    v  = Field(PS, (u,v)->v)
    return (maximum(uv) >= 1.0) || (minimum(v) <= 0.0)
end

function excision(bnd::NTuple{2, Field{S}})::Bool where {S}
    return any(map(isnan, bnd))
end

function computeUboundary(PS::ProductSpace{S1, S2})::NTuple{2, Field{S2}} where {S1, S2}
    if excision(PS) == true
        return map(x->extractUboundary(x, :incoming), singularity(PS))
    else
        return map(x->extractUboundary(x, :incoming), schwarzschild(PS))
    end
end

function computeVboundary(PS::ProductSpace{S1, S2})::NTuple{2, Field{S1}} where {S1, S2}
    if excision(PS) == true
        return map(x->extractVboundary(x, :incoming), singularity(PS))
    else
        return map(x->extractVboundary(x, :incoming), schwarzschild(PS))
    end
end

function computePatch(PS::ProductSpace{S1, S2}, ubnd::NTuple{2, Field{S2}},
                                                vbnd::NTuple{2, Field{S1}})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
    if excision(PS) || excision(ubnd) || excision(vbnd)
        return singularity(PS)
    else
        B = incomingboundary(PS)
        I = identity(PS)
        DU, DV = derivative(PS)
        ∂U = combineUVboundary(ubnd, vbnd, :incoming)

        function constr(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
            (a, η) = U
            C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η)
            C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η)
            C  = sqrt(C1*C1 + C2*C2)
            return (C, C)
        end

        function offaxis(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
            (a, η) = U
            F1 = DU*(DV*a) - (1/a)*(DU*a)*(DV*a) + (a/η)*(DU*(DV*η))
            F2 = DU*(DV*η) + (1/η)*(DU*η)*(DV*η) + (1/4)*(1/η)*(a^2)
            C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η)
            C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η)
            return (F1 , F2)
        end

        function boundary(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
            return U - ∂U
        end

        function residual(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
            F = offaxis(U) + λ * constr(U)
            enforcebc!(F, boundary(U))
            return F
        end

        function initialguess(PS::ProductSpace{S1, S2})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
            a = Field(PS, (u,v)->2)
            η = Field(PS, (u,v)->(v-u))
            return (a, η)
        end

        function f!(fvec::Array{T,1}, x::Array{T,1}) where {T}
            fvec[:] = reshapeFromTuple(residual(reshapeToTuple(PS, 2, x)))
        end

        s = nlsolve(f!, reshapeFromTuple(initialguess(PS)); method=:trust_region, autodiff=:forward,
                    show_trace=debug >= 3, ftol=1e-9, iterations=100)
        k = reshapeToTuple(PS, 2, s.zero)

        if debug > 1
            println("    Converged?            = ", converged(s))
        end

        if converged(s) == false
            return singularity(PS)
        else
            return k
        end
    end
end

function pconv(range::StepRange, nh::Int=1)
    println("Testing p convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        np = n_[index]
        params = Parameters((np, np), (nh,nh), urange, vrange, 2)
        q = distribute(params, computePatch, computeUboundary, computeVboundary)
        l_[index] = rmse(extract(map(constraint, q), 1))
        @printf("  n = %i, rmse = %e\n", n_[index], l_[index])
    end
    return (n_, l_)
end

function hconv(range::StepRange, np::Int=4)
    println("Testing h convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        nh = n_[index] = 2^n_[index]
        params = Parameters((np, np), (nh, nh), urange, vrange, 2)
        l_[index] = rmse(extract(map(constraint, distribute(params, computePatch, computeUboundary, computeVboundary)), 1))
        l0 = index[1] > 1 ? l_[index[1] - 1] : 1
        @printf("  n = %3i, rmse[2h] / rmse[h] = %e\n", n_[index], l0 / l_[index])
    end
    return (n_, l_)
end

function constraint(U::NTuple{2, Field{S}})::NTuple{1, Field{S}} where {S}
    (a, η) = U
    PS = first(U).space
    DU, DV = derivative(PS)
    C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η)
    C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η)
    C  = sqrt(C1*C1 + C2*C2)
    return (C,)
end

function error(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
    (a, η) = U
    PS = first(U).space
    if excision(PS) == false
        (as, ηs) = schwarzschild(PS)
    else
        (as, ηs) = singularity(PS)
    end
    return (a - as, η - ηs)
end

#-----------------------------------------
# Setup simulation grid parameters
#-----------------------------------------
npoints  = (10, 10)
npatches = (12, 12)
M        = 0.125
debug    = 1
λ        = 0.0

# urange   = (-3.0, 3.0)
# vrange   = ( 2.0, 8.0)
urange   = (-9.0, -3.0)
vrange   = ( 2.0,  3.0)

params   = Parameters(npoints, npatches, urange, vrange, 2)
s        = distribute(params, computePatch, computeUboundary, computeVboundary, debug)
@show minimum(extract(s, 2))
save("output/data/schwarzschild/schwarzschild_M_0p1_a", extract(s, 1))
save("output/data/schwarzschild/schwarzschild_M_0p1_eta", extract(s, 2))
save("output/data/schwarzschild/schwarzschild_M_0p1_constraints", extract(map(constraint, s), 1))

#-----------------------------------------
# Check convergence away from the singularity
# and then near the singularity
# NOTE: We only get convergence when we don't include
# the constraints
#-----------------------------------------
np, ep = pconv(2:2:12, 12)
nh, eh = hconv(1:1:7, 4)
save("output/data/schwarzschild/schwarzschild_M_0p1_pconv", np, ep)
save("output/data/schwarzschild/schwarzschild_M_0p1_hconv", nh, eh)


text:
"""
"""
function autodiff_array_gradient(a,i_to_x,j_to_i=IdentityVector(length(i_to_x)))

  i_to_xdual = lazy_map(i_to_x) do x
    cfg = ForwardDiff.GradientConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
    xdual = cfg.duals
    xdual
  end

  j_to_f = to_array_of_functions(a,i_to_xdual,j_to_i)
  j_to_x = lazy_map(Reindex(i_to_x),j_to_i)

  k = ForwardDiffGradientMap()
  lazy_map(k,j_to_f,j_to_x)

end

struct ForwardDiffGradientMap <: Map end

function return_cache(k::ForwardDiffGradientMap,f,x)
  cfg = ForwardDiff.GradientConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
  r = copy(x)
  (r, cfg)
end

@inline function evaluate!(cache,k::ForwardDiffGradientMap,f,x)
  r, cfg = cache
  @notimplementedif length(r) != length(x)
  ForwardDiff.gradient!(r,f,x,cfg)
  r
end

"""
"""
function autodiff_array_jacobian(a,i_to_x,j_to_i=IdentityVector(length(i_to_x)))

  i_to_xdual = lazy_map(i_to_x) do x
    cfg = ForwardDiff.JacobianConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
    xdual = cfg.duals
    xdual
  end

  j_to_f = to_array_of_functions(a,i_to_xdual,j_to_i)
  j_to_x = lazy_map(Reindex(i_to_x),j_to_i)

  k = ForwardDiffJacobianMap()
  lazy_map(k,j_to_f,j_to_x)

end

struct ForwardDiffJacobianMap <: Map end

function return_cache(k::ForwardDiffJacobianMap,f,x)
  cfg = ForwardDiff.JacobianConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
  n = length(x)
  j = zeros(eltype(x),n,n)
  (j, cfg)
end

@inline function evaluate!(cache,k::ForwardDiffJacobianMap,f,x)
  j, cfg = cache
  @notimplementedif size(j,1) != length(x)
  @notimplementedif size(j,2) != length(x)
  ForwardDiff.jacobian!(j,f,x,cfg)
  j
end

"""
"""
function autodiff_array_hessian(a,i_to_x,j_to_i=IdentityVector(length(i_to_x)))
   agrad = i_to_y -> autodiff_array_gradient(a,i_to_y,j_to_i)
   autodiff_array_jacobian(agrad,i_to_x,j_to_i)
end

function to_array_of_functions(a,x,ids=IdentityVector(length(x)))
  k = ArrayOfFunctionsMap(a,x)
  j = IdentityVector(length(ids))
  lazy_map(k,j)
end

struct ArrayOfFunctionsMap{A,X} <: Map
  a::A
  x::X
end

function return_cache(k::ArrayOfFunctionsMap,j)
  xi = testitem(k.x)
  l = length(k.x)
  x = MutableFill(xi,l)
  ax = k.a(x)
  axc = array_cache(ax)
  (ax, x, axc)
end

@inline function evaluate!(cache,k::ArrayOfFunctionsMap,j)
  ax, x, axc = cache
  @inline function f(xj)
    x.value = xj
    axj = getindex!(axc,ax,j)
  end
  f
end

mutable struct MutableFill{T} <: AbstractVector{T}
  value::T
  length::Int
end

Base.size(a::MutableFill) = (a.length,)

@inline Base.getindex(a::MutableFill,i::Integer) = a.value


text:
# Copyright 2021 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Computes a performance profile for a set of instances from the netlib
# infeasibility database.
# To run this script execute
# $ julia --project performance_netlib_experiment.jl <path_where_to_store_results> <path_to_directory_containing_mps_files>

import FirstOrderLp

include("primal_dual_hybrid_gradient.jl")
include("utils.jl")

function run_identifiability_profile(
  instance_names::Vector{String},
  accuracy::Float64,
)
  outputs = Vector{PdhgOutput}()
  for instance_name in instance_names
    params = PdhgParameters(
      nothing, # step_size (forces the solver to use a provably correct step size)
      1.0, # primal_weight
      1000, # printing frequency
      false, # verbose
      1000000, # iteration limit
      nothing, # two_dimensional_subspace
      accuracy, # eps_primal_infasible
      accuracy, # eps_dual_infeasible
      false, # store_iterates
      nothing, # reference_point
    )
    print("About to start solving " * instance_name * "\n")
    instance_path = joinpath(ARGS[2], (instance_name * ".mps"))
    lp = FirstOrderLp.qps_reader_to_standard_form(instance_path)
    push!(outputs, optimize(params, lp))
  end
  instance_path =
    joinpath(ARGS[1], "per_instance_detection" * string(accuracy) * ".csv")
  path_csv = joinpath(ARGS[1], "detected_fraction" * string(accuracy) * ".csv")
  save_fraction_of_infeasible_instances_detected_to_csv(
    outputs,
    instance_names,
    instance_path,
    path_csv,
  )
  generate_plot_fraction_infeasible_instances_detected(
    path_csv,
    joinpath(ARGS[1], "performance_" * string(accuracy) * ".pdf"),
  )
end

function main()
  instance_names = [
    "bgdbg1",
    "bgetam",
    "bgindy",
    "bgprtr",
    "box1",
    "ceria3d",
    "chemcom",
    "cplex1",
    "cplex2",
    "ex72a",
    "ex73a",
    "forest6",
    "galenet",
    "gosh",
    "gran",
    "greenbea",
    "itest2",
    "itest6",
    "klein1",
    "klein2",
    "klein3",
    "mondou2",
    "pang",
    "qual",
    "reactor",
    "refinery",
    "vol1",
    "woodinfe",
  ]
  run_identifiability_profile(instance_names, 1e-4)
  run_identifiability_profile(instance_names, 1e-8)
end

main()


text:
# Recurrent type
g = NEAT.Global(config)
g.cf.feedforward = false
g.cf.nn_activation = :sigm
ch1 = NEAT.create_unconnected(g)
NEAT.add_hidden_nodes!(g,ch1, 2, NEAT.Recurrent())
rec = NEAT.createPhenotype(ch1)
@test length(rec.neurons) == length(ch1.node_genes)
@test length(rec.synapses) == length(ch1.connection_genes)
@test rec.nntype == ch1.node_gene_type
@test rec.numInputs == ch1.inputCnt

@test typeof(NEAT.activate(rec.nntype, rec, [1.,1.,1.])) == Vector{Float64}
@test NEAT.activate(rec.nntype, rec, [1.,1.,1.])[1] > 0. # due to sigmoids

g = NEAT.Global(config)
g.cf.prob_addnode = 0.5
p = NEAT.Population(g)
p.evaluate = (chs)->map(ch->ch.fitness=max(0,randn()*0.3+ch.fitness*0.7), chs)
NEAT.epoch(g, p, 10, false)
rec = NEAT.createPhenotype(p.best_fitness[end])
@test typeof(NEAT.activate(rec.nntype, rec, [1.,1.,1.])) == Vector{Float64}
@test NEAT.activate(rec.nntype, rec, [1.,1.,1.])[1] > 0. # due to sigmoids

# FeedForward type
g = NEAT.Global(config)
g.cf.feedforward = true
ch2 = NEAT.create_unconnected(g)
NEAT.add_hidden_nodes!(g,ch2, 2,NEAT.FeedForward())
ff =  NEAT.createPhenotype(ch2)
@test length(ff.neurons) == length(ch2.node_genes)
@test length(ff.synapses) == length(ch2.connection_genes)
@test ff.nntype == ch2.node_gene_type
@test ff.numInputs == ch2.inputCnt

@test typeof(NEAT.activate(ff.nntype, ff, [1.,1.,1.])) == Vector{Float64}
@test NEAT.activate(ff.nntype, ff, [1.,1.,1.])[1] > 0. # due to sigmoids

g = NEAT.Global(config)
g.cf.prob_addnode = 0.5
p = NEAT.Population(g)
p.evaluate = (chs)->map(ch->ch.fitness=max(0,randn()*0.3+ch.fitness*0.7), chs)
NEAT.epoch(g, p, 10, false)
ff = NEAT.createPhenotype(p.best_fitness[end])
@test typeof(NEAT.activate(ff.nntype, rec, [1.,1.,1.])) == Vector{Float64}
@test NEAT.activate(rec.nntype, rec, [1.,1.,1.])[1] > 0. # due to sigmoids


text:
function uhlig01()
#
#  test polynomial used by F. Uhlig
#
    a = 0.01;
    p1 = [1;0;0;0;-a^4]
    p2 = reverse(poly([a;a;a;a]).a)
    p = conv(p1,p2)
    z = [-0.01000000000000 - 0.00000000000000im
         -0.00000000000000 + 0.01000000000000im
         -0.00000000000000 - 0.01000000000000im
         0.01000000000000 - 0.00000000000000im];
    z = [z [1;1;1;5]]
    p, PolyZeros(z)
end


text:
module TestRightIf

using DataTools
using Test
using Transducers: Map, Take

reduce_bs1(args...; kw...) = reduce(args...; basesize = 1, kw...)

@testset for fold in [foldl, reduce_bs1, reduce]
    foldl = nothing
    table = 43:100 |> Map(x -> (k = gcd(x, 42), v = x))
    @test fold(rightif(<), Map(x -> x.k), table) == 42
    @test fold(rightif(>), Map(x -> x.k), table) == 1
    @test fold(rightif(<, x -> x.k), table) == (k = 42, v = 84)
end

end  # module


text:
function checkord(A::Matrix{Ty}, tol=10) where {Ty<:Complex}
    n = size(A,1)
    ulp = eps(real(Ty))

    S = schur(A)
    T2 = copy(S.T)
    Z2 = copy(S.Z)
    select = fill(true,n)
    for i=1:(n>>1)
        select[i] = false
    end
    T2,Z2,v = invoke(LinearAlgebra._ordschur!,
           Tuple{StridedMatrix{T}, StridedMatrix{T},
                 Union{Vector{Bool},BitVector}} where T <: Complex,
        T2,Z2,select)

    # usual tests for Schur
    @test all(tril(T2,-1) .== 0)
    @test norm(Z2*T2*Z2' - A) / (n * norm(A) * ulp) < tol
    @test norm(I - Z2 * Z2') / (n * ulp) < tol

    # make sure we got the ones we asked for
    nwanted = count(select)
    wwanted = [S.values[j] for j=1:n if select[j]]
    errs = zeros(nwanted)
    for i=1:nwanted
        w = T2[i,i]
        errs[i] = minimum(abs.(wwanted .- w)) / (ulp + abs(w))
    end
    @test all(errs .< tol)
end

using LinearAlgebra.LAPACK: trsen!, BlasInt

function checkecond(A::Matrix{T}, nsub=3) where T
    n = size(A,1)
    select = fill(false,n)
    inds = rand(1:n,nsub)
    for i in inds
        select[i] = true
    end
    S = schur(A)

    @assert real(T) <: BlasFloat
    T1 = copy(S.T)
    Z1 = copy(S.Z)
    fselect = BlasInt.([x ? 1 : 0 for x in select])
    T1,Z1,w1,s1,sep1 = trsen!('B','V',fselect,T1,Z1)

    S2 = ordschur(S, select)
    s2 = eigvalscond(S2, count(select))
    sep2 = subspacesep(S2, count(select))

    @test s1 ≈ s2
    @test sep1 ≈ sep2

end

for Ty in [Complex{Float64}]
    @testset "ordschur $Ty" begin
        n = 32
        A = rand(Ty,n,n)
        checkord(A)
    end
    @testset "subspace condition $Ty" begin
        n = 32
        A = rand(Ty,n,n)
        checkecond(A)
    end
end

function checkord(A::Matrix{Ty}, tol=10) where {Ty<:Real}
    n = size(A,1)
    ulp = eps(real(Ty))

    S = schur(A)
    v = S.values
    # cases to try: 2 reals, 1 complex pair, 1 real & 1 pair, 2 pairs
    nr = [2,0,1,0]
    ni = [0,2,2,4]
    for icase in 1:4
        T2 = copy(S.T)
        Z2 = copy(S.Z)
        select = fill(false,n)
        i = n+1
        nr_needed = nr[icase]
        ni_needed = ni[icase]
        while (nr_needed + ni_needed > 0) && i>2
            i -= 1
            if nr_needed > 0
                if imag(v[i]) == 0
                    select[i] = true
                    nr_needed -= 1
                end
            end
            if ni_needed > 0
                ivi = imag(v[i])
                if (abs(ivi) > 1000ulp) && (imag(v[i-1]) == -ivi)
                    select[i] = true
                    select[i-1] = true
                    ni_needed -= 2
                    i -= 1
                end
            end
        end
        if nr_needed + ni_needed > 0
            @warn "failed to find acceptable subset (class $icase) for checking ordschur"
            continue
        end
        T2,Z2,v = invoke(LinearAlgebra._ordschur!,
                         Tuple{StridedMatrix{T}, StridedMatrix{T},
                               Union{Vector{Bool},BitVector}} where T <: Real,
                         T2,Z2,select)

        # usual tests for Schur
        @test all(tril(T2,-2) .== 0)
        @test norm(Z2*T2*Z2' - A) / (n * norm(A) * ulp) < tol
        @test norm(I - Z2 * Z2') / (n * ulp) < tol

        # make sure we got the ones we asked for
        nwanted = count(select)
        wwanted = [S.values[j] for j=1:n if select[j]]
        errs = zeros(nwanted)
        for i=1:nwanted
            w = T2[i,i]
            errs[i] = minimum(abs.(wwanted .- w)) / (ulp + abs(w))
        end
        @test all(errs .< tol)
    end
end

for Ty in [Float64, BigFloat]
    @testset "ordschur $Ty" begin
        # should be big enough to be nearly certain of meaningful selection
        n = 20
        A = rand(Ty,n,n)
        checkord(A)
    end
end


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

const _SemiSets{T} = Union{MOI.Semicontinuous{T},MOI.Semiinteger{T}}

"""
    SemiToBinaryBridge{T, S <: MOI.AbstractScalarSet}

The `SemiToBinaryBridge` replaces a Semicontinuous constraint:
``x \\in \\mathsf{Semicontinuous}(l, u)``
is replaced by:
``z \\in \\{0, 1\\}``,
``x \\leq z \\cdot u ``,
``x \\geq z \\cdot l ``.

The `SemiToBinaryBridge` replaces a Semiinteger constraint:
``x \\in Semiinteger(l, u)``
is replaced by:
``z \\in \\{0, 1\\}``,
``x \\in \\mathbb{Z}``,
``x \\leq z \\cdot u ``,
``x \\geq z \\cdot l ``.
"""
mutable struct SemiToBinaryBridge{T,S<:_SemiSets{T}} <: AbstractBridge
    semi_set::S
    variable::MOI.VariableIndex
    binary_variable::MOI.VariableIndex
    binary_constraint_index::MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}
    lower_bound_index::MOI.ConstraintIndex{
        MOI.ScalarAffineFunction{T},
        MOI.GreaterThan{T},
    }
    upper_bound_index::MOI.ConstraintIndex{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    }
    integer_index::Union{
        Nothing,
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer},
    }
end

const SemiToBinary{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SemiToBinaryBridge{T},OT}

function bridge_constraint(
    ::Type{SemiToBinaryBridge{T,S}},
    model::MOI.ModelLike,
    f::MOI.VariableIndex,
    s::S,
) where {T<:Real,S<:_SemiSets{T}}
    binary, binary_con = MOI.add_constrained_variable(model, MOI.ZeroOne())
    # var - LB * bin >= 0
    lb = MOI.Utilities.operate(*, T, -s.lower, binary)
    lb = MOI.Utilities.operate!(+, T, lb, f)
    lb_ci = MOI.add_constraint(model, lb, MOI.GreaterThan{T}(zero(T)))
    # var - UB * bin <= 0
    ub = MOI.Utilities.operate(*, T, -s.upper, binary)
    ub = MOI.Utilities.operate!(+, T, ub, f)
    ub_ci = MOI.add_constraint(model, ub, MOI.LessThan{T}(zero(T)))
    if s isa MOI.Semiinteger{T}
        int_ci = MOI.add_constraint(model, f, MOI.Integer())
    else
        int_ci = nothing
    end
    return SemiToBinaryBridge{T,S}(
        s,
        f,
        binary,
        binary_con,
        lb_ci,
        ub_ci,
        int_ci,
    )
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:SemiToBinaryBridge{T,S}},
) where {T,S}
    return Tuple{Type}[(MOI.ZeroOne,)]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:SemiToBinaryBridge{T,S}},
) where {T,S<:MOI.Semicontinuous{T}}
    return Tuple{Type,Type}[
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
    ]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:SemiToBinaryBridge{T,S}},
) where {T,S<:MOI.Semiinteger{T}}
    return Tuple{Type,Type}[
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
        (MOI.VariableIndex, MOI.Integer),
    ]
end

function concrete_bridge_type(
    ::Type{<:SemiToBinaryBridge{T}},
    ::Type{MOI.VariableIndex},
    ::Type{S},
) where {T,S<:_SemiSets}
    return SemiToBinaryBridge{T,S}
end

function MOI.supports_constraint(
    ::Type{<:SemiToBinaryBridge},
    ::Type{MOI.VariableIndex},
    ::Type{<:_SemiSets},
)
    return true
end

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, b::SemiToBinaryBridge)
    return b.semi_set
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::SemiToBinaryBridge{T,S},
    set::S,
) where {T,S}
    bridge.semi_set = set
    MOI.modify(
        model,
        bridge.upper_bound_index,
        MOI.ScalarCoefficientChange(bridge.binary_variable, -set.upper),
    )
    MOI.modify(
        model,
        bridge.lower_bound_index,
        MOI.ScalarCoefficientChange(bridge.binary_variable, -set.lower),
    )
    return
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    b::SemiToBinaryBridge{T},
) where {T}
    return b.variable
end

function MOI.delete(model::MOI.ModelLike, bridge::SemiToBinaryBridge)
    if bridge.integer_index !== nothing
        MOI.delete(model, bridge.integer_index)
    end
    MOI.delete(model, bridge.upper_bound_index)
    MOI.delete(model, bridge.lower_bound_index)
    MOI.delete(model, bridge.binary_constraint_index)
    MOI.delete(model, bridge.binary_variable)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::SemiToBinaryBridge,
)
    return MOI.get(
        model,
        MOI.VariablePrimal(attr.result_index),
        bridge.variable,
    )
end

function MOI.supports(
    model::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    ::Type{SemiToBinaryBridge{T,S}},
) where {T,S}
    ci_1 = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}}
    ci_2 = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.LessThan{T}}
    return MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex) &&
           MOI.supports(model, MOI.ConstraintPrimalStart(), ci_1) &&
           MOI.supports(model, MOI.ConstraintPrimalStart(), ci_2)
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    bridge::SemiToBinaryBridge,
)
    return MOI.get(model, MOI.VariablePrimalStart(), bridge.variable)
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    bridge::SemiToBinaryBridge{T},
    value,
) where {T}
    MOI.set(model, MOI.VariablePrimalStart(), bridge.variable, value)
    bin_value = iszero(value) ? 0.0 : 1.0
    MOI.set(model, MOI.VariablePrimalStart(), bridge.binary_variable, bin_value)
    MOI.set(
        model,
        MOI.ConstraintPrimalStart(),
        bridge.upper_bound_index,
        value - bridge.semi_set.upper * bin_value,
    )
    MOI.set(
        model,
        MOI.ConstraintPrimalStart(),
        bridge.lower_bound_index,
        value - bridge.semi_set.lower * bin_value,
    )
    return
end

MOI.get(::SemiToBinaryBridge, ::MOI.NumberOfVariables)::Int64 = 1

function MOI.get(b::SemiToBinaryBridge, ::MOI.ListOfVariableIndices)
    return [b.binary_variable]
end

function MOI.get(
    ::SemiToBinaryBridge{T,S},
    ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.ZeroOne},
)::Int64 where {T,S}
    return 1
end

function MOI.get(
    ::SemiToBinaryBridge{T,S},
    ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.Integer},
)::Int64 where {T,S<:MOI.Semiinteger}
    return 1
end

function MOI.get(
    ::SemiToBinaryBridge{T,S},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}},
)::Int64 where {T,S}
    return 1
end

function MOI.get(
    ::SemiToBinaryBridge{T,S},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
)::Int64 where {T,S}
    return 1
end

function MOI.get(
    b::SemiToBinaryBridge{T,S},
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne},
) where {T,S}
    return [b.binary_constraint_index]
end

function MOI.get(
    b::SemiToBinaryBridge{T,S},
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer},
) where {T,S<:MOI.Semiinteger}
    return [b.integer_index]
end

function MOI.get(
    b::SemiToBinaryBridge{T,S},
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
) where {T,S}
    return [b.upper_bound_index]
end

function MOI.get(
    b::SemiToBinaryBridge{T,S},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T},
        MOI.GreaterThan{T},
    },
) where {T,S}
    return [b.lower_bound_index]
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
typealias P2 FixedSizeArrays.Vec{2,Float64}
typealias P3 FixedSizeArrays.Vec{3,Float64}

nanpush!(a::AbstractVector{P2}, b) = (push!(a, P2(NaN,NaN)); push!(a, b))
nanappend!(a::AbstractVector{P2}, b) = (push!(a, P2(NaN,NaN)); append!(a, b))
nanpush!(a::AbstractVector{P3}, b) = (push!(a, P3(NaN,NaN,NaN)); push!(a, b))
nanappend!(a::AbstractVector{P3}, b) = (push!(a, P3(NaN,NaN,NaN)); append!(a, b))
compute_angle(v::P2) = (angle = atan2(v[2], v[1]); angle < 0 ? 2π - angle : angle)

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

immutable Shape
  # vertices::AVec
  x::AVec
  y::AVec
end

# Shape(x, y) = Shape(collect(zip(x, y)))
Shape(verts::AVec) = Shape(unzip(verts)...)

# get_xs(shape::Shape) = Float64[v[1] for v in shape.vertices]
# get_ys(shape::Shape) = Float64[v[2] for v in shape.vertices]
get_xs(shape::Shape) = shape.x
get_ys(shape::Shape) = shape.y
vertices(shape::Shape) = collect(zip(shape.x, shape.y))

function shape_coords(shape::Shape)
    # unzip(shape.vertices)
    shape.x, shape.y
end

function shape_coords(shapes::AVec{Shape})
    length(shapes) == 0 && return zeros(0), zeros(0)
    xs = map(get_xs, shapes)
    ys = map(get_ys, shapes)
    # x, y = shapes[1].x, shapes[1].y #unzip(shapes[1].vertices)
    x, y = map(copy, shape_coords(shapes[1]))
    for shape in shapes[2:end]
        # tmpx, tmpy = unzip(shape.vertices)
        nanappend!(x, shape.x)
        nanappend!(y, shape.y)
        # x = vcat(x, NaN, tmpx)
        # y = vcat(y, NaN, tmpy)
    end
    x, y
end

"get an array of tuples of points on a circle with radius `r`"
function partialcircle(start_θ, end_θ, n = 20, r=1)
  @compat(Tuple{Float64,Float64})[(r*cos(u),r*sin(u)) for u in linspace(start_θ, end_θ, n)]
end

"interleave 2 vectors into each other (like a zipper's teeth)"
function weave(x,y; ordering = Vector[x,y])
  ret = eltype(x)[]
  done = false
  while !done
    for o in ordering
      try
          push!(ret, shift!(o))
      end
      # try
      #     push!(ret, shift!(y))
      # end
    end
    done = isempty(x) && isempty(y)
  end
  ret
end

"create a star by weaving together points from an outer and inner circle.  `n` is the number of arms"
function makestar(n; offset = -0.5, radius = 1.0)
    z1 = offset * π
    z2 = z1 + π / (n)
    outercircle = partialcircle(z1, z1 + 2π, n+1, radius)
    innercircle = partialcircle(z2, z2 + 2π, n+1, 0.4radius)
    Shape(weave(outercircle, innercircle)[1:end-2])
end

"create a shape by picking points around the unit circle.  `n` is the number of point/sides, `offset` is the starting angle"
function makeshape(n; offset = -0.5, radius = 1.0)
    z = offset * π
    Shape(partialcircle(z, z + 2π, n+1, radius)[1:end-1])
end

function makecross(; offset = -0.5, radius = 1.0)
    z2 = offset * π
    z1 = z2 - π/8
    outercircle = partialcircle(z1, z1 + 2π, 9, radius)
    innercircle = partialcircle(z2, z2 + 2π, 5, 0.5radius)
    Shape(weave(outercircle, innercircle,
                ordering=Vector[outercircle,innercircle,outercircle])[1:end-2])
end

from_polar(angle, dist) = P2(dist*cos(angle), dist*sin(angle))

function makearrowhead(angle; h = 2.0, w = 0.4)
    tip = from_polar(angle, h)
    Shape(P2[(0,0), from_polar(angle - 0.5π, w) - tip,
        from_polar(angle + 0.5π, w) - tip, (0,0)])
end

const _shapes = KW(
    :ellipse    => makeshape(20),
    :rect       => makeshape(4, offset=-0.25),
    :diamond    => makeshape(4),
    :utriangle  => makeshape(3),
    :dtriangle  => makeshape(3, offset=0.5),
    :pentagon   => makeshape(5),
    :hexagon    => makeshape(6),
    :heptagon   => makeshape(7),
    :octagon    => makeshape(8),
    :cross      => makecross(offset=-0.25),
    :xcross     => makecross(),
    :vline      => Shape([(0,1),(0,-1)]),
    :hline      => Shape([(1,0),(-1,0)]),
  )

for n in [4,5,6,7,8]
  _shapes[symbol("star$n")] = makestar(n)
end

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

# center(shape::Shape) = (mean(shape.x), mean(shape.y))

# uses the centroid calculation from https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
function center(shape::Shape)
    x, y = shape_coords(shape)
    n = length(x)
    A, Cx, Cy = 0.0, 0.0, 0.0
    for i=1:n
        ip1 = i==n ? 1 : i+1
        A += x[i] * y[ip1] - x[ip1] * y[i]
    end
    A *= 0.5
    for i=1:n
        ip1 = i==n ? 1 : i+1
        m = (x[i] * y[ip1] - x[ip1] * y[i])
        Cx += (x[i] + x[ip1]) * m
        Cy += (y[i] + y[ip1]) * m
    end
    Cx / 6A, Cy / 6A
end

function Base.scale!(shape::Shape, x::Real, y::Real = x, c = center(shape))
    sx, sy = shape_coords(shape)
    cx, cy = c
    for i=1:length(sx)
        sx[i] = (sx[i] - cx) * x + cx
        sy[i] = (sy[i] - cy) * y + cy
    end
    shape
end

function Base.scale(shape::Shape, x::Real, y::Real = x, c = center(shape))
    shapecopy = deepcopy(shape)
    scale!(shape, x, y, c)
end

function translate!(shape::Shape, x::Real, y::Real = x)
    sx, sy = shape_coords(shape)
    for i=1:length(sx)
        sx[i] += x
        sy[i] += y
    end
    shape
end

function translate(shape::Shape, x::Real, y::Real = x)
    shapecopy = deepcopy(shape)
    translate!(shape, x, y)
end

function rotate_x(x::Real, y::Real, Θ::Real, centerx::Real, centery::Real)
    (x - centerx) * cos(Θ) - (y - centery) * sin(Θ) + centerx
end

function rotate_y(x::Real, y::Real, Θ::Real, centerx::Real, centery::Real)
    (y - centery) * cos(Θ) + (x - centerx) * sin(Θ) + centery
end

function rotate(x::Real, y::Real, θ::Real, c = center(shape))
    cx, cy = c
    rotate_x(x, y, Θ, cx, cy), rotate_y(x, y, Θ, cx, cy)
end

function rotate!(shape::Shape, Θ::Real, c = center(shape))
    x, y = shape_coords(shape)
    cx, cy = c
    for i=1:length(x)
        x[i] = rotate_x(x[i], y[i], Θ, cx, cy)
        y[i] = rotate_y(x[i], y[i], Θ, cx, cy)
    end
    shape
end

function rotate(shape::Shape, Θ::Real, c = center(shape))
    shapecopy = deepcopy(shape)
    rotate!(shapecopy, Θ, c)
end

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

immutable Font
  family::AbstractString
  pointsize::Int
  halign::Symbol
  valign::Symbol
  rotation::Float64
  color::Colorant
end

"Create a Font from a list of unordered features"
function font(args...)

  # defaults
  family = "Helvetica"
  pointsize = 14
  halign = :hcenter
  valign = :vcenter
  rotation = 0.0
  color = colorant"black"

  for arg in args
    T = typeof(arg)

    if arg == :center
      halign = :hcenter
      valign = :vcenter
    elseif arg in (:hcenter, :left, :right)
      halign = arg
    elseif arg in (:vcenter, :top, :bottom)
      valign = arg
    elseif T <: Colorant
      color = arg
    elseif T <: @compat Union{Symbol,AbstractString}
      try
        color = parse(Colorant, string(arg))
      catch
        family = string(arg)
      end
    elseif typeof(arg) <: Integer
      pointsize = arg
    elseif typeof(arg) <: Real
      rotation = convert(Float64, arg)
    else
      warn("Unused font arg: $arg ($(typeof(arg)))")
    end
  end

  Font(family, pointsize, halign, valign, rotation, color)
end

"Wrap a string with font info"
immutable PlotText
  str::@compat(AbstractString)
  font::Font
end
PlotText(str) = PlotText(string(str), font())

function text(str, args...)
  PlotText(string(str), font(args...))
end

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

immutable Stroke
  width
  color
  alpha
  style
end

function stroke(args...; alpha = nothing)
  # defaults
  # width = 1
  # color = colorant"black"
  # style = :solid
  width = nothing
  color = nothing
  style = nothing

  for arg in args
    T = typeof(arg)

    # if arg in _allStyles
    if allStyles(arg)
      style = arg
    elseif T <: Colorant
      color = arg
    elseif T <: @compat Union{Symbol,AbstractString}
      try
        color = parse(Colorant, string(arg))
      end
    # elseif trueOrAllTrue(a -> typeof(a) <: Real && a > 0 && a < 1, arg)
    elseif allAlphas(arg)
      alpha = arg
    # elseif typeof(arg) <: Real
    elseif allReals(arg)
      width = arg
    else
      warn("Unused stroke arg: $arg ($(typeof(arg)))")
    end
  end

  Stroke(width, color, alpha, style)
end

immutable Brush
  size  # fillrange, markersize, or any other sizey attribute
  color
  alpha
end

function brush(args...; alpha = nothing)
  # defaults
  # sz = 1
  # color = colorant"black"
  size = nothing
  color = nothing

  for arg in args
    T = typeof(arg)

    if T <: Colorant
      color = arg
    elseif T <: @compat Union{Symbol,AbstractString}
      try
        color = parse(Colorant, string(arg))
      end
    # elseif trueOrAllTrue(a -> typeof(a) <: Real && a > 0 && a < 1, arg)
    elseif allAlphas(arg)
      alpha = arg
    # elseif typeof(arg) <: Real
    elseif allReals(arg)
      size = arg
    else
      warn("Unused brush arg: $arg ($(typeof(arg)))")
    end
  end

  Brush(size, color, alpha)
end

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

"type which represents z-values for colors and sizes (and anything else that might come up)"
immutable ZValues
  values::Vector{Float64}
  zrange::Tuple{Float64,Float64}
end

function zvalues{T<:Real}(values::AVec{T}, zrange::Tuple{T,T} = (minimum(values), maximum(values)))
  ZValues(collect(float(values)), map(Float64, zrange))
end

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

abstract AbstractSurface

"represents a contour or surface mesh"
immutable Surface{M<:AMat} <: AbstractSurface
  # x::AVec
  # y::AVec
  surf::M
end

Surface(f::Function, x, y) = Surface(Float64[f(xi,yi) for xi in x, yi in y])

Base.Array(surf::Surface) = surf.surf

for f in (:length, :size)
  @eval Base.$f(surf::Surface, args...) = $f(surf.surf, args...)
end
Base.copy(surf::Surface) = Surface(copy(surf.surf))

"For the case of representing a surface as a function of x/y... can possibly avoid allocations."
immutable SurfaceFunction <: AbstractSurface
    f::Function
end

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

type OHLC{T<:Real}
  open::T
  high::T
  low::T
  close::T
end

# @require FixedSizeArrays begin

  type BezierCurve{T <: FixedSizeArrays.Vec}
      control_points::Vector{T}
  end

  function Base.call(bc::BezierCurve, t::Real)
      p = zero(P2)
      n = length(bc.control_points)-1
      for i in 0:n
          p += bc.control_points[i+1] * binomial(n, i) * (1-t)^(n-i) * t^i
      end
      p
  end

  Base.mean(x::Real, y::Real) = 0.5*(x+y)
  Base.mean{N,T<:Real}(ps::FixedSizeArrays.Vec{N,T}...) = sum(ps) / length(ps)

  curve_points(curve::BezierCurve, n::Integer = 30; range = [0,1]) = map(curve, linspace(range..., n))

  # build a BezierCurve which leaves point p vertically upwards and arrives point q vertically upwards.
  # may create a loop if necessary.  Assumes the view is [0,1]
  function directed_curve(p::P2, q::P2; xview = 0:1, yview = 0:1)
    mn = mean(p, q)
    diff = q - p

    minx, maxx = minimum(xview), maximum(xview)
    miny, maxy = minimum(yview), maximum(yview)
    diffpct = P2(diff[1] / (maxx - minx),
                 diff[2] / (maxy - miny))

    # these points give the initial/final "rise"
    # vertical_offset = P2(0, (maxy - miny) * max(0.03, min(abs(0.5diffpct[2]), 1.0)))
    vertical_offset = P2(0, max(0.15, 0.5norm(diff)))
    upper_control = p + vertical_offset
    lower_control = q - vertical_offset

    # try to figure out when to loop around vs just connecting straight
    # TODO: choose loop direction based on sign of p[1]??
    # x_close_together = abs(diffpct[1]) <= 0.05
    p_is_higher = diff[2] <= 0
    inside_control_points = if p_is_higher
      # add curve points which will create a loop
      sgn = mn[1] < 0.5 * (maxx + minx) ? -1 : 1
      inside_offset = P2(0.3 * (maxx - minx), 0)
      additional_offset = P2(sgn * diff[1], 0)  # make it even loopier
      [upper_control + sgn * (inside_offset + max(0,  additional_offset)),
       lower_control + sgn * (inside_offset + max(0, -additional_offset))]
    else
      []
    end

    BezierCurve([p, upper_control, inside_control_points..., lower_control, q])
  end

# end


text:
# _maxfloat(::Type{T}) where {T<:AbstractFloat} = prevfloat(typemax(T))
# _minfloat(::Type{T}) where {T<:AbstractFloat} = nextfloat(typemin(T))

Base.ndims(::Interval) = 1

function Base.rand(i::Interval{T,L,R}) where {T,L,R}
    # α = max(_minfloat(T), first(i))
    # β = min(_maxfloat(T), last(i))
    μ = maxintfloat(T, Int)
    α = max(-μ, first(i))
    β = min(μ, last(i))
    δ = β - α
    # if δ === Inf
    if δ > μ
        return rand(rand() < 0.5 ? α..zero(T) : zero(T)..β)
    else
        # r = α + exp10(log10(δ) * rand())
        r = α + δ * rand()
        r ∉ i && return rand(i)
        return r
    end
end

# TODO - Optimise the type of Intervals.Bound (currently abstract super type)
mutable struct IntervalsFold{T<:AbstractFloat,L<:Intervals.Bound,R<:Intervals.Bound}
    pattern::Interval{T,L,R}
    gap::T
    folds::Int
    current::Int
end

IntervalsFold(p, g, f, c=1) = IntervalsFold(p, g, f, c)

@forward IntervalsFold.pattern Base.isempty, Base.ndims

Base.lastindex(isf::IntervalsFold) = folds(isf)

Base.getindex(isf::IntervalsFold, key...) = map(k -> get_fold!(isf, k), key)

function pattern(isf::IntervalsFold)
    distortion = gap(isf) * (isf.current - 1)
    return isf.pattern + (-distortion)
end

function unfold(isf::IntervalsFold{T,L,R}) where {T,L,R}
    reset_pattern!(isf)
    x = first(pattern(isf))
    y = last(pattern(isf))
    g = gap(isf)
    f = folds(isf)
    return [Interval{T,L,R}(x + g * i, y + g * i) for i in 0:(f - 1)]
end

function set_fold!(isf::IntervalsFold, new_fold=isf.current + 1)
    if new_fold != isf.current && 0 < new_fold ≤ isf.folds
        distortion = gap(isf) * (new_fold - isf.current)
        isf.pattern += distortion
        isf.current = new_fold
    end
end

function get_fold!(isf::IntervalsFold, f)
    set_fold!(isf, f)
    return isf.pattern
end

function Base.iterate(iter::IntervalsFold)
    reset_pattern!(iter)
    return pattern(iter), 1
end

# "Induction" iterate method
function Base.iterate(iter::IntervalsFold, state::Int)
    state ≥ folds(iter) && return nothing
    set_fold!(iter)
    return iter.pattern, state + 1
end

# Reverse iterate method
function Base.iterate(
    r_iter::Base.Iterators.Reverse{<:IntervalsFold}, state::Int=length(r_iter.itr)
)
    state < 1 && return nothing
    iter = r_iter.itr
    next_state = state - 1
    set_fold!(iter, state)
    return iter.pattern, next_state
end

function Base.in(val, isf::IntervalsFold)
    reset_pattern!(isf)
    return any(i -> val ∈ i, isf)
end

Base.size(isf::IntervalsFold) = span(isf.pattern) * folds(isf)

Base.length(isf::IntervalsFold) = folds(isf)

Base.eltype(::Type{<:IntervalsFold{T,L,R}}) where {T,L,R} = Interval{T,L,R}

is_points(isf) = is_point(pattern(isf))

pattern_length(isf::IntervalsFold) = span(pattern(isf))

function Base.rand(isf::IntervalsFold)
    i = rand(1:folds(isf)) - 1
    p = pattern(isf)
    return rand(p) + i * gap(isf)
end

"""
    Base.rand(::Vector{IntervalsFold})
Extend the `Base.rand` function to `Vector{IntervalsFold}`.
"""
function Base.rand(v::V) where {V <: Set{<:IntervalsFold}}
    λ = map(size, collect(v))
    Λ = sum(λ)
    σ = [sum(λ[1:i]) / Λ for i in 1:length(v)]
    r = rand()
    for z in zip(σ, v)
        r < z[1] && return rand(z[2])
    end
end


text:
#################### Basis for Polynomial Space ####################3
struct PolynomialBasis{T}
  order::Int
  nodes::Vector{T}      #Quadrature nodes for numerical integration
  weights::Vector{T}    #Quadrature weights for numerical integration
  polynomials::Vector{Poly}
  φ::Matrix{T}     #Vandermonde matrix: basis polynomials evaluated on G-L nodes
  dφ::Matrix{T}
  invφ::Matrix{T}   #inverse of Gen. Vandermonde matrix
end

"""compute Legendre polynomials coefficients, normalized to be orthonormal"""
function poly_legendre{T<:Number}(n, ::Type{T}=Float64, var=:x)
    return poly_jacobi(n,0.0,0.0,T,var)
end

"""compute Jacobi polynomials coefficients, normalized to be orthonormal"""
function poly_jacobi{T<:Number}(n, a, b, ::Type{T}=Float64, var=:x)
    ox = one(T)
    zx = zero(T)
    #Compute initial P_0 and P_1
    γ0 = 2^(a+b+1)/(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+1);
    p0 = Poly{T}([one(T)/sqrt(γ0)], var)
    if n==0; return p0; end
    γ1 = (a+1)*(b+1)/(a+b+3)*γ0
    p1 = Poly{T}([(a-b)/2/sqrt(γ1), (a+b+2)/2/sqrt(γ1)], var)
    if n==1; return p1; end
    px = Poly{T}([zero(T), one(T)], var)
    aold = 2/(2+a+b)*sqrt((a+1)*(b+1)/(a+b+3))
    for i = 1:(n-1)
        h1 = 2*i+a+b;
        anew = 2/(h1+2)*sqrt((i+1)*(i+1+a+b)*(i+1+a)*(i+1+b)/(h1+1)/(h1+3))
        bnew = -(a^2-b^2)/h1/(h1+2);
        p2 = ox/anew*(-aold*p0 + (px-bnew)*p1);
        aold =anew;
        p0 = p1
        p1 = p2
    end
    return p1
end

function legendre_basis{T<:Number}(order, ::Type{T}=Float64)
  nodes, weights = gausslobatto(order+1)
  φ = zeros(T,order+1,order+1)
  dφ = zeros(T,order+1,order+1)
  polynomials = Vector{Poly}(order+1)
  for n = 0:order
    p = poly_legendre(n, T)
    dp = polyder(p)
    polynomials[n+1] = p
    # Eval interior nodes
    φ[:,n+1] = polyval(p, nodes)
    dφ[:,n+1] = polyval(dp, nodes)
  end
  invφ = inv(φ)
  PolynomialBasis{T}(order,nodes,weights,polynomials,φ,dφ,invφ)
end

"Maps reference coordinates (ξ ∈ [-1,1]) to interval coordinates (x)"
function reference_to_interval(ξ,a::Tuple)
   0.5*(a[2]-a[1])*ξ + 0.5*(a[2]+a[1])
end

"Project function f on polynomial space Vₕ"
function project_function(f, basis, interval::Tuple; component=1)
  nodes = reference_to_interval(basis.nodes, interval)
  f_val = zeros(nodes)
  for i in 1:size(nodes,1)
    f_val[i] = f(nodes[i])[component]
  end
  function model(x,p)
    result = zeros(x)
    for i in 1:size(p,1)
      result.+=p[i]*polyval(basis.polynomials[i], x)
    end
    result
  end
  p0 = zeros(eltype(basis.nodes),basis.order+1)
  curve_fit(model, basis.nodes, f_val, p0)
end

# using FEMBasis
# import FEMBasis: create_basis
# a = legendre_basis(2)
#
# code = create_basis(
# :LegSeg3,
# "3 node linear segment/line orthogonal legendre element",
# NTuple((
#  (-0.7745966692414834,), # N1
#  ( 0.0, ), # N2
#  ( 0.7745966692414834, ), # N3
# )),(
#  :(0.5 + 0.5),
#  :(1.0*u),
#  :(-0.5 + 1.5*u^2),)
# )
# eval(code)
#
# eval_basis!(LegSeg3(),zeros(3,1),(-0.7745966692414834))
# eval_dbasis!(LegSeg3(),zeros(1,3),(-0.7745966692414834))

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
# This test harness executes the specified BENCH_FILE several times, each with different
# JULIA_NUM_THREADS settings, to measure performance scaling of the benchmark with more
# threads. This is needed because currently Julia does not allow dynamically changing the
# number of available threads.
# To change which benchmark is being run, set the BENCH_FILE constant below.

using DataFrames
using Plots
using Plots.PlotMeasures
using Statistics

# ----------------------------------------------------------------
# Params for this test file -- EDIT THESE
const NUM_CORES = Sys.CPU_THREADS

function run_bench_file(file, nthreads, nqueries, query_size)
    withenv(
        "JULIA_NUM_THREADS" => nthreads,
        "NUM_QUERIES" => nqueries,
        "QUERY_SIZE" => query_size,
    ) do
         @assert Base.active_project() != nothing "Must run this benchmark w/ an active Julia Project.toml (and ideally Manifest.toml)!"
         lines = readlines(`$(Base.julia_cmd()) --project=$(Base.active_project()) $file`)
         return (;
             time_secs            = parse(Float64, lines[end-5]) / 1e9            ,
             allocs               = parse(Int,     lines[end-4])                  ,
             memory               = parse(Int,     lines[end-3])                  ,
             gctime_secs          = parse(Float64, lines[end-2]) / 1e9            ,
             thread_counts        = Vector{Int}(Meta.parse(lines[end-1]).args)    ,
             query_latencies_secs = Vector{Float64}(Meta.parse(lines[end]).args)  ,
         )
    end
end

function bench_nthreads_scaling(file, nqueries, query_size, all_nthreads, preprocess_results)
    println("NUM_QUERIES: $nqueries")
    println("WORK_SIZE: $query_size")
    results = DataFrame(nthreads=Int[], time_secs=Float64[], allocs=Int[], memory=Int[], gctime_secs=Float64[],
                                        thread_counts=Vector{Int}[], query_latencies_secs=Vector{Float64}[])
    for NTHREADS in all_nthreads
        println("$NTHREADS threads: ")
        benchresults = run_bench_file(file, NTHREADS, nqueries, query_size)
        push!(results, (; nthreads = NTHREADS, benchresults...))
        println("$(benchresults.time_secs) secs")
    end

    # ------------
    # Run provided preprocess_results()
    # ------------
    results = preprocess_results(results)

    println("Results (omitted printing latencies column):")
    println(results[:, 1:end-1])  # Pretty-print and hide latencies column (too long)

    #-------------
    # Plot the results
    #-------------

    df = DataFrame(nthreads = results.nthreads, time_secs = results.time_secs,
            speedup_factor = [results[1,:].time_secs / t for t in results.time_secs])

    utilization = df.speedup_factor ./ df.nthreads
    marginal_speedup = let out = Float64[]
         push!(out, 1.0)
         for i in 2:size(df, 1)
            push!(out, (df.speedup_factor[i] - df.speedup_factor[i-1]) ./ (df.nthreads[i] - df.nthreads[i-1]))
         end
         out
     end
     query_latencies_ms = sort.(results.query_latencies_secs .* 1e3)
     quantiles = (0.1, 0.5, 0.9, 0.99, 1)
     processed = DataFrame(nthreads = df.nthreads, time_secs = df.time_secs,
                           speedup_factor = df.speedup_factor,
                           marginal_speedup = marginal_speedup, utilization = utilization,
                           latency_quantiles_ms = [((q=>quantile(l, q) for q in quantiles)...,) for l in query_latencies_ms],
                           )

    println("Processed:")
    println(processed)  # Print the full processed results DataFrame
    return results, processed
end

plot_speedup(args...) = plot_speedup!(plot(), args...)
function plot_speedup!(p, processed, label)
    Plots.plot!(p, processed.nthreads, processed.speedup_factor,
        xticks = processed.nthreads,
        yticks = processed.nthreads,
        ylims = (1, min(Sys.CPU_THREADS, maximum(processed.nthreads))),
        xlabel = "JULIA_NUM_THREADS",
        ylabel = "speedup factor",
        label = label,
        left_margin=50px,
        bottom_margin=50px
    )
end
plot_abs_time(args...) = plot_abs_time!(plot(), args...)
function plot_abs_time!(p, processed, label)
    yticks = range(minimum(processed.time_secs), maximum(processed.time_secs), length=6)
    Plots.plot!(p, processed.nthreads, processed.time_secs,
        xticks = processed.nthreads,
        yticks = (yticks, ["$(Int(round(t))) secs" for t in yticks]),
        #yscale = :log10,
        xlabel = "JULIA_NUM_THREADS",
        #ylabel = "absolute time (secs)",
        label = label,
        left_margin=50px,
        bottom_margin=50px
    )
end
scalevec(v) = v ./ sum(v)
function plot_thread_counts(results, label)
    numrows = size(results, 1)
    range = numrows÷2 : numrows÷2 + 2
    normalized_tc = scalevec.(results.thread_counts[range])
    Plots.bar(normalized_tc,
        xlabel = "Thread ID",
        ylabel = "Num Queries",
        fillalpha = 0.2,
        title = label,
        labels = ["$(results.nthreads[r]) threads" for r in range],
        left_margin=50px,
        bottom_margin=50px
    )
end
function plot_latencies(results, label)
    nqueries = length(results.query_latencies_secs[1])
    xticks_labels = 0:0.10:1
    Plots.plot(sort.(results.query_latencies_secs),
        xticks = (xticks_labels * nqueries, ["$(Int(round(x*100)))%" for x in xticks_labels]),
        xlabel = "Percentile",
        yscale = :log10,
        labels = ["$t threads" for t in results.nthreads],
        title = "Query latencies (secs, log): \"$label\"",
        left_margin=50px,
        bottom_margin=50px
    )
end

function perform_scaling_experiment(;
        bench_file, nqueries, num_ops, plot_series_name, num_datapoints,
        # Test from 1 ... NUM_CORES, because julia doesn't let you go over that maximum
        nthreads_to_test = [1, Int.(round.(range(0, stop=NUM_CORES, length=num_datapoints)))[2:end]...],
        preprocess_results = identity,
        )

    println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
    println("Running benchmark for $bench_file")
    println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

    println("------ TEST PARAMETERS ------")
    @show nthreads_to_test
    @show nqueries
    @show num_ops

    println("warmup")
    bench_nthreads_scaling(bench_file, 1,  10, [5,6], # Warmup w/ nthreads > 1 (some tests require it)
                           preprocess_results);

    println("run benchmark")
    global results,processed = bench_nthreads_scaling(bench_file, nqueries,  num_ops,
                                                      nthreads_to_test, preprocess_results)

    plot_basename = splitext(basename(bench_file))[1]

    println("plot results")
    println("-- saving figures in $(pwd()) --")
    pp = plot_speedup(processed, plot_series_name)
    savefig(pp, "$plot_basename-speedup_plot.png")

    ap = plot_abs_time(processed, plot_series_name)
    savefig(ap, "$plot_basename-abs_time_plot.png")

    tp = plot_thread_counts(results, plot_series_name)
    savefig(tp, "$plot_basename-thread_counts.png")

    lp = plot_latencies(results, plot_series_name)
    savefig(lp, "$plot_basename-query_latencies.png")
end


text:
using ARCTools
using Test

@testset "ARCTools.jl" begin
    q = QSub(;
        resources = Resources(
            walltime = "00:05:00",
            select = 2,
            ncpus = 3,
            mem = 5,
        ),
        modules = [
            "gcc/9.1.0",  # for git/2.21.0
            "git/2.21.0", # for julia 1.6+
            "python/3.7.3",
            "cuda/10.0.130",
        ],
        account = "st-arausch-1",
        jobname = "test-job",
    )

    jl = JuliaCmd(;
        bindir = Sys.BINDIR,
        threads = q.resources.ompthreads,
        project = ARCTools.copy_project(;
            project = joinpath(@__DIR__, "TestProj"),
            dest = mktempdir(
                "/scratch/st-arausch-1/jcd1994";
                prefix = "TestProj_",
            ),
        ),
        script = "./scripts/test_script.jl",
        args = ["Alice", "Bob", "Carol"],
        env = Dict{String, Any}(
            "JL_SECRET_VAL1" => "SECRET1",
            "JL_SECRET_VAL2" => "SECRET2",
            "JL_PUBLIC_VAL" => "PUBLIC",
        ),
        secrets = String[
            "JL_SECRET_VAL1",
            "JL_SECRET_VAL2",
        ],
    )

    run(qsub(q, jl))
end


text:
# Stencil: applies a space-invariant, linear, symmetric filter (stencil)
# to a square grid.
#
# Include this file after setting `Context` to one of the `ClimaComms.jl`
# backends' context types.

using CUDA
using Logging
using OffsetArrays
using Printf

using ClimaComms

# defaults for grid dimension, stencil radius, etc.
const default_n = 256
const default_radius = 2
const default_niterations = 50

# floating point precisions to test with
const default_float_types = (Float64, Float32)

# array types with which to test
#const default_array_types = CUDA.has_cuda() ? (CuArray, Array) : (Array,)
const default_array_types = (Array,)

puts(s...) = ccall(:puts, Cint, (Cstring,), string(s...))

epsilon(::Type{Float64}) = 1.0e-8
epsilon(::Type{Float32}) = 0.0001f0
coefx(::Type{Float64}) = 1.0
coefx(::Type{Float32}) = 1.0f0
coefy(::Type{Float64}) = 1.0
coefy(::Type{Float32}) = 1.0f0

# find 2 factors of `n` closest to `√n`
function factor(n::Int)
    f1 = floor(Int, sqrt(n + 1.0))
    f2 = 0
    while f1 > 0
        if n % f1 == 0
            f2 = convert(Int, n / f1)
            break
        end
        f1 -= 1
    end
    return f1, f2
end

# partition `n` into `nr` parts and return the size and the offsets
# of pid `rdim`s part
function compute_dim(n::Int, nr::Int, rdim::Int)
    dim, dimleft = divrem(n, nr)
    if rdim < dimleft
        ostart = (dim + 1) * rdim + 1
        oend = ostart + dim
    else
        ostart = (dim + 1) * dimleft + dim * (rdim - dimleft) + 1
        oend = ostart + dim - 1
    end
    dim = oend - ostart + 1
    return dim, ostart, oend
end

# set up stencil weights to reflect a discrete divergence operator
function setup_stencil(AT, FT, radius)
    weight = zeros(FT, 2radius + 1, 2radius + 1)
    oweight = OffsetArray(weight, (-radius):radius, (-radius):radius)
    for i in 1:radius
        oweight[0, i] = oweight[i, 0] = 1.0 / (2.0 * i * radius)
        oweight[0, -i] = oweight[-i, 0] = -1.0 / (2.0 * i * radius)
    end
    return OffsetArray(AT(weight), (-radius):radius, (-radius):radius)
end

# set up input array with ghost elements
function setup_input(AT, FT, height, cstart, cend, width, rstart, rend, radius)
    input = zeros(FT, height + 2radius, width + 2radius)
    oinput = OffsetArray(
        input,
        (rstart - radius):(rend + radius),
        (cstart - radius):(cend + radius),
    )
    @inbounds for c in cstart:cend
        @simd for r in rstart:rend
            oinput[r, c] = coefx(FT) * (c - 1) + coefy(FT) * (r - 1)
        end
    end
    return OffsetArray(
        AT(input),
        (rstart - radius):(rend + radius),
        (cstart - radius):(cend + radius),
    )
end

# fill the send stage buffer for the specified neighbor from `input`
function fill_send_buf(nbr, input, rrange, crange)
    if ClimaComms.pid(nbr) > 0
        buf = ClimaComms.send_stage(nbr)
        buf .= view(input, rrange, crange)
    end
end

# copy the receive stage buffer for the specified neighbor into `input`
function empty_recv_buf(nbr, input, rrange, crange)
    if ClimaComms.pid(nbr) > 0
        buf = ClimaComms.recv_stage(nbr)
        view(input, rrange, crange) .= buf
    end
end

# apply the stencil operator
# TODO: ensure this uses FMAs
function stencil(out, weight, radius, input, r, c)
    @inbounds begin
        @simd for ii in (-radius):radius
            out[r, c] += weight[ii, 0] * input[r + ii, c]
        end
        @simd for jj in (-radius):-1
            out[r, c] += weight[0, jj] * input[r, c + jj]
        end
        @simd for jj in 1:radius
            out[r, c] += weight[0, jj] * input[r, c + jj]
        end
    end
end

function stencil_test(
    Context::Type{<:ClimaComms.AbstractCommsContext};
    array_types = default_array_types,
    float_types = default_float_types,
    n = default_n,
    radius = default_radius,
    niterations = default_niterations,
)
    stencil_size = 4radius + 1

    # initialize and get processor info
    pid, nprocs = ClimaComms.init(Context)

    # log output only from pid 0
    logger_stream = ClimaComms.iamroot(Context) ? stderr : devnull
    prev_logger = global_logger(ConsoleLogger(logger_stream, Logging.Info))
    atexit() do
        global_logger(prev_logger)
    end

    # determine 2D grid of pids (closest to square)
    nprocsx, nprocsy = factor(nprocs)
    pidy, pidx = divrem(pid - 1, nprocsx)

    # determine neighbors
    left_nbr = pidx > 0 ? pid - 1 : -1
    right_nbr = pidx < nprocsx - 1 ? pid + 1 : -1
    top_nbr = pidy < nprocsy - 1 ? pid + nprocsx : -1
    bottom_nbr = pidy > 0 ? pid - nprocsx : -1

    # determine dimensions of input and output arrays
    height, rstart, rend = compute_dim(n, nprocsx, pidx)
    width, cstart, cend = compute_dim(n, nprocsy, pidy)
    active_points = (n - 2radius) * (n - 2radius)

    # start tests
    for AT in array_types, FT in float_types
        weight = setup_stencil(AT, FT, radius)
        input = setup_input(
            AT,
            FT,
            height,
            cstart,
            cend,
            width,
            rstart,
            rend,
            radius,
        )
        outa = AT{FT}(undef, height, width)
        fill!(outa, zero(FT))
        out = OffsetArray(outa, rstart:rend, cstart:cend)

        # set up communication buffers for neighbors
        make_neighbor(pid, send_dims, recv_dims = send_dims) =
            ClimaComms.Neighbor(Context, pid, AT, FT, send_dims, recv_dims)
        neighbors = (;
            left = make_neighbor(left_nbr, (radius, width)),
            right = make_neighbor(right_nbr, (radius, width)),
            top = make_neighbor(top_nbr, (height, radius)),
            bottom = make_neighbor(bottom_nbr, (height, radius)),
        )

        # create ClimaComms context
        comms_ctx = Context(
            filter(nbr -> ClimaComms.pid(nbr) > 0, [nbr for nbr in neighbors]),
        )

        # compute loop
        local_stencil_time = 0
        for iter in 0:niterations
            # start timing after a warmup iteration
            if iter == 1
                ClimaComms.barrier(comms_ctx)
                local_stencil_time = time()
            end

            # fill send buffers
            fill_send_buf(
                neighbors.left,
                input,
                rstart:(rstart + radius - 1),
                cstart:cend,
            )
            fill_send_buf(
                neighbors.right,
                input,
                (rend - radius + 1):rend,
                cstart:cend,
            )
            fill_send_buf(
                neighbors.top,
                input,
                rstart:rend,
                (cend - radius + 1):cend,
            )
            fill_send_buf(
                neighbors.bottom,
                input,
                rstart:rend,
                cstart:(cstart + radius - 1),
            )

            # initiate communication
            ClimaComms.start(comms_ctx)

            # apply the stencil operator to the interior while waiting
            # for communication
            @inbounds for c in (cstart + radius):(cend - radius)
                for r in (rstart + radius):(rend - radius)
                    stencil(out, weight, radius, input, r, c)
                end
                ClimaComms.progress(comms_ctx)
            end

            # complete communication
            ClimaComms.finish(comms_ctx)

            # move receive buffers to the input array
            empty_recv_buf(
                neighbors.left,
                input,
                (rstart - radius):(rstart - 1),
                cstart:cend,
            )
            empty_recv_buf(
                neighbors.right,
                input,
                (rend + 1):(rend + radius),
                cstart:cend,
            )
            empty_recv_buf(
                neighbors.top,
                input,
                rstart:rend,
                (cend + 1):(cend + radius),
            )
            empty_recv_buf(
                neighbors.bottom,
                input,
                rstart:rend,
                (cstart - radius):(cstart - 1),
            )

            # now we can apply the stencil operator to the exterior
            @inbounds for c in
                          max(cstart, radius + 1):min(cend - radius, n - radius)
                for r in max(rstart, radius + 1):(rstart + radius - 1)
                    stencil(out, weight, radius, input, r, c)
                end
                for r in (rend - radius + 1):min(rend, n - radius)
                    stencil(out, weight, radius, input, r, c)
                end
            end
            @inbounds for c in max(cstart, radius + 1):(cstart + radius - 1)
                for r in max(rstart, radius + 1):(rend - radius)
                    stencil(out, weight, radius, input, r, c)
                end
            end
            @inbounds for c in (cend - radius + 1):min(cend, n - radius)
                for r in max(rstart, radius + 1):(rend - radius)
                    stencil(out, weight, radius, input, r, c)
                end
            end

            # add constant to solution to force refresh of neighbor data, if any
            view(input, rstart:rend, cstart:cend) .+= 1.0
        end

        # get time
        local_stencil_time = (time() - local_stencil_time)
        stencil_time = ClimaComms.reduce(comms_ctx, local_stencil_time, max)

        # compute L1 norm
        local_norm = zero(FT)
        for c in max(cstart, radius + 1):min(cend, n - radius)
            for r in max(rstart, radius + 1):min(rend, n - radius)
                local_norm += abs(out[r, c])
            end
        end
        norm = ClimaComms.reduce(comms_ctx, local_norm, +)

        # verify correctness
        if ClimaComms.iamroot(Context)
            norm /= active_points
            refnorm = (niterations + 1) * (coefx(FT) + coefy(FT))
            if abs(norm - refnorm) > epsilon(FT)
                @error "$AT{$FT}: norm = $(norm), reference = $(refnorm)\n"
                # TODO: uncomment this when ClimaCommsSA is fixed
                # see: https://github.com/CliMA/ClimaComms.jl/issues/6
                #ClimaComms.abort(comms_ctx, -1)
            end

            # flops/stencil: 2 flops (fma) for each point in the stencil
            # plus one flop for the update of the input array
            flops = (2 * stencil_size + 1) * active_points
            avgtime = stencil_time / niterations
            @info @sprintf "%f MFlops/s\n" 1.0e-6 * flops / avgtime
        end
    end

    return nothing
end


text:
export
    MultiScalarRule,
    dnaSize,
    resultSize,
    initialize!,
    decode,
    asValue

immutable MultiScalarRule <: SynthRule
    count::Int
    min::Float32
    max::Float32
    round::Bool
    variationCount::Int
    dnaSize::Int

    function MultiScalarRule(;count = 1, min = 0f0, max = 1f0, round = false, variationCount = 2)
        dnaSize = variationCount * count * 2

        new(count, min, max, round, variationCount, dnaSize)
    end
end

dnaSize(rule::MultiScalarRule) = rule.dnaSize

resultSize(rule::MultiScalarRule) = rule.count

function initialize!(rule::MultiScalarRule, ctx, state)
    state[:s] = range(ctx, UInt32, [rule.count])
end

function decode(rule::MultiScalarRule, pars)
    scope!(pars.ctx) do this
        const depth = rule.variationCount
        const blockSize = rule.count * depth
        values = moddims(pars.dnaFragment[seq(pars.ctx, 0, blockSize - 1)], rule.count, depth)
        dist = moddims(pars.dnaFragment[seq(pars.ctx, blockSize, blockSize * 2 - 1)], rule.count, depth)
        rnd = randu(pars.ctx, Float32, rule.count, depth)
        test = rnd ./ maxOf(dist, 0.000001f0)
        max, idx = imax(test, 1)
        idx2 = idx .* rule.count .+ pars.state[:s]
        result = values[idx2]

        if rule.min == 0f0 && rule.max == 1f0
            pars.result[pars.resultSeq] = rule.round ? round(result) : result
        else
            const d = rule.max - rule.min
            result = result .* d .+ rule.min
            pars.result[pars.resultSeq] = rule.round ? round(result) : result
        end
    end
end

function asValues(rule::MultiScalarRule, values, startIndex)
    result = Vector(resultSize(rule))
    for i in 1:resultSize(rule)
        result[i] = values[startIndex - 1 + i]
    end
    result
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# Knots

struct Knots
    vector::Vector{Float64}
    @doc raw"""
    Construct knot vector from given array.
    ```math
    k=(k_1,\dots,k_l)
    ```

    # Examples
    ```jldoctest
    julia> k = Knots([1,2,3])
    Knots([1.0, 2.0, 3.0])

    julia> k = Knots(1:3)
    Knots([1.0, 2.0, 3.0])
    ```
    """
    function Knots(vector::AbstractVector{<:Real})
        return new(sort(convert(Vector{Float64}, vector)))
    end
end

function Knots(vector::AbstractVector)
    if isempty(vector)
        @warn "use Knots() or Knots(Float64[]), instead of Knots([])."
        return Knots(Float64[])
    else
        @warn "the type $(typeof(vector)) is deprecated."
        return Knots(convert(Vector{Float64}, vector))
    end
end

@doc raw"""
Construct knot vector from given real numbers.

# Examples
```jldoctest
julia> k = Knots(1,2,3)
Knots([1.0, 2.0, 3.0])

julia> k = Knots()
Knots([])
```
"""
function Knots(knot::Real...)
    return Knots(collect(knot))
end

function Base.show(io::IO, k::Knots)
    if k.vector == Float64[]
        print(io, "Knots([])")
    else
        print(io, "Knots($(k.vector))")
    end
end

Base.zero(::Type{Knots}) = Knots(Float64[])
Base.:(==)(k₁::Knots, k₂::Knots) = (k₁.vector == k₂.vector)

@doc raw"""
Sum of knot vectors

```math
\begin{aligned}
k^{(1)}+k^{(2)}
&=(k^{(1)}_1, \dots, k^{(1)}_{l^{(1)}}) + (k^{(2)}_1, \dots, k^{(2)}_{l^{(2)}}) \\
&=(\text{sort of union of} \  k^{(1)} \ \text{and} \  k^{(2)} \text{)}
\end{aligned}
```

For example, ``(1,2,3,5)+(4,5,8)=(1,2,3,4,5,5,8)``.

# Examples
```jldoctest
julia> k1 = Knots(1,2,3,5);

julia> k2 = Knots(4,5,8);

julia> k1 + k2
Knots([1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 8.0])
```
"""
Base.:+(k₁::Knots, k₂::Knots) = Knots(sort([k₁.vector..., k₂.vector...]))

@doc raw"""
Product of integer and knot vector

```math
\begin{aligned}
m\cdot k&=\underbrace{k+\cdots+k}_{m}
\end{aligned}
```

For example, ``2\cdot (1,2,2,5)=(1,1,2,2,2,2,5,5)``.

# Examples
```jldoctest
julia> k = Knots(1,2,2,5);

julia> 2 * k
Knots([1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 5.0, 5.0])
```
"""
function Base.:*(p₊::Integer, k::Knots)
    if p₊ == 0
        return zero(Knots)
    elseif p₊ > 0
        return sum(k for _ in 1:p₊)
    else
        throw(DomainError(p₊, "additional degree of polynominal must be non-negative"))
    end
end

Base.in(r::Real, k::Knots) = in(r, k.vector)
Base.getindex(k::Knots, i::Integer) = k.vector[i]
Base.getindex(k::Knots, v::AbstractVector{<:Integer}) = Knots(k.vector[v])

@doc raw"""
Length of knot vector

# Examples
```jldoctest
julia> k = Knots(1,2,3,5);

julia> length(k)
4
```
"""
Base.length(k::Knots) = length(k.vector)

function ♯(k::Knots)
    @warn "BasicBSpline.♯ is deprecated. Use `length` instead."
    return length(k::Knots)
end
Base.firstindex(k::Knots) = 1
Base.lastindex(k::Knots) = length(k)

@doc raw"""
Unique elements of knot vector.

```math
\begin{aligned}
\widehat{k}
&=(\text{remove duplicates of} \  k) \\
\end{aligned}
```

For example, ``\widehat{(1,2,2,3)}=(1,2,3)``.

# Examples
```jldoctest
julia> k = Knots([1,2,2,3]);

julia> unique(k)
Knots([1.0, 2.0, 3.0])
```
"""
Base.unique(k::Knots) = Knots(unique(k.vector))
Base.iterate(k::Knots) = iterate(k.vector)
Base.iterate(k::Knots, i::Integer) = iterate(k.vector, i)
Base.searchsortedfirst(k::Knots,t) = searchsortedfirst(k.vector,t)
Base.searchsortedlast(k::Knots,t) = searchsortedlast(k.vector,t)
Base.searchsorted(k::Knots,t) = searchsorted(k.vector,t)

@doc raw"""
Check a inclusive relation ship ``k\subset k'``.
```math
(1,2,3) \subseteq (1,2,3,4)
(1,2,5) \nsubseteq (1,2,3,4)
```
"""
function Base.issubset(k::Knots, k′::Knots)
    v = k′.vector
    l = length(v)
    i = 0
    for kᵢ in k
        i = findnext(==(kᵢ), v, i+1)
        if isnothing(i)
            return false
        end
    end
    return true
end

@doc raw"""
For Given knot vector ``k``, the following function ``\mathfrak{n}_k:\mathbb{R}\to\mathbb{Z}`` represents the number of knots that duplicate the knot vector ``k``.

```math
\mathfrak{n}_k(t) = \sharp\{i \mid k_i=t \}
```
For example, if ``k=(1,2,2,3)``, then ``\mathfrak{n}_k(0.3)=0``, ``\mathfrak{n}_k(1)=1``, ``\mathfrak{n}_k(2)=2``.

```jldoctest
julia> k = Knots([1,2,2,3]);

julia> 𝔫(k,0.3)
0

julia> 𝔫(k,1.0)
1

julia> 𝔫(k,2.0)
2
```
"""
function 𝔫(k::Knots, t::Real)
    # for small case, this is faster
    # return count(==(t), k.vector)

    # for large case, this is faster
    return length(searchsorted(k,t))
end

"""
Find an index ``i`` such that ``k_{i} ≤ t < k_{i+1}``.
"""
function _knotindex₊₀(k::Union{Knots, AbstractVector{<:Real}}, t::Real)
    return searchsortedlast(k, t)
end

"""
Find an index ``i`` such that ``k_{i} < t ≤ k_{i+1}``.
"""
function _knotindex₋₀(k::Union{Knots, AbstractVector{<:Real}}, t::Real)
    return searchsortedfirst(k, t) - 1
end

"""
Find an index ``i`` such that ``k_{i} ≤ t < k_{i+1}`` or ``k_{i} < t = k_{i+1} = k_{\text{end}}``.
"""
function _knotindex(k::Union{Knots, AbstractVector{<:Real}}, t::Real)
    j = searchsortedlast(k, t)
    if j < length(k)
        return j
    else
        return searchsortedfirst(k, t) - 1
    end
end