text:
module ProblemBenchmarks

# This module contains benchmarks that test against general problem cases taken
# from real-world examples (as opposed to microbenchmarks, language-agnostic
# benchmark suites, or benchmarks that stress specific implementation details).

# A lot of the benchmarks here originated from JuliaLang/julia/test/perf/kernel,
# where much of the code was naively translated to Julia from other languages,
# and thus is written non-idiomatically.

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

const SUITE = BenchmarkGroup(["example", "kernel"])
const PROBLEM_DATA_DIR = joinpath(dirname(@__FILE__), "data")

#######################################
# IMDB Actor Centrality (Issue #1163) #
#######################################

include("IMDBGraphs.jl")

g = addgroup!(SUITE, "imdb", ["graph", "centrality"])

g["centrality"] = @benchmarkable IMDBGraphs.perf_imdb_centrality(50)

###############
# Monte Carlo #
###############

include("MonteCarlo.jl")

g = addgroup!(SUITE, "monte carlo", ["finance", "vectorization", "random"])

g["euro_option_devec"] = @benchmarkable MonteCarlo.perf_euro_option_devec(10^4)
g["euro_option_vec"] = @benchmarkable MonteCarlo.perf_euro_option_vec(10^4)

###################################
# Laplacian (Issues #1168, #4707) #
###################################

include("Laplacian.jl")

g = addgroup!(SUITE, "laplacian", ["iterative", "sparse", "vectorization", "subarray", "array"])

g["laplace_sparse_matvec"] = @benchmarkable Laplacian.perf_laplace_sparse_matvec(8^5)
g["laplace_iter_devec"] = @benchmarkable Laplacian.perf_laplace_iter_devec(8^2)
g["laplace_iter_vec"] = @benchmarkable Laplacian.perf_laplace_iter_vec(8^2)
g["laplace_iter_sub"] = @benchmarkable Laplacian.perf_laplace_iter_sub(8^2)

###################################################
# Grigoriadis Khachiyan Matrix Games (Issue #950) #
###################################################

include("GrigoriadisKhachiyan.jl")

g = addgroup!(SUITE, "grigoriadis khachiyan", ["grigoriadis", "khachiyan", "game"])

g["grigoriadis_khachiyan"] = @benchmarkable GrigoriadisKhachiyan.perf_gk(350, [0.1])

####################################
# Go Game Simulation (Issue #1169) #
####################################

include("GoGame.jl")

g = addgroup!(SUITE, "go", ["game"])

g["go_game"] = @benchmarkable GoGame.perf_go_game(10)

################
# JSON Parsing #
################

include("JSONParse.jl")

g = addgroup!(SUITE, "json", ["parse", "closure"])

jstr = read(joinpath(PROBLEM_DATA_DIR, "test.json"), String)

g["parse_json"] = @benchmarkable JSONParse.perf_parse_json($(jstr))

############################
# Raytracing (Issue #3811) #
############################

include("Raytracer.jl")

g = addgroup!(SUITE, "raytrace")

g["raytrace"] = @benchmarkable Raytracer.perf_raytrace(5, 256, 4)

############################################
# Correlated Asset Simulation (Issue #445) #
############################################

include("StockCorr.jl")

g = addgroup!(SUITE, "stockcorr", ["finance"])

g["stockcorr"] = @benchmarkable StockCorr.perf_stockcorr()

#########################
# Simplex (Issue #3142) #
#########################

include("Simplex.jl")

g = addgroup!(SUITE, "simplex")

g["simplex"] = @benchmarkable Simplex.perf_simplex()

####################################################
# Ziggurat Gaussian Number Generator (Issue #1211) #
####################################################

include("Ziggurat.jl")

g = addgroup!(SUITE, "ziggurat")

g["ziggurat"] = @benchmarkable Ziggurat.perf_ziggurat(10^6)

######################
# Seismic Simulation #
######################

include("SeismicSimulation.jl")

g = addgroup!(SUITE, "seismic", ["simulation", "simd"])

for T in (Float32, Float64)
    g["seismic", string(T)] = @benchmarkable SeismicSimulation.perf_seismic_sim($T)
end

############################
# Sparse FEM (Issue #9668) #
############################

include("SparseFEM.jl")

g = addgroup!(SUITE, "fem", ["sparse"])

g["sparse_fem"] = @benchmarkable SparseFEM.perf_sparse_fem(256)

###############
# Spell Check #
###############

include("SpellCheck.jl")

g = addgroup!(SUITE, "spellcheck", ["spell", "string"])

g["spellcheck"] = @benchmarkable SpellCheck.perf_spellcheck()

#########
# Chaos #
#########

include("ChaosGame.jl")
g = addgroup!(SUITE, "chaosgame")
g["chaosgame"] = @benchmarkable ChaosGameFractals.perf_chaos()

end # module


text:
module IoTests

using Test

@testset "Interfaces" begin include("IoInterfacesTests.jl") end

@testset "Json" begin include("JsonTests.jl") end

@testset "Bson" begin include("BsonTests.jl") end

@testset "JLD2" begin include("JLD2Tests.jl") end

end # module

ext: jl


text:
using DynamicGrids, DynamicGridsGtk, Gtk, Test, Colors, ColorSchemes

# life glider sims

init =  Bool[0 0 0 0 0 0
             0 0 0 0 0 0
             0 0 0 0 0 0
             0 0 0 1 1 1
             0 0 0 0 0 1
             0 0 0 0 1 0]

test =  Bool[0 0 0 0 0 0
             0 0 0 0 0 0
             0 0 0 0 1 1
             0 0 0 1 0 1
             0 0 0 0 0 1
             0 0 0 0 0 0]

test2 = Bool[0 0 0 0 0 0
             0 0 0 0 0 0
             1 0 0 0 1 1
             1 0 0 0 0 0
             0 0 0 0 0 1
             0 0 0 0 0 0]

@testset "Simulation" begin

    ruleset = Ruleset(Life(); boundary=Wrap())

  @testset "GtkOutput stored simulation matches glider behavior" begin
        # TODO test Gtk canvas image
        # g0 = RGB24(0)
        # g1 = RGB24(1)
        # grey2 = [g0 g0 g0 g0 g0 g0
        #          g0 g0 g0 g0 g0 g0
        #          g1 g0 g0 g0 g1 g1
        #          g1 g0 g0 g0 g0 g0
        #          g0 g0 g0 g0 g0 g1
        #          g0 g0 g0 g0 g0 g0]

        # l0 = get(ColorSchemes.leonardo, 0)
        # l1 = get(ColorSchemes.leonardo, 1)

        # leonardo2 = [l0 l0 l0 l0 l0 l0
        #              l0 l0 l0 l0 l0 l0
        #              l1 l0 l0 l0 l1 l1
        #              l1 l0 l0 l0 l0 l0
        #              l0 l0 l0 l0 l0 l1
        #              l0 l0 l0 l0 l0 l0]

        output = GtkOutput(init; tspan=1:2, store=true)
        sim!(output, ruleset)
        resume!(output, ruleset; tstop=5)
        @test output[3] == test
        @test output[5] == test2
        # TODO @test the canvaas images == leonardo2
        destroy(output.window)
        @testset "display" begin
            @test DynamicGridsGtk._isalive(output) == false
            display(output)
            @test DynamicGridsGtk._isalive(output) == true
            destroy(output.window)
            @test DynamicGridsGtk._isalive(output) == false
        end

    end
end

@testset "Float output" begin

    flt = [0.0 0.0 0.0 0.1 0.0 0.0
           0.0 0.3 0.0 0.0 0.6 0.0
           0.2 0.0 0.2 0.1 0.0 0.6
           0.0 0.0 0.0 1.0 1.0 1.0
           0.0 0.3 0.3 0.7 0.8 1.0
           0.0 0.0 0.0 0.0 1.0 0.6]

    int = [0 0 0 0 0 0
           0 0 0 0 0 0
           0 0 0 0 0 0
           0 0 0 1 1 1
           0 0 0 0 0 1
           0 0 0 0 1 0]

    ruleset = Ruleset(Life(); boundary=Wrap())

    @testset "GtkOutput works" begin
        output = GtkOutput(int; tspan=1:1)
        data = DynamicGrids.SimData(output.extent, ruleset)
        DynamicGrids.showframe(output, data)
        Gtk.destroy(output.window)
        output = GtkOutput(flt; tspan=1:1)
        DynamicGrids.showframe(output, data)
        Gtk.destroy(output.window)
    end

end


text:
# Functions to calculate an ethnicity estimate.

using CSV
using DataFrames
using DelimitedFiles
using Printf

const DEBUG = true

const CHROMOSOMES = 22

# Length of the chromosome segments that are used for mapping.
const SEGMENT_LENGTH = 1000000

# Matching segments that are larger than MAX_MATCH_LENGTH are
# excluded from the evaluation because they indicate close family
# members. We do not want them to obfuscate the results..
const MAX_MATCH_LENGTH = 20000000

# Sizes of the different chromosomes. We use larger numbers here
# because chromosome sizes can change with different genome reference versions.
const CHROMOSOME_SIZES = [250000000, 250000000, 200000000, 200000000, 190000000,
                   180000000, 160000000, 150000000, 150000000, 140000000,
                   140000000, 140000000, 120000000, 110000000, 110000000,
                   100000000,  90000000,  90000000,  60000000,  70000000,
                    50000000,  60000000]

"""
    ethnicity(matches, segments, [parent_matches,] birth_country="", excludes=["USA", "Canada", "Australia"])

Calculate an ethnicity estimate from a MyHeritage "DNA matches
list" and the corresponding "shared DNA segments" list.
If a birth country is provided the calculation will be more accurate.
You can specify a number of countries that are excluded from the calculation
because of massive recent migrations.

Important: The CSV files must be in MyHeritage's English file format!

# Examples
```
julia> ethnicity("matches.csv", "shared_segments.csv")
Total: 1311 segments
        Germany  465  35.5%
        Denmark  336  25.6%
    Netherlands  206  15.7%

julia> ethnicity("matches.csv", "shared_segments.csv", birth_country="Germany")
Total: 1124 segments
        Germany  710  63.2%
        Denmark  164  14.6%
    Netherlands  106   9.4%
```
"""
function ethnicity end

function ethnicity(matches_file::AbstractString,
                   segments_file::AbstractString;
                   birth_country = "",
                   excludes=["USA", "Canada", "Australia"])
    local matches
    local segments

    # Read files.
    try
        matches = open(matches_file) do f
            readdlm(f, ',', String)
        end
        segments = open(segments_file) do f
            readdlm(f, ',', String)
        end
    catch ex
        error("Could not read file. ", ex)
    end

    # Convert matches and segments to DataFrames.
    try
        matches = matches_to_dataframe(matches)
        segments = segments_to_dataframe(segments)
    catch ex
        error("Could not interpret files. Wrong format? ", ex)
    end

    # Map matches' segments onto DNA.
    #matches_segments = innerjoin(matches, segments, on = :Name)
    matches_segments = innerjoin(matches, segments, on = :Match_ID)

    dna = DNA()
    map_countries!(dna, matches_segments, birth_country, excludes)

    # Evaluate ethnicities.
    eths, segments_total = ethnicities(dna)
    print_ethnicities(eths, segments_total)
end

function ethnicity(matches_file::AbstractString,
                   segments_file::AbstractString,
                   parent_matches_file::AbstractString;
                   birth_country = "",
                   excludes=["USA", "Canada", "Australia"])
    local matches
    local segments
    local parent

    # Read files.
    try
        matches = open(matches_file) do f
            readdlm(f, ',', String)
        end
        segments = open(segments_file) do f
            readdlm(f, ',', String)
        end
        parent = open(parent_matches_file) do f
            readdlm(f, ',', String)
        end
    catch ex
        error("Could not read file. ", ex)
    end

    # Convert matches and segments to DataFrames.
    try
        matches = matches_to_dataframe(matches)
        segments = segments_to_dataframe(segments)
        parent = matches_to_dataframe(parent)
    catch ex
        error("Could not interpret files. Wrong format? ", ex)
    end

    # Phase matches.
    # Matches that match child and parent.
    child_parent = innerjoin(matches, parent, on = :Match_ID, makeunique=true)
    # Matches that only match the second parent.
    child_parent2 = antijoin(matches, child_parent, on = :Match_ID, makeunique=true)
    # Remove matches that match both parents from the child-parent matches.
    both, child_parent = remove_false_matches(child_parent, parent)

    # Add segments to the matches lists.
    child_parent_segments = innerjoin(child_parent, segments, on = :Match_ID)
    child_parent2_segments = innerjoin(child_parent2, segments, on = :Match_ID)
    both_segments = innerjoin(both, segments, on = :Match_ID)

    # Map matches' segments onto DNA.
    parent_dna = DNA()
    parent2_dna = DNA()
    both_dna = DNA()
    map_countries!(parent_dna, child_parent_segments, birth_country, excludes)
    map_countries!(parent2_dna, child_parent2_segments, birth_country, excludes)
    map_countries!(both_dna, both_segments, birth_country, excludes)

    # Evaluate ethnicities.
    parent_eths, segments_total_1 = ethnicities(parent_dna)
    parent2_eths, segments_total_2 = ethnicities(parent2_dna)
    both_eths, segments_total_3 = ethnicities(both_dna)

    println("\nEthnicities you inherited from your parent.")
    print_ethnicities(parent_eths, segments_total_1)

    println("\nEthnicities you inherited from your other parent.")
    print_ethnicities(parent2_eths, segments_total_2)

    println("\nEthnicities that contributed to you and both of your parents.")
    print_ethnicities(both_eths, segments_total_3)
