text:
lp_solvers = Any[]
ip_solvers = Any[]
ip_dual_solvers = Any[]
semi_solvers = Any[]
sos_solvers = Any[]
lazy_solvers, lazy_soc_solvers, lazylocal_solvers, cut_solvers, cutlocal_solvers, heur_solvers, info_solvers = Any[], Any[], Any[], Any[], Any[], Any[], Any[]
quad_solvers = Any[]
quad_soc_solvers = Any[]
rsoc_solvers = Any[]
nlp_solvers = Any[]
convex_nlp_solvers = Any[]
minlp_solvers = Any[]
sdp_solvers = Any[]

lp_solvers = ["baron" "knitro"]

lp_exclude = Dict("knitro" => ["linear8c", "linear8b"])

quad_solvers = [("baron", "nodual"),
                ("knitro", ""),
               ]
quad_exclude = Dict("baron" => ["qcp1"], # SLOW_PROGRESS
                    )
# mosek fails?
global_solvers = [("knitro", "")]

# fails with "'gmsbaxnx.exe' is not recognized as an internal or external command, operable program or batch file"
if !Sys.iswindows()
    push!(global_solvers, ("baron", ""))
end

push!(nlp_solvers, "baron")
push!(nlp_solvers, "knitro")

append!(convex_nlp_solvers, nlp_solvers)

push!(minlp_solvers, "")

# Default configuration.
const config_mosek = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                               optimal_status=MOI.OPTIMAL,
                               query=false,
                               infeas_certificates=false, # Do not ask for infeasibility certificates.
                               modify_lhs=false)

const config_mosek_noduals = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                                       optimal_status=MOI.OPTIMAL,
                                       query=false,
                                       duals=false,
                                       infeas_certificates=false,
                                       modify_lhs=false)

# Default configuration.
const config_knitro = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                               optimal_status=MOI.LOCALLY_SOLVED,
                               query=false,
                               infeas_certificates=false, # Do not ask for infeasibility certificates.
                               modify_lhs=false)

const config_knitro_noduals = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                                       optimal_status=MOI.LOCALLY_SOLVED,
                                       query=false,
                                       duals=false,
                                       infeas_certificates=false,
                                       modify_lhs=false)

config_solver = Dict(
"knitro" => config_knitro
)

config_solver_noduals = Dict(
"knitro" => config_knitro_noduals
)


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# contributed by David Campbell
# based on the Go version
# modified by Jarret Revels, Alex Arslan, Yichao Yu
#
# Bit-twiddle optimizations added by Kristoffer Carlsson

using Printf

struct KNucleotides
    i::UInt64
end
Base.hash(kn::KNucleotides, h::UInt64) = hash(kn.i, h)
Base.isequal(kn1::KNucleotides, kn2::KNucleotides) = kn1.i == kn2.i
KNucleotides() = KNucleotides(0)
Base.show(io::IO, kn::KNucleotides) = print(io, '[', string(kn), ']')

const NucleotideLUT = zeros(UInt8, 256)
NucleotideLUT['A'%UInt8] = 0
NucleotideLUT['C'%UInt8] = 1
NucleotideLUT['G'%UInt8] = 2
NucleotideLUT['T'%UInt8] = 3

length_bit(l) = (UInt64(1) << 2l)

# Should be called after the nucelotides are added!
set_length(kn::KNucleotides, l::Integer) = KNucleotides(kn.i | length_bit(l))
add_nucleotide(kn::KNucleotides, c::UInt8) = KNucleotides(@inbounds (kn.i << 2) | NucleotideLUT[c])
Base.length(kn::KNucleotides) = (64 - leading_zeros(kn.i) - 1) ÷ 2
function KNucleotides(str::String, n=length(str), offset=0)
    # @assert isascii(str) && n <= 29
    kn = KNucleotides()
    @inbounds for i in 1:n
        kn = add_nucleotide(kn, codeunit(str, i + offset))
    end
    kn = set_length(kn, n)
    return kn
end

function Base.string(kn::KNucleotides)
    l = length(kn)
    i = kn.i - length_bit(l)
    sprint() do io
        for j in 1:l
            iiii = i
            mask = 3 << (2(l-j))
            z = (i & mask) >> 2(l-j)
            write(io,
                z == 0 ? 'A' :
                z == 1 ? 'C' :
                z == 2 ? 'G' :
                z == 3 ? 'T' :
                error())
        end
    end
end

function count_data(data::String, n::Int)
    counts = Dict{KNucleotides, Int}()
    top = length(data) - n
    @inbounds for offset = 0:top
        kn = KNucleotides(data, n, offset)
        token = Base.ht_keyindex2!(counts, kn)
        if token > 0
            counts.vals[token] += 1
        else
            Base._setindex!(counts, 1, kn, -token)
        end
    end
    return counts
end

function count_one(data::String, s::String)
    k = KNucleotides(s)
    d = count_data(data, length(s))
    return haskey(d, k) ? d[k] : 0
end

struct KNuc
    name::String
    count::Int
end

# sort down
function Base.isless(x::KNuc, y::KNuc)
    if x.count == y.count
        return x.name > y.name
    end
    x.count > y.count
end

function sorted_array(m)
    kn = Vector{KNuc}(undef, length(m))
    i = 1
    for (k, v) in m
        kn[i] = KNuc(string(k), v)
        i += 1
    end
    sort!(kn)
end

function print_knucs(a::Array{KNuc, 1})
    sum = 0
    for kn in a
        sum += kn.count
    end
    for kn in a
        @printf("%s %.3f\n", kn.name, 100.0kn.count/sum)
    end
    println()
end

do_work(str::String, i::Int) = sorted_array(count_data(str, i))
do_work(str::String, i::String) = count_one(str, i)

function perf_k_nucleotide(io = stdin)
    three = ">THREE "
    while true
        line = readline(io)
        if startswith(line, three)
            break
        end
    end
    data = read(io, String)
    str = uppercase(data)
    str = filter(!isequal('\n'), str)

    vs = [1, 2, "GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"]
    results = Vector{Any}(undef, length(vs))
    Threads.@threads for i in 1:length(vs)
        results[i] = do_work(str, vs[i])
    end

    for (v, result) in zip(vs, results)
        if result isa Array
            print_knucs(result)
        end
        if result isa Int
            @printf("%d\t%s\n", result, v)
        end
    end
end

perf_k_nucleotide()

ext: jl


text:
module RandomXorshiro

export AbstractFastRNG, AbstractXoshroRNG,
       srand, rand, rand!, jump

import Random: AbstractRNG, srand, rand, rand!

abstract type AbstractFastRNG   <: AbstractRNG      end
abstract type AbstractXoshroRNG <: AbstractFastRNG  end

abstract type AbstractInternalState end

const IUI32 = Union{Int32, UInt32}
const IUI64 = Union{Int64, UInt64}

include("rotatebits.jl")
include("states.jl")
include("seed.jl")
include("pseudorandom.jl")
include("randvecs.jl")

end # RandomXorshiro


text:
using DifferentialEquations
using Plots

# Parameters
const f = 2.0 # depensation in piscivores (B)
const m = 0.01 # yr⁻¹ natural piscivore mortality
const a = 0.37 # B recruitment assymptote (shape parameter)
const r = 8.7 # mg⋅m⁻²⋅d⁻¹ maximum recycling of P
const q = 8 # Or 2, steepness coefficient of R
const k = 90 # mg⋅m⁻² half saturation constant of P recycling
const s = 0.001 # d⁻¹ flushing rate
const b = 0.5 # adult density when recruitment = a/2 (B recruitment shape parameter)
const α = 0.8 # mg⋅m⁻²⋅d⁻¹ maximum growth rate of A per P input rate
const v = 0.07 # m⋅d⁻¹ sinking rate per algal cell
const g = 0.1 # d⁻¹ C grazing coefficient

# Drivers
# H yr⁻¹ harvest rate
H = 0.1 # Transition above 0.367
# DOC mg/m2 dissolved organic carbon
DOC = 6.0 # or 18
# L mg⋅m⁻²⋅d⁻¹ P loading
L = 2.0 # [0-6]

# State variables
# B g/m2 piscivore population biomass (pike)
## W g/m2 planktivore population biomass (bream)
#W = @. 3.2*exp(-3.6926*B)
## C mm zooplankton community length
#C = W < 2.58 ? 0.5 : 0.35 # mm
## Z m depth of thermocline
z = 13.1826/(1 + DOC)^0.61 # log(z) = 1.12 - 0.61log(DOC+1)
## P mg/m2 phosphorus
## R mg/m2 recycled phosphorus
## A mg/m2 chlorophyll

function beisnerd!(du,u,h,p,t)
    A, P, B = u
    α, v, z, g, s, r, q, a, f, b, m, H, L, tau = p

    R = (r*(A^q))/(k^q + A^q)
    W = 3.2*exp(-3.6926*B)
    C = W < 2.58 ? 0.5 : 0.35 # mm

    Blag3 = h(p, t-tau)[1]

    du[1] = dA = α*P*A - ((v/z) + g*C + s)*A
    du[2] = dP = L + R - α*P*A - s*P
    du[3] = dB = B*exp(a*Blag3^f/(b^f+Blag3^f))-m*B-H*B
end

u0 = [1.0;0.0;1.5]
h(p,t) = u0[3] # assume B0 for initial delay values
tau = 3
lags = [tau]
p = (α, v, z, g, s, r, q, a, f, b, m, H, L, tau)
tspan = (0.0,10.0)
prob = DDEProblem(beisnerd!,u0,h,tspan,p;constant_lags=lags)
sol = solve(prob, MethodOfSteps(Tsit5()))

plot(sol.t, [a[1] for a in sol.u])
plot!(sol.t, [a[2] for a in sol.u])


text:
using MPCDocs
using Base.Test


text:
module Perftests

export get_perf_groups, run_perf_groups, @perf, meta

using Benchmarks
using Compat

# This is where we expect to get our performance regression code
const perfdir = abspath(joinpath(dirname(@__FILE__),"../benchmarks/"))

# This is the default group, set by run_perf_groups() below
default_group = ""

# This is the directory where we'll plop out our .csv files
const resultsdir = joinpath(perfdir, "../test/results-$(Base.GIT_VERSION_INFO.commit)")

# This object represents metadata about a perf test; group/name/variant is a hierarchical
# organization of the perf tests where name is the only required member of the triplet; group
# defaults to the name of the directory holding the test, variant defaults to the empty string.
# Description is a human-readable description of the test, and issue is the JuliaLang github repo
# issue (if any) related to this test.
immutable meta
    group::AbstractString
    name::AbstractString
    variant::AbstractString

    description::AbstractString
    issue::UInt32

    meta(name; issue=0) = new(default_group, name, "", "", issue)
    meta(name, desc; issue=0) = new(default_group, name, "", desc, issue)
    meta(name, variant, desc; issue=0) = new(default_group, name, variant, desc, issue)
    meta(group, name, variant, desc; issue=0) = new(group, name, variant, desc, issue)
end

# This macro takes in a test expression and a PerfMetadata object
macro perf(ex, outer_meta)
    quote
        meta = $(esc(outer_meta))
        if contains(meta.group, "-") || contains(meta.name, "-") || contains(meta.variant, "-")
            throw(ArgumentError("Benchmark group/name/variant cannot contain '-'!"))
        end
        result = @benchmark $(esc(ex))
        stats = Benchmarks.SummaryStatistics(result)
        pts = Benchmarks.pretty_time_string
        time_str = "$(pts(stats.elapsed_time_center))"
        if !isnull(stats.elapsed_time_lower) && !isnull(stats.elapsed_time_upper)
            lower = stats.elapsed_time_lower.value
            upper = stats.elapsed_time_upper.value
            time_str *= " [$(pts(lower)), $(pts(upper))]"
        end
        if length(meta.variant) > 0
            csvpath = "$resultsdir/$(meta.group)-$(meta.name)-$(meta.variant).csv"
            println("  $(meta.group)/$(meta.name)/$(meta.variant) done in $time_str")
        else
            csvpath = "$resultsdir/$(meta.group)-$(meta.name).csv"
            println("  $(meta.group)/$(meta.name) done in $time_str")
        end
        writecsv(joinpath(perfdir, csvpath), result.samples)
    end
end

# Let's find all of our test groups; that's any directory inside benchmarks/
# us that has a perf.jl in it that we can call:
function get_perf_groups()
    return filter(x -> isfile(joinpath(perfdir,x,"perf.jl")), readdir(perfdir))
end

function run_perf_groups()
    return run_perf_groups(get_perf_groups())
end

function run_perf_groups(perf_groups)
    global default_group
    # Write out environment to $resultsdir/env.csv
    cd(perfdir)
    try mkdir(resultsdir) end
    writecsv(joinpath(perfdir, resultsdir, "env.csv"), Benchmarks.Environment())

    # Iterate, my friend.  Iterate.
    for dir in perf_groups
        # Set the default group first
        default_group = dir
        start_time = time()
        println("Running $dir/perf.jl...")
        include(joinpath(perfdir,dir,"perf.jl"))

        ns_elapsed = (time() - start_time)*1e9
        println("Finished $dir tests in $(Benchmarks.pretty_time_string(ns_elapsed))")
    end
end

end # module


text:
export reductions, Reductions, is_first_visit, IsFirstVisit, countmap

import Base:iterate, length, size, eltype
import StatsBase:countmap

"extend the countmap in StatsBase to support general iterator"
function countmap(iter)
    res = Dict{eltype(iter), Int}()
    for x in iter
        if haskey(res, x)
            res[x] += 1
        else
            res[x] = 1
        end
    end
    res
end

"""
    struct Reductions{T<:Union{NamedTuple{(:init,)}, NamedTuple{()}}, F, I}
        op::F
        itr::I
        kw::T
    end

Construct an iterator to get all the intermediate values while calling `reduce(op, itr;kw...)`
"""
struct Reductions{T<:Union{NamedTuple{(:init,)}, NamedTuple{()}}, F, I}
    op::F
    itr::I
    kw::T
end

Reductions(op, itr) = Reductions(op, itr, NamedTuple())

"""
    reductions(op, iter; init)

Return an Iterator of the intermediate values of the reduction
(as per reduce) of `iter` by `op`.

!!! note
    You can not apply `Iterators.reverse` to `Reductions` (due to time complexity).
    If you really want to get the reversed `Reductions`, consider `collect` first and then call `reverse!`

# Example
```julia
julia> reductions(+, [2,3,4])
3-element Array{Int64,1}:
 2
 5
 9

julia> reductions(+, [2,3,4], init=3)
Reductions{NamedTuple{(:init,),Tuple{Int64}},typeof(+),Array{Int64,1}}(+, [2, 3, 4], (init = 3,))

julia> collect(reductions(+, [2,3,4], init=3))
4-element Array{Any,1}:
  3
  5
  8
 12
```
"""
reductions(op, itr; kw...) = Reductions(op, itr, kw.data)

function iterate(r::Reductions{<:NamedTuple{()}})
    y = iterate(r.itr)
    if y === nothing
        throw(ArgumentError("reductions over an empty collection without init value is not allowed"))
    end
    (v, s) = y
    (v, (v, s))
end

function iterate(r::Reductions{<:NamedTuple{(:init,)}})
    v = r.kw.init
    (v, (v,))
end

function iterate(r::Reductions, (v,)::Tuple{Any})
    y = iterate(r.itr)
    y === nothing && return nothing
    v = r.op(v, y[1])
    (v, (v, y[2]))
end

function iterate(r::Reductions, (v, s)::Tuple{Any, Any})
    y = iterate(r.itr, s)
    y === nothing && return nothing
    v = r.op(v, y[1])
    (v, (v, y[2]))
end

length(r::Reductions{<:NamedTuple{()}}) = length(r.itr)
length(r::Reductions{<:NamedTuple{(:init,)}}) = length(r.itr) + 1
size(r::Reductions{<:NamedTuple{()}}) = size(r.itr)
Base.IteratorEltype(::Reductions{<:NamedTuple{()}}) = Base.EltypeUnknown()
Base.IteratorSize(r::Reductions{<:NamedTuple{()}}) = Base.IteratorSize(r.itr)

# TODO: handle multi-dim Arrays along different dims?
reductions(op, A::AbstractArray) = Base.collect_similar(A, Reductions(op, A))

"""
    IsFirstVisit(itr)

Return an iterator which signifies whether each element in `itr` occurs for the first time.

# Example
```julia
julia> s = [1,2,3,1,4,2,5];

julia> is_first_visit(s)
7-element Array{Bool,1}:
  true
  true
  true
 false
  true
 false
  true

julia> Iterators.reverse(IsFirstVisit(s))
Base.Iterators.Reverse{IsFirstVisit{Array{Int64,1}}}(IsFirstVisit{Array{Int64,1}}([1, 2, 3, 1, 4, 2, 5]))

julia> collect(Iterators.reverse(IsFirstVisit(s)))
7-element Array{Any,1}:
  true
 false
  true
 false
  true
  true
  true
```

!!! warning
    Although `IsFirstVisit` supports `Iterators.reverse`, you should still take care of the memory usage.
    Internally we will walk through the `itr` first and calculate the count of each unique element.
"""
struct IsFirstVisit{I}
    itr::I
end

is_first_visit(itr) = IsFirstVisit(itr)

function iterate(itr::IsFirstVisit)
    y = iterate(itr.itr)
    y === nothing && return nothing
    (v, s) = y
    (true, (push!(Set(), v), s))
end

function iterate(itr::IsFirstVisit, (ss, s))
    y = iterate(itr.itr, s)
    y === nothing && return nothing
    (v, s) = y
    (!(v in ss), (push!(ss, v), s))
end

# TODO: handle multi-dim Arrays along different dims?
is_first_visit(A::AbstractArray) = Base.collect_similar(A, IsFirstVisit(A))
length(itr::IsFirstVisit) = length(itr.itr)
size(itr::IsFirstVisit) = size(itr.itr)
eltype(itr::IsFirstVisit) = Bool
Base.IteratorSize(itr::IsFirstVisit) = Base.IteratorSize(itr.itr)

function iterate(itr::Iterators.Reverse{<:IsFirstVisit})
    y = iterate(Iterators.reverse(itr.itr.itr))
    y === nothing && return nothing
    counts = countmap(itr.itr.itr)
    (v, s) = y
    counts[v] -= 1
    counts[v] == 0 && delete!(counts, v)
    (!haskey(counts, v), (counts, s))
end

function iterate(itr::Iterators.Reverse{<:IsFirstVisit}, (counts, s))
    y = iterate(Iterators.reverse(itr.itr.itr), s)
    y === nothing && return nothing
    (v, s) = y
    counts[v] -= 1
    counts[v] == 0 && delete!(counts, v)
    (!haskey(counts, v), (counts, s))
end


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

const EMPTYSTRING = ""

# Implementation of MOI for AbstractModel
abstract type AbstractModelLike{T} <: MOI.ModelLike end
abstract type AbstractOptimizer{T} <: MOI.AbstractOptimizer end
const AbstractModel{T} = Union{AbstractModelLike{T},AbstractOptimizer{T}}

# Variables
function MOI.get(model::AbstractModel, attr::MOI.NumberOfVariables)::Int64
    return MOI.get(model.variables, attr)
end

"""
    function _add_variable end

This is called by `AbstractModel` to inform the `constraints` field that a
variable has been added. This is similar to
[`MathOptInterface.add_variable`](@ref) except that it should return `nothing`.
"""
function _add_variable end

function _add_variable(::Nothing) end
function _add_variables(::Nothing, ::Int64) end

function MOI.add_variable(model::AbstractModel)
    x = MOI.add_variable(model.variables)
    _add_variable(model.constraints)
    return x
end

function MOI.add_variables(model::AbstractModel, n::Integer)
    return [MOI.add_variable(model) for i in 1:n]
end

"""
    remove_variable(f::MOI.AbstractFunction, s::MOI.AbstractSet, vi::MOI.VariableIndex)

Return a tuple `(g, t)` representing the constraint `f`-in-`s` with the
variable `vi` removed. That is, the terms containing the variable `vi` in the
function `f` are removed and the dimension of the set `s` is updated if
needed (e.g. when `f` is a `VectorOfVariables` with `vi` being one of the
variables).
"""
remove_variable(f, s, vi::VI) = remove_variable(f, vi), s
function remove_variable(f::MOI.VectorOfVariables, s, vi::VI)
    g = remove_variable(f, vi)
    if length(g.variables) != length(f.variables)
        t = MOI.update_dimension(s, length(g.variables))
    else
        t = s
    end
    return g, t
end

function filter_variables(keep::F, f, s) where {F<:Function}
    return filter_variables(keep, f), s
end
function filter_variables(
    keep::F,
    f::MOI.VectorOfVariables,
    s,
) where {F<:Function}
    g = filter_variables(keep, f)
    if length(g.variables) != length(f.variables)
        t = MOI.update_dimension(s, length(g.variables))
    else
        t = s
    end
    return g, t
end

function _delete_variable(
    model::AbstractModel{T},
    vi::MOI.VariableIndex,
) where {T}
    MOI.delete(model.variables, vi)
    delete!(model.var_to_name, vi)
    model.name_to_var = nothing
    model.name_to_con = nothing
    return
end

function MOI.delete(model::AbstractModel, vi::MOI.VariableIndex)
    _throw_if_cannot_delete(model.constraints, [vi], vi)
    _delete_variable(model, vi)
    _deleted_constraints(model.constraints, vi) do ci
        return delete!(model.con_to_name, ci)
    end
    MOI.delete(model.objective, vi)
    model.name_to_con = nothing
    return
end

function MOI.delete(model::AbstractModel, vis::Vector{MOI.VariableIndex})
    if isempty(vis)
        return
    end
    _throw_if_cannot_delete(model.constraints, vis, Set(vis))
    _deleted_constraints(model.constraints, vis) do ci
        return delete!(model.con_to_name, ci)
    end
    for vi in vis
        _delete_variable(model, vi)
    end
    MOI.delete(model.objective, vis)
    model.name_to_con = nothing
    return
end

# `AbstractModel` fields `.variables` and `.constraints` act like a
# `StructOfConstraints` where `.variables` contains the `VariableIndex`-in-`S`
# constraints and `.constraints` contains the other constraints.
function constraints(
    model::AbstractModel,
    ci::MOI.ConstraintIndex{MOI.VariableIndex},
)
    return model.variables
end
function constraints(model::AbstractModel, ci::MOI.ConstraintIndex)
    return model.constraints
end

function MOI.is_valid(model::AbstractModel, ci::MOI.ConstraintIndex)
    return MOI.is_valid(constraints(model, ci), ci)
end

function MOI.is_valid(model::AbstractModel, x::MOI.VariableIndex)
    return MOI.is_valid(model.variables, x)
end

function MOI.get(model::AbstractModel, attr::MOI.ListOfVariableIndices)
    return MOI.get(model.variables, attr)
end

# Names
MOI.supports(::AbstractModel, ::MOI.Name) = true
function MOI.set(model::AbstractModel, ::MOI.Name, name::String)
    return model.name = name
end
MOI.get(model::AbstractModel, ::MOI.Name) = model.name

MOI.supports(::AbstractModel, ::MOI.VariableName, vi::Type{VI}) = true
function MOI.set(model::AbstractModel, ::MOI.VariableName, vi::VI, name::String)
    model.var_to_name[vi] = name
    model.name_to_var = nothing # Invalidate the name map.
    return
end

function MOI.get(model::AbstractModel, ::MOI.VariableName, vi::VI)
    return get(model.var_to_name, vi, EMPTYSTRING)
end

"""
    build_name_to_var_map(con_to_name::Dict{MOI.VariableIndex, String})

Create and return a reverse map from name to variable index, given a map from
variable index to name. The special value `MOI.VariableIndex(0)` is used to
indicate that multiple variables have the same name.
"""
function build_name_to_var_map(var_to_name::Dict{VI,String})
    name_to_var = Dict{String,VI}()
    for (var, var_name) in var_to_name
        if haskey(name_to_var, var_name)
            # 0 is a special value that means this string does not map to
            # a unique variable name.
            name_to_var[var_name] = VI(0)
        else
            name_to_var[var_name] = var
        end
    end
    return name_to_var
end

function throw_multiple_name_error(::Type{MOI.VariableIndex}, name::String)
    return error("Multiple variables have the name $name.")
end
function throw_multiple_name_error(::Type{<:MOI.ConstraintIndex}, name::String)
    return error("Multiple constraints have the name $name.")
end
function throw_if_multiple_with_name(::Nothing, ::String) end
function throw_if_multiple_with_name(index::MOI.Index, name::String)
    if iszero(index.value)
        throw_multiple_name_error(typeof(index), name)
    end
end

function MOI.get(model::AbstractModel, ::Type{VI}, name::String)
    if model.name_to_var === nothing
        # Rebuild the map.
        model.name_to_var = build_name_to_var_map(model.var_to_name)
    end
    result = get(model.name_to_var, name, nothing)
    throw_if_multiple_with_name(result, name)
    return result
end

function MOI.get(
    model::AbstractModel,
    ::MOI.ListOfVariableAttributesSet,
)::Vector{MOI.AbstractVariableAttribute}
    return isempty(model.var_to_name) ? [] : [MOI.VariableName()]
end

MOI.supports(model::AbstractModel, ::MOI.ConstraintName, ::Type{<:CI}) = true

function MOI.set(
    model::AbstractModel,
    ::MOI.ConstraintName,
    ci::CI,
    name::String,
)
    model.con_to_name[ci] = name
    model.name_to_con = nothing # Invalidate the name map.
    return
end

function MOI.supports(
    ::AbstractModel,
    ::MOI.ConstraintName,
    ::Type{<:MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.AbstractScalarSet}},
)
    return throw(MOI.VariableIndexConstraintNameError())
end

function MOI.set(
    ::AbstractModel,
    ::MOI.ConstraintName,
    ::MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.AbstractScalarSet},
    ::String,
)
    return throw(MOI.VariableIndexConstraintNameError())
end

function MOI.get(model::AbstractModel, ::MOI.ConstraintName, ci::CI)
    return get(model.con_to_name, ci, EMPTYSTRING)
end

"""
    build_name_to_con_map(con_to_name::Dict{MOI.ConstraintIndex, String})

Create and return a reverse map from name to constraint index, given a map from
constraint index to name. The special value
`MOI.ConstraintIndex{Nothing, Nothing}(0)` is used to indicate that multiple
constraints have the same name.
"""
function build_name_to_con_map(con_to_name::Dict{CI,String})
    name_to_con = Dict{String,CI}()
    for (con, con_name) in con_to_name
        if haskey(name_to_con, con_name)
            name_to_con[con_name] = CI{Nothing,Nothing}(0)
        else
            name_to_con[con_name] = con
        end
    end
    return name_to_con
end

function MOI.get(model::AbstractModel, ConType::Type{<:CI}, name::String)
    if model.name_to_con === nothing
        # Rebuild the map.
        model.name_to_con = build_name_to_con_map(model.con_to_name)
    end
    ci = get(model.name_to_con, name, nothing)
    throw_if_multiple_with_name(ci, name)
    return ci isa ConType ? ci : nothing
end

function MOI.get(
    model::AbstractModel,
    ::MOI.ListOfConstraintAttributesSet{F,S},
) where {F,S}
    if any(k -> k isa MOI.ConstraintIndex{F,S}, keys(model.con_to_name))
        return MOI.AbstractConstraintAttribute[MOI.ConstraintName()]
    end
    return MOI.AbstractConstraintAttribute[]
end

# Objective

function MOI.get(
    model::AbstractModel,
    attr::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction,
        MOI.ObjectiveFunctionType,
    },
)
    return MOI.get(model.objective, attr)
end

function MOI.supports(
    model::AbstractModel,
    attr::Union{MOI.ObjectiveSense,MOI.ObjectiveFunction},
)
    return MOI.supports(model.objective, attr)