end

"""
DNA is modeled as an array of chromosomes.
Each chromosome consists of a number of segments represented
by a dictionary that maps country names to the number of matches
from each country.
"""
mutable struct DNA
    chromosomes::Array{Array{Dict{String, Integer}}}
end

function DNA()
    chromosomes = Array{Array{Dict{String, Integer}}}(undef, CHROMOSOMES)
    for i in 1:CHROMOSOMES
        nSegs = div(CHROMOSOME_SIZES[i] + SEGMENT_LENGTH - 1, SEGMENT_LENGTH)
        segments = Array{Dict{String, Integer}}(undef, nSegs)
        for j in 1:nSegs
            segments[j] = Dict{String, Integer}()
        end
        chromosomes[i] = segments
    end
    DNA(chromosomes)
end

"""
    map_countries!(dna::DNA, df::DataFrame)

Map country information from matches onto the DNA chromosomes.
The data frame "segments" must contain the following columns:
:Country, :Chromosome, :StartLocation, :EndLocation.
excludes contains a list of countries that should be excluded from mapping.
"""
function map_countries!(dna::DNA, segments::DataFrame, birth_country::String, excludes::Array{String})
    # Map matches' segments onto the DNA.
    nrows, _ = size(segments)
    for i in 1: nrows
        # Ignore close relatives.
        if segments[i, :EndLocation] - segments[i, :StartLocation] > MAX_MATCH_LENGTH
            continue
        end

        # Map matching segment onto the DNA.
        a = segments[i, :StartLocation] + SEGMENT_LENGTH
        b = segments[i, :EndLocation]
        c = segments[i, :Chromosome]
        country = segments[i, :Country]
        # Ignore countries that should be excluded.
        if in(country, excludes)
            continue
        end
        for j in a:SEGMENT_LENGTH:b
            seg = div(j, SEGMENT_LENGTH)
            if haskey(dna.chromosomes[c][seg], country)
                dna.chromosomes[c][seg][country] += 1
            else
                dna.chromosomes[c][seg][country] = 1
            end
        end
    end
    # Map birth country onto DNA segments that are already populated.
    if birth_country == "" || in(birth_country, excludes)
        return
    end
    for i in eachindex(dna.chromosomes)
        for j in eachindex(dna.chromosomes[i])
            # Check if DNA segment is populated.
            if length(keys(dna.chromosomes[i][j])) == 0
                continue
            end
            # Add birth country.
            if haskey(dna.chromosomes[i][j], birth_country)
                dna.chromosomes[i][j][birth_country] += 1
            else
                dna.chromosomes[i][j][birth_country] = 1
            end
        end
    end
end

"""
    ethnicities(dna::DNA)

Calculate ethnicities for a given DNA sample.

Returns an Array{Array{String}} and the total number of
segments that were evaluated. For each chromosome an array
exists that consists of segment entries. A segment entry is
a String that holds the most likely country origin of the
segment.
"""
function ethnicities(dna::DNA)
    # Create an empty country entry for each DNA segment.
    ethSegments = Array{Array{String}}(undef, CHROMOSOMES)
    for i in 1:CHROMOSOMES
        nSegs = length(dna.chromosomes[i])
        chromosome = Array{String}(undef, nSegs)
        for j in 1:nSegs
            chromosome[j] = ""
        end
        ethSegments[i] = chromosome
    end

    # Determine ethnicity for each segment.
    for i in 1:CHROMOSOMES
        for j in eachindex(dna.chromosomes[i])
            ethSegments[i][j] = most_common_ethnicity(dna.chromosomes[i][j])
        end
    end

    # Count segment ethnicities.
    ethnicities = Dict{String, Integer}()
    # Total number of evaluated segments.
    total = 0
    for i in 1:CHROMOSOMES
        for j in eachindex(ethSegments[i])
            if ethSegments[i][j] != ""
                total += 1
                if haskey(ethnicities, ethSegments[i][j])
                    ethnicities[ethSegments[i][j]] += 1
                else
                    ethnicities[ethSegments[i][j]] = 1
                end
            end
        end
    end
    return ethnicities, total
end

"""
    print_ethnicities(eth::Dict{String, Integer}, segments_total::Integer)

Print the result of the ethnicities function.
"""
function print_ethnicities(eth::Dict{String, Integer}, segments_total::Integer)
    # Sort results.
    results = []
    for (country, value) in eth
        percentage = value / segments_total * 100
        push!(results, (country, value, percentage))
    end
    sort!(results, by = x -> x[2], rev = true)

    # Output results.
    println("Total: ", segments_total, " segments")
    for i in results
        @printf("%15s %4d %5.1f%%\n", i[1], i[2], i[3])
    end
end

"""
    matches_to_dataframe(matches::Array{String, 2})

Converts matches to a valid DataFrame. MyHeritage CSV files
sometimes contain an invalid CSV format and can not be read by
the standard Julia CSV package. So we must do this manually.
"""
function matches_to_dataframe(matches::Array{String, 2})
    df = DataFrame(Match_ID = String[], Match_Name = String[], Country = String[], Total_cM_shared = Float64[])

    # Copy matches rows into DataFrame.
    match_names = Dict{String, Int64}()
    rows, _ = size(matches)
    for i = 2:rows
        match_id = matches[i, 1]
        country = matches[i, 4]
        match_name = matches[i, 2]
        matches_cM = matches[i, 10]
        if country == "" || match_name == "" || matches_cM == ""
            continue
        end

        # Convert cM to Float.
        cM = parse(Float64, replace(matches_cM, "," => ""))
        push!(df, [match_id, match_name, country, cM])
    end
    return df
end

"""
    segments_to_dataframe(segments::Array{String, 2})

Converts matches to a valid DataFrame. MyHeritage CSV files
sometimes contain an invalid CSV format and no valid key column.
They can not always be read by the standard Julia CSV package.
So we must do this manually.
"""
function segments_to_dataframe(segments::Array{String, 2})
    df = DataFrame(Match_ID = String[], Chromosome = Integer[], StartLocation = Integer[], EndLocation = Integer[])

    # Copy segments rows into DataFrame.
    rows, _ = size(segments)
    for i = 2:rows
        match_id = segments[i, 1]
        chromosome = segments[i, 4]
        startloc = segments[i, 5]
        endloc = segments[i, 6]
        if  chromosome == "" || startloc == "" || endloc == ""
            continue
        end

        chromosome = parse(Int64, chromosome)
        startloc = parse(Int64, startloc)
        endloc = parse(Int64, endloc)
        push!(df, [match_id, chromosome, startloc, endloc])
    end
    return df
end

"""
    remove_false_matches(child_matches, parent_matches)

Remove matches from the child_matches if corresponding parent_matches
share less cM with the match than the child.

Return the matches where both parents have contributed and
the matches that are from one parent only.
"""
function remove_false_matches(child_matches::DataFrame, parent_matches::DataFrame)
    both = similar(child_matches, 0)
    parent_cleaned = similar(child_matches, 0)
    # Create a dictionary of parent matches cM.
    parentCM = Dict{String, Float64}()
    for row in eachrow(parent_matches)
        name = row[:Name]
        cM = row[:Total_cM_shared]
        parentCM[name] = cM
    end
    # From the matches that were inherited from the parent
    # remove thosw where child shared cM > parent shared cM.
    for row in eachrow(child_matches)
        name = row[:Name]
        childCM = row[:Total_cM_shared]
        if childCM > parentCM[name]
            push!(both, row)
        else
            push!(parent_cleaned, row)
        end
    end
    return both, parent_cleaned
end

"""
    most_common_ethnicity(d::Dict{String, Integer})::String

Eetermine the ethnicity based on a map that contains a count
of matches from each country.

d is a dictionary that contains the number of matches for each country.

Returns the name of the country that harbors the most matches
or an empty string if there is no clear result.
"""
function most_common_ethnicity(d::Dict{String, Integer})::String
    if length(keys(d)) == 0
        return ""
    end

    # Determine country with the most cousins.
    country = ""
    top = 0
    second = 0
    for (key, value) in d
        DEBUG && print(key, " ", value, "  ")
        if value > second && value <= top
            second = value
        end
        if value > top
            country = key
            second = top
            top = value
        end
    end
    DEBUG && println()

    # Look if the most common country is valid.
    # Experiments with different conditions are allowed here.

    # All matches from one country.
    if second == 0
        return country
    end

    # Condition for the top most matches to be valid.
    if top / second >= 1.5
        return country
    else
        return ""
    end
end


text:
using RobustAndOptimalControl, ControlSystems

G = ssrand(1,1,3, proper=true)
Gd = add_low_frequency_disturbance(G)
@test Gd.nx == 4
@test rank(obsv(Gd)) == 4
@test rank(ctrb(Gd)) == 3
@test any(isapprox(0, atol=eps()), pole(Gd))

Gd = add_low_frequency_disturbance(G, measurement=true)
Gd.C[end] == 1
@test Gd.nx == 4
@test rank(obsv(Gd)) == 4
@test rank(ctrb(Gd)) == 3
@test any(isapprox(0, atol=eps()), pole(Gd))

G = ssrand(1,1,3, proper=true)
Gd = add_resonant_disturbance(G, 1, 0, 3)
@test Gd.nx == 5
@test rank(obsv(Gd)) == 5
@test rank(ctrb(Gd)) == 3
@test any(isapprox(1, atol=eps()), imag.(pole(Gd)))

Gd = add_resonant_disturbance(G, 1, 0, 1, measurement=true)
@test Gd.nx == 5
@test rank(obsv(Gd)) == 5
@test rank(ctrb(Gd)) == 3
@test any(isapprox(1, atol=eps()), imag.(pole(Gd)))

##

G = ssrand(1,1,3, proper=true, Ts=1)
GD = ssrand(1,1,3, proper=false, Ts=1)

## Diff
Gd = add_output_differentiator(G)
Gd2 = [tf(1,1); tf([1, -1], [1], 1)]*tf(G)
@test tf(Gd) ≈ Gd2
# @test hinfnorm(Gd-Gd2)[1] < 1e-10 hinfnorm not robust

## Int
Gd = add_output_integrator(G)
Gd2 = [tf(1,1); tf(1, [1, -1], 1)]*G
@test Gd ≈ Gd2
@test sminreal(Gd[1,1]) == G # Exact equivalence should hold here
@test Gd.nx == 4 # To guard agains changes in realization of tf as ss

Gd = add_input_integrator(G)
@test sminreal(Gd[1,1]) == G # Exact equivalence should hold here
@test Gd.nx == 4 # To guard agains changes in realization of tf as ss
@test tf(sminreal(Gd[2,1])) == tf(1, [1,-1], 1)

G = ssrand(2,1,3, proper=true, Ts=1)
Gd = add_input_integrator(G)
@test sminreal(Gd[1:2,1]) == G # Exact equivalence should hold here
@test Gd.nx == 4 # To guard agains changes in realization of tf as ss
@test tf(sminreal(Gd[3,1])) == tf(1, [1,-1], 1)

G = ssrand(1,2,3, proper=true, Ts=1)
Gd = add_input_integrator(G)
@test sminreal(Gd[1,1:2]) == G # Exact equivalence should hold here
@test Gd.nx == 4 # To guard agains changes in realization of tf as ss
@test tf(sminreal(Gd[2,1])) == tf(1, [1,-1], 1)
@test tf(sminreal(Gd[2,2])) == tf(0, 1)

Gd = add_input_integrator(G, 2) # other input
@test sminreal(Gd[1,1:2]) == G # Exact equivalence should hold here
@test Gd.nx == 4 # To guard agains changes in realization of tf as ss
@test tf(sminreal(Gd[2,2])) == tf(1, [1,-1], 1)
@test tf(sminreal(Gd[2,1])) == tf(0, 1)