end

function MOI.set(
    model::AbstractModel,
    attr::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    MOI.set(model.objective, attr, sense)
    return
end

function MOI.set(
    model::AbstractModel,
    attr::MOI.ObjectiveFunction{F},
    f::F,
) where {F<:MOI.AbstractFunction}
    if !MOI.supports(model, attr)
        throw(MOI.UnsupportedAttribute(attr))
    end
    MOI.set(model.objective, attr, f)
    return
end

function MOI.modify(
    model::AbstractModel,
    attr::MOI.ObjectiveFunction,
    change::MOI.AbstractFunctionModification,
)
    MOI.modify(model.objective, attr, change)
    return
end

function MOI.get(::AbstractModel, ::MOI.ListOfOptimizerAttributesSet)
    return MOI.AbstractOptimizerAttribute[]
end

function MOI.get(
    model::AbstractModel,
    attr::MOI.ListOfModelAttributesSet,
)::Vector{MOI.AbstractModelAttribute}
    ret = MOI.get(model.objective, attr)
    if !isempty(model.name)
        push!(ret, MOI.Name())
    end
    return ret
end

# Constraints

function MOI.supports_constraint(
    model::AbstractModel,
    ::Type{MOI.VariableIndex},
    ::Type{S},
) where {S<:MOI.AbstractScalarSet}
    return MOI.supports_constraint(model.variables, MOI.VariableIndex, S)
end
function MOI.supports_constraint(
    model::AbstractModel,
    ::Type{F},
    ::Type{S},
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    return MOI.supports_constraint(model.constraints, F, S)
end

function MOI.add_constraint(
    model::AbstractModel,
    f::MOI.VariableIndex,
    s::MOI.AbstractScalarSet,
)
    return MOI.add_constraint(model.variables, f, s)
end

function MOI.add_constraint(
    model::AbstractModel,
    func::MOI.AbstractFunction,
    set::MOI.AbstractSet,
)
    return MOI.add_constraint(model.constraints, func, set)
end

function MOI.get(
    model::AbstractModel,
    attr::Union{MOI.AbstractFunction,MOI.AbstractSet},
    ci::MOI.ConstraintIndex,
)
    return MOI.get(model.constraints, attr, ci)
end

function MOI.delete(model::AbstractModel, ci::MOI.ConstraintIndex)
    MOI.delete(constraints(model, ci), ci)
    model.name_to_con = nothing
    delete!(model.con_to_name, ci)
    return
end

function MOI.modify(
    model::AbstractModel,
    ci::MOI.ConstraintIndex,
    change::MOI.AbstractFunctionModification,
)
    MOI.modify(model.constraints, ci, change)
    return
end

function MOI.set(
    model::AbstractModel,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
    value,
)
    MOI.set(constraints(model, ci), attr, ci, value)
    return
end

function MOI.get(
    model::AbstractModel,
    attr::MOI.NumberOfConstraints{MOI.VariableIndex,S},
) where {S}
    return MOI.get(model.variables, attr)
end

function MOI.get(
    model::AbstractModel,
    noc::MOI.NumberOfConstraints{F,S},
)::Int64 where {F,S}
    return MOI.get(model.constraints, noc)
end

function MOI.get(
    model::AbstractModel{T},
    attr::MOI.ListOfConstraintTypesPresent,
) where {T}
    return vcat(
        MOI.get(model.constraints, attr)::Vector{Tuple{Type,Type}},
        MOI.get(model.variables, attr)::Vector{Tuple{Type,Type}},
    )
end

function MOI.get(
    model::AbstractModel,
    attr::MOI.ListOfConstraintIndices{MOI.VariableIndex,S},
) where {S}
    return MOI.get(model.variables, attr)
end

function MOI.get(
    model::AbstractModel,
    loc::MOI.ListOfConstraintIndices{F,S},
) where {F,S}
    return MOI.get(model.constraints, loc)
end

function MOI.get(
    model::AbstractModel,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
)
    return MOI.get(constraints(model, ci), attr, ci)
end

function MOI.is_empty(model::AbstractModel)
    return isempty(model.name) &&
           MOI.is_empty(model.objective) &&
           MOI.is_empty(model.constraints) &&
           MOI.is_empty(model.variables)
end

function MOI.empty!(model::AbstractModel{T}) where {T}
    model.name = ""
    MOI.empty!(model.objective)
    MOI.empty!(model.variables)
    empty!(model.var_to_name)
    model.name_to_var = nothing
    empty!(model.con_to_name)
    model.name_to_con = nothing
    MOI.empty!(model.constraints)
    return
end

function pass_nonvariable_constraints(
    dest::AbstractModel,
    src::MOI.ModelLike,
    idxmap::IndexMap,
    constraint_types,
)
    return pass_nonvariable_constraints(
        dest.constraints,
        src,
        idxmap,
        constraint_types,
    )
end

function MOI.copy_to(dest::AbstractModel, src::MOI.ModelLike)
    return default_copy_to(dest, src)
end

MOI.supports_incremental_interface(::AbstractModel) = true
function final_touch(model::AbstractModel, index_map)
    return final_touch(model.constraints, index_map)
end

# Macro to generate Model

function _struct_of_constraints_type(name, subtypes, parametrized_type)
    if length(subtypes) == 1
        # Only one type, no need for a `StructOfConstraints`.
        return subtypes[1]
    else
        T = esc(:T)
        t = :($name{$T})
        if parametrized_type
            append!(t.args, subtypes)
        end
        return t
    end
end

# This macro is for expert/internal use only. Prefer the concrete Model type
# instantiated below.
"""
    macro model(
        model_name,
        scalar_sets,
        typed_scalar_sets,
        vector_sets,
        typed_vector_sets,
        scalar_functions,
        typed_scalar_functions,
        vector_functions,
        typed_vector_functions,
        is_optimizer = false
    )

Creates a type `model_name` implementing the MOI model interface and containing
`scalar_sets` scalar sets `typed_scalar_sets` typed scalar sets, `vector_sets`
vector sets, `typed_vector_sets` typed vector sets, `scalar_functions` scalar
functions, `typed_scalar_functions` typed scalar functions, `vector_functions`
vector functions and `typed_vector_functions` typed vector functions.
To give no set/function, write `()`, to give one set `S`, write `(S,)`.

The function [`MathOptInterface.VariableIndex`](@ref) should not be given in
`scalar_functions`. The model supports [`MathOptInterface.VariableIndex`](@ref)-in-`S`
constraints where `S` is [`MathOptInterface.EqualTo`](@ref),
[`MathOptInterface.GreaterThan`](@ref), [`MathOptInterface.LessThan`](@ref),
[`MathOptInterface.Interval`](@ref), [`MathOptInterface.Integer`](@ref),
[`MathOptInterface.ZeroOne`](@ref), [`MathOptInterface.Semicontinuous`](@ref)
or [`MathOptInterface.Semiinteger`](@ref). The sets supported
with the [`MathOptInterface.VariableIndex`](@ref) cannot be controlled from the
macro, use the [`UniversalFallback`](@ref) to support more sets.

This macro creates a model specialized for specific types of constraint,
by defining specialized structures and methods. To create a model that,
in addition to be optimized for specific constraints, also support arbitrary
constraints and attributes, use [`UniversalFallback`](@ref).

If `is_optimizer = true`, the resulting struct is a
of [`GenericOptimizer`](@ref), which is a subtype of
[`MathOptInterface.AbstractOptimizer`](@ref), otherwise, it is a
[`GenericModel`](@ref), which is a subtype of
[`MathOptInterface.ModelLike`](@ref).

### Examples

The model describing an linear program would be:
```julia
@model(LPModel,                                                   # Name of model
      (),                                                         # untyped scalar sets
      (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval), #   typed scalar sets
      (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives),            # untyped vector sets
      (),                                                         #   typed vector sets
      (),                                                         # untyped scalar functions
      (MOI.ScalarAffineFunction,),                                #   typed scalar functions
      (MOI.VectorOfVariables,),                                   # untyped vector functions
      (MOI.VectorAffineFunction,),                                #   typed vector functions
      false
    )
```

Let `MOI` denote `MathOptInterface`, `MOIU` denote `MOI.Utilities`.
The macro would create the following types with
[`struct_of_constraint_code`](@ref):
```julia
struct LPModelScalarConstraints{T, C1, C2, C3, C4} <: MOIU.StructOfConstraints
    moi_equalto::C1
    moi_greaterthan::C2
    moi_lessthan::C3
    moi_interval::C4
end
struct LPModelVectorConstraints{T, C1, C2, C3} <: MOIU.StructOfConstraints
    moi_zeros::C1
    moi_nonnegatives::C2
    moi_nonpositives::C3
end
struct LPModelFunctionConstraints{T} <: MOIU.StructOfConstraints
    moi_scalaraffinefunction::LPModelScalarConstraints{
        T,
        MOIU.VectorOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}},
        MOIU.VectorOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}},
        MOIU.VectorOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}},
        MOIU.VectorOfConstraints{MOI.ScalarAffineFunction{T}, MOI.Interval{T}}
    }
    moi_vectorofvariables::LPModelVectorConstraints{
        T,
        MOIU.VectorOfConstraints{MOI.VectorOfVariables, MOI.Zeros},
        MOIU.VectorOfConstraints{MOI.VectorOfVariables, MOI.Nonnegatives},
        MOIU.VectorOfConstraints{MOI.VectorOfVariables, MOI.Nonpositives}
    }
    moi_vectoraffinefunction::LPModelVectorConstraints{
        T,
        MOIU.VectorOfConstraints{MOI.VectorAffineFunction{T}, MOI.Zeros},
        MOIU.VectorOfConstraints{MOI.VectorAffineFunction{T}, MOI.Nonnegatives},
        MOIU.VectorOfConstraints{MOI.VectorAffineFunction{T}, MOI.Nonpositives}
    }
end
const LPModel{T} = MOIU.GenericModel{T,MOIU.ObjectiveContainer{T},MOIU.VariablesContainer{T},LPModelFunctionConstraints{T}}
```
The type `LPModel` implements the MathOptInterface API except methods specific
to optimizers like `optimize!` or `get` with `VariablePrimal`.
"""
macro model(
    model_name,
    ss,
    sst,
    vs,
    vst,
    sf,
    sft,
    vf,
    vft,
    is_optimizer = false,
)
    scalar_sets = [SymbolSet.(ss.args, false); SymbolSet.(sst.args, true)]
    vector_sets = [SymbolSet.(vs.args, false); SymbolSet.(vst.args, true)]

    scname = esc(Symbol(string(model_name) * "ScalarConstraints"))
    vcname = esc(Symbol(string(model_name) * "VectorConstraints"))

    esc_model_name = esc(model_name)
    # TODO if there is only one function or one set, remove the layer

    scalar_funs = [
        SymbolFun.(sf.args, false)
        SymbolFun.(sft.args, true)
    ]
    vector_funs = [
        SymbolFun.(vf.args, false)
        SymbolFun.(vft.args, true)
    ]
    funs = [scalar_funs; vector_funs]
    set_struct_types = map(eachindex(funs)) do i
        if i <= length(scalar_funs)
            cname = scname
            sets = scalar_sets
        else
            cname = vcname
            sets = vector_sets
        end
        voc = map(sets) do set
            return :(VectorOfConstraints{$(_typed(funs[i])),$(_typed(set))})
        end
        return _struct_of_constraints_type(cname, voc, true)
    end
    func_name = esc(Symbol(string(model_name) * "FunctionConstraints"))
    func_typed = _struct_of_constraints_type(func_name, set_struct_types, false)
    T = esc(:T)
    generic = if is_optimizer
        :(GenericOptimizer{
            $T,
            ObjectiveContainer{$T},
            VariablesContainer{$T},
            $func_typed,
        })
    else
        :(GenericModel{
            $T,
            ObjectiveContainer{$T},
            VariablesContainer{$T},
            $func_typed,
        })
    end
    model_code = :(const $esc_model_name{$T} = $generic)
    expr = Expr(:block)
    if length(scalar_sets) >= 2
        push!(expr.args, struct_of_constraint_code(scname, scalar_sets))
    end
    if length(vector_sets) >= 2
        push!(expr.args, struct_of_constraint_code(vcname, vector_sets))
    end
    if length(funs) != 1
        push!(
            expr.args,
            struct_of_constraint_code(func_name, funs, set_struct_types),
        )
    end
    push!(expr.args, model_code)
    return expr
end

for (loop_name, loop_super_type) in [
    (:GenericModel, :AbstractModelLike),
    (:GenericOptimizer, :AbstractOptimizer),
]
    global name = loop_name
    global super_type = loop_super_type
    @eval begin
        """
            mutable struct $name{T,O,V,C} <: $super_type{T}

        Implements a model supporting coefficients of type `T` and:

         * An objective function stored in `.objective::O`
         * Variables and `VariableIndex` constraints stored in `.variable_bounds::V`
         * `F`-in-`S` constraints (excluding `VariableIndex` constraints)
           stored in `.constraints::C`

        All interactions should take place via the MOI interface, so the types
        `O`, `V`, and `C` should implement the API as needed for their
        functionality.
        """
        mutable struct $name{T,O,V,C} <: $super_type{T}
            name::String
            objective::O
            variables::V
            constraints::C
            var_to_name::Dict{MOI.VariableIndex,String}
            # If `nothing`, the dictionary hasn't been constructed yet.
            name_to_var::Union{Dict{String,MOI.VariableIndex},Nothing}
            con_to_name::Dict{MOI.ConstraintIndex,String}
            name_to_con::Union{Dict{String,MOI.ConstraintIndex},Nothing}
            # A useful dictionary for extensions to store things. These are
            # _not_ copied between models!
            ext::Dict{Symbol,Any}
            function $name{T,O,V,C}() where {T,O,V,C}
                return new{T,O,V,C}(
                    EMPTYSTRING,
                    O(),
                    V(),
                    C(),
                    Dict{MOI.VariableIndex,String}(),
                    nothing,
                    Dict{MOI.ConstraintIndex,String}(),
                    nothing,
                    Dict{Symbol,Any}(),
                )
            end
        end
    end
end

const LessThanIndicatorOne{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.LessThan{T}}
const LessThanIndicatorZero{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ZERO,MOI.LessThan{T}}

@model(
    Model,
    (MOI.ZeroOne, MOI.Integer),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        MOI.Semicontinuous,
        MOI.Semiinteger,
    ),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.Complements,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.PositiveSemidefiniteConeSquare,
        MOI.RootDetConeTriangle,
        MOI.RootDetConeSquare,
        MOI.LogDetConeTriangle,
        MOI.LogDetConeSquare,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        MOI.SOS1,
        MOI.SOS2,
        LessThanIndicatorOne,
        LessThanIndicatorZero,
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

@doc raw"""

An implementation of `ModelLike` that supports all functions and sets defined
in MOI. It is parameterized by the coefficient type.

# Examples
```jl
model = Model{Float64}()
x = add_variable(model)
```
"""
Model

# This export makes the type be printed as:
# ```julia
# julia> Base.show(IOContext(stdout, :compact => true), MathOptInterface.Utilities.Model)
# Model{T} where T
# julia> print(MathOptInterface.Utilities.Model)
# MathOptInterface.Utilities.Model{T} where T
# julia> MathOptInterface.Utilities.Model
# MathOptInterface.Utilities.Model{T} where T (alias for MathOptInterface.Utilities.GenericModel{T, MathOptInterface.Utilities.ObjectiveContainer{T}, MathOptInterface.Utilities.VariablesContainer{T}, MathOptInterface.Utilities.ModelFunctionConstraints{T}} where T)
# ```
# As MOI is not doing `using .Utilities` and is not exporting `Model`, the user
# still needs to do `MOI.Utilities.Model` unless he does
# `using MathOptInterface.Utilities`.
export Model

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
include("Utils.jl")
install("EntityComponentSystem")
install("Gtk")

using EntityComponentSystem, Gtk, Graphics

world = World()

includedir("components")
includedir("systems")
includedir("entities")

canvas = GtkCanvas()
win = GtkWindow(canvas, "Pong")
@guarded draw(canvas) do widget
    ctx = getgc(canvas)
    h = height(canvas)
    w = width(canvas)
    # Paint red rectangle
    rectangle(ctx, 0, 0, w, h/2)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)
    # Paint blue rectangle
    rectangle(ctx, 0, 3h/4, w, h/4)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)
end
show(canvas)

exit_condition = Condition()
signal_connect(win, :destroy) do widget
	notify(exit_condition)
end
wait(exit_condition)


text:
module SaddlePointSystems

using LinearAlgebra
using LinearMaps
using IterativeSolvers
#using ..Fields

import Base: *, \
import LinearAlgebra: ldiv!

export SaddleSystem

"""
    SaddleSystem((u,f),(A⁻¹,B₁ᵀ,B₂);[tol=1e-3],[issymmetric=false],[isposdef=false],[conditioner=Identity],[store=false])

Construct the computational operators for a saddle-point system of the form
\$[A B₁ᵀ; B₂ 0][u;f]\$. Note that the constituent operators are passed in as a
tuple in the order seen here. Each of these operators could act on its corresponding
data type in a function-like way, e.g. `A⁻¹(u)`, or in a matrix-like way, e.g.,
`A⁻¹*u`.

The optional argument `tol` sets the tolerance for iterative solution (if
  applicable). Its default is 1e-3.

The optional argument `conditioner` can be used to supply a function that acts
upon the result `f` to 'condition' it (e.g. filter it). It is, by default, set
to the identity.

The optional Boolean argument `store` will compute and store the Schur
complement matrix's factorization. This makes the inversion faster, though
it comes at the expense of memory and overhead time for pre-computing it. The
resulting solution is somewhat noiser, too.

# Arguments

- `u` : example of state vector data.
- `f` : example of constraint force vector data. This data must be of
        AbstractVector supertype.
- `A⁻¹` : operator evaluating the inverse of `A` on data of type `u`, return type `u`
- `B₁ᵀ` : operator evaluating the influence of constraint force,
            acting on `f` and returning type `u`
- `B₂` : operator evaluating the influence of state vector on constraints,
            acting on `u` and returning type `f`
"""
struct SaddleSystem{FA,FAB,FBA,FP,TU,TF,T,N,Storage}
    A⁻¹B₁ᵀf :: TU
    B₂A⁻¹r₁ :: TF
    tmpvec :: Vector{T}
    tmpvecout :: Vector{T}
    A⁻¹ :: FA
    A⁻¹B₁ᵀ :: FAB
    B₂A⁻¹ :: FBA
    P :: FP
    S  :: LinearMap
    S⁻¹ :: Union{Factorization{T},Nothing}
    tol :: Float64
    _issymmetric :: Bool
    _isposdef :: Bool
end

function (::Type{SaddleSystem})(state::Tuple{TU,TF},sys::Tuple{FA,FB1,FB2};
                                tol::Float64=1e-3,
                                conditioner::FP=x->x,
                                issymmetric::Bool=false,
                                isposdef::Bool=false,
                                store::Bool=false,
                                precompile::Bool=true) where {TU,TF,FA,FB1,FB2,FP}
    u,f = state

    optypes = (TU,TF,TU)
    opnames = ("A⁻¹","B₁ᵀ","B₂")
    ops = []
    T = eltype(TU)

    # check for methods
    for (i,typ) in enumerate(optypes)
      if hasmethod(sys[i],Tuple{typ})
        push!(ops,sys[i])
      elseif hasmethod(*,Tuple{typeof(sys[i]),typ})
        # generate a method that acts on TU
        push!(ops,x->sys[i]*x)
      else
        error("No valid operator for $(opnames[i]) supplied")
      end
    end

    A⁻¹, B₁ᵀ, B₂ = ops
    ubuffer = deepcopy(u)
    fbuffer = deepcopy(f)
    N = length(f)
    tmpvec = zeros(T,N)
    tmpvecout = zeros(T,N)

    # Schur complement
    function Schur!(fout::AbstractVector{T},fin::AbstractVector{T}) where {T<:Number}
       fbuffer .= fin
       fbuffer .= (B₂∘A⁻¹∘B₁ᵀ)(fbuffer)
       fout .= -fbuffer
       return fout
    end
    S = LinearMap{T}(Schur!,N;ismutating=true,issymmetric=issymmetric,isposdef=isposdef)

    if store && N > 0
      Smat = zeros(T,N,N)
      fill!(f,0.0)
      for i = 1:N
        f[i] = 1.0
        fbuffer .= S*f
        Smat[1:N,i] .= fbuffer
        f[i] = 0.0
      end
      #S⁻¹ = Nullable(factorize(Smat))
      S⁻¹ = factorize(Smat)
    else
      #S⁻¹ = Nullable()
      S⁻¹ = nothing
    end

    A⁻¹B₁ᵀ(f::TF) = (A⁻¹∘B₁ᵀ)(f)

    B₂A⁻¹(w::TU) = (B₂∘A⁻¹)(w)

    saddlesys = SaddleSystem{typeof(A⁻¹),typeof(A⁻¹B₁ᵀ),typeof(B₂A⁻¹),typeof(conditioner),TU,TF,T,N,store}(
                                ubuffer,fbuffer,tmpvec,tmpvecout,
                                A⁻¹,A⁻¹B₁ᵀ,B₂A⁻¹,conditioner,S,S⁻¹,
                                tol,issymmetric,isposdef)
    # run once in order to precompile it
    if precompile
      saddlesys\(u,f)
    end

    return saddlesys

end

(::Type{SaddleSystem})(state::Tuple{TU,TF},sys::Array{Any,2};kwarg...) where {TU,TF} =
            SaddleSystem(state,(sys[1,1],sys[1,2],sys[2,1]);kwarg...)

function Base.show(io::IO, S::SaddleSystem{FA,FAB,FBA,FP,TU,TF,T,N,Storage}) where {FA,FAB,FBA,FP,TU,TF,T,N,Storage}
    println(io, "Saddle system with $N constraints and")
    println(io, "   State of type $TU")
    println(io, "   Force of type $TF")
end

"""
    ldiv!(state,sys::SaddleSystem,rhs)

Solve a saddle-point system. `rhs` is a tuple of the right-hand side `(ru,rf)`.
Output `state`, a tuple (u,f), is updated. Note that `sys` is also mutated:
its scratch space `sys.B₂A⁻¹r₁` and `sys.A⁻¹B₁ᵀf` hold the intermediate results
of the solution.

A shorthand can be used for this operation: `state = sys\\rhs`
"""
function ldiv!(state::Tuple{TU,TF},
                    sys::SaddleSystem{FA,FAB,FBA,FP,TU,TF,T,N,false},
                    rhs::Tuple{TU,TF}) where {TU,TF,T,FA,FAB,FBA,FP,N}
  # non-stored matrix

  ru, rf = rhs
  u, f = state
  sys.B₂A⁻¹r₁ .= sys.B₂A⁻¹(ru)
  sys.tmpvec .= rf
  sys.tmpvec .-= sys.B₂A⁻¹r₁
  if N > 0
    sys.tmpvecout .= f
    cg!(sys.tmpvecout,sys.S,sys.tmpvec,tol=sys.tol)
    f .= sys.tmpvecout
    f .= sys.P(f)
  end
  u .= sys.A⁻¹(ru)
  sys.A⁻¹B₁ᵀf .= sys.A⁻¹B₁ᵀ(f)
  u .-= sys.A⁻¹B₁ᵀf
  state = u, f
end

# stored matrix
function ldiv!(state::Tuple{TU,TF},
                    sys::SaddleSystem{FA,FAB,FBA,FP,TU,TF,T,N,true},
                    rhs::Tuple{TU,TF}) where {TU,TF,T,FA,FAB,FBA,FP,N}

  ru, rf = rhs
  u, f = state
  sys.B₂A⁻¹r₁ .= sys.B₂A⁻¹(ru)
  sys.tmpvec .= rf
  sys.tmpvec .-= sys.B₂A⁻¹r₁
  if N > 0
    ldiv!(sys.S⁻¹,sys.tmpvec)
    f .= sys.tmpvec
    f .= sys.P(f)
  end
  u .= sys.A⁻¹(ru)
  sys.A⁻¹B₁ᵀf .= sys.A⁻¹B₁ᵀ(f)
  u .-= sys.A⁻¹B₁ᵀf
  state = u, f
end

\(sys::SaddleSystem{FA,FAB,FBA,FP,TU,TF,T,N,Storage},rhs::Tuple{TU,TF}) where {TU,TF,T,FA,FAB,FBA,FP,N,Storage} =
      ldiv!(similar.(rhs),sys,rhs)

# solving tuples of systems
function ldiv!(state,sys::NTuple{M,SaddleSystem},rhs) where {M}
  for (i,sysi) in enumerate(sys)
    ldiv!(state[i],sysi,rhs[i])
  end
  state
end

\(sys::NTuple{M,SaddleSystem},rhs) where {M} =
      ldiv!(deepcopy.(rhs),sys,rhs)

end


text:
using Elefridge
using ZfpCompression
using PyPlot
using PyCall
using Blosc
xr = pyimport("xarray")
ccrs = pyimport("cartopy.crs")

path = "/Users/milan/cams/gridded/"
filelist = filter(x->endswith(x,"_go3.grib"),readdir(path))
grib = xr.open_dataarray(joinpath(path,filelist[end]),engine="cfgrib")
X = grib.data
lat = grib.latitude.data
lon = grib.longitude.data