# Both
Gd = add_output_differentiator(G)
Gd = add_output_integrator(Gd, 1)
Gd2 = [tf(1,1); tf([1, -1], [1], 1); tf(1, [1, -1], 1)]*tf(G)
w = exp10.(LinRange(-2, 2, 100))

# These are harder to compare, the freqresp method seems most reliable
f1 = freqresp(Gd, w)
f2 = freqresp(Gd2, w)
@test f1 ≈ f2
# bodeplot([Gd, Gd2])

# Input diff
nx,nu,ny = G.nx, G.nu, G.ny
Gd = add_input_differentiator(G, 1:G.nu)
@test Gd.A[nx+1:end, nx+1:end] == 0I
@test Gd.B[nx+1:end, :] == I
@test Gd.C[ny+1:end, nx+1:end] == -I
@test Gd.D[ny+1:end, :] == I


text:
using Test, BFloat16s, Printf

@testset "comparisons" begin
	@test BFloat16(1)   <  BFloat16(2)
	@test BFloat16(1f0) <  BFloat16(2f0)
	@test BFloat16(1.0) <  BFloat16(2.0)
	@test BFloat16(1)   <= BFloat16(2)
	@test BFloat16(1f0) <= BFloat16(2f0)
	@test BFloat16(1.0) <= BFloat16(2.0)
	@test BFloat16(2)   >  BFloat16(1)
	@test BFloat16(2f0) >  BFloat16(1f0)
	@test BFloat16(2.0) >  BFloat16(1.0)
	@test BFloat16(2)   >= BFloat16(1)
	@test BFloat16(2f0) >= BFloat16(1f0)
	@test BFloat16(2.0) >= BFloat16(1.0)
	@test BFloat16(2)   != BFloat16(1)
	@test BFloat16(2f0) != BFloat16(1f0)
	@test BFloat16(2.0) != BFloat16(1.0)
  @test iszero(BFloat16(0)) == true
  @test iszero(BFloat16(3.45)) == false
end

@testset "conversions" begin
	@test Float32(BFloat16(10)) == 1f1
	@test Float64(BFloat16(10)) == 10.0
	@test Int32(BFloat16(10)) == Int32(10)
	@test Int64(BFloat16(10)) == Int64(10)
end

@testset "functions" begin
	@test abs(BFloat16(-10)) == BFloat16(10)
	@test BFloat16(2) ^ BFloat16(4) == BFloat16(16)
	@test eps(BFloat16) == BFloat16(0.0078125)
	@test sqrt(BFloat16(4f0)) == BFloat16(2f0)
	@test round(BFloat16(10.4), RoundUp) == BFloat16(11.0)
	@test round(BFloat16(10.6), RoundDown) == BFloat16(10.0)
	@test round(BFloat16(3.2), RoundNearest) == BFloat16(3.0)
	@test round(BFloat16(4.8), RoundNearest) == BFloat16(5.0)
end

@testset "arithmetic" begin
  @test BFloat16(0.2) * BFloat16(5.0) == BFloat16(1.0)
  @test BFloat16(1.0) / BFloat16(5.0) == BFloat16(0.2)
  @test inv(BFloat16(5.0)) == BFloat16(0.2)
  @test zero(BFloat16) == BFloat16(0.0f0)
  @test one(BFloat16) == BFloat16(1.0)
end

@testset "printf" begin
  for (fmt, val) in (("%7.2f", "   1.23"),
                     ("%-7.2f", "1.23   "),
                     ("%07.2f", "0001.23"),
                     ("%.0f", "1"),
                     ("%#.0f", "1."),
                     ("%.3e", "1.234e+00"),
                     ("%.3E", "1.234E+00"),
                     ("%.2a", "0x1.3cp+0"),
                     ("%.2A", "0X1.3CP+0")),
      num in (BFloat16(1.234),)
      @test @eval(@sprintf($fmt, $num) == $val)
  end
  @test (@sprintf "%f" BFloat16(Inf)) == "Inf"
  @test (@sprintf "%f" BFloat16(NaN)) == "NaN"
  @test (@sprintf "%.0e" BFloat16(3e2)) == "3e+02"
  @test (@sprintf "%#.0e" BFloat16(3e2)) == "3.e+02"

  for (fmt, val) in (("%10.5g", "     123.5"),
                     ("%+10.5g", "    +123.5"),
                     ("% 10.5g","     123.5"),
                     ("%#10.5g", "    123.50"),
                     ("%-10.5g", "123.5     "),
                     ("%-+10.5g", "+123.5    "),
                     ("%010.5g", "00000123.5")),
      num in (BFloat16(123.5),)
      @test @eval(@sprintf($fmt, $num) == $val)
  end
  @test( @sprintf( "%10.5g", BFloat16(-123.5) ) == "    -123.5")
  @test( @sprintf( "%010.5g", BFloat16(-123.5) ) == "-0000123.5")
  @test (@sprintf "%a" BFloat16(1.5)) == "0x1.8p+0"
end


text:
using ThreadingUtilities
using Test

import InteractiveUtils


text:
mutable struct JobQueue
    func
    itr
    state::Union{Integer, Nothing}

    JobQueue(f, itr) = new(f, itr, nothing)
end

# Returns true if there was a return value.
function process(q::JobQueue)
    y = q.state === nothing ? iterate(q.itr) : iterate(q.itr, q.state)
    y === nothing && return false
    q.func(y[1])
    q.state = y[2]::Integer

    return true
end


text:
using CovarianceMatrices
using Test
using CSV
using LinearAlgebra
using Statistics
using Random
using GLM
using DataFrames
using JSON

const CM = CovarianceMatrices

datapath = joinpath(@__DIR__)

@testset "HAC - Basic checks..............................." begin
    X = [0 0; 3 7; 4 8; 5 9]
    X_demean = [-3  -6; 0   1; 1   2; 2   3]
    cache = CM.HACCache(X, TruncatedKernel(prewhiten=true));
    CM.demean!(cache, X, Val{true})
    @test all(cache.μ .== [3 6])
    @test all(cache.X_demean .== X_demean)
    CM.fit_var!(cache)
    @test all(cache.D   .≈ [-9/5 -4; 1 2])
    @test all(cache.XX  .≈ [3/5 1; 0 0; 9/5 3])
end

@testset "HAC - Optimal Bandwidth Calculations............." begin
    X = [1, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7,  .8,
         .3, .4, .5, .6, .7, .8, .9, .10, .11, .12]
    X = reshape(X, 8, 3)

    k = ParzenKernel(prewhiten=true)
    cache = CovarianceMatrices.HACCache(X, k)
    CovarianceMatrices.demean!(cache, X, Val{true})

    Xd = reshape([ -2.75,    -1.75,   -0.75,     0.25,
                   1.25,     2.25,    3.25,    -1.75,
                   -0.312500000000002,   0.6875,   1.6875,  2.6875,
                   3.6875,  -2.5125,  -3.0125, -2.9125,
                   0.0212499999999998,  0.12125,  0.22125, 0.32125,
                   0.42125, -0.37875, -0.36875, -0.35875], 8, 3)
    @test all(Xd .≈ cache.X_demean)

    CovarianceMatrices.prewhiten!(cache)

    D = reshape([-0.409214603090833,   -0.546926852567603, -0.0793279159360713,
                 7.17805137100394,    -1.09268085748226,  -0.152674132870786,
                 -56.8116543172651,    12.0360552297997,    1.52457939755027], 3,3)'

    @test all(D .≈ cache.D)

    XX = reshape([0.575048549180823,     0.487377212994233,   0.399705876807639,
                  0.312034540621058,      0.22436320443446,   0.688172853937608,
                  0.254529685703059,     -1.41377778615736,  0.0222244009125362,
                  1.45822658798243,       2.89422877505233,   -2.86976903787778,
                  0.0313806821395586,   0.0116065536680327,   -0.17700974754426,
                  0.00253436150757022,     0.1820784705594,    0.36162257961123,
                  -0.35883331133694,   0.00357849884047592, 0.00132355436565176], 7, 3)
    @test all(XX .≈ cache.XX)

    CovarianceMatrices.fit_ar!(cache)
    ρ = [-0.411560904449045,  -0.202815142161935, -0.202567957968907]
    σ⁴= [ 0.00043730337156934, 8.49515993964198,   0.0020740008652368]
    @test all(ρ .≈ cache.ρ)
    @test all(σ⁴ .≈ cache.σ⁴)

    CovarianceMatrices.fit_ar!(cache)
    if isempty(k.weights)
        for j in 1:size(X,2)
            push!(k.weights, 1.0)
        end
    end
    a1, a2 = CovarianceMatrices.getalpha!(cache, k.weights)
    @test a1 ≈ 0.1789771071933
    @test a2 ≈ 0.07861018601427
    @test CovarianceMatrices.bw_andrews(k, a1, a2, 7) ≈ 2.361704327253
    k.bw .= CovarianceMatrices.bw_andrews(k, a1, a2, 7)
    bw = k.bw[1]
    fill!(cache.V, zero(eltype(cache.V)))
    mul!(cache.V, cache.XX', cache.XX)
    V = reshape([0.707026904730492,   0.0322404268557179, 0.00367653990459955,
                 0.0322404268557179, 10.3694701634646,    1.29616765021965,
                 0.00367653990459955, 1.29616765021965,   0.162019118007504], 3, 3)*2

    @test all(V ≈ cache.V)
    triu!(cache.V)
    @test CovarianceMatrices.kernel(k, 1/bw) ≈  0.379763254768776
    @test CovarianceMatrices.kernel(k, 2/bw) ≈  0.00718479751373071

    @test floor(Int, bw) == 2

    for j in -floor(Int, bw):-1
        k_j = CovarianceMatrices.kernel(k, j/bw)
        LinearAlgebra.axpy!(k_j, CovarianceMatrices.Γ!(cache, j), cache.V)
    end

    for j in 1:floor(Int, bw)
        k_j = CovarianceMatrices.kernel(k, j/bw)
        LinearAlgebra.axpy!(k_j, CovarianceMatrices.Γ!(cache, j), cache.V)
    end

    LinearAlgebra.copytri!(cache.V, 'U')

    V = reshape([ 2.18378130064673,  -0.126247404742945, -0.0168728162431914,
                  -0.126247404742945, 17.4806466521012,    2.18514395521502,
                  -0.0168728162431914, 2.18514395521502,   0.273151430809999], 3, 3)

    @test all(V .≈ cache.V)

    fill!(cache.Q, zero(eltype(cache.Q)))
    for i = 1:size(cache.Q, 2)
        cache.Q[i,i] = one(eltype(cache.Q))
    end
    v = ldiv!(qr(I-cache.D'), cache.Q)
    cache.V .= v*cache.V*v'

    V = reshape([10.3056593241458,   7.62922768968777,  0.731955791380014,
                 7.62922768968773,  8.49755610886395,  0.859802385984295,
                 0.731955791380006, 0.859802385984297, 0.0874594553955033], 3, 3)

    @test all(V .≈ cache.V)
end

andrews_kernels = [:TruncatedKernel, :ParzenKernel, :TukeyHanningKernel, :QuadraticSpectralKernel, :BartlettKernel]
neweywest_kernels = [:ParzenKernel, :QuadraticSpectralKernel, :BartlettKernel]

@testset "HAC Univariate Long Run Variance................." begin
    univariate = JSON.parse(read("testdata/univariate.json", String))

    global x = [0.164199142438412, -0.22231320001538, -2.29596418288347,
         -0.562717239408665, 1.11832433510705, 0.259810317798567,
         0.647885100553029, 2.53438209392891, -0.419561292138475,
         1.19138801574376, 2.52661839907567, -0.443382492040113,
         -0.137169008509379, -0.967782699857501, 0.150507152028812,
         -1.27098181862663, 1.38639734998711, 0.231229342441316,
         -0.943827510026301, -1.11211457440442]

    function fopt!(u)
        global da = Dict{String, Any}()
        global dn = Dict{String, Any}()
        for pre in (:false, :true)
            da["bwtype"] = "auto"
            da["prewhite"] = pre == :true ? 1 : 0
            dn["bwtype"] = "auto"
            dn["prewhite"] = pre == :true ? 1 : 0
            for k in andrews_kernels
                eval(quote
                     tmp = CovarianceMatrices.covariance(x[1:end,:], ($k)(prewhiten=$pre))
                     da[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     if Symbol($k) in neweywest_kernels
                     tmp = CovarianceMatrices.covariance(x[1:end,:], ($k)(prewhiten=$pre))
                     dn[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     end
                     end)
            end
            push!(u, Dict("andrews" => da, "neweywest" => dn))
            da = Dict{String, Any}()
            dn = Dict{String, Any}()
        end
    end

    function ffix!(u)
        global da = Dict{String, Any}()
        global dn = Dict{String, Any}()
        for pre in (:false, :true)
            da["bwtype"] = "auto"
            da["prewhite"] = pre == :true ? 1 : 0
            dn["bwtype"] = "auto"
            dn["prewhite"] = pre == :true ? 1 : 0
            for k in andrews_kernels
                eval(quote
                     tmp = CovarianceMatrices.covariance(x[1:end,:], ($k)(1.5, prewhiten=$pre))
                     da[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     if Symbol($k) in neweywest_kernels
                         tmp = CovarianceMatrices.covariance(x[1:end,:], ($k)(1.5, prewhiten=$pre))
                         dn[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     end
                end)
            end
            push!(u, Dict("andrews" => da, "neweywest" => dn))
            da = Dict{String, Any}()
            dn = Dict{String, Any}()
        end
    end

    u = Any[]
    fopt!(u)
    ffix!(u)
    for j in 1:4, h in ("andrews",), k in ("Truncated", "Bartlett", "Tukey-Hanning", "Quadratic Spectral")
        @test univariate[j][h][k]["V"] ≈ u[j][h][k]["V"]
        @test univariate[j][h][k]["bw"] ≈ u[j][h][k]["bw"]
    end

end

@testset "HAC Multivariate Long Run Variance..............." begin
    multivariate = JSON.parse(read("testdata/multivariate.json", String))

    global X = reshape([-0.392343481403554, 0.369702177262769, -0.283239954622452, -1.71955159119934,
                 -0.196311837779395, 0.567935573633728, 0.675556050609011, -0.59266740997335,
                 0.433501752110958, -0.108134518170885, 0.686698254275096,
                 0.905166380758529, 0.997306068740448, -0.498167979332402, -0.52789448547848,
                 0.100420689892425, 1.49495246674266, -0.601579368927597, -0.166121087348934,
                 0.545348413754815, 0.294595939676804, -1.4327640501377, 0.719755536433108,
                 -1.2570984424727, -1.5357895599007, 1.0027040264714, 1.08593183743541,
                 -0.711342086638968, -0.772154611498373, 1.30668004120162, 2.89166460485415,
                 0.614941151434828, -1.59157215625112, -0.51709522643053, 1.97543651877893,
                 1.92940239570577, 0.889679654045449, -0.471409046820902, -1.30435924644088,
                 0.424385975482526, -1.63421078193559, -0.562676024339334, -1.9773880104471,
                 -1.13903619779294, 0.586555209142932, -1.60335681344097, -1.19189517293108,
                 2.13456541115109, -1.42078068631219, -0.207019328929601, -0.600736667890819,
                 -1.41872438508684, -0.608094864262906, 1.07318477908557, -0.477244503174433,
                 0.567645883530209, -0.149728929145769, 1.41918460373266, 0.462399751718563,
                 -0.132320093478005, 1.27839375300926, -0.480093346919616, -0.0428876036353262,
                 -1.56004471316687, -0.134394994431227, 2.04942053641657, -1.8022396128532,
                 -1.72537103051051, -0.657108341488561, 1.34392729540088, 1.90019159830845,
                 0.126395933995686, -0.770308826277006, 0.457784471260252, -0.160271362465108,
                 0.0598471594810112, 1.75048422563306, -0.737657566749125, -0.462941254536989,
                 0.699940101308202, 1.2689352813542, -0.336296482224711, 1.42781944149594,
                 -1.69159952993731, -1.15816200645816, 0.83309270822555, -1.34670872662577,
                 2.24540522326547, 1.14409536415596, -0.959417691691381], 30, 3)

    function fopt!(u)
        global da = Dict{String, Any}()
        global dn = Dict{String, Any}()
        for pre in (:false, :true)
            da["bwtype"] = "auto"
            da["prewhite"] = pre == :true ? 1 : 0
            dn["bwtype"] = "auto"
            dn["prewhite"] = pre == :true ? 1 : 0
            for k in andrews_kernels
                eval(quote
                     tmp = CovarianceMatrices.covariance(X, ($k)(prewhiten=$pre))
                     da[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     if Symbol($k) in neweywest_kernels
                     tmp = CovarianceMatrices.covariance(X, ($k)(prewhiten=$pre))
                     dn[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     end
                     end)
            end
            push!(u, Dict("andrews" => da, "neweywest" => dn))
            da = Dict{String, Any}()
            dn = Dict{String, Any}()
        end
    end

    function ffix!(u)
        global da = Dict{String, Any}()
        global dn = Dict{String, Any}()
        for pre in (:false, :true)
            da["bwtype"] = "auto"
            da["prewhite"] = pre == :true ? 1 : 0
            dn["bwtype"] = "auto"
            dn["prewhite"] = pre == :true ? 1 : 0
            for k in andrews_kernels
                eval(quote
                     tmp = CovarianceMatrices.covariance(X, ($k)(1.5, prewhiten=$pre))
                     da[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     if Symbol($k) in neweywest_kernels
                         tmp = CovarianceMatrices.covariance(X, ($k)(1.5, prewhiten=$pre))
                         dn[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     end
                     end)
            end
            push!(u, Dict("andrews" => da, "neweywest" => dn))
            da = Dict{String, Any}()
            dn = Dict{String, Any}()
        end
    end

    u = Any[]
    fopt!(u)
    ffix!(u)
    for j in 1:4, h in ("andrews",), k in ("Truncated", "Bartlett", "Tukey-Hanning", "Quadratic Spectral")
        @test hcat(multivariate[j][h][k]["V"]...) ≈ u[j][h][k]["V"]
        @test multivariate[j][h][k]["bw"] ≈ u[j][h][k]["bw"]
    end

end

@testset "HAC OLS VCOV....................................." begin
    reg = JSON.parse(read("testdata/regression.json", String))
    global df = CSV.read("testdata/ols_df.csv")

    function fopt!(u)
        global da = Dict{String, Any}()
        global dn = Dict{String, Any}()
        for pre in (:false, :true)
            da["bwtype"] = "auto"
            da["prewhite"] = pre == :true ? 1 : 0
            dn["bwtype"] = "auto"
            dn["prewhite"] = pre == :true ? 1 : 0
            for k in andrews_kernels
                eval(quote
                     ols = glm(@formula(y~x1+x2+x3), df, Normal(), IdentityLink())
                     tmp = CovarianceMatrices.vcov(ols, ($k)(prewhiten=$pre), dof_adjustment = false)
                     da[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     if Symbol($k) in neweywest_kernels
                     tmp = CovarianceMatrices.vcov(ols, ($k)(prewhiten=$pre), dof_adjustment = false)
                     dn[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     end
                     end)
            end
            push!(u, Dict("andrews" => da, "neweywest" => dn))
            da = Dict{String, Any}()
            dn = Dict{String, Any}()
        end
    end

    function ffix!(u)
        global da = Dict{String, Any}()
        global dn = Dict{String, Any}()
        for pre in (:false, :true)
            da["bwtype"] = "auto"
            da["prewhite"] = pre == :true ? 1 : 0
            dn["bwtype"] = "auto"
            dn["prewhite"] = pre == :true ? 1 : 0
            for k in andrews_kernels
                eval(quote
                     ols = glm(@formula(y~x1+x2+x3), df, Normal(), IdentityLink())
                     tmp = CovarianceMatrices.vcov(ols, ($k)(1.5, prewhiten=$pre), dof_adjustment = false)
                     da[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     if Symbol($k) in neweywest_kernels
                         tmp = CovarianceMatrices.vcov(ols, ($k)(1.5, prewhiten=$pre), dof_adjustment = false)
                         dn[String($k)] = Dict{String, Any}("bw" => tmp.K.bw, "V" => tmp.V)
                     end
                     end)
            end
            push!(u, Dict("andrews" => da, "neweywest" => dn))
            da = Dict{String, Any}()
            dn = Dict{String, Any}()
        end
    end

    u = Any[]
    fopt!(u)
    ffix!(u)
    for j in 1:4, h in ("andrews",), k in ("Truncated", "Bartlett", "Tukey-Hanning", "Quadratic Spectral")
        @test hcat(reg[j][h][k]["V"]...) ≈ u[j][h][k]["V"]
        @test reg[j][h][k]["bw"] ≈ u[j][h][k]["bw"]
    end

end

@testset "Accessor functions .............................." begin
    Random.seed!(1)
    y = randn(100)
    X = [ones(100) randn(100,4)]
    df = DataFrame(y = y, x1 = X[:,2], x2 = X[:,3], x3 = X[:,4], x4 = X[:,5])

    OLS1 = fit(GeneralizedLinearModel, X, y, Normal(), IdentityLink())
    OLS2 = fit(LinearModel, X, y)
    OLS3 = glm(@formula(y~x1+x2+x3+x4), df, Normal(), IdentityLink())
    OLS4 = lm(X,y)

    @test CovarianceMatrices.modelmatrix(OLS1) == X
    @test CovarianceMatrices.modelmatrix(OLS2) == X
    @test CovarianceMatrices.modelmatrix(OLS3) == X
    @test CovarianceMatrices.modelmatrix(OLS4) == X

    @test CovarianceMatrices.unweighted_nobs(OLS1) == 100
    @test CovarianceMatrices.unweighted_nobs(OLS1) == CovarianceMatrices.unweighted_nobs(OLS2)
    @test CovarianceMatrices.unweighted_nobs(OLS2) == CovarianceMatrices.unweighted_nobs(OLS3)
    @test CovarianceMatrices.unweighted_nobs(OLS3) == CovarianceMatrices.unweighted_nobs(OLS4)

    @test CovarianceMatrices.pseudohessian(OLS1) ≈ inv(X'X)
    @test CovarianceMatrices.pseudohessian(OLS2) ≈ inv(X'X)
    @test CovarianceMatrices.pseudohessian(OLS3) ≈ inv(X'X)
    @test CovarianceMatrices.pseudohessian(OLS4) ≈ inv(X'X)

    res = y - X*coef(OLS1)
    @test CovarianceMatrices.residuals(OLS1) ≈ res
    @test CovarianceMatrices.residuals(OLS2) ≈ res
    @test CovarianceMatrices.residuals(OLS3) ≈ res
    @test CovarianceMatrices.residuals(OLS4) ≈ res

end

@testset "HAC - GLM........................................" begin
    clotting = DataFrame(
        u    = log.([5,10,15,20,30,40,60,80,100]),
        lot1 = [118,58,42,35,27,25,21,19,18],
        lot2 = [69,35,26,21,18,16,13,12,12],
        w    = 9.0*[1/8, 1/9, 1/25, 1/6, 1/14, 1/25, 1/15, 1/13, 0.3022039]
    )

    GAMMA = glm(@formula(lot1~u), clotting, Gamma(),InverseLink(), wts = convert(Array, clotting[:w]))
    V = vcov(GAMMA, ParzenKernel())
    Vp = [5.48898e-7 -2.60409e-7; -2.60409e-7 1.4226e-7]
    @test V ≈ Vp atol = 1e-08

    GAMMA = glm(@formula(lot1~u), clotting, Gamma(),InverseLink())
    V = vcov(GAMMA, ParzenKernel())
    Vp = [5.81672e-7 -2.24162e-7; -2.24162e-7 1.09657e-7]
    @test V ≈ Vp atol = 1e-08
end

@testset "HC..............................................." begin

    # A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
    clotting = DataFrame(
        u    = log.([5,10,15,20,30,40,60,80,100]),
        lot1 = [118,58,42,35,27,25,21,19,18],
        lot2 = [69,35,26,21,18,16,13,12,12],
        w    = 9.0*[1/8, 1/9, 1/25, 1/6, 1/14, 1/25, 1/15, 1/13, 0.3022039]
    )

    ## Unweighted OLS though GLM interface
    OLS = fit(GeneralizedLinearModel, @formula(lot1~u),clotting, Normal(), IdentityLink())
    mf = ModelFrame(@formula(lot1~u),clotting)
    X = ModelMatrix(mf).m
    y = clotting[:lot1]
    GL  = fit(GeneralizedLinearModel, X,y, Normal(), IdentityLink())
    LM  = lm(X,y)

    S0 = vcov(OLS, HC0())
    S1 = vcov(OLS, HC1())
    S2 = vcov(OLS, HC2())
    S3 = vcov(OLS, HC3())
    S4 = vcov(OLS, HC4())
    S4m = vcov(OLS, HC4m())
    S5 = vcov(OLS, HC5())

    St0 = [720.621306411 -190.064512543; -190.064512543 51.163333742]
    St1 = [926.513108242 -244.368658984; -244.368658984 65.781429097]
    St2 = [1300.895673284 -343.330672699; -343.330672699  91.997186033]
    St3 = [2384.50393347 -628.97499232; -628.97499232 167.78976878]
    St4 = [2538.74635384 -667.95972319; -667.95972319 177.26308957]
    St4m = [3221.09520169 -849.64802585; -849.64802585 226.17046981]
    St5 = [1334.670541439 -351.751377823; -351.751377823  93.949230276]

    @test S0 ≈ St0
    @test S1 ≈ St1
    @test S2 ≈ St2
    @test S3 ≈ St3
    @test S4 ≈ St4
    @test S4m ≈ St4m
    @test S5 ≈ St5

    S0 = vcov(GL, HC0())
    S1 = vcov(GL, HC1())
    S2 = vcov(GL, HC2())
    S3 = vcov(GL, HC3())
    S4 = vcov(GL, HC4())
    S4m = vcov(GL, HC4m())
    S5 = vcov(GL, HC5())

    @test S0 ≈ St0
    @test S1 ≈ St1
    @test S2 ≈ St2
    @test S3 ≈ St3
    @test S4 ≈ St4
    @test S4m ≈ St4m
    @test S5 ≈ St5

    ## Weighted OLS though GLM interface
    wOLS = fit(GeneralizedLinearModel, @formula(lot1~u), clotting, Normal(),
               IdentityLink(), wts = Vector{Float64}(clotting[:w]))

    wts = Vector{Float64}(clotting[:w])
    X = [fill(1,size(clotting[:u])) clotting[:u]]
    y = clotting[:lot1]
    wLM = lm(X, y)
    wGL = fit(GeneralizedLinearModel, X, y, Normal(),
              IdentityLink(), wts = wts)

    S0 = vcov(wOLS, HC0())
    S1 = vcov(wOLS, HC1())
    S2 = vcov(wOLS, HC2())
    S3 = vcov(wOLS, HC3())
    S4 = vcov(wOLS, HC4())
    S4m= vcov(wOLS, HC4m())
    S5 = vcov(wOLS, HC5())

    St0 = [717.736178076 -178.404274981; -178.404274981   45.822730697]
    St1 = [922.803657527 -229.376924975; -229.376924975 58.914939468]
    St2 = [1412.940497584 -361.329969345; -361.329969345 95.912520696]
    St3 = [2869.53068690 -756.29761027; -756.29761027 208.23437869]
    St4 = [3969.9130263 -1131.3577578; -1131.3577578 342.2858663]
    St4m= [4111.62611908 -1103.17362711; -1103.17362711   310.19430896]
    St5 = [1597.40932634 -420.66907485; -420.66907485 115.99180777]

    @test S0  ≈ St0
    @test S1  ≈ St1
    @test S2  ≈ St2
    @test S3  ≈ St3
    @test S4  ≈ St4
    @test S4m ≈ St4m
    @test S5  ≈ St5

    ## Unweighted GLM - Gamma
    GAMMA = glm(@formula(lot1~u), clotting, Gamma(),InverseLink())

    S0 = vcov(GAMMA, HC0())
    S1 = vcov(GAMMA, HC1())
    S2 = vcov(GAMMA, HC2())
    S3 = vcov(GAMMA, HC3())
    S4 = vcov(GAMMA, HC4())
    S4m = vcov(GAMMA, HC4m())
    S5 = vcov(GAMMA, HC5())

    St0 = [4.504287921232951e-07 -1.700020601541489e-07;
           -1.700020601541490e-07  8.203697048568913e-08]

    St1 = [5.791227327299548e-07 -2.185740773410504e-07;
       -2.185740773410510e-07  1.054761049101728e-07]

    St2 = [3.192633083111232e-06 -9.942484630848573e-07;
           -9.942484630848578e-07  3.329973305723091e-07]

    St3 = [2.982697811926944e-05 -8.948137019946751e-06;
           -8.948137019946738e-06  2.712024459305714e-06]

    St4 = [0.002840158946368653 -0.0008474436578800430;
           -0.000847443657880045  0.0002528819761961959]

    St4m= [9.2891282926e-05  -2.7759505159e-05;
           -2.7759505159e-05  8.3203461732e-06]

    St5 = [2.9781374021e-05  -8.9232514073e-06
           -8.9232514073e-06  2.6952175350e-06]

    @test S0  ≈ St0 atol = 1e-06
    @test S1  ≈ St1 atol = 1e-06
    @test S2  ≈ St2 atol = 1e-06
    @test S3  ≈ St3 atol = 1e-06
    @test S4  ≈ St4 atol = 1e-05
    @test S4m ≈ St4m atol = 1e-06
    @test S5  ≈ St5 atol = 1e-06

    ## Weighted Gamma

    GAMMA = glm(@formula(lot1~u), clotting, Gamma(),InverseLink(), wts = convert(Array, clotting[:w]))

    S0 = vcov(GAMMA, HC0())
    S1 = vcov(GAMMA, HC1())
    S2 = vcov(GAMMA, HC2())
    S3 = vcov(GAMMA, HC3())
    S4 = vcov(GAMMA, HC4())
    S4m = vcov(GAMMA, HC4m())
    S5 = vcov(GAMMA, HC5())

    St0 = [4.015104e-07 -1.615094e-07;
           -1.615094e-07  8.378363e-08]

    St1 = [5.162277e-07 -2.076549e-07;
           -2.076549e-07  1.077218e-07]

    St2 = [2.720127e-06 -8.490977e-07;
           -8.490977e-07  2.963563e-07]

    St3 = [2.638128e-05 -7.639883e-06;
           -7.639883e-06  2.259590e-06]

    St4 = [0.0029025754 -0.0008275858;
           -0.0008275858  0.0002360053]

    St4m = [8.493064e-05 -2.436180e-05;
            -2.436180e-05  7.042101e-06]

    St5 = [2.6206554518e-05 -7.5421496876e-06
          -7.5421496876e-06  2.2017813312e-06]

    @test S0  ≈ St0 atol = 1e-06
    @test S1  ≈ St1 atol = 1e-06
    @test S2  ≈ St2 atol = 1e-06
    @test S3  ≈ St3 atol = 1e-06
    @test S4  ≈ St4 atol = 1e-05
    @test S4m ≈ St4m atol = 1e-06
    @test S5  ≈ St5 atol = 1e-06
end

@testset "CRHC............................................." begin

    df = CSV.read("testdata/wols_test.csv")
    df_sorted = sort!(copy(df), :cl)

    St1 = [.0374668, .0497666, .0472636, .0437952, .0513613, .0435369]

    OLS = fit(GeneralizedLinearModel, @formula(Y~X1+X2+X3+X4+X5), df, Normal(), IdentityLink())
    cl = convert(Array, df[:cl])
    k0 = CRHC0(cl)
    k1 = CRHC1(cl)
    k2 = CRHC2(cl)
    k3 = CRHC3(cl)

    V0 = vcov(OLS, k0)
    V1 = vcov(OLS, k1)
    V2 = vcov(OLS, k2)
    V3 = vcov(OLS, k3)

    @test V1 ≈ [0.00140376 0.000215526 -5.99768e-5 0.000296271 0.000460622 -0.000139741;
                0.000215526 0.00247671 -0.000270429 0.000218622 0.000610127 7.23345e-5;
                -5.99768e-5 -0.000270429 0.00223385 -0.000145166 -0.00018859 -0.000903561;
                0.000296271 0.000218622 -0.000145166 0.00191802 -0.000444364 -0.000420563;
                0.000460622 0.000610127 -0.00018859 -0.000444364 0.00263798 0.000736363;
                -0.000139741 7.23345e-5 -0.000903561 -0.000420563 0.000736363 0.00189546] atol = 1e-08

    ## Note sandwich in R has HC3 without G/(G-1) and CRHC2 is problematic

    @test stderror(OLS, k0, sorted = false) == sqrt.(diag(V0))

    OLS_sorted = fit(GeneralizedLinearModel, @formula(Y~X1+X2+X3+X4+X5), df_sorted, Normal(), IdentityLink())
    V0s = vcov(OLS_sorted, k0, sorted = true)
    V1s = vcov(OLS_sorted, k1, sorted = true)
    V2s = vcov(OLS_sorted, k2, sorted = true)
    V3s = vcov(OLS_sorted, k3, sorted = true)

    @test V0s == V0
    @test V1s == V1
    @test V2s == V2
    @test V3s == V3

    wOLS = fit(GeneralizedLinearModel, @formula(Y~X1+X2+X3+X4+X5), df,
               Normal(), IdentityLink(), wts = convert(Array{Float64}, df[:w]))

    cl = convert(Array, df[:cl])
    k0 = CRHC0(cl)
    k1 = CRHC1(cl)
    k2 = CRHC2(cl)
    k3 = CRHC3(cl)

    V0 = vcov(wOLS, k0)
    V1 = vcov(wOLS, k1)
    V2 = vcov(wOLS, k2)
    V3 = vcov(wOLS, k3)

    @test V1 ≈ [0.00183525 0.000137208 -0.00038971 0.000389943 0.000619903 0.00019496;
                0.000137208 0.00242781 -0.000272316 0.000462353 2.99597e-5 0.000133303;
                -0.00038971 -0.000272316 0.00273479 -0.000113765 -7.26396e-5 -0.000998524;
                0.000389943 0.000462353 -0.000113765 0.00171538 -0.00067357 -0.000416268;
                0.000619903 2.99597e-5 -7.26396e-5 -0.00067357 0.00225446 0.00106796;
                0.00019496 0.000133303 -0.000998524 -0.000416268 0.00106796 0.00226444] atol = 1e-07

    wOLS_sorted = fit(GeneralizedLinearModel, @formula(Y~X1+X2+X3+X4+X5), df_sorted,
                      Normal(), IdentityLink(), wts = convert(Array{Float64}, df[:w]))

    V0s = vcov(wOLS_sorted, k0, sorted = true)
    V1s = vcov(wOLS_sorted, k1, sorted = true)
    V2s = vcov(wOLS_sorted, k2, sorted = true)
    V3s = vcov(wOLS_sorted, k3, sorted = true)

    @test V1 ≈ [0.00183525 0.000137208 -0.00038971 0.000389943 0.000619903 0.00019496;
                0.000137208 0.00242781 -0.000272316 0.000462353 2.99597e-5 0.000133303;
                -0.00038971 -0.000272316 0.00273479 -0.000113765 -7.26396e-5 -0.000998524;
                0.000389943 0.000462353 -0.000113765 0.00171538 -0.00067357 -0.000416268;
                0.000619903 2.99597e-5 -7.26396e-5 -0.00067357 0.00225446 0.00106796;
                0.00019496 0.000133303 -0.000998524 -0.000416268 0.00106796 0.00226444] atol = 1e-07

    innovation = CSV.read("testdata/InstInnovation.csv", allowmissing=:none)

    innovation[:capemp] = log.(innovation[:capital]./innovation[:employment])
    innovation[:lsales] = log.(innovation[:sales])
    innovation[:year] = categorical(innovation[:year])
    innovation[:industry] = categorical(innovation[:industry])
    #innovation[:company] = categorical(innovation[:company])
    pois = glm(@formula(cites ~ institutions + capemp + lsales + industry + year), innovation, Poisson(), LogLink())

    Vt = [0.817387, 5.7907e-6, 0.0184833, 0.00172419]
    @test diag(vcov(pois, CRHC0(innovation[:company]), Matrix))[1:4] ≈ Vt atol = 1e-5

end

@testset "CovarianceMatrices Methods......................." begin
    df = CSV.read("testdata/wols_test.csv")
    df_sorted = sort!(copy(df), :cl)
    cl = convert(Array, df[:cl])
    wOLS = fit(GeneralizedLinearModel, @formula(Y~X1+X2+X3+X4+X5), df,
               Normal(), IdentityLink(), wts = convert(Array{Float64}, df[:w]))

    V0 = vcov(wOLS, HC0())
    V1 = vcov(wOLS, ParzenKernel())
    V2 = vcov(wOLS, CRHC0(cl))

    V0s = vcov(wOLS, HC0(), CovarianceMatrix, SVD)
    V1s = vcov(wOLS, ParzenKernel(), CovarianceMatrix, SVD)
    V2s = vcov(wOLS, CRHC0(cl), CovarianceMatrix, SVD)

    V0m = vcov(wOLS, HC0(), Matrix)
    V1m = vcov(wOLS, ParzenKernel(), Matrix)
    V2m = vcov(wOLS, CRHC0(cl), Matrix)

    @test V0 == V0m
    @test V1 == V1m
    @test V2 == V2m

    @test V0.F == cholesky(Symmetric(V0m))
    @test V1.F == cholesky(Symmetric(V1m))
    @test V2.F == cholesky(Symmetric(V2m))

    @test V0s.F == svd(V0m)
    @test V1s.F == svd(V1m)
    @test V2s.F == svd(V2m)

end


text:
## Batch result figures: compare estimators
using EstimateSDE
using Serialization
using Pipe: @pipe
using Plots
pgfplotsx()
using DataFrames
using Statistics
using StatsPlots
using REPL.TerminalMenus

## Settings
name = "experiment01"
ssize = (300,220)

## Load settings and aggregated results
results_dir = settings_resultsdir(name)
S = deserialize(joinpath(results_dir,"settings.jld"))
R = deserialize(joinpath(results_dir,"aggregate.jld")).df

figure_dir = joinpath(results_dir,"figures")
if !isdir(figure_dir) mkpath(figure_dir) end

## Preprocessing for plotting
R.proclabel = string.(R.proclabel)
R.proclabel = map(x->split(String(x),"_")[end],R.proclabel)
em = Dict(SdeRootsEngine=>"root",ElimQ=>"coef")
R.param = [em[x] for x in R.engine]
mmap = Dict("MLE EXPAND"=>"expand","MLE init ALL_RANDOM"=>"random init")
R.cmethod = map(eachrow(R)) do x
  m = get(mmap,x.method,x.method)
  m*" "*x.param
end

## KLD of MLE by process and computation method
dfn = @pipe(
  filter(R) do x
    x.cmethod ∈ ["MLE root","MLE coef","random init root"] &&
    x.proclabel∈["A","B","C"] &&
    x.kld_label=="regular" && x.kld_T==0.2
  end |>
  select(_,[:proclabel,:cmethod,:kld,:kld_std]) |>
  groupby(_,[:proclabel]) |>
  combine(_) do df
    kld_norm = filter(x->x.cmethod=="MLE root",df).kld[1]
    dfc = copy(df)
    dfc[:,:kld] = df.kld./kld_norm
    dfc[:,:kld_std] = df.kld_std./kld_norm
    dfc
  end
)

## Bar plot of relative KLD per engine, per process label
plt = @df dfn groupedbar(:proclabel,:kld,group=:cmethod,bar_width=0.8,yerror=:kld_std,
  linecolor=:black) #,markerstrokecolor=:auto) #lw=0 removes line around box
plot!(legend=:topleft,fg_legend=nothing,bg_legend=nothing, size=ssize,
  xlabel="process",ylabel="relative KLD",ylim=(0,5))
display(plt)

## Save figure
for ext in ["tex","pdf"]
  f = "cmp-estimators.$ext"
  savefig(plt,joinpath(figure_dir,f))
  @info "Saved $f"
end


text:
@testset "Vector Constraints - Nonnegatives - Parameter in affine part" begin
    """
        min x + y
            x - t + 1 >= 0
            y - t + 2 >= 0

        opt
            x* = t-1
            y* = t-2
            obj = 2*t-3
    """
    cached = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        ECOS.Optimizer(),
    )
    model = POI.Optimizer(cached)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    t, ct = MOI.add_constrained_variable(model, POI.Parameter(5))

    A = [1.0 0 -1; 0 1 -1]
    b = [1.0; 2]
    terms =
        MOI.VectorAffineTerm.(
            1:2,
            MOI.ScalarAffineTerm.(A, reshape([x, y, t], 1, 3)),
        )
    f = MOI.VectorAffineFunction(vec(terms), b)
    set = MOI.Nonnegatives(2)

    cnn = MOI.add_constraint(model, f, MOI.Nonnegatives(2))

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [y, x]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 4 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 3 atol = ATOL
    @test MOI.get(model, MOI.ConstraintPrimal(), cnn) ≈ [0.0, 0.0] atol = ATOL
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 7 atol = ATOL
    @test MOI.get(model, MOI.DualObjectiveValue()) ≈ 7 atol = ATOL
    @test MOI.get(model, MOI.ConstraintDual(), cnn) ≈ [1.0, 1.0] atol = ATOL

    MOI.set(model, POI.ParameterValue(), t, 6)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 5 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 4 atol = ATOL
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 9 atol = ATOL
end

@testset "JuMP direct model - Vector Constraints - Nonnegatives - Parameter in affine part" begin
    """
        min x + y
            x - t + 1 >= 0
            y - t + 2 >= 0

        opt
            x* = t-1
            y* = t-2
            obj = 2*t-3
    """

    cached = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        ECOS.Optimizer(),
    )
    optimizer = POI.Optimizer(cached)

    model = direct_model(optimizer)
    @variable(model, x)
    @variable(model, y)
    @variable(model, t in POI.Parameter(5))
    @constraint(model, [(x - t + 1), (y - t + 2)...] in MOI.Nonnegatives(2))
    @objective(model, Min, x + y)
    optimize!(model)

    @test isapprox.(value(x), 4.0, atol = ATOL)
    @test isapprox.(value(y), 3.0, atol = ATOL)

    MOI.set(model, POI.ParameterValue(), t, 6)
    optimize!(model)

    @test isapprox.(value(x), 5.0, atol = ATOL)
    @test isapprox.(value(y), 4.0, atol = ATOL)
end

@testset "Vector Constraints - Nonpositives - Parameter in affine part" begin
    """
        min x + y
            - x + t - 1 ≤ 0
            - y + t - 2 ≤ 0

        opt
            x* = t-1
            y* = t-2
            obj = 2*t-3
    """
    cached = MOI.Bridges.full_bridge_optimizer(
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{Float64}()),
            ECOS.Optimizer(),
        ),
        Float64,
    )
    model = POI.Optimizer(cached)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    t, ct = MOI.add_constrained_variable(model, POI.Parameter(5))

    A = [-1.0 0 1; 0 -1 1]
    b = [-1.0; -2]
    terms =
        MOI.VectorAffineTerm.(
            1:2,
            MOI.ScalarAffineTerm.(A, reshape([x, y, t], 1, 3)),
        )
    f = MOI.VectorAffineFunction(vec(terms), b)
    set = MOI.Nonnegatives(2)

    cnn = MOI.add_constraint(model, f, MOI.Nonpositives(2))

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [y, x]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 4 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 3 atol = ATOL
    @test MOI.get(model, MOI.ConstraintPrimal(), cnn) ≈ [0.0, 0.0] atol = ATOL
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 7 atol = ATOL
    @test MOI.get(model, MOI.DualObjectiveValue()) ≈ 7 atol = ATOL
    @test MOI.get(model, MOI.ConstraintDual(), cnn) ≈ [-1.0, -1.0] atol = ATOL

    MOI.set(model, POI.ParameterValue(), t, 6)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 5 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 4 atol = ATOL
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 9 atol = ATOL
end

@testset "JuMP direct model - Vector Constraints - Nonnegatives - Parameter in affine part" begin
    """
        min x + y
            - x + t - 1 ≤ 0
            - y + t - 2 ≤ 0

        opt
            x* = t-1
            y* = t-2
            obj = 2*t-3
    """

    cached = MOI.Bridges.full_bridge_optimizer(
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{Float64}()),
            ECOS.Optimizer(),
        ),
        Float64,
    )
    optimizer = POI.Optimizer(cached)

    model = direct_model(optimizer)
    @variable(model, x)
    @variable(model, y)
    @variable(model, t in POI.Parameter(5))
    @constraint(model, [(-x + t - 1), (-y + t - 2)...] in MOI.Nonpositives(2))
    @objective(model, Min, x + y)
    optimize!(model)

    @test isapprox.(value(x), 4.0, atol = ATOL)
    @test isapprox.(value(y), 3.0, atol = ATOL)

    MOI.set(model, POI.ParameterValue(), t, 6)
    optimize!(model)

    @test isapprox.(value(x), 5.0, atol = ATOL)
    @test isapprox.(value(y), 4.0, atol = ATOL)
end

@testset "Vector Constraints - SOC - Parameter in affine part" begin
    """
        Problem SOC2 from MOI

        min  x
        s.t. y ≥ 1/√2
            (x-p)² + y² ≤ 1

        in conic form:

        min  x
        s.t.  -1/√2 + y ∈ R₊
            1 - t ∈ {0}
            (t, x-p ,y) ∈ SOC₃

        opt
            x* = p - 1/√2
            y* = 1/√2
    """
    cached = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        ECOS.Optimizer(),
    )
    model = POI.Optimizer(cached)

    x, y, t = MOI.add_variables(model, 3)
    p, cp = MOI.add_constrained_variable(model, POI.Parameter(0))

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    cnon = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y))],
            [-1 / √2],
        ),
        MOI.Nonnegatives(1),
    )

    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, t))],
            [1.0],
        ),
        MOI.Zeros(1),
    )

    A = [
        1.0 0.0 0.0 0.0
        0.0 1.0 0.0 -1
        0.0 0.0 1.0 0.0
    ]
    f = MOI.VectorAffineFunction(
        vec(
            MOI.VectorAffineTerm.(
                1:3,
                MOI.ScalarAffineTerm.(A, reshape([t, x, y, p], 1, 4)),
            ),
        ),
        zeros(3),
    )
    csoc = MOI.add_constraint(model, f, MOI.SecondOrderCone(3))

    f_error = MOI.VectorOfVariables([t, p, y])
    @test_throws ErrorException MOI.add_constraint(
        model,
        f_error,
        MOI.SecondOrderCone(3),
    )

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ -1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), t) ≈ 1 atol = ATOL

    MOI.set(model, POI.ParameterValue(), p, 1)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 - 1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1 - 1 / √2 atol = ATOL