for level in [5:5:137]
    o3 = X[level,:,:]
    o3 = copy(o3')

    ## compression
    rbits = [10,7,5,3,1]
    cfs = fill(0.0,length(rbits))     # compression factors
    Blosc.set_compressor("lz4hc")

    for (i,r) in enumerate(rbits)
        o3r = round(o3,r)
        o3c = compress(o3r)
        cfs[i] = sizeof(o3)/sizeof(o3c)
    end

    lat_div(n::Integer) = Array(-90:180/(n-1):90)
    lon_div(n::Integer) = Array(0:360/n:360)

    ## PLOT
    ioff()
    vmin,vmax = extrema(o3)
    fig,ax1 = plt.subplots(1,1,figsize=(10,6),sharex=true,sharey=true,
                        subplot_kw=Dict("projection"=>ccrs.Robinson(central_longitude=180)))

    lond = lon_div(length(rbits)+1)

    # Float32
    ax1.pcolormesh(lon[lon .<= lond[2]],lat,o3[lon .<= lond[2],:]',vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree())

    for i in 3:length(lond)
        lonwhere = (lon .<= lond[i]) .& (lon .> lond[i-1]-0.5)
        ax1.pcolormesh(lon[lonwhere],lat,round(o3,rbits[i-2])[lonwhere,:]',vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree())
    end

    ax1.coastlines(color="k",lw=0.5)

    for longi in lond[2:end-1]
        ax1.plot(repeat([longi,],31),lat_div(31),"w",transform=ccrs.PlateCarree(),lw=0.25)
        # ax1.plot(repeat([longi,],6),lat_div(31)[end-5:end],"w",transform=ccrs.PlateCarree(),lw=1.5)
    end

    ax1.text(0.127,0.105,"Compression\nfactor",transform=ax1.transAxes,fontweight="bold",color="white")

    for (i,x) in enumerate([0.23,0.33,0.43,0.54,0.64,0.76])
        c = i == 1 ? 1 : Int(round(cfs[i-1]))
        ax1.text(x,0.03,"$(c)x",transform=ax1.transAxes,fontweight="bold",fontsize=14,
                                color="white")

    end

    ax1.set_title(L"O$_3$ compression: Round + LZ4HC lossless",fontweight="bold")

    tight_layout()
    savefig("/Users/milan/git/Elefridge.jl/maps/o3/round_o3_$level.png",dpi=200)
    close(fig)
end


text:
using JUtils
using Test

@testset "JUtils.jl" begin
    # Write your tests here.
end


text:
__precompile__()
module Giac

using Libdl
using SpecialFunctions

import Base: string, show, write, collect
import Base: diff, sum, zeros, length, size, getindex
import Base: +, -, (*), /, ^, ==, >, <, >=, <=
import Base: convert, one, oneunit, zero, inv
import Base: round, floor, ceil, trunc
import Base: real, imag, conj, adjoint, abs, sign
import Base: sqrt, exp, log, log10, sin, cos, tan, sec, csc, cot
import Base: sinh, cosh, tanh, asin, acos, atan, acot, asec, acsc
import Base: asinh, acosh, atanh, expm1, log1p
import Base: factorial, binomial, numerator, denominator
import SpecialFunctions: gamma, beta, zeta, besselj, bessely, erfc, erf
import SpecialFunctions: airyai, airybi, sinint, cosint, polygamma, digamma

export @giac, @vars, giac, giac_identifier, type, subtype, isgiacnumber, isgiacreal
export evaluate, evaluatef, evalf, simplify, expand, to_julia, set!, purge!, giac_vars
export unapply, latex, prettyprint, head, args, lambdify
export infinity, plus_inf, minus_inf

export partfrac, subst, left, right
export ⩦, equal
export Ei, Psi
export integrate, limit, series, curl, grad, divergence, hessian
export preval, sum_riemann, taylor
export solve, cSolve, cZeros, fSolve, deSolve, linsolve
export logcollect, powexpand,texpand, exp2pow, pow2exp, exp2trig, expexpand
export asin2acos, asin2atan, acos2asin, acos2atan
export sin2costan, cos2sintan
export atan2asin, atan2acos, tan2sincos, halftan
export trigsin, trigcos, trigtan, atrig2log, tlin, tcollect, trigexpand
export trig2exp
export gbasis, greduce, factor

mutable struct giac
    g::Ptr{Nothing}
    function giac(g::Ptr{Nothing})
        gg = new(g)
        finalizer(_delete, gg)
        gg
    end
end

const libgiac_c = Ref{Ptr{Cvoid}}(0)
const context_ptr = Ref{Ptr{Cvoid}}()
const giac_undef = Ref{giac}()
const giac_one = Ref{giac}()
const giac_zero = Ref{giac}()
const giac_pi = Ref{giac}()
const giac_e = Ref{giac}()

function __init__()
    libgiac_c[] = dlopen(joinpath(dirname(@__FILE__), "..", "deps", "lib",
                         string("libgiac_c.", dlext)))
    context_ptr[] = ccall(dlsym(libgiac_c[], "giac_context_ptr"), Ptr{Nothing}, () )
    giac_undef[] = giac("undef")
    global infinity = giac("infinity")
    global plus_inf = giac("plus_inf")
    global minus_inf = giac("minus_inf")
    giac_e[] = giac("e")
    giac_pi[] = giac("pi")
    giac_one[] = giac(1)
    giac_zero[] = giac(0)
end

const config_vars = (:Digits, :epsilon, :approx_mode, :complex_mode, :complex_variables)

const giac_julia_dictionay = Dict(
    :ln      => :log,
    :Gamma   => :gamma,
    :Beta    => :beta,
    :BesselJ => :besselj,
    :BesselY => :bessely,
    :Si      => :sinint,
    :Ci      => :cosint,
    :re      => :real,
    :im      => :imag,
    )

function type(g::giac)
    t = unsafe_load(Ptr{UInt8}(g.g), 1) & 31
    if t<0 || t>21
        return :UNKNOWN
    end
    [
    # from giac/dispatch.h:
    :INT,        #  0 # int val
    :DOUBLE,     #  1 # double _DOUBLE_val
    :ZINT,       #  2 # mpz_t * _ZINTptr
    :REAL,       #  3 # mpf_t * _REALptr
    :CPLX,       #  4 # gen * _CPLXptr
    :POLY,       #  5 # polynome * _POLYptr
    :IDNT,       #  6 # identificateur * _IDNTptr
    :VECT,       #  7 # vecteur * _VECTptr
    :SYMB,       #  8 # symbolic * _SYMBptr
    :SPOL1,      #  9 # sparse_poly1 * _SPOL1ptr
    :FRAC,       # 10 # fraction * _FRACptr
    :EXT,        # 11 # gen * _EXTptr
    :STRNG,      # 12 # string * _STRNGptr
    :FUNC,       # 13 # unary_fonction_ptr * _FUNCptr
    :ROOT,       # 14 # real_complex_rootof *_ROOTptr
    :MOD,        # 15 # gen * _MODptr
    :USER,       # 16 # giac_user * _USERptr
    :MAP,        # 17 # map<gen.gen> * _MAPptr
    :EQW,        # 18 # eqwdata * _EQWptr
    :GROB,       # 19 # grob * _GROBptr
    :POINTER,    # 20 # void * _POINTER_val
    :FLOAT,      # 21 # immediate, _FLOAT_val
    ][t+1]
end

isgiacreal(g::giac) = type(g) in [:INT, :DOUBLE, :ZINT, :REAL, :FLOAT]
isgiacnumber(g::giac)  = type(g) in [:INT, :DOUBLE, :ZINT, :REAL, :FLOAT, :CPLX]

subtype(g::giac) = unsafe_load(Ptr{UInt8}(g.g), 2) & 31

function change_subtype(g::giac, subtype::Integer)
    ccall( dlsym(libgiac_c[], "giac_change_subtype"), Nothing, (Ptr{Nothing},Cint), g.g, subtype)
end

function _delete(g::giac)
    ccall( dlsym(libgiac_c[], "giac_delete"), Nothing, (Ptr{Nothing},), g.g)
end

giac(x) = giac_undef[]

function giac(g::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_copy_gen"), Ptr{Nothing}, (Ptr{Nothing},), g.g))
end

function giac(val::Bool)
    g = giac(ccall(dlsym(libgiac_c[], "giac_new_int"), Ptr{Nothing}, (Cint,), val ? 1 : 0))
    change_subtype(g, 6)
    g
end

function giac(val::Cint)
    giac(ccall(dlsym(libgiac_c[], "giac_new_int"), Ptr{Nothing}, (Cint,), val))
end

function giac(val::Int64)
    giac(ccall(dlsym(libgiac_c[], "giac_new_int64_t"), Ptr{Nothing}, (Int64,), val))
end

function giac(val::BigInt)
    giac(ccall(dlsym(libgiac_c[], "giac_new_bigint"), Ptr{Nothing}, (Ref{BigInt},), val))
end

function giac(val::Rational)
    giac(ccall(dlsym(libgiac_c[], "giac_new_rational"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), giac(val.num).g, giac(val.den).g))
end

function giac(val::Int128)
    giac(ccall(dlsym(libgiac_c[], "giac_new_int128_t"), Ptr{Nothing}, (Int128,), val))
end

function giac(val::Cdouble)
    giac(ccall(dlsym(libgiac_c[], "giac_new_double"), Ptr{Nothing}, (Cdouble,), val))
end

function giac(val::BigFloat)
    giac(ccall(dlsym(libgiac_c[], "giac_new_bigfloat"), Ptr{Nothing}, (Ref{BigFloat},), val))
end

function giac(val::Complex)
    giac(ccall(dlsym(libgiac_c[], "giac_new_complex"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), giac(real(val)).g, giac(imag(val)).g))
end

function giac(val::Complex{Cint})
    giac(ccall(dlsym(libgiac_c[], "giac_new_complex_int"), Ptr{Nothing}, (Cint, Cint), real(val), real(val)))
end

function giac(val::Complex{Cdouble})
    giac(ccall(dlsym(libgiac_c[], "giac_new_complex_double"), Ptr{Nothing}, (Cdouble, Cdouble), real(val), real(val)))
end

function giac_identifier(s::AbstractString)
    giac(ccall(dlsym(libgiac_c[], "giac_new_ident"), Ptr{Nothing}, (Cstring,), s))
end

function giac(s::AbstractString)
    giac(ccall(dlsym(libgiac_c[], "giac_new_symbolic"), Ptr{Nothing}, (Cstring,Ptr{Nothing}), s, context_ptr[]))
end

function giac(v::Array{T,1}; subtype::Integer=0) where T
    v1 = Ptr{Nothing}[giac(i).g for i in v]
    giac(ccall(dlsym(libgiac_c[], "giac_new_vector"), Ptr{Nothing},
              (Ptr{Ptr{Nothing}},Cint,Cshort),
               v1, length(v1), subtype))
end

function giac(A::Array{T,2}) where T
    giac([reshape(A[i,:], size(A,2)) for i in 1:size(A,1)])
end

giac(x::AbstractRange{T}) where T = giac(collect(x))

giac(::Irrational{:ℯ}) = giac_e[]
giac(::Irrational{:π}) = giac_pi[]

one(::Type{giac}) = giac_one[]
one(::giac) = giac_one[]
oneunit(::Type{giac}) = giac_one[]
oneunit(::giac) = giac_one[]
zero(::Type{giac}) = giac_zero[]
zero(::giac) = giac_zero[]

const _r_giac_julia = Regex(string("\\b(",join(keys(giac_julia_dictionay),"|"),")\\("))
_f_giac_julia(x) = string(giac_julia_dictionay[Symbol(x[1:end-1])], x[end])

function string(g::giac)
   cs = ccall(dlsym(libgiac_c[], "giac_to_string"), Ptr{UInt8}, (Ptr{Nothing},Ptr{Nothing}), g.g, context_ptr[])
   s = unsafe_string(cs)
   ccall(dlsym(libgiac_c[], "giac_free"), Nothing, (Ptr{Nothing},), cs)
   replace(s, _r_giac_julia => _f_giac_julia)
end

function latex(g::giac)
   cs = ccall(dlsym(libgiac_c[], "giac_to_latex"), Ptr{UInt8}, (Ptr{Nothing},Ptr{Nothing}), g.g, context_ptr[])
   s = unsafe_string(cs)
   ccall(dlsym(libgiac_c[], "giac_free"), Nothing, (Ptr{Nothing},), cs)
   s
end

prettyprint(ex::giac) = display("text/latex", string("\$\$",latex(ex),"\$\$"))

#show(io::IO, g::giac) = print(io, string(g))
const giac_color_on = Base.text_colors[:light_blue]
const giac_color_off = Base.text_colors[:normal]
show(io::IO, g::giac) = print(io, giac_color_on, string(g), giac_color_off)
#show(io::IO, g::giac) = print(io, "giac(\"", string(g), "\")")

#function writemime(io::IO, mime::MIME"text/latex", ex::giac)
#    write(io, "\$\$", latex(ex), "\$\$")
#end

function +(a::giac, b::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_plus"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g))
end
+(a::giac, b::Number) = a+giac(b)
+(a::Number, b::giac) = giac(a)+b
+(a::giac) = a

function -(a::giac, b::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_minus"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g))
end
-(a::giac, b::Number) = a-giac(b)
-(a::Number, b::giac) = giac(a)-b

function -(a::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_uminus"), Ptr{Nothing}, (Ptr{Nothing},), a.g))
end

function *(a::giac, b::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_times"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g))
end
*(a::giac, b::Number) = a*giac(b)
*(a::Number, b::giac) = giac(a)*b

function /(a::giac, b::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_rdiv"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g))
end
/(a::giac, b::Number) = a/giac(b)
/(a::Number, b::giac) = giac(a)/b
inv(a::giac) = giac_one[]/a

function ^(a::giac, b::giac)
    giac(ccall(dlsym(libgiac_c[], "giac_pow"), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing}), a.g, b.g, context_ptr[]))
end
^(a::giac, b::Integer) = a^giac(b)  # to ovverride ^(::Any, Integer) which is already defined
^(a::giac, b::Number) = a^giac(b)
^(a::Number, b::giac) = giac(a)^b

function ==(a::giac, b::giac)
    ccall(dlsym(libgiac_c[], "giac_equal_bool"), Cint, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g)!=0
end
==(a::giac, b::Number) = a==giac(b)
==(a::Number, b::giac) = giac(a)==b

function >(a::giac, b::giac)
    ccall(dlsym(libgiac_c[], "giac_greater_than"), Cint, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g) !=0
end

function >=(a::giac, b::giac)
    ccall(dlsym(libgiac_c[], "giac_greater_equal"), Cint, (Ptr{Nothing},Ptr{Nothing}), a.g, b.g) !=0
end

>(a::giac, b::Real) = a>giac(b)
>(a::Real, b::giac) = giac(a)>b
>=(a::giac, b::Real) = a>=giac(b)
>=(a::Real, b::giac) = giac(a)>=b
<(a::giac, b::giac) = b>a
<(a::Real, b::giac) = giac(a)<b
<(a::giac, b::Real) = a<giac(b)
<=(a::giac, b::giac) = b>=a
<=(a::giac, b::Real) = a<=giac(b)
<=(a::Real, b::giac) = giac(a)<=b

function size(g::giac)
    ccall(dlsym(libgiac_c[], "giac_size1"), Cint, (Ptr{Nothing},), g.g)
end

length(g::giac) = size(g) #

function getindex(g::giac, i)
    giac(ccall(dlsym(libgiac_c[], "giac_getindex"), Ptr{Nothing}, (Ptr{Nothing},Cint), g.g, i-1))
end

function (g::giac)(x)
    giac(ccall(dlsym(libgiac_c[], "giac_call"), Ptr{Nothing},
    (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing}), g.g, giac(x).g, context_ptr[]))
end

function (g::giac)(x...)
    giac(ccall(dlsym(libgiac_c[], "giac_call"), Ptr{Nothing},
    (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing}), g.g, giac([x...],subtype=1).g, context_ptr[]))
end

# unary functions with context_ptr:
for F in (:real, :imag, :conj, :abs,
          :sqrt, :exp, :log, :sin, :cos, :tan,
          :sinh, :cosh, :tanh, :asin, :acos, :atan,
          :asinh, :acosh, :atanh)
   @eval begin
       function ($F)(a::giac)
           giac(ccall(dlsym(libgiac_c[], $(string("giac_",F))), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), a.g, context_ptr[]))
       end
   end
end

adjoint(x::giac) = conj(x)

function giac(f::Symbol, arg)
    giac(ccall(dlsym(libgiac_c[], string("giac_", f)), Ptr{Nothing}, (Ptr{Nothing},Ptr{Nothing}), giac(arg).g, context_ptr[]))
end

giac(f::Symbol, args...) = giac(f, giac(Any[args...], subtype=1))

include("library.jl")

unapply(ex::giac,var) = giac(:unapply, ex, var)
set!(var::giac, val) = giac(:sto, val, var)
purge!(var::giac) = giac(:purge, var)
giac_vars() = [Pair(left(x), right(x)) for x in to_julia(giac(:VARS, 1))]

#This magic code inspired by SymPy, cf.
#https://github.com/jverzani/SymPy.jl/blob/master/src/utils.jl
macro giac(x...)
    q=Expr(:block)
    r=Expr(:tuple)
    if length(x) == 1 && isa(x[1],Expr)
        if (x[1].head === :tuple)
            x = x[1].args
        end
    end
    for s in x
        if isa(s,Symbol)
            push!(q.args, Expr(:(=), s, Expr(:call, :giac, Expr(:quote, string(s)))))
            # Use here :giac (:giac_identifier should also work, but this leads to
            # strange behavior...
            if length(x)>1
                push!(r.args, s)
            end
        elseif isa(s, Expr)&&s.head==:(=)&&isa(s.args[1],Symbol)
            push!(q.args, Expr(:(=), s.args[1], Expr(:call, :giac,
                  Expr(:quote, string(s.args[1])))))
            if s.args[1] in config_vars
                push!(q.args, Expr(:call, s.args[1], s.args[2]))
            else
                push!(q.args, Expr(:call, :set!, s.args[1], s.args[2]))
            end
            if length(x)>1
                push!(r.args, s.args[2])
            end
        elseif (isa(s, Expr)&&s.head==:(=)&&isa(s.args[1], Expr)
              &&s.args[1].head==:(call)&&isa(s.args[1].args[1], Symbol))
            push!(q.args, Expr(:(=), s.args[1].args[1],
                 Expr(:call, :giac, Expr(:quote, string(s.args[1].args[1])))))
            if length(s.args[1].args)>2
                push!(q.args, Expr(:call, :set!, s.args[1].args[1], Expr(:call, :unapply, s.args[2],
                      Expr(:vect, s.args[1].args[2:end]...))))
            else
                push!(q.args, Expr(:call, :set!, s.args[1].args[1], Expr(:call, :unapply, s.args[2],
                      s.args[1].args[2])))
            end
            if length(x)>1
                push!(r.args, s.args[1].args[2])
            end
        else
            @assert false "@giac expected a list of symbols or var=val expressions"
        end
    end
    if length(x)>1
        push!(q.args, r)
    end
    Core.eval(Main, q)
end

@eval const $(Symbol("@vars")) = $(Symbol("@giac"))

to_julia(x) = x

function to_julia(g::giac)
  T = type(g)
  if T == :INT

   z = ccall(dlsym(libgiac_c[], "giac_get_int"), Cint, (Ptr{Nothing},), g.g)
   if subtype(g)==6 #Bool
       return z==1
   else
       return z
   end

  elseif T == :DOUBLE || T == :FLOAT
    return ccall(dlsym(libgiac_c[], "giac_get_double"), Cdouble, (Ptr{Nothing},), g.g)

  elseif T == :ZINT
    z = BigInt()
    m = ccall(dlsym(libgiac_c[], "giac_get_bigint"), Ptr{BigInt}, (Ptr{Nothing},), g.g)
    ccall((:__gmpz_set,:libgmp), Nothing, (Ref{BigInt}, Ptr{BigInt}), z, m)
    return z

  elseif T == :REAL
    z = BigFloat()
    m = ccall(dlsym(libgiac_c[], "giac_get_bigfloat"), Ptr{BigFloat}, (Ptr{Nothing},), g.g)
    ccall((:mpfr_set, :libmpfr), Int32, (Ref{BigFloat}, Ptr{BigFloat}, Int32),
          z, m, 0)
    return z

  elseif T == :CPLX
    return complex(to_julia(real(g)), to_julia(imag(g)))

  elseif T == :FRAC
    return to_julia(numerator(g))//to_julia(denominator(g))

  elseif T == :VECT
    return [to_julia(g[i]) for i=1:length(g)]
  else
    return g
  end
end

convert(::Type{giac}, x) = giac(x)
convert(::Type{giac}, x::giac) = x # to avoid some ambiguities

function convert(::Type{T}, ex::giac) where T<:Number
    if isgiacnumber(ex)
        return convert(T,to_julia(ex))
    end
    ex = evaluate(ex)
    if isgiacnumber(ex)
       return convert(T,to_julia(ex))
    end
    ex = evaluatef(ex)
    if isgiacnumber(ex)
       return convert(T,to_julia(ex))
    end
    error("Could not convert giac to Number")
end

_head(ex::giac) = Symbol(string(head(ex))[2:end-1])

function _args(ex)
    a = args(ex)
    if isa(a, giac) && type(a)==:VECT && subtype(a)==1
        return to_julia(a)
    else
        return Any[to_julia(a)]
    end
end

_expr(x) = x

function _expr(ex::giac)
  T = type(ex)
  if T == :IDNT
    return Symbol(string(ex))
  elseif  T == :SYMB
    h = _head(ex)
    if h==:id
        return _expr(args(ex))
    end
    if h==:of
        ex = evaluate(ex)
        h = _head(ex)
        @assert h!=:of "could not evaluate user-defined function"
    end
    h = get(giac_julia_dictionay, h, h)
    return Expr(:call, h, [_expr(arg) for arg in _args(ex)]...)
  else
    return ex
  end
end

"""
    lambdify(ex, vars...)

Take a symbolic expression and return an anonymous `Julia` function

"""
function lambdify(ex::giac, var::giac, vars...)
    @assert type(ex)==:SYMB && type(var)==:IDNT
    if length(vars) == 0
        eval(Expr(:->, _expr(var), _expr(evaluate(ex))))
    else
        eval(Expr(:->, Expr(:tuple, _expr(var),
            [_expr(v) for v in vars]...), _expr(evaluate(ex))))
    end
end

function convert(::Type{Function}, f::giac)
    if type(f)==:SYMB || type(f)==:IDNT
         a = args(evaluate(f))
         return lambdify(to_julia(a[3]),to_julia(a[1])...)
    else
         error("Could not convert giac to Function")
    end
end

end # module


text:
using LsqFit, Distributions, DSP, DelimitedFiles

struct input_data
    t::Array{Float64,1}
	DTOF::Array{Float64,1}
	IRF::Array{Float64,1}
	DTOFerrors::Array{Float64,1}
	ρ::Float64
	nmed::Float64
	ndet::Float64
	lambda::Float64
	normpeaks::Bool
end

struct model_params
	#yerrors::Array{Float64,1}          #uncertainties in dependent variable
    model::Function                    #model to be fitted against
	initparams::Array{Float64,1}
	lb::Array{Float64,1}
	ub::Array{Float64,1}
	risefactor::Float64  			#percent of peak to fit on rising edge
	tailfactor::Float64				#perfecnt of peak to fit on failing tail
	dof::Int64
end

struct fitresult
    #Fit inputs
    xdata::Array{Float64,1}     #independent variable
    ydata::Array{Float64,1}     #dependent variable
    yerrors::Array{Float64,1}   #uncertainties in dependent variable
	model::Function             #model to be fitted against

    #Fit results
    residuals::Array{Float64,1}  	 #weighted residuals
    params::Array{Float64,1}    		#best fit parameters for model
	marginerror::Array{Float64,1}        #product of standard error and the critical value of each parameter at a certain significance level (default is 5%) from t-distribution.
    perrors::Array{Float64,1}            #uncertainties in best fit parameters
    dof::Int64                  		#number of degrees of freedom in fit
    chisq::Float64             		 #chi^2 value of fit
    probability::Float64      		  #probability of chi^2 value
end

function DT_model(t,β,ρ,nmed,ndet)
	#diffusion theory model for homogenous semi-infinite layer
	μa = β[1]
	μsp = β[2]

	if (nmed == ndet)
		Afac = 1
	elseif nmed > ndet
		costhetac = sqrt(1 - (ndet/nmed).^2)
		R0 = ((nmed/ndet-1)/(nmed/ndet +1)).^2
		Afac = (2/(1-R0) -1 + abs(costhetac.^3))./(1-abs(costhetac.^2))
	else nmed < ndet
		R0 = ((nmed/ndet-1)/(nmed/ndet +1)).^2
		Afac = 2/(1-R0) - 1
	end

	z0 = 1/(μa+μsp)
	D = 1/(3*(μa + μsp))
	zb = 2*Afac*D

	v = 29.9792345/nmed # speed of light in medium

	Rt1 = @. v*exp(-μa*v*t)
	Rt1 = @. Rt1/((4*π*D*v*t)^1.5)
	Rt1 = @. Rt1*(exp(-(z0^2 + ρ^2)/(4*D*v*t)) - exp(-((2*zb + z0)^2 + ρ^2)/(4*D*v*t)))

	Rt2 = @. 3*exp(-μa*v*t)/(2*((4*π*D*v)^1.5)*(t^2.5))
	Rt2 = @. Rt2*(z0*exp(-((z0^2 + ρ^2)/(4*D*v*t))) + (2*zb + z0)*exp(-((2*zb+z0)^2 + ρ^2)/(4*D*v*t)))

	Rt = @. abs(Rt1)+abs(Rt2)
	replace!(Rt, NaN => 0)
	m = findmax(Rt)
	Rt = Rt./m[1]
end

function conv_DT(t, β, data)

	RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)

	convDT = conv(data.IRF, RtDT)
	convDT = convDT./maximum(convDT)

	tidx = findfirst(x -> x == t[1], data.t)

	if data.normpeaks == true
		peakindDTOF, peakindconvDT = findmax(data.DTOF)[2], findmax(convDT)[2]
		if peakindDTOF > peakindconvDT
			convDT = [zeros(peakindDTOF-peakindconvDT); convDT]
		elseif peakindDTOF < peakindconvDT
			convDT = convDT[peakindconvDT-peakindDTOF:end]
		end
	end

	convDT = convDT[tidx:tidx+length(t)-1]
	return log.(convDT)
end

function getfit(input_data, model_params; alpha=0.1)

	#get fitting window
	counts = input_data.DTOF
	maxvalue, maxindex = findmax(counts)

	ind1 = findfirst(x -> x > maxvalue*model_params.risefactor, counts)
	ind2 =  findlast(x -> x >maxvalue*model_params.tailfactor, counts)

	weights = 1 ./ input_data.DTOFerrors

	fit = curve_fit((t, β) -> conv_DT(t, β, input_data), input_data.t[ind1:ind2], log.(input_data.DTOF[ind1:ind2]), model_params.initparams, lower=model_params.lb, upper=model_params.ub)

	perrors = estimate_errors(fit,alpha)
	marginerror = margin_error(fit,alpha)

	chisq = sum(fit.resid.^2)
	chisqdist = Distributions.Chisq(2)
	probability = Distributions.ccdf(chisqdist,chisq)
	return fitresult(input_data.t[ind1:ind2],log.(input_data.DTOF[ind1:ind2]),input_data.DTOFerrors, DT_model, fit.resid,fit.param,marginerror,perrors,2,chisq,probability)
end

function plotfit(fit::fitresult1, data::input_data)
	scatter(fit.xdata, fit.ydata, color="black", label="DTOF", markersize=3, alpha=0.8)
	plot!(fit.xdata, conv_DT(fit.xdata, fit.params, data), color="red", lw=3, label="DT-fit", alpha=0.8)

end

function loadData(filename)
    #loads data from asc
    data_array = readdlm(filename,skipstart=10)
    t = data_array[1:end-1,1]
    Rt = data_array[1:end-1,2]
    (time = convert(Array{Float64},t), counts = convert(Array{Float64},Rt))
end

max_stars_repo_name: heltonmc/load-asc-data

max_issues_repo_name: heltonmc/load-asc-data

max_forks_repo_name: heltonmc/load-asc-data


text:
using Test
using HMatrices
using LinearAlgebra
using Random
using StaticArrays

using HMatrices: RkMatrix

include(joinpath(HMatrices.PROJECT_ROOT,"test","testutils.jl"))

Random.seed!(1)

m = 1000
T = Float64
X    = points_on_sphere(m)
Y    = X
struct ExponentialKernel <: AbstractMatrix{Float64}
    X::Vector{SVector{3,Float64}}
    Y::Vector{SVector{3,Float64}}
end
function Base.getindex(K::ExponentialKernel,i::Int,j::Int)
    x,y = K.X[i], K.Y[j]
    exp(-norm(x-y))
end
Base.size(K::ExponentialKernel) = length(K.X), length(K.Y)
K = ExponentialKernel(X,X)

splitter  = CardinalitySplitter(nmax=50)
Xclt      = ClusterTree(X,splitter)
Yclt      = ClusterTree(Y,splitter)
H = assemble_hmat(K,Xclt,Yclt;threads=false,distributed=false)
H_full      = Matrix(H;global_index=false)

@testset "ldiv!" begin
    B = rand(m,2)
    R = RkMatrix(rand(m,3),rand(m,3))
    R_full = Matrix(R)

    ## 3.1
    exact  = ldiv!(UnitLowerTriangular(H_full),copy(B))
    approx = ldiv!(UnitLowerTriangular(H),copy(B))
    @test exact ≈ approx

    exact  = ldiv!(UpperTriangular(H_full),copy(B))
    approx = ldiv!(UpperTriangular(H),copy(B))
    @test exact ≈ approx

    ## 3.2
    exact  = ldiv!(UnitLowerTriangular(H_full),copy(R_full))
    approx = ldiv!(UnitLowerTriangular(H),copy(R))
    @test exact ≈ Matrix(approx)

    ## 3.3
    compressor = PartialACA(;atol=1e-8)
    exact      = ldiv!(UnitLowerTriangular(H_full),copy(H_full))
    approx     = ldiv!(UnitLowerTriangular(H),deepcopy(H),compressor)
    @test exact ≈ Matrix(approx;global_index=false)
end

@testset "rdiv!" begin
    B = rand(2,m)
    R = RkMatrix(rand(m,3),rand(m,3))
    R_full = Matrix(R)
    # 3.1
    exact  = rdiv!(copy(B),UpperTriangular(H_full))
    approx = rdiv!(copy(B),UpperTriangular(H))
    @test exact ≈ approx

    ## 3.2
    exact  = rdiv!(copy(R_full),UpperTriangular(H_full))
    approx = rdiv!(copy(R),UpperTriangular(H))
    @test exact ≈ approx

    ## 3.3
    compressor = PartialACA(;atol=1e-10)
    exact  = rdiv!(deepcopy(H_full),UpperTriangular(H_full))
    approx = rdiv!(deepcopy(H),UpperTriangular(H),compressor)
    @test exact ≈ Matrix(approx;global_index=false)
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: https://julialang.org/license

const options = [
    "--cflags",
    "--ldflags",
    "--ldlibs"
];

threadingOn() = ccall(:jl_threading_enabled, Cint, ()) != 0

function shell_escape(str)
    str = replace(str, "'", "'\''")
    return "'$str'"
end

function imagePath()
    opts = Base.JLOptions()
    return unsafe_string(opts.image_file)
end

function libDir()
    return if ccall(:jl_is_debugbuild, Cint, ()) != 0
        dirname(abspath(Libdl.dlpath("libjulia-debug")))
    else
        dirname(abspath(Libdl.dlpath("libjulia")))
    end
end

private_libDir() = abspath(JULIA_HOME, Base.PRIVATE_LIBDIR)

function includeDir()
    return abspath(JULIA_HOME, Base.INCLUDEDIR, "julia")
end

function ldflags()
    fl = "-L$(shell_escape(libDir()))"
    if Sys.iswindows()
        fl = fl * " -Wl,--stack,8388608"
    elseif Sys.islinux()
        fl = fl * " -Wl,--export-dynamic"
    end
    return fl
end

function ldlibs()
    libname = if ccall(:jl_is_debugbuild, Cint, ()) != 0
        "julia-debug"
    else
        "julia"
    end
    if Sys.isunix()
        return "-Wl,-rpath,$(shell_escape(libDir())) -Wl,-rpath,$(shell_escape(private_libDir())) -l$libname"
    else
        return "-l$libname -lopenlibm"
    end
end

function cflags()
    flags = IOBuffer()
    print(flags, "-std=gnu99")
    include = shell_escape(includeDir())
    print(flags, " -I", include)
    if threadingOn()
        print(flags, " -DJULIA_ENABLE_THREADING=1")
    end
    if Sys.isunix()
        print(flags, " -fPIC")
    end
    return String(take!(flags))
end

function check_args(args)
    checked = intersect(args, options)
    if length(checked) == 0 || length(checked) != length(args)
        println(STDERR, "Usage: julia-config [", join(options, " | "), "]")
        exit(1)
    end
end

function main()
    check_args(ARGS)
    for args in ARGS
        if args == "--ldflags"
            println(ldflags())
        elseif args == "--cflags"
            println(cflags())
        elseif args == "--ldlibs"
            println(ldlibs())
        end
    end
end

main()

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
import DataStructures: DefaultDictBase

@testset "DefaultDict" begin

    @testset "DefaultDictBase" begin
        #construction
        @test_throws ArgumentError DefaultDictBase()
        @test isa(DefaultDictBase(0.0), DefaultDictBase{Any, Any, Float64, Dict{Any, Any}})
        @test isa(DefaultDictBase(0.0, 1 => 1.0), DefaultDictBase{Int, Float64, Float64, Dict{Int, Float64}})
        @test isa(DefaultDictBase(0.0, [(1, 1.0)]), DefaultDictBase{Int, Float64, Float64, Dict{Int, Float64}})
        #@test isa(DefaultDictBase(0.0, Dict()), DefaultDictBase{Any, Any, Float64, Dict{Any, Any}}))

        ddb = DefaultDictBase{Int, Float64}(0.0)
        @test isa(ddb, DefaultDictBase{Int, Float64, Float64, Dict{Int,Float64}})
        @test isa(DefaultDictBase(1.0, ddb), DefaultDictBase{Int, Float64, Float64, Dict{Int,Float64}})

        @test isa(DefaultDictBase{Int, String}(String), DefaultDictBase{Int, String, typeof(String), Dict{Int, String}})
        @test isa(DefaultDictBase{Int, String}(String; passkey=true), DefaultDictBase{Int, String, typeof(String), Dict{Int, String}})
        @test isa(DefaultDictBase{Int, String}(String; passkey=false), DefaultDictBase{Int, String, typeof(String), Dict{Int, String}})
    end

    @testset "DefaultDict" begin
        @testset "construction" begin
            @test_throws ArgumentError DefaultDict()
            @test_throws ArgumentError DefaultDict(AbstractString, Int)
            @test_throws ArgumentError DefaultDict{AbstractString, Int}()

            @test isa(DefaultDict(0.0, 1 => 1.0), DefaultDict{Int, Float64, Float64})
        end

        @testset "Core Functionality" begin
            for F in [1, ()->1]
                # empty dictionary
                d = DefaultDict{Char, Int}(F)
                @test length(d) == 0
                @test isempty(d)
                @test d['c'] == 1
                @test setindex!(d, 10, 'd') === d
                @test !isempty(d)
                @test empty!(d) === d
                @test isempty(d)

                # access, modification
                @test (d['a'] += 1) == 2
                @test 'a' in keys(d)
                @test haskey(d, 'a')
                @test get(d, 'b', 0) == 0
                @test !('b' in keys(d))
                @test !haskey(d, 'b')
                @test pop!(d, 'a') == 2
                @test isempty(d)

                # pushing multiple pairs
                @test push!(d, 'c'=>3, 'd'=>4, 'e'=>5, 'f'=>6) === d
                empty!(d)

                # pushing multiple tuples
                @test push!(d, ('c',3), ('d',4), ('e',5), ('f',6)) === d
                empty!(d)

                @test sizehint!(d, 26) === d

                for c in 'a':'z'
                    d[c] = c-'a'+1
                end

                @test d['z'] == 26
                @test d['@'] == 1
                @test length(d) == 27
                @test delete!(d, '@') === d
                @test length(d) == 26

                for (k,v) in d
                    @test v == k-'a'+1
                end

                s = empty(d)
                @test typeof(s) === typeof(d)
                @test s.d.default == d.d.default

                @test sort(collect(keys(d))) == collect('a':'z')
                @test sort(collect(values(d))) == collect(1:26)
            end
        end

        @testset "From Dict" begin
            # Starting from an existing dictionary
            # Note: dictionary is copied upon construction
            e = Dict([('a',1), ('b',3), ('c',5)])
            f = DefaultDict(0, e)
            @test f['d'] == 0
            @test_throws KeyError e['d']
            e['e'] = 9
            @test e['e'] == 9
            @test f['e'] == 0
        end

        @testset "passkey" begin
            calls = 0
            g = DefaultDict{String, Int}(passkey=true) do key
                calls += 1
                return length(key)
            end
            @test g["foobar"] == 6
            @test calls == 1
            @test length(g) == 1
            @test g["baz"] == 3
            @test calls == 2
            @test length(g) == 2
            @test g["foobar"] == 6
            @test calls == 2
            @test length(g) == 2

            @testset "Incorrect Usage" begin
                bad_dds = [
                    DefaultDict{Int, Int}(k -> k, passkey=false),
                    DefaultDict{Int, Int}(() -> 3, passkey=true),
                ]
                for bad_dd in bad_dds
                    bad_dd[3] = 10
                    @test bad_dd[3] == 10
                    @test_throws MethodError bad_dd[234]
                end
            end
        end

        # Alternate constructor
        @test isa(DefaultDict(0.0), DefaultDict{Any, Any, Float64})
        @test isa(DefaultDict(0.0, [(1, 1.0)]), DefaultDict{Int, Float64, Float64})
    end

    @testset "DefaultOrderedDict" begin
        @testset "construction" begin
            @test_throws ArgumentError DefaultOrderedDict()
            @test_throws ArgumentError DefaultOrderedDict{AbstractString, Int}()
        end

        @testset "Core Functionality" begin
            for F in [1, ()->1]
                # empty dictionary
                d = DefaultOrderedDict{Char, Int}(F)
                @test length(d) == 0
                @test isempty(d)
                @test d['c'] == 1
                @test setindex!(d, 10, 'd') === d
                @test !isempty(d)
                @test empty!(d) === d
                @test isempty(d)

                # access, modification
                @test (d['a'] += 1) == 2
                @test 'a' in keys(d)
                @test haskey(d, 'a')
                @test get(d, 'b', 0) == 0
                @test !('b' in keys(d))
                @test !haskey(d, 'b')
                @test pop!(d, 'a') == 2
                @test isempty(d)

                # pushing multiple pairs
                @test push!(d, 'c'=>3, 'd'=>4, 'e'=>5, 'f'=>6) === d
                empty!(d)

                # pushing multiple tuples
                @test push!(d, ('c',3), ('d',4), ('e',5), ('f',6)) === d
                empty!(d)

                @test sizehint!(d, 26) === d

                for c in 'a':'z'
                    d[c] = c-'a'+1
                end

                @test d['z'] == 26
                @test d['@'] == 1
                @test length(d) == 27
                @test delete!(d, '@') === d
                @test length(d) == 26

                for (k,v) in d
                    @test v == k-'a'+1
                end

                @test collect(keys(d)) == collect('a':'z')
                @test collect(values(d)) == collect(1:26)

                s = empty(d)
                @test typeof(s) === typeof(d)
                @test s.d.default == d.d.default
            end
        end

        @testset "Alternate constructor" begin
            @test isa(DefaultOrderedDict(0.0), DefaultOrderedDict{Any, Any, Float64})
            @test isa(DefaultOrderedDict(0.0, [(1, 1.0)]), DefaultOrderedDict{Int, Float64, Float64})
        end

        @testset "issue #216" begin
            @test DataStructures.isordered(DefaultOrderedDict{Int, String})
        end

    end

end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Diban
using Thrift

using Parquet, ParquetFiles,

load("c:/scratch/python-test.parquet")

p = ParFile("c:/scratch/python-test.parquet")
p

mutable struct Ghi2
    a::Vector{UInt8}
end

rc = RecCursor(p, 1:nrows(p), colnames(p), JuliaBuilder(p, Ghi2))

@which RecCursor(p, 1:nrows(p), colnames(p), JuliaBuilder(p, Ghi2))

using Debugger
Debugger.@enter ColCursor(par, 1:2, "a", 2)
@which ColCursor(par, 1:2, "a", 2)
@which RowCursor(par, 1:3, "a", 1)

rowgroups(par)
rowgroups(par, "a", 1:3)

rc.colcursors

rc.colstates
rc.colcursors[1]

record_state = iterate(rc)

while record_state != nothing
    global record_state
    record = record_state[1]
    state = record_state[2]
    println(record)
    record_state = iterate(rc, state)
end

rci, state = iterate(rc)

rci, state = iterate(rc, state)

rci, state = iterate(rc, state)

state

rci

a =  iterate(rc, state)

a == nothing

iterate(rc, state)

T_row_name = Symbol("RCType$(String(gensym())[3:end])")

@which schema(JuliaConverter(ParquetFiles), p, T_row_name)

conv = JuliaConverter(ParquetFiles)
par = p
schema_name = T_row_name
@which schema(conv, par.schema, schema_name)

sch = par.schema
@which Parquet.schema_to_julia_types(conv.to, sch, schema_name)

io = IOBuffer()
@which Parquet.schema_to_julia_types(io, sch, schema_name)

@which Parquet.schema_to_julia_types(io, sch.schema, schema_name)

sch.schema[1]
sch.schema[2]

typeof(sch.schema[1])

Parquet.isfilled(sch.schema[1], :_type)
Parquet.isfilled(sch.schema[2], :_type)

sch.schema[1].name

fieldnames(typeof(sch.schema[1]))

typestr = "begin\n" * String(take!(io)) * "\nend"
parsedtypes = Meta.parse(typestr)
Core.eval(mod, parsedtypes)

T_row = eval(:haha)

col_names = [:a]
col_types = [i <: Vector{UInt8} ? String : i for i in [String]]

T = NamedTuple{(col_names...,),Tuple{col_types...}}

rc = RecCursor(p, 1:nrows(p), colnames(p), JuliaBuilder(p, T_row))

it = ParquetNamedTupleIterator{T,T_row}(rc, nrows(p))

return it

using Debugger

using ParquetFiles

Debugger.@enter(load("c:/scratch/python-test.parquet"))

io = open("c:/scratch/python-test.parquet", "r")
io = open("C:/data/parquet-test/fannie_mae_perf_small/fannie_mae_perf_small.parquet")
io = open("c:/scratch/test.parquet")
#b = read(io)
#close(io)

using ParquetWriter

metadata("c:/scratch/test.parquet")

par = ParFile("c:/scratch/test.parquet")
Parquet.pages(par, filemetadata.row_groups[1].columns[1])

filemetadata

Parquet.PLAIN_JTYPES

const TYPES = (Bool, Int32, Int64, Int128, Float32, Float64, String, UInt8)

filemetadata.schema[1]._type
filemetadata.schema[2]._type

filemetadata

nrows(par)

pgs = Parquet.pages(par, filemetadata.row_groups[1].columns[1])

vals = values(par, pgs[1])

@which values(par, pgs[1])

@which bytes(pgs[1])

page =  pgs[1]
rawbytes = Parquet.bytes(page)
num_values = Parquet.page_num_values(page)

io = IOBuffer(rawbytes)

ctype = Parquet.coltype(page.colchunk)

ok = Parquet.read_plain_values(io, 131078, ctype)

data = page.data
codec = page.colchunk.meta_data.codec
uncompressed=true

@which Parquet.page_num_values(page)
@which Parquet.page_num_values(page.hdr)

Thrift.isfilled(page.hdr, :data_page_header_v2)
Thrift.isfilled(page.hdr, :data_page_header)
Thrift.isfilled(page.hdr, :dictionary_page_header)

@which Parquet.page_num_values(page.hdr.dictionary_page_header)

page.hdr.dictionary_page_header

codec != Parquet.CompressionCodec.UNCOMPRESSED
using Snappy
data = Snappy.uncompress(data)

SZ_PAR_MAGIC = 4#
SZ_FOOTER = 4

sz = filesize(io)
seek(io, sz -8)
len = read(io, Int32)

data_size = sz -  len -12
SZ_PAR_MAGIC + data_size
seek(io, SZ_PAR_MAGIC + data_size)
filemetadata = read_thrift(io, PAR2.FileMetaData)

sz = filesize(io)

seek(io, sz - SZ_PAR_MAGIC - SZ_FOOTER)

# read footer size as little endian signed Int32

pg = Parquet.pages(par, a.row_groups[1].columns[1])
values(par, pg)
par

load_parquet

a.row_groups[1].file_offset

seek(io, 113)

read_thrift(io, ParquetThrift.ColumnMetaData)

seek(io, 100)
position(io)
read_thrift(io, ParquetThrift.PageHeader )

ParquetThrift.DataPageHeader

a.row_groups[1].columns

a.row_groups

using Parquet

ParFile("c:/scratch/python-test.parquet")

using Test
@test String(b[1:4]) == "PAR1"
@test String(b[end-3:end]) == "PAR1"
footer_length = reinterpret(Int32, b[end-7:end-4])[1]

data = b[5:5+length(b)-4-4-4-footer_length-1]

footer = b[end-8-footer_length+1:end-8]

Int.(footer[1:4])
fileversion = reinterpret(Int32, footer[1:4])[1]

# next is the schema

# schema starts with a type
Int.(footer[5:8])
footer[5:8]
reinterpret(Int32, footer[5:8])[1]
593 613
findall(footer .== UInt('y'))

footer[594]
footer[614]

for i = 1:length(footer)-7
    if reinterpret(Int64, footer[i:i+7])[1] == 1
        println(i)
    end
end


text:
# Formulation K2: (if regul==:classic, adds additional regularization parmeters -ρ (top left) and δ (bottom right))
# [-Q - D     A' ] [x] = rhs
# [ A         0  ] [y]
export K2LDLParams

"""
Type to use the K2 formulation with a LDLᵀ factorization, using the package
[`LDLFactorizations.jl`](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl).
The outer constructor

    sp = K2LDLParams(; regul :: Symbol = :classic)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
`regul = :dynamic` uses a dynamic regularization (the regularization is only added if the LDLᵀ factorization
encounters a pivot that has a small magnitude).
`regul = :none` uses no regularization (not recommended).
"""
struct K2LDLParams <: SolverParams
    regul :: Symbol
end

function K2LDLParams(; regul :: Symbol = :classic)
    regul == :classic || regul == :dynamic || regul == :none || error("regul should be :classic or :dynamic or :none")
    return K2LDLParams(regul)
end

mutable struct PreallocatedData_K2{T<:Real} <: PreallocatedData{T}
    D                :: Vector{T}                                        # temporary top-left diagonal
    regu             :: Regularization{T}
    diag_Q           :: SparseVector{T,Int} # Q diagonal
    K                :: SparseMatrixCSC{T,Int} # augmented matrix
    K_fact           :: LDLFactorizations.LDLFactorization{T,Int,Int,Int} # factorized matrix
    fact_fail        :: Bool # true if factorization failed
    diagind_K        :: Vector{Int} # diagonal indices of J
end

# outer constructor
function PreallocatedData(sp :: K2LDLParams, fd :: QM_FloatData{T}, id :: QM_IntData,
                          iconf :: InputConfig{Tconf}) where {T<:Real, Tconf<:Real}

    # init Regularization values
    if iconf.mode == :mono
        regu = Regularization(T(sqrt(eps())*1e5), T(sqrt(eps())*1e5), 1e-5*sqrt(eps(T)), 1e0*sqrt(eps(T)), sp.regul)
        D = -T(1.0e0)/2 .* ones(T, id.nvar)
    else
        regu = Regularization(T(sqrt(eps())*1e5), T(sqrt(eps())*1e5), T(sqrt(eps(T))*1e0), T(sqrt(eps(T))*1e0), sp.regul)
        D = -T(1.0e-2) .* ones(T, id.nvar)
    end
    diag_Q = get_diag_Q(fd.Q.colptr, fd.Q.rowval, fd.Q.nzval, id.nvar)
    K = create_K2(id, D, fd.Q, fd.AT, diag_Q, regu)

    diagind_K = get_diag_sparseCSC(K.colptr, id.ncon+id.nvar)
    K_fact = ldl_analyze(Symmetric(K, :U))
    if regu.regul == :dynamic
        Amax = @views norm(K.nzval[diagind_K], Inf)
        regu.ρ, regu.δ = -T(eps(T)^(3/4)), T(eps(T)^(0.45))
        K_fact.r1, K_fact.r2 = regu.ρ, regu.δ
        K_fact.tol = Amax*T(eps(T))
        K_fact.n_d = id.nvar
    elseif regu.regul == :none
        regu.ρ, regu.δ = zero(T), zero(T)
    end
    K_fact = ldl_factorize!(Symmetric(K, :U), K_fact)
    K_fact.__factorized = true

    return PreallocatedData_K2(D,
                               regu,
                               diag_Q, #diag_Q
                               K, #K
                               K_fact, #K_fact
                               false,
                               diagind_K #diagind_K
                               )
end

function convertpad(::Type{<:PreallocatedData{T}}, pad :: PreallocatedData_K2{T_old}, T0 :: DataType) where {T<:Real, T_old<:Real}
    pad = PreallocatedData_K2(convert(Array{T}, pad.D),
                              convert(Regularization{T}, pad.regu),
                              convert(SparseVector{T,Int}, pad.diag_Q),
                              convert(SparseMatrixCSC{T,Int}, pad.K),
                              convertldl(T, pad.K_fact),
                              pad.fact_fail,
                              pad.diagind_K
                              )

    if pad.regu.regul == :classic
        if T == Float64 && T0 == Float64
            pad.regu.ρ_min, pad.regu.δ_min = T(sqrt(eps())*1e-5), T(sqrt(eps())*1e0)
        else
            pad.regu.ρ_min, pad.regu.δ_min = T(sqrt(eps(T))*1e1), T(sqrt(eps(T))*1e1)
        end
    elseif pad.regu.regul == :dynamic
        pad.regu.ρ, pad.regu.δ = -T(eps(T)^(3/4)), T(eps(T)^(0.45))
        pad.K_fact.r1, pad.K_fact.r2 = pad.regu.ρ, pad.regu.δ
    end

    return pad
end

# function used to solve problems
# solver LDLFactorization
function solver!(pad :: PreallocatedData_K2{T}, dda :: DescentDirectionAllocs{T}, pt :: Point{T}, itd :: IterData{T},
                 fd :: Abstract_QM_FloatData{T}, id :: QM_IntData, res :: Residuals{T}, cnts :: Counters, T0 :: DataType,
                 step :: Symbol) where {T<:Real}
    # erase dda.Δxy_aff only for affine predictor step with PC method
    step == :aff ? ldiv!(pad.K_fact, dda.Δxy_aff) : ldiv!(pad.K_fact, itd.Δxy)
    return 0
end

function update_pad!(pad :: PreallocatedData_K2{T}, dda :: DescentDirectionAllocs{T}, pt :: Point{T}, itd :: IterData{T},
                     fd :: Abstract_QM_FloatData{T}, id :: QM_IntData, res :: Residuals{T}, cnts :: Counters, T0 :: DataType) where {T<:Real}

    if pad.regu.regul == :classic && cnts.k != 0 # update ρ and δ values, check K diag magnitude
        out = update_regu_diagK2!(pad.regu, pad.K.nzval, pad.diagind_K, id.nvar, itd.pdd,
                                    itd.l_pdd, itd.mean_pdd, cnts, T, T0)
        out == 1 && return out
    end

    out = factorize_K2!(pad.K, pad.K_fact, pad.D, pad.diag_Q, pad.diagind_K, pad.regu,
                        pt.s_l, pt.s_u, itd.x_m_lvar, itd.uvar_m_x, id.ilow, id.iupp,
                        id.ncon, id.nvar, cnts, itd.qp, T, T0) # update D and factorize K

    if out == 1
        pad.fact_fail = true
        return out
    end

    return 0
end

# Init functions for the K2 system
function fill_K2!(K_colptr, K_rowval, K_nzval, D, Q_colptr, Q_rowval, Q_nzval,
                  AT_colptr, AT_rowval, AT_nzval, diag_Q_nzind, δ, ncon, nvar, regul)

    added_coeffs_diag = 0 # we add coefficients that do not appear in Q in position i,i if Q[i,i] = 0
    n_nz = length(diag_Q_nzind)
    c_nz = n_nz > 0 ? 1 : 0
    @inbounds for j=1:nvar  # Q coeffs, tmp diag coefs.
        K_colptr[j+1] = Q_colptr[j+1] + added_coeffs_diag
        for k=Q_colptr[j]:(Q_colptr[j+1]-2)
            nz_idx = k + added_coeffs_diag
            K_rowval[nz_idx] = Q_rowval[k]
            K_nzval[nz_idx] = -Q_nzval[k]
        end
        if c_nz == 0 || c_nz > n_nz || diag_Q_nzind[c_nz] != j
            added_coeffs_diag += 1
            K_colptr[j+1] += 1
            nz_idx = K_colptr[j+1] - 1
            K_rowval[nz_idx] = j
            K_nzval[nz_idx] = D[j]
        else
            if c_nz != 0
                c_nz += 1
            end
            nz_idx = K_colptr[j+1] - 1
            K_rowval[nz_idx] = j
            K_nzval[nz_idx] = D[j] - Q_nzval[Q_colptr[j+1]-1]
        end
    end

    countsum = K_colptr[nvar+1] # current value of K_colptr[Q.n+j+1]
    nnz_top_left = countsum # number of coefficients + 1 already added
    @inbounds for j=1:ncon
        countsum += AT_colptr[j+1] - AT_colptr[j]
        if regul == :classic
            countsum += 1
        end
        K_colptr[nvar+j+1] = countsum
        for k=AT_colptr[j]:(AT_colptr[j+1]-1)
            nz_idx = regul == :classic ? k + nnz_top_left + j - 2 : k + nnz_top_left - 1
            K_rowval[nz_idx] = AT_rowval[k]
            K_nzval[nz_idx] = AT_nzval[k]
        end
        if regul == :classic
            nz_idx = K_colptr[nvar+j+1] - 1
            K_rowval[nz_idx] = nvar + j
            K_nzval[nz_idx] = δ
        end
    end
end

function create_K2(id, D, Q, AT, diag_Q, regu)
    # for classic regul only
    n_nz = length(D) - length(diag_Q.nzind) + length(AT.nzval) + length(Q.nzval)
    T = eltype(D)
    if regu.regul == :classic
        n_nz += id.ncon
    end
    K_colptr = Vector{Int}(undef, id.ncon+id.nvar+1)
    K_colptr[1] = 1
    K_rowval = Vector{Int}(undef, n_nz)
    K_nzval = Vector{T}(undef, n_nz)
    # [-Q -D    AT]
    # [0        δI]

    fill_K2!(K_colptr, K_rowval, K_nzval, D, Q.colptr, Q.rowval, Q.nzval,
             AT.colptr, AT.rowval, AT.nzval, diag_Q.nzind, regu.δ, id.ncon, id.nvar, regu.regul)

    return SparseMatrixCSC(id.ncon+id.nvar, id.ncon+id.nvar,
                           K_colptr, K_rowval, K_nzval)
end

# iteration functions for the K2 system
function factorize_K2!(K, K_fact, D, diag_Q , diagind_K, regu, s_l, s_u, x_m_lvar, uvar_m_x,
                       ilow, iupp, ncon, nvar, cnts, qp, T, T0)

    if regu.regul == :classic
        D .= -regu.ρ
        K.nzval[view(diagind_K, nvar+1:ncon+nvar)] .= regu.δ
    else
        D .= zero(T)
    end
    D[ilow] .-= s_l ./ x_m_lvar
    D[iupp] .-= s_u ./ uvar_m_x
    D[diag_Q.nzind] .-= diag_Q.nzval
    K.nzval[view(diagind_K,1:nvar)] = D

    if regu.regul == :dynamic
        Amax = @views norm(K.nzval[diagind_K], Inf)
        if Amax > T(1e6) / K_fact.r2 && cnts.c_pdd < 8
            if T == Float32
                return one(Int) # update to Float64
            elseif qp || cnts.c_pdd < 4
                cnts.c_pdd += 1
                regu.δ /= 10
                K_fact.r2 = regu.δ
            end
        end
        K_fact.tol = Amax * T(eps(T))
        ldl_factorize!(Symmetric(K, :U), K_fact)

    elseif regu.regul == :classic
        ldl_factorize!(Symmetric(K, :U), K_fact)
        while !factorized(K_fact)
            out = update_regu_trycatch!(regu, cnts, T, T0)
            out == 1 && return out
            cnts.c_catch += 1
            cnts.c_catch >= 4 && return 1
            D .= -regu.ρ
            D[ilow] .-= s_l ./ x_m_lvar
            D[iupp] .-= s_u ./ uvar_m_x
            D[diag_Q.nzind] .-= diag_Q.nzval
            K.nzval[view(diagind_K,1:nvar)] = D
            K.nzval[view(diagind_K, nvar+1:ncon+nvar)] .= regu.δ
            ldl_factorize!(Symmetric(K, :U), K_fact)
        end

    else # no Regularization
        ldl_factorize!(Symmetric(K, :U), K_fact)
    end

    return 0 # factorization succeeded
end