end

@testset "Vector Constraints - SOC - No parameters" begin
    """
        Problem SOC2 from MOI

        min  x
        s.t. y ≥ 1/√2
            x² + y² ≤ 1

        in conic form:

        min  x
        s.t.  -1/√2 + y ∈ R₊
            1 - t ∈ {0}
            (t, x ,y) ∈ SOC₃

        opt
            x* = 1/√2
            y* = 1/√2
    """
    cached = MOI.Bridges.full_bridge_optimizer(
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{Float64}()),
            ECOS.Optimizer(),
        ),
        Float64,
    )
    model = POI.Optimizer(cached)

    x, y, t = MOI.add_variables(model, 3)

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    cnon = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y))],
            [-1 / √2],
        ),
        MOI.Nonnegatives(1),
    )

    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, t))],
            [1.0],
        ),
        MOI.Zeros(1),
    )

    f = MOI.VectorOfVariables([t, x, y])
    csoc = MOI.add_constraint(model, f, MOI.SecondOrderCone(3))

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ -1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1 / √2 atol = ATOL
    @test MOI.get(model, MOI.VariablePrimal(), t) ≈ 1 atol = ATOL
end

@testset "JuMP direct model - Vector Constraints - SOC - Parameter in affine part" begin
    """
        Problem SOC2 from MOI

        min  x
        s.t. y ≥ 1/√2
            (x-p)² + y² ≤ 1

        in conic form:

        min  x
        s.t.  -1/√2 + y ∈ R₊
            1 - t ∈ {0}
            (t, x-p ,y) ∈ SOC₃

        opt
            x* = p - 1/√2
            y* = 1/√2
    """

    cached = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        ECOS.Optimizer(),
    )
    optimizer = POI.Optimizer(cached)

    model = direct_model(optimizer)
    @variable(model, x)
    @variable(model, y)
    @variable(model, t)
    @variable(model, p in POI.Parameter(0))

    @constraint(model, [y - 1 / √2] in MOI.Nonnegatives(1))
    @constraint(model, [t - 1] in MOI.Zeros(1))
    @constraint(model, [t, (x - p), y...] in SecondOrderCone())
    @objective(model, Min, 1.0 * x)
    optimize!(model)

    @test objective_value(model) ≈ -1 / √2 atol = ATOL
    @test value(x) ≈ -1 / √2 atol = ATOL
    @test value(y) ≈ 1 / √2 atol = ATOL
    @test value(t) ≈ 1 atol = ATOL

    MOI.set(model, POI.ParameterValue(), p, 1)
    optimize!(model)

    @test objective_value(model) ≈ 1 - 1 / √2 atol = ATOL
    @test value(x) ≈ 1 - 1 / √2 atol = ATOL
end


text:
struct softSAFTParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type softSAFTModel <: SAFTModel end
@newmodel softSAFT softSAFTModel softSAFTParam

export softSAFT
function softSAFT(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params,sites = getparams(components, ["SAFT/softSAFT"]; userlocations=userlocations, verbose=verbose)

    segment = params["m"]
    k = params["k"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]

    packagedparams = softSAFTParam(segment, sigma, epsilon, epsilon_assoc, bondvol)
    references = ["todo"]

    model = softSAFT(packagedparams, sites, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::softSAFTModel, V, T, z)
    return @f(a_LJ) + @f(a_chain) + @f(a_assoc)
end

function a_LJ(model::softSAFTModel, V, T, z)
    m = model.params.segment.values
    m̄ = ∑(z/∑(z) .* m)
    ϵ̄ = @f(ϵ_m)
    σ̄ = @f(σ_m)
    ρ̄ = @f(ρ_S)*σ̄^3
    T̄ = T/ϵ̄
    γ = 3
    F = exp(-γ*ρ̄^2)
    x = softSAFTconsts.x
    a = [
         x[1]*T̄+x[2]*√(T̄)+x[3]+x[4]/T̄+x[5]/T̄^2,
         x[6]*T̄+x[7]+x[8]/T̄+x[9]/T̄^2,
         x[11]+x[10]*T̄+x[12]/T̄,
         x[13],
         x[14]/T̄+x[15]/T̄^2,
         x[16]/T̄,
         x[17]/T̄+x[18]/T̄^2,
         x[19]/T̄^2,
        ]
    b = [
         x[20]/T̄^2+x[21]/T̄^3,
         x[22]/T̄^2+x[23]/T̄^4,
         x[24]/T̄^2+x[25]/T̄^3,
         x[26]/T̄^2+x[27]/T̄^4,
         x[28]/T̄^2+x[29]/T̄^3,
         x[30]/T̄^2+x[31]/T̄^3+x[32]/T̄^4,
        ]
    _0 = zero(V+T+first(z))
    G = [_0,_0,_0,_0,_0,_0]
    G[1] = (1-F)/(2γ)
    for i in 2:6
        k = 2*(i-1)
        G[i] = -(F*ρ̄^k - k*G[i-1]) / 2γ
    end
    return m̄*(∑(a[i]*ρ̄^i/i for i ∈ 1:8)+∑(b[i]*G[i] for i ∈ 1:6))/T̄
end

function ϵ_m(model::softSAFTModel, V, T, z)
    comps = @comps
    x = z/∑(z)
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    m = model.params.segment.values
    return ∑(m[i]*m[j]*x[i]*x[j]*σ[i,j]^3*ϵ[i,j] for i ∈ comps for j ∈ comps)/∑(m[i]*m[j]*x[i]*x[j]*σ[i,j]^3 for i ∈ comps for j ∈ comps)
end

function σ_m(model::softSAFTModel, V, T, z)
    comps = @comps
    x = z/∑(z)
    σ = model.params.sigma.values
    m = model.params.segment.values
    return (∑(m[i]*m[j]*x[i]*x[j]*σ[i,j]^3 for i ∈ comps for j ∈ comps)/∑(m[i]*m[j]*x[i]*x[j] for i ∈ comps for j ∈ comps))^(1/3)
end

function ρ_S(model::softSAFTModel, V, T, z)
    ∑z = ∑(z)
    N = N_A*∑z
    x = z/∑z
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return N/V*m̄
end

function a_chain(model::softSAFTModel, V, T, z)
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return -log(@f(y_LJ))*(m̄-1)
end

function y_LJ(model::softSAFTModel, V, T, z)
    ϵ̄ = @f(ϵ_m)
    gLJ = @f(g_LJ)
    return gLJ*exp(-ϵ̄/T)
end

function g_LJ(model::softSAFTModel, V, T, z)
    ϵ̄ = @f(ϵ_m)
    σ̄ = @f(σ_m)
    T̄ = T/ϵ̄
    ρ̄ = @f(ρ_S)*σ̄^3
    a = LJSAFTconsts.a
    return 1+∑(a[i,j]*ρ̄^i*T̄^(1-j) for i ∈ 1:5 for j ∈ 1:5)
end

function a_assoc(model::softSAFTModel, V, T, z)
    x = z/∑(z)
    X_ = @f(X)
    n = model.sites.n_sites
    return ∑(x[i]*∑(n[i][a]*(log(X_[i][a])+(1-X_[i][a])/2) for a ∈ @sites(i)) for i ∈ @comps)
end

function X(model::softSAFTModel, V, T, z)
    _1 = one(V+T+first(z))
    ∑z = ∑(z)
    x = z/∑z
    ρ = N_A*∑z/V
    n = model.sites.n_sites
    itermax = 500
    dampingfactor = 0.5
    error = 1.
    tol = model.absolutetolerance
    iter = 1
    X_ = [[_1 for a ∈ @sites(i)] for i ∈ @comps]
    X_old = deepcopy(X_)
    while error > tol
        iter > itermax && error("X has failed to converge after $itermax iterations")
        for i ∈ @comps, a ∈ @sites(i)
            rhs = 1/(1+∑(ρ*x[j]*∑(n[j][b]*X_old[j][b]*@f(Δ,i,j,a,b) for b ∈ @sites(j)) for j ∈ @comps))
            X_[i][a] = (1-dampingfactor)*X_old[i][a] + dampingfactor*rhs
        end
        error = sqrt(∑(∑((X_[i][a] - X_old[i][a])^2 for a ∈ @sites(i)) for i ∈ @comps))
        for i = 1:length(X_)
            X_old[i] .= X_[i]
        end
        iter += 1
    end
    return X_
end

function Δ(model::softSAFTModel, V, T, z, i, j, a, b)
    ϵ̄ = @f(ϵ_m)
    σ̄ = @f(σ_m)
    T̄ = T/ϵ̄
    ρ̄ = @f(ρ_S)*σ̄^3
    ϵ_assoc = model.params.epsilon_assoc.values
    κ = model.params.bondvol.values
    b_ = LJSAFTconsts.b

    I = ∑(b_[i+1,j+1]*ρ̄^i*T̄^j for i ∈ 0:4 for j ∈ 0:4)/3.84/1e4
    return 4π*(exp(ϵ_assoc[i,j][a,b]/T)-1)*κ[i,j][a,b]*I
end

const softSAFTconsts =
(
 x = [ 0.8623085097507421, 2.976218765822098,-8.402230115796038,0.1054136629203555,
      -0.8564583828174598, 1.582759470107601,0.7639421948305453, 1.753173414312048,
      2.798291772190376e3,-4.8394220260857657e-2,0.9963265197721935,-3.698000291272493e1,
      2.084012299434647e1, 8.305402124717285e1,-9.574799715203068e2,-1.477746229234994e2,
      6.398607852471505e1, 1.603993673294834e1, 6.805916615864377e1,-2.791293578795945e3,
      -6.245128304568454, -8.116836104958410e3, 1.488735559561229e1,-1.059346754655084e4,
      -1.131607632802822e2,-8.867771540418822e3,-3.986982844450543e1,-4.689270299917261e3,
      2.593535277438717e2,-2.694523589434903e3,-7.218487631550215e2, 1.721802063863269e2],

 a = [0.49304346593882 2.1528349894745 -15.955682329017 24.035999666294 -8.6437958513990;
      -0.47031983115362 1.1471647487376  37.889828024211 -84.667121491179 39.643914108411;
      5.0325486243620 -25.915399226419 -18.862251310090 107.63707381726 -66.602649735720;
      -7.3633150434385  51.553565337453 -40.519369256098 -38.796692647218 44.605139198378;
      2.9043607296043 -24.478812869291  31.500186765040 -5.3368920371407 -9.5183440180133],

 b = [-0.03915181 0.08450471 0.06889053 -0.01034279  0.5728662e-3;
      -0.5915018  0.9838141 -0.4862279   0.1029708 -0.6919154e-2;
      1.908368  -3.415721   2.124052  -0.4298159  0.02798384;
      -0.7957312  0.7187330 -0.9678804   0.2431675 -0.01644710;
      -0.9399577   2.314054 -0.4877045  0.03932058 -0.1600850e-2],
)


text:
struct IRDump
    id::Int
    var"function"::Union{String,Nothing}
    pass::Union{String,Nothing}
    code::String
end

struct IRCollection <: AbstractVector{IRDump}
    var"function"::Union{String,Nothing}
    dumps::Vector{IRDump}
end

const FunctionCollectionDict = Dict{Union{String,Nothing},IRCollection}

struct FunctionCollection <: AbstractDict{String,IRCollection}
    dumps::FunctionCollectionDict
    nirs::typeof(Ref(0))
end

FunctionCollection() = FunctionCollection(FunctionCollectionDict(), Ref(0))

Base.getindex(irs::IRCollection, i::Int) = irs.dumps[i]
Base.size(irs::IRCollection) = size(irs.dumps)

Base.length(fns::FunctionCollection) = length(fns.dumps)
Base.keys(fns::FunctionCollection) = keys(fns.dumps)
Base.values(fns::FunctionCollection) = values(fns.dumps)
Base.getindex(fns::FunctionCollection, name::AbstractString) = fns.dumps[name]

function Base.getindex(fns::FunctionCollection, pattern::AbstractPattern)
    dumps = FunctionCollectionDict()
    nirs = 0
    for irs in values(fns)
        irs.function === nothing && continue
        if match(pattern, irs.function) !== nothing
            dumps[irs.function] = irs
            nirs += length(irs)
        end
    end
    return FunctionCollection(dumps, Ref(nirs))
end

function Base.push!(fns::FunctionCollection, ir::IRDump)
    irs = get!(fns.dumps, ir.function) do
        IRCollection(ir.function, IRDump[])
    end::IRCollection
    push!(irs.dumps, ir)
    fns.nirs[] += 1
    return fns
end

function Base.iterate(it::Union{FunctionCollection,IRCollection}, state = iterate(it.dumps))
    state === nothing && return nothing
    item, istate = state
    return (item, iterate(it.dumps, istate))
end

# function Base.summary(io::IO, irs::IRCollection)
#     print(io, "Collection of ", length(irs.dumps), " IRs")
# end

# function Base.show(io::IO, ::MIME"text/plain", irs::IRCollection)
#     summary(io, irs)
# end

# function Base.summary(io::IO, irs::FunctionCollection)
#     print(io, "Collection of ", length(irs.dumps), " functions and ", irs.nirs[], " IRs")
# end

# function Base.show(io::IO, ::MIME"text/plain", irs::FunctionCollection)
#     summary(io, irs)
# end

function Base.summary(io::IO, ir::IRDump)
    print(io, "*** IR Dump of ")
    printstyled(io, something(ir.function, "UNKNOWN FUNCTION"); color = :blue, bold = true)
    print(io, " After ")
    printstyled(io, something(ir.pass, "UNKNOWN PASS"); color = :blue, bold = true)
    print(io, " ***")
end

Base.print(ir::IRDump; options...) = print(stdout, ir; options...)
function Base.print(io::IO, ir::IRDump; raw::Bool = true, dump_module::Bool = true)
    summary(io, ir)
    println(io)
    code = stripoff_passname(ir.code)
    if !dump_module
        code = extract_function(code)
    end
    if !raw
        code = cleanup_ir(code)
    end
    print_llvm(io, code)
end

function Base.show(io::IO, ::MIME"text/plain", ir::IRDump)
    print(io, ir; raw = false, dump_module = false)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libevdev"
version = v"1.11.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://gitlab.freedesktop.org/libevdev/libevdev", "8855f1ac59a69d5bdff51e3f8980697f0127c270"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libevdev/
atomic_patch ../patches/link_rt.patch
mkdir build
cd build

meson -D tests=disabled -D documentation=disabled ../ --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libevdev", :libevdev),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Julia wrapper for header: /usr/include/scip/pub_sepa.h
# Automatically generated using Clang.jl wrap_c

function SCIPsepaComp(elem1, elem2)
    ccall((:SCIPsepaComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPsepaCompName(elem1, elem2)
    ccall((:SCIPsepaCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPsepaGetData(sepa)
    ccall((:SCIPsepaGetData, libscip), Ptr{SCIP_SEPADATA}, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaSetData(sepa, sepadata)
    ccall((:SCIPsepaSetData, libscip), Cvoid, (Ptr{SCIP_SEPA}, Ptr{SCIP_SEPADATA}), sepa, sepadata)
end

function SCIPsepaGetName(sepa)
    ccall((:SCIPsepaGetName, libscip), Cstring, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetDesc(sepa)
    ccall((:SCIPsepaGetDesc, libscip), Cstring, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetPriority(sepa)
    ccall((:SCIPsepaGetPriority, libscip), Cint, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetFreq(sepa)
    ccall((:SCIPsepaGetFreq, libscip), Cint, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaSetFreq(sepa, freq)
    ccall((:SCIPsepaSetFreq, libscip), Cvoid, (Ptr{SCIP_SEPA}, Cint), sepa, freq)
end

function SCIPsepaGetMaxbounddist(sepa)
    ccall((:SCIPsepaGetMaxbounddist, libscip), Cdouble, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaUsesSubscip(sepa)
    ccall((:SCIPsepaUsesSubscip, libscip), UInt32, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetSetupTime(sepa)
    ccall((:SCIPsepaGetSetupTime, libscip), Cdouble, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetTime(sepa)
    ccall((:SCIPsepaGetTime, libscip), Cdouble, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNCalls(sepa)
    ccall((:SCIPsepaGetNCalls, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNCallsAtNode(sepa)
    ccall((:SCIPsepaGetNCallsAtNode, libscip), Cint, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNCutoffs(sepa)
    ccall((:SCIPsepaGetNCutoffs, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNCutsFound(sepa)
    ccall((:SCIPsepaGetNCutsFound, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNCutsApplied(sepa)
    ccall((:SCIPsepaGetNCutsApplied, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNCutsFoundAtNode(sepa)
    ccall((:SCIPsepaGetNCutsFoundAtNode, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNConssFound(sepa)
    ccall((:SCIPsepaGetNConssFound, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaGetNDomredsFound(sepa)
    ccall((:SCIPsepaGetNDomredsFound, libscip), Clonglong, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaIsDelayed(sepa)
    ccall((:SCIPsepaIsDelayed, libscip), UInt32, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaWasLPDelayed(sepa)
    ccall((:SCIPsepaWasLPDelayed, libscip), UInt32, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaWasSolDelayed(sepa)
    ccall((:SCIPsepaWasSolDelayed, libscip), UInt32, (Ptr{SCIP_SEPA},), sepa)
end

function SCIPsepaIsInitialized(sepa)
    ccall((:SCIPsepaIsInitialized, libscip), UInt32, (Ptr{SCIP_SEPA},), sepa)
end


text:
module ROC

using Infinity
using Missings # to use Missings.nonmissingtype
using RecipesBase # for creating a Plots.jl recipe

export	ROCData,
		roc,
		AUC,
		PPV,
		cutoffs

include("rocdata.jl")
include("roc_main.jl")
include("rocplot.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 4655bdc2-1ed7-11ec-022d-7df0b8fdb907
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 661cc8fe-44f3-422a-91f5-96956b748847
begin
	using SqState
	using Plots
	using QuantumStateBase
	using Flux
	using LinearAlgebra
end

# ╔═╡ c55d7af1-0cd3-4b0a-bbd5-0048f15ea851
m = SqState.get_model("model_args_sqth_th")

# ╔═╡ 8f82e42f-8c2c-4bd0-ab20-95b0930670dc
# args = (0.9,  4.17646, 0.4);
# args = SqState.rand_arg_sqth_th((0, 1), (0, 2π), (0, 1))
args = (1., π/3, 0.4, 0.07, 0.8, 0.2)

# ╔═╡ f0cfd6f3-1e55-46cf-bd58-937e52cb7daa
begin
	dim = 100
	state = SqState.construct_state_sqth_th(args..., 1000)
	d = Float32.(rand(state, 4096, IsGaussian))
	argŝ = m(reshape(d[2, :], :, 1, 1))
	statê = SqState.construct_state_sqth_th(argŝ..., 1000)
	ρ = statê.𝛒
	argŝ
end

# ╔═╡ 57e75729-4bf4-4e0c-9ef2-bdaf6fb4e22b
function fidelity(ρ1, ρ2)
	sqrt_ρ1 = sqrt(ρ1)
	return tr(sqrt(sqrt_ρ1*ρ2*sqrt_ρ1))^2
end

# ╔═╡ 6739ac26-7590-46a8-a75a-2160a1d7c11d
fidelity(ρ, state.𝛒)

# ╔═╡ 00c09086-1eef-4b3f-b764-fe7201e1a354
scatter(d[1, :], d[2, :])

# ╔═╡ be08410c-4fce-405c-9b27-ababf152402e
begin
	d̂ = rand(statê, 4096, IsGaussian)
	scatter(d̂[1, :], d̂[2, :])
end

# ╔═╡ 32c40e60-0989-44f6-b518-9f6203aa6209

# ╔═╡ Cell order:
# ╠═4655bdc2-1ed7-11ec-022d-7df0b8fdb907
# ╠═661cc8fe-44f3-422a-91f5-96956b748847
# ╠═c55d7af1-0cd3-4b0a-bbd5-0048f15ea851
# ╠═8f82e42f-8c2c-4bd0-ab20-95b0930670dc
# ╠═f0cfd6f3-1e55-46cf-bd58-937e52cb7daa
# ╟─57e75729-4bf4-4e0c-9ef2-bdaf6fb4e22b
# ╟─6739ac26-7590-46a8-a75a-2160a1d7c11d
# ╠═00c09086-1eef-4b3f-b764-fe7201e1a354
# ╠═be08410c-4fce-405c-9b27-ababf152402e
# ╠═32c40e60-0989-44f6-b518-9f6203aa6209


text:
"""
    modelflops(model::Chain)::Float64

Calculates the approximate number of Floating Point Operations that the model will require

```julia
weight = rand(Float64, 3, 3, 5)
bias = zeros(Float64, 5)
Conv(weight = weight,
    bias = bias,
    σ = sigmoid)
```
"""
function modelflops(model::Chain,inputsize::Tuple,samplesize::Float64,batchsize::Float64)
    x = 0; y = 0; Fm = 1
    if length(inputsize) == 3
        x,y,Fm = inputsize
    elseif length(inputsize) == 2
        x,y = inputsize
    end
    modellayers = collect(model)
    nelayers = Array{Any,1}()
    layeroutput = Array{Tuple,1}()
    outsizes = Array{Tuple,1}()
    lossandgradient = 1
    for ml in modellayers
        if islayer(ml)
            push!(nelayers,ml)
        end
    end
    output = (0,0)
    for mli in 1:length(nelayers)
        if mli == 1
            noflops, output = layerflops(nelayers[mli],inputsize)
            layeroutput = vcat(layeroutput, noflops)
        else
            noflops, output = layerflops(nelayers[mli],output)
            layeroutput = vcat(layeroutput, noflops)
        end
    end
    numberoflayers = length(layeroutput)
    return sum(layeroutput) * samplesize + lossandgradient*batchsize
end


text:
#=
  Run this script as follows:
  $ cd /path/to/Schemata.jl
  $ /path/to/julia scripts/compare.jl /path/to/config.toml /path/to/inputdata sorted_by_primarykey
  The 3rd argument, sorted_by_primarykey is either "true" or "false".
  If "true" the compare function assumes that your table is sorted by its primary key,
  which enables a faster comparison to the schema to be made.
=#

configfile      = ARGS[1]
input_data_file = ARGS[2]
sorted_by_primarykey = false
if ARGS[3] == "true"
    sorted_by_primarykey = true
end

using Pkg
Pkg.activate(".")

using Dates
using Schemata

# Construct the TableSchema
ts     = nothing
schema = readschema(configfile)
if schema isa Schema
    if isnothing(tablename)
        length(schema.tables) != 1 && error("The schema has more than 1 table. Please specify a table name as a 4th command line argument.")
        for (tablename, tableschema) in schema.tables
            ts = tableschema
            break
        end
    else
        ts = schema.tables[tablename]
    end
elseif schema isa TableSchema
    ts = schema
else
    error("The schema is neither a Schema or a TableSchema.")
end

# Compare data to schema
println("$(now()) Starting comparison.")
compare(ts, input_data_file; sorted_by_primarykey=sorted_by_primarykey)
println("$(now()) A table of transformed data, input issues and output issues have been stored at $(dirname(input_data_file)).")


text:
"""
    ForwardDiffBackend <: AbstractDiffBackend

Differentiation backend based on the ForwardDiff.jl package.
"""
struct ForwardDiffBackend <: AbstractDiffBackend end

function Manifolds._derivative(f, p, ::ForwardDiffBackend)
    return ForwardDiff.derivative(f, p)
end

function _derivative!(f, X, t, ::ForwardDiffBackend)
    return ForwardDiff.derivative!(X, f, t)
end

function _gradient(f, p, ::ForwardDiffBackend)
    return ForwardDiff.gradient(f, p)
end

function _gradient!(f, X, t, ::ForwardDiffBackend)
    return ForwardDiff.gradient!(X, f, t)
end

function _jacobian(f, p, ::ForwardDiffBackend)
    return ForwardDiff.jacobian(f, p)
end

if default_differential_backend() === NoneDiffBackend()
    set_default_differential_backend!(ForwardDiffBackend())
end


text:
using QMatrices
using LinearAlgebra: LinearAlgebra

@testset "Basic algebraic props" begin
    for m in (X, Y, Z, H)
        @test isapprox(m^2, LinearAlgebra.I) # TODO: (slow) norm is unnecessary
    end
end

@testset "control" begin
    toffoli = [1 0 0 0 0 0 0 0;
               0 1 0 0 0 0 0 0;
               0 0 1 0 0 0 0 0;
               0 0 0 1 0 0 0 0;
               0 0 0 0 1 0 0 0;
               0 0 0 0 0 1 0 0;
               0 0 0 0 0 0 0 1;
               0 0 0 0 0 0 1 0]

    @test control(X, 2) == toffoli
end

@testset "Bell state" begin
    k = zeroket(2)
    @test CX * (H ⊗ I2) * k == [0.7071067811865475, 0.0, 0.0, 0.7071067811865475]
end


text:
module JuMPTutorials

using Weave

const repo_directory = joinpath(@__DIR__,"..")

"""
    weave_file(folder::AbstractString, file::AbstractString)

Checks if the file has been modified since the last Weave and updates the notebook and tests accordingly.

* `folder` = Name of the folder the tutorial is in
* `file` = Name of the tutorial file
"""
function weave_file(folder,file)
    cd(joinpath(repo_directory, "script", folder))

    filename = split(file, ".")[1]
    srcpath = joinpath(repo_directory, "script", folder, file)
    testpath = joinpath(repo_directory, "test", folder, file)
    notebookpath = joinpath(repo_directory, "notebook", folder)

    if mtime(srcpath) > mtime(testpath) || mtime(testpath)==0
        @warn "Updating tests for $filename as it has been updated since the last weave."
        tangle(srcpath, out_path=testpath)
    else
        @warn "Skipping tests for $filename as it has not been updated."
    end

    if mtime(srcpath) > mtime(notebookpath) || mtime(notebookpath)==0
        @warn "Weaving $filename to Jupyter Notebook as it has been updated since the last weave."
        notebook(srcpath, out_path=notebookpath, timeout=-1, nbconvert_options="--allow-errors")
    else
        @warn "Skipping Jupyter Notebook for $filename as it has not been updated."
    end

    cd(joinpath(repo_directory,"src"))
end

"""
    weave_file(folder::AbstractString)

Checks the files present in the specified folder for modifications and updates the notebook and tests accordingly.

* `folder` = Name of the folder to check
"""
function weave_folder(folder)
    for file in readdir(joinpath(repo_directory, "script", folder))
        if endswith(file, ".jl")
            println("")
            println("Building $(joinpath(folder,file))")
            try
                weave_file(folder,file)
            catch e
                @warn("Error weaving $folder/$file:\n$e")
            end
            println("")
        end
    end
end

"""
    weave_all()

Checks every tutorial for modifications and updates the notebook and tests accordingly.
"""
function weave_all()
    for folder in readdir(joinpath(repo_directory, "script"))
        weave_folder(folder)
    end
end

"""
    weave_file_f(folder::AbstractString, file::AbstractString)

Use Weave to convert the given file irrespective of whether it has been modified or not.

* `folder` = Name of the folder the tutorial is in
* `file` = Name of the tutorial file
"""
function weave_file_f(folder,file)
    cd(joinpath(repo_directory,"script",folder))

    srcpath = joinpath(repo_directory, "script", folder, file)
    testpath = joinpath(repo_directory, "test", folder, file)
    notebookpath = joinpath(repo_directory, "notebook", folder)

    tangle(srcpath, out_path = testpath)
    notebook(srcpath, out_path=notebookpath, timeout=-1, nbconvert_options="--allow-errors")

    cd(joinpath(repo_directory,"src"))
end

"""
`weave_folder_f(folder::AbstractString, file::AbstractString)`

Use Weave to convert every file in the specified folder irrespective of whether it has been modified or not.

* `folder` = Name of the folder the tutorial is in
"""
function weave_folder_f(folder)
    for file in readdir(joinpath(repo_directory, "script", folder))
        if endswith(file, ".jl")
            println("")
            println("Building $(joinpath(folder,file))")
            try
                weave_file_f(folder,file)
            catch e
                @warn("Error weaving $(folder)/$(file):\n$(e)")
            end
            println("")
        end
    end
end

"""
`weave_all_f()`

Use Weave to convert every tutorial irrespective of whether it has been modified or not.
"""
function weave_all_f()
    for folder in readdir(joinpath(repo_directory, "script"))
        weave_folder_f(folder)
    end
end

end