text:
module RandomSeparations

using Statistics, Random, StatsBase

export run

struct Options
    n_subfeat::Int
    n_thresholds::Int
    zero_score::Float64
    pos_subsample::Float64
    ignore_samples::Set{Int}
end

struct Split
    n_separated::Int
    is_left::Bool
    feature::Int
    threshold::Float32
end

Split() = Split(0, true, 0, 0.0f0)

function make_split(V, labels, n_pos_samples, n_unl_samples, feature, threshold)
    n_left_pos, n_left_unl = 0, 0
    for (vv, ll) in zip(V, labels)
        if vv < threshold
            if ll
                n_left_pos += 1
            else
                n_left_unl += 1
            end
        end
    end
    n_right_pos = n_pos_samples - n_left_pos
    n_right_unl = n_unl_samples - n_left_unl

    if n_left_pos == n_pos_samples
        return Split(n_right_unl, true, feature, threshold)
    end

    if n_right_pos == n_pos_samples
        return Split(n_left_unl, false, feature, threshold)
    end

    return Split()
end

function score(X, Y, opt, seed)
    rng = MersenneTwister(seed)

    scores = Int[]
    iterations = zeros(Int, length(Y))

    features = collect(1:size(X, 2))

    pos_samples = shuffle!(rng, [ss for ss in 1:length(Y) if (!(ss in opt.ignore_samples) && Y[ss])])
    pos_subsample_set = Set(pos_samples[1:round(Int, length(pos_samples) * opt.pos_subsample)])

    samples = [ss for ss in 1:length(Y) if (!(ss in opt.ignore_samples) && (Y[ss] ? (ss in pos_subsample_set) : true))]

    # iterations
    while true

        labels = Y[samples]
        n_pos_samples = sum(labels)
        n_unl_samples = length(labels) - n_pos_samples

        if n_pos_samples == 0 || n_unl_samples == 0
            break
        end

        best_split = Split()
        first_usable_feat = 1
        last_nonconst_feat = length(features)
        V = Vector{Float32}(undef, length(samples))

        mtry = 1
        while (mtry <= opt.n_subfeat) && (first_usable_feat <= last_nonconst_feat)

            ff = rand(rng, first_usable_feat:last_nonconst_feat)
            feature = features[ff]
            features[ff], features[first_usable_feat] = features[first_usable_feat], features[ff]
            first_usable_feat += 1

            minv, maxv = Inf32, -Inf32
            for (ii, ss) in enumerate(samples)
                V[ii] = X[ss, feature]
                if minv > V[ii]
                    minv = V[ii]
                end
                if maxv < V[ii]
                    maxv = V[ii]
                end
            end

            if minv == maxv
                first_usable_feat -= 1
                features[first_usable_feat], features[last_nonconst_feat] = features[last_nonconst_feat], features[first_usable_feat]
                last_nonconst_feat -= 1
                continue
            end

            for _ in 1:opt.n_thresholds
                threshold = minv + (maxv - minv) * rand(rng, Float32)
                split = make_split(V, labels, n_pos_samples, n_unl_samples, feature, threshold)
                if split.n_separated > best_split.n_separated
                    best_split = split
                end
            end

            mtry += 1
        end

        push!(scores, best_split.n_separated)
        interation = length(scores)

        if best_split.n_separated == 0
            for (ss, ll) in zip(samples, labels)
                if !ll
                    iterations[ss] = interation
                end
            end
            break
        end

        nn = 1
        new_samples = Vector{Int}(undef, length(samples) - best_split.n_separated)
        for ss in samples
            if best_split.is_left == (X[ss, best_split.feature] < best_split.threshold)
                new_samples[nn] = ss
                nn += 1
            else
                iterations[ss] = interation
            end
        end

        samples = new_samples
        features = features[1:last_nonconst_feat]
    end

    return [ii > 0 ? (scores[ii] > 0 ? 1.0 / scores[ii] : opt.zero_score) : 0 for ii in iterations]

end

function run(
            X               ::Matrix{Float32},
            Y               ::AbstractArray{Bool};
            n_repetitions   ::Int=100,
            n_subfeat       ::Int=0,
            n_thresholds    ::Int=8,
            zero_score      ::Float64=1.0,
            pos_subsample   ::Float64=1.0, #0.67
            ignore_samples  ::Set{Int}=Set{Int}(),
            seed            ::Int=1234)

    @assert any(Y)
    @assert !all(Y)
    @assert size(X, 1) == length(Y)
    @assert zero_score >= 1.0

    if !(1 <= n_subfeat <= size(X, 2))
        n_subfeat = round(Int, sqrt(size(X, 2)))
    end

    opt = Options(n_subfeat, n_thresholds, zero_score, pos_subsample, ignore_samples)

    repetitions = Vector{Vector{Float64}}(undef, n_repetitions)
    seeds = abs.(rand(MersenneTwister(seed), Int, n_repetitions))
    Threads.@threads for rr in 1:n_repetitions
        repetitions[rr] = score(X, Y, opt, seeds[rr])
    end

    # return reduce((x, y) -> min.(x, y), repetitions)
    return reduce((x, y) -> (x .+= y), repetitions) / n_repetitions

end

end #module


text:
ts = Dates.today()

println("Today's date is:")
println("\t$ts")
println("\t", Dates.format(ts, "E, U dd, yyyy"))

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Tests based on the Raster API tutorial found at http://gdal.org/gdal_tutorial.html

# Opening the File
dataset = GDAL.open("data/utmsmall.tif", GDAL.GA_ReadOnly) # file from GDAL SVN

# Getting Dataset Information
driver = GDAL.getdriverbyname("GTiff")
@test GDAL.getdrivershortname(driver) == "GTiff"
@test GDAL.getdriverlongname(driver) == "GeoTIFF"

@test GDAL.getrasterxsize(dataset) == 100
@test GDAL.getrasterysize(dataset) == 100
@test GDAL.getrastercount(dataset) == 1
nad27_prefix = "PROJCS[\"NAD27 / UTM zone 11N\",GEOGCS[\"NAD27\",DATUM[\"North_American_Datum_1927\","
@test startswith(GDAL.getprojectionref(dataset), nad27_prefix)

geotransform = zeros(6)
GDAL.getgeotransform(dataset, geotransform)
@test geotransform == [440720.0,60.0,0.0,3.75132e6,0.0,-60.0]

# Fetching a Raster Band
band = GDAL.getrasterband(dataset, 1)

blockxsize, blockysize = Ref(Cint(-1)), Ref(Cint(-1))
GDAL.getblocksize(band, blockxsize, blockysize)
@test blockxsize[] == 100
@test blockysize[] == 81

@test GDAL.getdatatypename(GDAL.getrasterdatatype(band)) == "Byte"
@test GDAL.getcolorinterpretationname(GDAL.getrastercolorinterpretation(band)) == "Gray"

gotmin, gotmax = Ref(Cint(-1)), Ref(Cint(-1))
@test GDAL.getrasterminimum(band, gotmin) == 0.0
@test GDAL.getrastermaximum(band, gotmax) == 255.0
@test gotmin[] == gotmax[] == 0

@test GDAL.getoverviewcount(band) == 0
# note that if there is no raster color table, no error (C_None) is set by GDAL
# so there is no GDALError thrown here
@test GDAL.getrastercolortable(band) == C_NULL

# Reading Raster Data
xsize = GDAL.getrasterbandxsize(band)
@test xsize == 100
scanline = Array{Float32}(xsize)
GDAL.rasterio(band, GDAL.GF_Read, 0, 0, xsize, 1,
              scanline, xsize, 1, GDAL.GDT_Float32, 0, 0)
@test scanline == Float32[107.0f0,123.0f0,132.0f0,115.0f0,132.0f0,132.0f0,140.0f0,132.0f0,132.0f0,132.0f0,107.0f0,132.0f0,107.0f0,132.0f0,132.0f0,107.0f0,123.0f0,115.0f0,156.0f0,148.0f0,107.0f0,132.0f0,107.0f0,115.0f0,99.0f0,123.0f0,99.0f0,74.0f0,115.0f0,82.0f0,115.0f0,115.0f0,107.0f0,123.0f0,123.0f0,99.0f0,123.0f0,123.0f0,115.0f0,115.0f0,107.0f0,90.0f0,99.0f0,107.0f0,107.0f0,99.0f0,123.0f0,107.0f0,140.0f0,123.0f0,123.0f0,115.0f0,99.0f0,132.0f0,123.0f0,115.0f0,115.0f0,123.0f0,132.0f0,115.0f0,123.0f0,132.0f0,214.0f0,156.0f0,165.0f0,148.0f0,115.0f0,148.0f0,156.0f0,148.0f0,140.0f0,165.0f0,156.0f0,197.0f0,156.0f0,197.0f0,140.0f0,173.0f0,156.0f0,165.0f0,148.0f0,156.0f0,206.0f0,214.0f0,181.0f0,206.0f0,173.0f0,222.0f0,206.0f0,255.0f0,214.0f0,173.0f0,214.0f0,255.0f0,214.0f0,247.0f0,255.0f0,230.0f0,206.0f0,197.0f0]

# Closing the Dataset
GDAL.close(dataset)

# Techniques for Creating Files
@test GDAL.getmetadataitem(driver, "DCAP_CREATE", "") == "YES"
@test GDAL.getmetadataitem(driver, "DCAP_CREATECOPY", "") == "YES"

# Using CreateCopy()
ds_src = GDAL.open("data/utmsmall.tif", GDAL.GA_ReadOnly)

progressfunc = convert(Ptr{GDAL.GDALProgressFunc}, C_NULL)

ds_dst = GDAL.createcopy(driver, "tmp/utmsmall.tif", ds_src, 0,
                         C_NULL, progressfunc, C_NULL)
GDAL.close(ds_dst)
GDAL.close(ds_src)

rm("tmp/utmsmall.tif")

# rest of the tutorial is skipped, perhaps we should wrap "cpl_string.h"
# to be able to build papszOptions, or write it in julia
# #include "cpl_string.h"
# ...
#     char **papszOptions = NULL;
#
#     papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
#     papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
#     hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE,
#                              papszOptions, GDALTermProgres, NULL );
#     /* Once we're done, close properly the dataset */
#     if( hDstDS != NULL )
#         GDALClose( hDstDS );
#     CSLDestroy( papszOptions );


text:
using JSON, HTTP
include("es.jl")

const BlockNames = Dict(
    :MustQueryBlock => "query",
    :HighlightBlock => "highlight",
    :SourceBlock => "_source",
    :FromBlock => "from",
    :SizeBlock => "size",
    :SortBlock => "sort"
)

abstract type FilterBody end

struct RangeBody <: FilterBody
    field::String
    range::Dict{String,Any}
end
JSON.lower(rb::RangeBody) = Dict("range" => Dict(rb.field => rb.range))

struct MultiMatchBody <: FilterBody
    fields::Vector{String}
    query::String
end
JSON.lower(mm::MultiMatchBody) =
    Dict("multi_match" => Dict("fields" => mm.fields, "query" => mm.query))

struct TermsBody <: FilterBody
    field::String
    values::Vector
    TermsBody(field, values::Vector) = new(field, values)
    TermsBody(field, value::Union{String, Number}) = new(field, [value])
end
JSON.lower(tb::TermsBody) = Dict("terms"=>Dict(tb.field=>tb.values))

abstract type Block end

struct MustQueryBlock <: Block
    bodys::Vector{FilterBody}
    MustQueryBlock(bodys::Vector{FilterBody}) = new(bodys)
    MustQueryBlock(body::FilterBody) = new([body])
end
JSON.lower(mq::MustQueryBlock) = Dict("bool" => Dict("must" => mq.bodys))

struct HighlightBlock <: Block
    fields::Vector{String}
end
JSON.lower(hb::HighlightBlock) =
    Dict("fields" => Dict(field => Dict() for field in hb.fields))

struct SourceBlock <: Block
    fields::Vector{String}
end
JSON.lower(sb::SourceBlock) = sb.fields

struct FromBlock <: Block
    from::Integer
end
JSON.lower(fb::FromBlock) = fb.from

struct SizeBlock <: Block
    size::Integer
end
JSON.lower(sb::SizeBlock) = sb.size

struct SortBody
    field::String
    order::String
    SortBody(field, ascending::Bool = true) =
        new(field, ascending ? "asc" : "desc")
end
JSON.lower(sb::SortBody) = Dict(sb.field => Dict("order" => sb.order))
struct SortBlock <: Block
    bodys::Vector{Union{SortBody,String}}
end
JSON.lower(sb::SortBlock) = sb.bodys

struct Query
    blocks::Vector{Block}
    Query(blocks::Vector{Block}) = new(blocks)
    Query(block::Block) = new([blocks])
end
JSON.lower(q::Query) = begin
    Dict(
        get(BlockNames, block |> typeof |> nameof |> Symbol, "error") =>
            block for block in q.blocks
    )
end

function search(es::ES, index::String, query::String)
    res = HTTP.get(
        "$(es.url)/$index/_search",
        headers = Dict("content-type" => "application/json"),
        body = query,
    )
    String(res.body) |> JSON.parse
end

function search(es::ES, index::String, query::Union{Query,Dict})
    search(es, index, JSON.json(query))
end

function search(query_func::Function; es::ES, index::String)
    search(es, index, query_func())
end


text:
# This file is part of Kpax3. License is MIT.

function test_state_list_constructor()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  x = Kpax3.AminoAcidData(settings)

  (m, n) = size(x.data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(x.data, settings.γ, settings.r)

  slist = Kpax3.AminoAcidStateList(x.data, [3; 3; 3; 1; 1; 2], priorR, priorC, settings)

  @test isa(slist.state, Vector{Kpax3.AminoAcidState})
  @test isa(slist.logpp, Vector{Float64})
  @test isa(slist.rank, Vector{Int})
  @test all(slist.logpp .< 0.0)
  @test length(unique(slist.logpp)) > 1
  @test all(0 .< slist.rank .< settings.popsize + 1)
  @test length(unique(slist.rank)) == settings.popsize
  @test all(diff(slist.logpp[slist.rank]) .<= 0)

  @test slist.state[1].R == [3; 3; 3; 1; 1; 2]

  for i in 1:settings.popsize
    t = Kpax3.AminoAcidState(x.data, slist.state[i].R, priorR, priorC, settings)

    l = t.cl[1:t.k]

    @test slist.state[i].R == t.R
    @test slist.state[i].C == t.C
    @test slist.state[i].emptycluster == t.emptycluster
    @test slist.state[i].cl == t.cl
    @test slist.state[i].k == t.k
    @test slist.state[i].v == t.v
    @test slist.state[i].n1s == t.n1s
    for g in l
      @test slist.state[i].unit[g][1:slist.state[i].v[g]]==t.unit[g][1:t.v[g]]
    end
    @test slist.state[i].logpR == t.logpR
    @test slist.state[i].logpC == t.logpC
    @test slist.state[i].loglik == t.loglik
    @test slist.state[i].logpp == t.logpp

    @test isapprox(slist.logpp[i], slist.state[i].logpp, atol=ε)
  end

  state = Kpax3.AminoAcidState(x.data, [3; 3; 3; 1; 1; 2], priorR, priorC, settings)
  slist = Kpax3.AminoAcidStateList(settings.popsize, state)

  @test isa(slist.state, Vector{Kpax3.AminoAcidState})
  @test isa(slist.logpp, Vector{Float64})
  @test isa(slist.rank, Vector{Int})
  @test all(slist.logpp .< 0.0)
  @test length(unique(slist.logpp)) == 1
  @test all(0 .< slist.rank .< settings.popsize + 1)
  @test length(unique(slist.rank)) == settings.popsize
  @test all(diff(slist.logpp[slist.rank]) .<= 0)
  @test slist.rank == Int[i for i in 1:settings.popsize]

  for i in 1:settings.popsize
    @test slist.state[i].R == state.R
    @test slist.state[i].C == state.C
    @test slist.state[i].emptycluster == state.emptycluster
    @test slist.state[i].cl == state.cl
    @test slist.state[i].k == state.k
    @test slist.state[i].v == state.v
    @test slist.state[i].n1s == state.n1s
    for g in 1:3
      @test (slist.state[i].unit[g][1:slist.state[i].v[g]] ==
             state.unit[g][1:state.v[g]])
    end
    @test slist.state[i].logpR == state.logpR
    @test slist.state[i].logpC == state.logpC
    @test slist.state[i].loglik == state.loglik
    @test slist.state[i].logpp == state.logpp
    @test slist.logpp[i] == state.logpp
  end

  nothing
end

test_state_list_constructor()

function test_state_list_copy_basic()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  x = Kpax3.AminoAcidData(settings)

  (m, n) = size(x.data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(x.data, settings.γ, settings.r)

  s1 = Kpax3.AminoAcidStateList(x.data, [3; 3; 3; 1; 1; 2], priorR, priorC, settings)
  s2 = Kpax3.AminoAcidStateList(x.data, [1; 1; 1; 1; 1; 1], priorR, priorC, settings)

  Kpax3.copystatelist!(s2, s1, settings.popsize)

  for i in 1:settings.popsize
    j = s1.rank[i]

    r = 1:s1.state[j].k
    l = s1.state[j].cl[r]

    @test s2.state[i].R == s1.state[j].R
    @test s2.state[i].C[l, :] == s1.state[j].C[l, :]
    @test s2.state[i].emptycluster[r] == s1.state[j].emptycluster[r]
    @test s2.state[i].cl[r] == s1.state[j].cl[r]
    @test s2.state[i].k == s1.state[j].k
    @test s2.state[i].v[l] == s1.state[j].v[l]
    @test s2.state[i].n1s[l, :] == s1.state[j].n1s[l, :]
    for g in l
      @test (s2.state[i].unit[g][1:s2.state[i].v[g]] ==
             s1.state[j].unit[g][1:s1.state[j].v[g]])
    end
    @test s2.state[i].logpR == s1.state[j].logpR
    @test s2.state[i].logpC == s1.state[j].logpC
    @test s2.state[i].loglik == s1.state[j].loglik
    @test s2.state[i].logpp == s1.state[j].logpp

    @test s2.logpp[i] == s1.logpp[j]

    @test s2.rank[i] == i
  end

  nothing
end

test_state_list_copy_basic()


text:
const DUMMY_PRECISION = 1e-12

Base.isfinite(a::AbstractVec) = all(isfinite, a)
Base.isinf(a::AbstractVec) = any(isinf, a)
Base.isnan(a::AbstractVec) = any(isnan, a)

function StaticArrays.similar_type(::Type{V}, ::Type{F}, size::Size{N}) where {V<:AbstractVec, F<:AbstractFloat, N}
    if size == Size(V) && eltype(V) == F
        return V
    else # convert back to SArray
        return SArray{Tuple{N...},F,length(size),prod(size)}
    end
end

Base.abs(a::AbstractVec) = error("abs(v::AbstractVec) has been removed. Use norm(v) to get the norm; abs.(v) to get the element-wise absolute value.")


text:
using Test

using DiffieHellman
using Sockets
using CryptoGroups
using CryptoSignatures
using Random
using Pkg.TOML

function rngint(len::Integer)
    max_n = ( BigInt(1) << len ) - 1
    if len > 2
        min_n = BigInt(1) << (len - 1)
        return rand(min_n:max_n)
    end
    return rand(1:max_n)
end

hash(x::AbstractString) = BigInt(Base.hash(x))

function chash(envelope1::Vector{UInt8},envelope2::Vector{UInt8},key::BigInt)
    str = "$(String(copy(envelope1))) $(String(copy(envelope2))) $key"
    inthash = hash(str)
    strhash = string(inthash,base=16)
    return Vector{UInt8}(strhash)
end

id(s) = hash("$(s.pubkey)")

function wrap(value::BigInt,signer::Signer)
    signature = DSASignature(hash("$value"),signer)
    signaturedict = Dict(signature)
    dict = Dict("value"=>string(value,base=16),"signature"=>signaturedict)
    io = IOBuffer()
    TOML.print(io,dict)
    return take!(io)
end

function unwrap(envelope::Vector{UInt8})
    dict = TOML.parse(String(copy(envelope)))
    value = parse(BigInt,dict["value"],base=16)
    signature = DSASignature{BigInt}(dict["signature"])
    @assert verify(signature,G) && signature.hash==hash("$value")
    return value, id(signature)
end

G = CryptoGroups.Scep256k1Group()
#G = CryptoGroups.MODP160Group()

# master
server = Signer(G)
serversign(data) = DSASignature(hash(data),server)
serverid = id(server)

# slave
slave = Signer(G)
slavesign(data) = DSASignature(hash(data),slave)
slaveid = id(slave)

value = BigInt(234234324)
value2, signerid2 = unwrap(wrap(value,server))
@test value==value2
@test id(server)==signerid2

### Let's assume that maintainer had contacted the server

@sync begin
    server = listen(2000)
    @async global serversocket = accept(server)
    global slavesocket = connect(2000)
end

### Let's say that we want slave to connect only to a true master. Thus he owns a certificate and from master had obtained a valid master id.
keyserver = @async diffiehellman(serversocket,value->wrap(value,server),unwrap,G,chash,rngint(100))
keyslave = @async diffiehellman(slavesocket,value->wrap(value,slave),unwrap,G,chash,rngint(100)) # x==serverid

keyserv,idserv = fetch(keyserver)
keyslav,idslav = fetch(keyslave)

@test idserv==slaveid
@test idslav==serverid
@test keyserv==keyslav


text:
a= Tokensgraphql("""
#
query {
  Region(name:"The North") {
      NobleHouse(name:"Stark"){
        castle{
          name
        }
        members{
          name
          alias
      }
    }
  }
}
""")

s=""
next = iterate(a)
while next !== nothing
  global next
  global s
    (i, state) = next
   s*= " $(i.kind) $(i.val) ($(i.startpos) , $(i.endpos) ) \n"
   next = iterate(a, state)
end
@test s == " NAME query ((2, 1) , (2, 5) ) \n LBRACE { ((2, 7) , (2, 7) ) \n NAME Region ((3, 3) , (3, 8) ) \n LPAREN ( ((3, 9) , (3, 9) ) \n NAME name ((3, 10) , (3, 13) ) \n COLON : ((3, 14) , (3, 14) ) \n STRING The North ((3, 15) , (3, 25) ) \n RPAREN ) ((3, 26) , (3, 26) ) \n LBRACE { ((3, 28) , (3, 28) ) \n NAME NobleHouse ((4, 7) , (4, 16) ) \n LPAREN ( ((4, 17) , (4, 17) ) \n NAME name ((4, 18) , (4, 21) ) \n COLON : ((4, 22) , (4, 22) ) \n STRING Stark ((4, 23) , (4, 29) ) \n RPAREN ) ((4, 30) , (4, 30) ) \n LBRACE { ((4, 31) , (4, 31) ) \n NAME castle ((5, 9) , (5, 14) ) \n LBRACE { ((5, 15) , (5, 15) ) \n NAME name ((6, 11) , (6, 14) ) \n RBRACE } ((7, 9) , (7, 9) ) \n NAME members ((8, 9) , (8, 15) ) \n LBRACE { ((8, 16) , (8, 16) ) \n NAME name ((9, 11) , (9, 14) ) \n NAME alias ((10, 11) , (10, 15) ) \n RBRACE } ((11, 7) , (11, 7) ) \n RBRACE } ((12, 5) , (12, 5) ) \n RBRACE } ((13, 3) , (13, 3) ) \n RBRACE } ((14, 1) , (14, 1) ) \n ENDMARKER eof ((15, 1) , (15, 1) ) \n"

b= Tokensgraphql(
"""
{
  leftComparison: hero(episode: EMPIRE) {
    ...comparisonFields
  }
  rightComparison: hero(episode: JEDI) {
    ...comparisonFields
  }
}

fragment comparisonFields on Character {
  name
  appearsIn
  friends {
    name
  }
}
""")

s=""
next = iterate(b)
while next !== nothing
  global next
  global s
    (i, state) = next
    s*= " $(i.kind) $(i.val) ($(i.startpos) , $(i.endpos) ) \n"
   next = iterate(b, state)
end
@test s == " LBRACE { ((1, 1) , (1, 1) ) \n NAME leftComparison ((2, 3) , (2, 16) ) \n COLON : ((2, 17) , (2, 17) ) \n NAME hero ((2, 19) , (2, 22) ) \n LPAREN ( ((2, 23) , (2, 23) ) \n NAME episode ((2, 24) , (2, 30) ) \n COLON : ((2, 31) , (2, 31) ) \n NAME EMPIRE ((2, 33) , (2, 38) ) \n RPAREN ) ((2, 39) , (2, 39) ) \n LBRACE { ((2, 41) , (2, 41) ) \n SPREAD ... ((3, 5) , (3, 5) ) \n NAME comparisonFields ((3, 6) , (3, 21) ) \n RBRACE } ((4, 3) , (4, 3) ) \n NAME rightComparison ((5, 3) , (5, 17) ) \n COLON : ((5, 18) , (5, 18) ) \n NAME hero ((5, 20) , (5, 23) ) \n LPAREN ( ((5, 24) , (5, 24) ) \n NAME episode ((5, 25) , (5, 31) ) \n COLON : ((5, 32) , (5, 32) ) \n NAME JEDI ((5, 34) , (5, 37) ) \n RPAREN ) ((5, 38) , (5, 38) ) \n LBRACE { ((5, 40) , (5, 40) ) \n SPREAD ... ((6, 5) , (6, 5) ) \n NAME comparisonFields ((6, 6) , (6, 21) ) \n RBRACE } ((7, 3) , (7, 3) ) \n RBRACE } ((8, 1) , (8, 1) ) \n NAME fragment ((10, 1) , (10, 8) ) \n NAME comparisonFields ((10, 10) , (10, 25) ) \n NAME on ((10, 27) , (10, 28) ) \n NAME Character ((10, 30) , (10, 38) ) \n LBRACE { ((10, 40) , (10, 40) ) \n NAME name ((11, 3) , (11, 6) ) \n NAME appearsIn ((12, 3) , (12, 11) ) \n NAME friends ((13, 3) , (13, 9) ) \n LBRACE { ((13, 11) , (13, 11) ) \n NAME name ((14, 5) , (14, 8) ) \n RBRACE } ((15, 3) , (15, 3) ) \n RBRACE } ((16, 1) , (16, 1) ) \n ENDMARKER eof ((17, 1) , (17, 1) ) \n"

b= Tokensgraphql(
"""
mutation {
  sendEmail(message: \"\"\"
    Hello,
      World!

    Yours,
      GraphQL.
  \"\"\")
}
""")

s=""
next = iterate(b)
while next !== nothing
  global next
  global s
    (i, state) = next
    s*= " $(i.kind) $(i.val) ($(i.startpos) , $(i.endpos) ) \n"
   next = iterate(b, state)
end
@test s ==  " NAME mutation ((1, 1) , (1, 8) ) \n LBRACE { ((1, 10) , (1, 10) ) \n NAME sendEmail ((2, 3) , (2, 11) ) \n LPAREN ( ((2, 12) , (2, 12) ) \n NAME message ((2, 13) , (2, 19) ) \n COLON : ((2, 20) , (2, 20) ) \n STRING \n    Hello,\n      World!\n\n    Yours,\n      GraphQL.\n   ((2, 22) , (8, 6) ) \n RPAREN ) ((8, 7) , (8, 7) ) \n RBRACE } ((9, 1) , (9, 1) ) \n ENDMARKER eof ((10, 1) , (10, 1) ) \n"

b= Tokensgraphql(
"""
type Query {
  dog: Dog
}

enum DogCommand { SIT, DOWN, HEEL }

type Dog implements Pet {
  name: String!
  nickname: String
  barkVolume: Int
  doesKnowCommand(dogCommand: DogCommand!): Boolean!
  isHousetrained(atOtherHomes: Boolean): Boolean!
  owner: Human
}

interface Sentient {
  name: String!
}

interface Pet {
  name: String!
}

type Alien implements Sentient {
  name: String!
  homePlanet: String
}

type Human implements Sentient {
  name: String!
}

enum CatCommand { JUMP }

type Cat implements Pet {
  name: String!
  nickname: String
  doesKnowCommand(catCommand: CatCommand!): Boolean!
  meowVolume: Int
}

union CatOrDog = Cat | Dog
union DogOrHuman = Dog | Human
union HumanOrAlien = Human | Alien
""")

s=""
next = iterate(b)
while next !== nothing
  global next
  global s
    (i, state) = next
    s*= " $(i.kind) $(i.val) ($(i.startpos) , $(i.endpos) ) \n"
   next = iterate(b, state)
end
@test s == " NAME type ((1, 1) , (1, 4) ) \n NAME Query ((1, 6) , (1, 10) ) \n LBRACE { ((1, 12) , (1, 12) ) \n NAME dog ((2, 3) , (2, 5) ) \n COLON : ((2, 6) , (2, 6) ) \n NAME Dog ((2, 8) , (2, 10) ) \n RBRACE } ((3, 1) , (3, 1) ) \n NAME enum ((5, 1) , (5, 4) ) \n NAME DogCommand ((5, 6) , (5, 15) ) \n LBRACE { ((5, 17) , (5, 17) ) \n NAME SIT ((5, 19) , (5, 21) ) \n NAME DOWN ((5, 24) , (5, 27) ) \n NAME HEEL ((5, 30) , (5, 33) ) \n RBRACE } ((5, 35) , (5, 35) ) \n NAME type ((7, 1) , (7, 4) ) \n NAME Dog ((7, 6) , (7, 8) ) \n NAME implements ((7, 10) , (7, 19) ) \n NAME Pet ((7, 21) , (7, 23) ) \n LBRACE { ((7, 25) , (7, 25) ) \n NAME name ((8, 3) , (8, 6) ) \n COLON : ((8, 7) , (8, 7) ) \n NAME String ((8, 9) , (8, 14) ) \n BANG ! ((8, 15) , (8, 15) ) \n NAME nickname ((9, 3) , (9, 10) ) \n COLON : ((9, 11) , (9, 11) ) \n NAME String ((9, 13) , (9, 18) ) \n NAME barkVolume ((10, 3) , (10, 12) ) \n COLON : ((10, 13) , (10, 13) ) \n NAME Int ((10, 15) , (10, 17) ) \n NAME doesKnowCommand ((11, 3) , (11, 17) ) \n LPAREN ( ((11, 18) , (11, 18) ) \n NAME dogCommand ((11, 19) , (11, 28) ) \n COLON : ((11, 29) , (11, 29) ) \n NAME DogCommand ((11, 31) , (11, 40) ) \n BANG ! ((11, 41) , (11, 41) ) \n RPAREN ) ((11, 42) , (11, 42) ) \n COLON : ((11, 43) , (11, 43) ) \n NAME Boolean ((11, 45) , (11, 51) ) \n BANG ! ((11, 52) , (11, 52) ) \n NAME isHousetrained ((12, 3) , (12, 16) ) \n LPAREN ( ((12, 17) , (12, 17) ) \n NAME atOtherHomes ((12, 18) , (12, 29) ) \n COLON : ((12, 30) , (12, 30) ) \n NAME Boolean ((12, 32) , (12, 38) ) \n RPAREN ) ((12, 39) , (12, 39) ) \n COLON : ((12, 40) , (12, 40) ) \n NAME Boolean ((12, 42) , (12, 48) ) \n BANG ! ((12, 49) , (12, 49) ) \n NAME owner ((13, 3) , (13, 7) ) \n COLON : ((13, 8) , (13, 8) ) \n NAME Human ((13, 10) , (13, 14) ) \n RBRACE } ((14, 1) , (14, 1) ) \n NAME interface ((16, 1) , (16, 9) ) \n NAME Sentient ((16, 11) , (16, 18) ) \n LBRACE { ((16, 20) , (16, 20) ) \n NAME name ((17, 3) , (17, 6) ) \n COLON : ((17, 7) , (17, 7) ) \n NAME String ((17, 9) , (17, 14) ) \n BANG ! ((17, 15) , (17, 15) ) \n RBRACE } ((18, 1) , (18, 1) ) \n NAME interface ((20, 1) , (20, 9) ) \n NAME Pet ((20, 11) , (20, 13) ) \n LBRACE { ((20, 15) , (20, 15) ) \n NAME name ((21, 3) , (21, 6) ) \n COLON : ((21, 7) , (21, 7) ) \n NAME String ((21, 9) , (21, 14) ) \n BANG ! ((21, 15) , (21, 15) ) \n RBRACE } ((22, 1) , (22, 1) ) \n NAME type ((24, 1) , (24, 4) ) \n NAME Alien ((24, 6) , (24, 10) ) \n NAME implements ((24, 12) , (24, 21) ) \n NAME Sentient ((24, 23) , (24, 30) ) \n LBRACE { ((24, 32) , (24, 32) ) \n NAME name ((25, 3) , (25, 6) ) \n COLON : ((25, 7) , (25, 7) ) \n NAME String ((25, 9) , (25, 14) ) \n BANG ! ((25, 15) , (25, 15) ) \n NAME homePlanet ((26, 3) , (26, 12) ) \n COLON : ((26, 13) , (26, 13) ) \n NAME String ((26, 15) , (26, 20) ) \n RBRACE } ((27, 1) , (27, 1) ) \n NAME type ((29, 1) , (29, 4) ) \n NAME Human ((29, 6) , (29, 10) ) \n NAME implements ((29, 12) , (29, 21) ) \n NAME Sentient ((29, 23) , (29, 30) ) \n LBRACE { ((29, 32) , (29, 32) ) \n NAME name ((30, 3) , (30, 6) ) \n COLON : ((30, 7) , (30, 7) ) \n NAME String ((30, 9) , (30, 14) ) \n BANG ! ((30, 15) , (30, 15) ) \n RBRACE } ((31, 1) , (31, 1) ) \n NAME enum ((33, 1) , (33, 4) ) \n NAME CatCommand ((33, 6) , (33, 15) ) \n LBRACE { ((33, 17) , (33, 17) ) \n NAME JUMP ((33, 19) , (33, 22) ) \n RBRACE } ((33, 24) , (33, 24) ) \n NAME type ((35, 1) , (35, 4) ) \n NAME Cat ((35, 6) , (35, 8) ) \n NAME implements ((35, 10) , (35, 19) ) \n NAME Pet ((35, 21) , (35, 23) ) \n LBRACE { ((35, 25) , (35, 25) ) \n NAME name ((36, 3) , (36, 6) ) \n COLON : ((36, 7) , (36, 7) ) \n NAME String ((36, 9) , (36, 14) ) \n BANG ! ((36, 15) , (36, 15) ) \n NAME nickname ((37, 3) , (37, 10) ) \n COLON : ((37, 11) , (37, 11) ) \n NAME String ((37, 13) , (37, 18) ) \n NAME doesKnowCommand ((38, 3) , (38, 17) ) \n LPAREN ( ((38, 18) , (38, 18) ) \n NAME catCommand ((38, 19) , (38, 28) ) \n COLON : ((38, 29) , (38, 29) ) \n NAME CatCommand ((38, 31) , (38, 40) ) \n BANG ! ((38, 41) , (38, 41) ) \n RPAREN ) ((38, 42) , (38, 42) ) \n COLON : ((38, 43) , (38, 43) ) \n NAME Boolean ((38, 45) , (38, 51) ) \n BANG ! ((38, 52) , (38, 52) ) \n NAME meowVolume ((39, 3) , (39, 12) ) \n COLON : ((39, 13) , (39, 13) ) \n NAME Int ((39, 15) , (39, 17) ) \n RBRACE } ((40, 1) , (40, 1) ) \n NAME union ((42, 1) , (42, 5) ) \n NAME CatOrDog ((42, 7) , (42, 14) ) \n EQUALS = ((42, 16) , (42, 16) ) \n NAME Cat ((42, 18) , (42, 20) ) \n PIPE | ((42, 22) , (42, 22) ) \n NAME Dog ((42, 24) , (42, 26) ) \n NAME union ((43, 1) , (43, 5) ) \n NAME DogOrHuman ((43, 7) , (43, 16) ) \n EQUALS = ((43, 18) , (43, 18) ) \n NAME Dog ((43, 20) , (43, 22) ) \n PIPE | ((43, 24) , (43, 24) ) \n NAME Human ((43, 26) , (43, 30) ) \n NAME union ((44, 1) , (44, 5) ) \n NAME HumanOrAlien ((44, 7) , (44, 18) ) \n EQUALS = ((44, 20) , (44, 20) ) \n NAME Human ((44, 22) , (44, 26) ) \n PIPE | ((44, 28) , (44, 28) ) \n NAME Alien ((44, 30) , (44, 34) ) \n ENDMARKER eof ((45, 1) , (45, 1) ) \n"

b= Tokensgraphql(
"""
{
  user(id: 4) {
    id
    name
    profilePic(width: 100, height: -50)
  }
}
""")

s=""
next = iterate(b)
while next !== nothing
  global next
  global s
    (i, state) = next
    s*= " $(i.kind) $(i.val) ($(i.startpos) , $(i.endpos) ) \n"
   next = iterate(b, state)
end
@test s == " LBRACE { ((1, 1) , (1, 1) ) \n NAME user ((2, 3) , (2, 6) ) \n LPAREN ( ((2, 7) , (2, 7) ) \n NAME id ((2, 8) , (2, 9) ) \n COLON : ((2, 10) , (2, 10) ) \n INT 4 ((2, 12) , (2, 12) ) \n RPAREN ) ((2, 13) , (2, 13) ) \n LBRACE { ((2, 15) , (2, 15) ) \n NAME id ((3, 5) , (3, 6) ) \n NAME name ((4, 5) , (4, 8) ) \n NAME profilePic ((5, 5) , (5, 14) ) \n LPAREN ( ((5, 15) , (5, 15) ) \n NAME width ((5, 16) , (5, 20) ) \n COLON : ((5, 21) , (5, 21) ) \n INT 100 ((5, 23) , (5, 25) ) \n NAME height ((5, 28) , (5, 33) ) \n COLON : ((5, 34) , (5, 34) ) \n INT -50 ((5, 36) , (5, 38) ) \n RPAREN ) ((5, 39) , (5, 39) ) \n RBRACE } ((6, 3) , (6, 3) ) \n RBRACE } ((7, 1) , (7, 1) ) \n ENDMARKER eof ((8, 1) , (8, 1) ) \n"

b= Tokensgraphql(
"""
{
  user(id: 4.0) {
    id
    name
    profilePic(width: 1e50, height: -6.0221413e23)
  }
}
""")

s=""
next = iterate(b)
while next !== nothing
  global next
  global s
    (i, state) = next
    s*= " $(i.kind) $(i.val) ($(i.startpos) , $(i.endpos) ) \n"
   next = iterate(b, state)
end
@test s ==" LBRACE { ((1, 1) , (1, 1) ) \n NAME user ((2, 3) , (2, 6) ) \n LPAREN ( ((2, 7) , (2, 7) ) \n NAME id ((2, 8) , (2, 9) ) \n COLON : ((2, 10) , (2, 10) ) \n FLOAT 4.0 ((2, 12) , (2, 14) ) \n RPAREN ) ((2, 15) , (2, 15) ) \n LBRACE { ((2, 17) , (2, 17) ) \n NAME id ((3, 5) , (3, 6) ) \n NAME name ((4, 5) , (4, 8) ) \n NAME profilePic ((5, 5) , (5, 14) ) \n LPAREN ( ((5, 15) , (5, 15) ) \n NAME width ((5, 16) , (5, 20) ) \n COLON : ((5, 21) , (5, 21) ) \n FLOAT 1.0e50 ((5, 23) , (5, 26) ) \n NAME height ((5, 29) , (5, 34) ) \n COLON : ((5, 35) , (5, 35) ) \n FLOAT -6.0221413e23 ((5, 37) , (5, 49) ) \n RPAREN ) ((5, 50) , (5, 50) ) \n RBRACE } ((6, 3) , (6, 3) ) \n RBRACE } ((7, 1) , (7, 1) ) \n ENDMARKER eof ((8, 1) , (8, 1) ) \n"

try
b= Tokensgraphql(
"""
{
  user(id: 4.) {
    id
    name
    profilePic(width: 1e50, height: -6.0221413e23)
  }
}
""")
catch e
  if e isa Diana.GraphQLError
        @test string(e) ==  "Diana.GraphQLError(\"{\\\"errors\\\":[{\\\"locations\\\": [{\\\"column\\\": 14,\\\"line\\\": 2}],\\\"message\\\": \\\"Syntax Error GraphQL request (2:14) Unexpected character ) \\\"}]}\")"
  else
   	 @test string(e) == ""
  end
end

try
b= Tokensgraphql(
"""
{
  user(id: 4.0) {
    id
    name
    profilePic(width: 1e, height: -6.0221413e23)
  }
}
""")
catch e
  if e isa Diana.GraphQLError
        @test string(e) ==  "Diana.GraphQLError(\"{\\\"errors\\\":[{\\\"locations\\\": [{\\\"column\\\": 25,\\\"line\\\": 5}],\\\"message\\\": \\\"Syntax Error GraphQL request (5:25) Unexpected character , \\\"}]}\")"
  else
   	 @test string(e) == ""
  end
end

try
b= Tokensgraphql(
"""
{
  user(id: 4.0) {
    id
    name
    profilePic(width: 1e50, height: -6.e23)
  }
}
""")
catch e
  if e isa Diana.GraphQLError
        @test string(e) ==  "Diana.GraphQLError(\"{\\\"errors\\\":[{\\\"locations\\\": [{\\\"column\\\": 40,\\\"line\\\": 5}],\\\"message\\\": \\\"Syntax Error GraphQL request (5:40) Unexpected character e \\\"}]}\")"
  else
   	 @test string(e) == ""
  end
end

try
b= Tokensgraphql(
"""
{
  user(id: 4.0) {
    id
    name
    profilePic(width: 1e50, height: -6.0e
  }
}
""")
catch e
  if e isa Diana.GraphQLError
        @test string(e) == "Diana.GraphQLError(\"{\\\"errors\\\":[{\\\"locations\\\": [{\\\"column\\\": 42,\\\"line\\\": 5}],\\\"message\\\": \\\"Syntax Error GraphQL request (5:42) Unexpected character \\\\n \\\"}]}\")"
   else
   	 @test string(e) == ""
  end
end


text:
# Surrogate functions
hinge(x, ϑ::Real = 1) = max(zero(x), 1 + ϑ * x)
quadratic(x, ϑ::Real = 1) = max(zero(x), 1 + ϑ * x)^2

# Objective functions
struct FNRate <: Objective end

Base.show(io::IO, ::FNRate) = print(io, "false-negative rate")

function objective(::FNRate, y, s::AbstractArray{T}, t, surrogate) where {T<:Real}
    inds = zero(T) .+ (y .== 1)
    return sum(surrogate.(t .-  s) .* inds; dims = 2) ./ sum(inds; dims = 2)
end

struct FPRate <: Objective end

Base.show(io::IO, ::FPRate) = print(io, "false-positive rate")

function objective(::FPRate, y, s::AbstractArray{T}, t, surrogate) where {T<:Real}
    inds = zero(T) .+ (y .== 1)
    return sum(surrogate.(s .-  t) .* inds; dims = 2) ./ sum(inds; dims = 2)
end

struct FNFPRate <: Objective
    α::Real
end

function Base.show(io::IO, o::FNFPRate)
    return print(io, "$(o.α)⋅false-negative + $(1-o.α)⋅false-positive rate")
end

function objective(o::FNFPRate, y, s, t, surrogate)
    return o.α .* objective(FNRate(), y, s, t, surrogate) .+
           (1 - o.α) .* objective(FPRate(), y, s, t, surrogate)
end

# Accuracy at Top formulation
struct AccAtTop
    threshold_type::Threshold
    objective_type::Objective
end

function Base.show(io::IO, m::AccAtTop)
    println(io, "Accuracy at the top:")
    println(io, " - threshold: $(m.threshold_type)")
    print(io, " - objective function: $(m.objective_type)")
    return
end

function objective(
    m::AccAtTop,
    y::AbstractMatrix,
    s::AbstractMatrix;
    surrogate = hinge,
    agg = mean,
    update_buffer = true,
)

    if size(y) != size(s)
        throw(DimensionMismatch("dimensions must match: y has dims $(size(y)), s has dims $(size(s))"))
    end

    ts = threshold(m.threshold_type, y, s; update_buffer)
    return agg(objective(m.objective_type, y, s, ts, surrogate))
end

function predict(
    m::AccAtTop,
    y::AbstractMatrix,
    s::AbstractMatrix;
    ts = threshold(m.threshold_type, y, s; update_buffer = false),
)

    if size(y) != size(s)
        throw(DimensionMismatch("dimensions must match: y has dims $(size(y)), s has dims $(size(s))"))
    end

    if length(ts) != size(s, 1)
        throw(DimensionMismatch("length of thresholds $(length(ts)) does not match first dim of scores $(size(s, 1))"))
    end
    return s .>= ts
end

# Basic models
DeepTopPush() = AccAtTop(Maximum(Neg), FNRate())
DeepTopPushK(K) = AccAtTop(Kth(K, Neg), FNRate())
PatMat(τ::Real) = AccAtTop(Quantile(τ, All; rev = true), FNRate())
PatMat(sampler) = AccAtTop(SampledQuantile(sampler, All; rev = true), FNRate())
PatMatNP(τ::Real) = AccAtTop(Quantile(τ, Neg; rev = true), FNRate())
PatMatNP(sampler) = AccAtTop(SampledQuantile(sampler, Neg; rev = true), FNRate())


text:
# Media response object
mutable struct Media <: TwitterType
    media_id::Union{Int, Nothing}
    media_id_string::Union{String, Nothing}
    size::Union{Int, Nothing}
    expires_after_secs::Union{Int, Nothing}
    image_type::Union{String, Nothing}
    w::Union{Int, Nothing}
    h::Union{Int, Nothing}
end

# Methods to convert JSON/Dict to Twitter types
Media(object::Dict) =
    Media(get(object, "media_id", nothing),
            get(object, "media_id_string", nothing),
            get(object, "size", nothing),
            get(object, "expires_after_secs", nothing),
            get(object, "image_type", nothing),
            get(object, "w", nothing),
            get(object, "h", nothing)
    )

Media(response::Vector) = Media[Media(x) for x in response]


text:
# Utility helpers for testing operators
export run_testset, run_proxyshowcheck, @ts

using Rocket
using Test

import Rocket: test_on_source
import Rocket: @ts

const TESTSET_MAXIMUM_TIMEOUT = 5000 # 5 sec
const TESTSET_MAXIMUM_COMPILATION_TIMEOUT = 5000 # 5 sec

# ------------------------------------------------- #
# Helpers for running custom testsets for operators
# ------------------------------------------------- #

function testset_parameters(testset::NamedTuple{V, T}) where { V, T <: Tuple }
    parameters = Dict([
        (:source, nothing),
        (:values, nothing),
        (:wait_for, TESTSET_MAXIMUM_TIMEOUT),
        (:throws, nothing),
        (:source_type, Any)
    ])

    for key in V
        parameters[key] = testset[key]
    end

    return parameters
end

function test(testset, check_timings)
    parameters = testset_parameters(testset)

    # Force JIT compilation of all statements before actual test execution
    # Iffy approach, TODO
    try
        test_on_source(parameters[:source], parameters[:values]; maximum_wait = convert(Float64, TESTSET_MAXIMUM_COMPILATION_TIMEOUT), check_timings = false)
    catch _
    end

    if parameters[:source] === nothing
        error("Missing :source field in the testset")
    elseif parameters[:throws] !== nothing
        @test_throws parameters[:throws] test_on_source(parameters[:source], parameters[:values]; maximum_wait = parameters[:wait_for], check_timings = check_timings)
    else
        @test test_on_source(parameters[:source], parameters[:values]; maximum_wait = convert(Float64, parameters[:wait_for]), check_timings = check_timings)
    end

    if parameters[:source_type] !== Any
        @test subscribable_extract_type(parameters[:source]) === parameters[:source_type]
    end

    return true
end

function run_testset(testsets; check_timings = true)
    for ts in testsets
        @testset begin
            test(ts, check_timings)
        end
    end
end

# --------------------------------------------------- #
# Helpers for checking Base.show methods for operators
# --------------------------------------------------- #

function run_proxyshowcheck(name::String, operator; args::Union{NamedTuple{V, T}, Nothing} = nothing) where { V, T <: Tuple }

    parameters = Dict([
        (:custom_source, never(Any)),
        (:custom_actor,  void(Any)),
        (:check_operator, true),
        (:check_proxy, true),
        (:check_actor, :auto),
        (:check_observable, :auto),
        (:check_subscription, false)
    ])

    if args !== nothing
        for key in V
            parameters[key] = args[key]
        end
    end

    withoperator = parameters[:custom_source] |> operator

    if !(withoperator isa ProxyObservable)
        error("run_proxyshowcheck(...) is a helper function only for proxy-based operators")
    end

    proxytrait = Rocket.as_proxy(typeof(withoperator))

    @testset "Check Base.show for $name operator" begin
        io = IOBuffer()

        if parameters[:check_operator] === true
            show(io, operator)
            printed = String(take!(io))
            if !occursin("$(name)Operator", printed)
                error("Check for Base.show(io, operator::$(typeof(operator)) has failed")
            end
        end

        if parameters[:check_proxy] === true
            show(io, withoperator.proxy)
            printed = String(take!(io))
            if !occursin("$(name)Proxy", printed)
                error("Check for Base.show(io, proxy::$(typeof(withoperator.proxy)) has failed")
            end
        end

        if parameters[:check_actor] === true || (parameters[:check_actor] === :auto && (proxytrait isa Rocket.ValidActorProxy || proxytrait isa Rocket.ValidActorSourceProxy))
            actor = Rocket.actor_proxy!(withoperator.proxy, parameters[:custom_actor])
            show(io, actor)
            printed = String(take!(io))
            if !occursin("$(name)Actor", printed)
                error("Check for Base.show(io, actor::$(typeof(actor)) has failed")
            end
        end

        if parameters[:check_observable] === true || (parameters[:check_observable] === :auto && (proxytrait isa Rocket.ValidSourceProxy || proxytrait isa Rocket.ValidActorSourceProxy))
            source = Rocket.source_proxy!(withoperator.proxy, parameters[:custom_source])
            show(io, source)
            printed = String(take!(io))
            if !occursin("$(name)Observable", printed)
                error("Check for Base.show(io, observable::$(typeof(source)) has failed")
            end
        end

        if parameters[:check_subscription] === true
            subscription = subscribe!(withoperator, parameters[:custom_actor])
            show(io, subscription)
            printed = String(take!(io))
            if !occursin("$(name)Subscription", printed)
                error("Check for Base.show(io, subscription::$(typeof(subscription)) has failed")
            end
            unsubscribe!(subscription)
        else
            subscription = subscribe!(withoperator, parameters[:custom_actor])
            show(io, subscription)
            printed = String(take!(io))
            if occursin("$(name)", printed)
                error("Check for Base.show(io, subscription::$(typeof(subscription)) has failed")
            end
            unsubscribe!(subscription)
        end

        @test true
    end

end


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

@everywhere module TestDecisionObjective

using StochasticPrograms
using Test

function test_SingleDecision_objective(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Min, x)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Min, y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionRef
    @test JuMP.objective_function(sp, 1) == x
    @test x == @inferred JuMP.objective_function(sp, 1, DecisionRef)
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 2, DecisionRef)
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2, JuMP.objective_function(sp))
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 2, DecisionRef)
        @test JuMP.objective_function(sp) == x
        @test x == @inferred JuMP.objective_function(sp, DecisionRef)
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 + y2)
        @test JuMP.isequal_canonical(x1 + y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(x2 + y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
    @first_stage sp = begin
        @decision(sp, x)
        @objective(sp, Max, x)
    end
    @second_stage sp = begin
        @recourse(sp, y)
        @objective(sp, Max, y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionRef
    @test JuMP.objective_function(sp, 1) == x
    @test x == @inferred JuMP.objective_function(sp, 1, DecisionRef)
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 2, DecisionRef)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), x + 0.5*y1 + 0.5*y2)
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionRef
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionRef
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 1) == y2
        @test y1 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test y2 == @inferred JuMP.objective_function(sp, 2, 1, DecisionRef)
        @test JuMP.objective_function(sp) == x
        @test x == @inferred JuMP.objective_function(sp, DecisionRef)
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
        @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 + y2)
        @test JuMP.isequal_canonical(x1 + y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(x2 + y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_DecisionAffExpr_objective(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Min, 2x)
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Min, a*y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionAffExpr{Float64}
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), 2x)
    @test JuMP.isequal_canonical(
        2x, @inferred JuMP.objective_function(sp, 1, DecisionAffExpr{Float64}))
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
    @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        @test JuMP.isequal_canonical(
            2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x + y1 + 2y2)
        @test JuMP.isequal_canonical(2x + y1 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        @test JuMP.isequal_canonical(
            2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x)
        @test JuMP.isequal_canonical(2x,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x1 + 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x2 + 4y2)
        @test JuMP.isequal_canonical(
            2x1 + 2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            2x2 + 4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Max, x + 3x + 1)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Max, x + a*x + 1 + y + 3y + 1)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionAffExpr{Float64}
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), 4x + 1)
    @test JuMP.isequal_canonical(
        4x + 1, @inferred JuMP.objective_function(sp, 1, DecisionAffExpr{Float64}))
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test JuMP.objective_function_type(sp, 2, 1) == DecisionAffExpr{Float64}
    @test JuMP.objective_function_type(sp, 2, 2) == DecisionAffExpr{Float64}
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 3x + 4y1 + 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 5x + 4y2 + 2)
        @test JuMP.isequal_canonical(
            3x + 4y1 + 2, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            5x + 4y2 + 2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 8x + 2y1 + 2y2 + 3)
        @test JuMP.isequal_canonical(8x + 2y1 + 2y2 + 3,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 3x + 4y1 + 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 5x + 4y2 + 2)
        @test JuMP.isequal_canonical(
            3x + 4y1 + 2, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            5x + 4y2 + 2, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 4x + 1)
        @test JuMP.isequal_canonical(4x + 1,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 7x1 + 4y1 + 3)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 9x2 + 4y2 + 3)
        @test JuMP.isequal_canonical(
            7x1 + 4y1 + 3, @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(
            9x2 + 4y2 + 3, @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_DecisionQuadExpr_objective(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Min, x^2 + 2x)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Min, y^2 + a*y + x^2 + 2x)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    @test JuMP.objective_function_type(sp, 1) == DecisionQuadExpr{Float64}
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), x^2 + 2x)
    @test JuMP.isequal_canonical(
        x^2 + 2x, @inferred JuMP.objective_function(sp, 1, DecisionQuadExpr{Float64}))
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test JuMP.objective_function_type(sp, 2, 1) == DecisionQuadExpr{Float64}
    # Structure specific
    x = DecisionRef(sp[1,:x])
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), y1^2 + 2y1 + x^2 + 2x)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), y2^2 + 4y2 + x^2 + 2x)
        @test JuMP.isequal_canonical(
            y1^2 + 2y1 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(
            y2^2 + 4y2 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x^2 + 4x + 0.5*y1^2 + y1 + 0.5*y2^2 + 2y2)
        @test JuMP.isequal_canonical(2x^2 + 4x + 0.5*y1^2 + y1 + 0.5*y2^2 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), y1^2 + 2y1 + x^2 + 2x)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), y2^2 + 4y2 + x^2 + 2x)
        @test JuMP.isequal_canonical(
            y1^2 + 2y1 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(
            y2^2 + 4y2 + x^2 + 2x, @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), x^2 + 2x)
        @test JuMP.isequal_canonical(x^2 + 2x,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x1^2 + 4x1 + y1^2 + 2y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x2^2 + 4x2 + y2^2 + 4y2)
        @test JuMP.isequal_canonical(
            2x1^2 + 4x1 + y1^2 + 2y1, @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(
            2x2^2 + 4x2 + y2^2 + 4y2, @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_decision_objective_coefficient_modifiction(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Min, x)
    end
    @second_stage sp = begin
        @known(sp, x)
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Min, y)
    end
    # First-stage
    x = sp[1,:x]
    JuMP.set_objective_coefficient(sp, x, 1, 4.0)
    x = DecisionRef(x)
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), 4x)
    # Second-stage
    x = sp[1,:x]
    y = sp[2,:y]
    JuMP.set_objective_coefficient(sp, y, 2, 1, 4.0)
    JuMP.set_objective_coefficient(sp, y, 2, 2, 4.0)
    y1 = DecisionRef(y, 1)
    y2 = DecisionRef(y, 2)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 7x + 2y1 + 2y2)
        @test JuMP.isequal_canonical(7x + 2y1 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4y2)
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 4x)
        @test JuMP.isequal_canonical(4x,
                                     @inferred JuMP.objective_function(sp, DecisionAffExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x = DecisionRef(sp[1,:x], 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 4x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 4x + 4y2)
        @test JuMP.isequal_canonical(4x + 4y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionAffExpr{Float64}))
        @test JuMP.isequal_canonical(4x + 4y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionAffExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Min, x^2 + x)
    end
    @second_stage sp = begin
        @known(sp, x)
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Min, x^2 + y^2 + y)
    end
    # First-stage
    x = sp[1,:x]
    JuMP.set_objective_coefficient(sp, x, 1, 4.0)
    x = DecisionRef(x)
    @test JuMP.isequal_canonical(JuMP.objective_function(sp, 1), x^2 + 4x)
    # Second-stage
    x = sp[1,:x]
    y = sp[2,:y]
    JuMP.set_objective_coefficient(sp, y, 2, 1, 4.0)
    JuMP.set_objective_coefficient(sp, y, 2, 2, 4.0)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        x = DecisionRef(x, 2, 1)
        y1 = DecisionRef(y, 1)
        y2 = DecisionRef(y, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4y2)
        x = sp[1,:x]
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(x, 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), 2x^2 + 7x + 0.5*y1^2 + 2y1 + 0.5*y2^2 + 2y2)
        @test JuMP.isequal_canonical(2x^2 + 7x + 0.5*y1^2 + 2y1 + 0.5*y2^2 + 2y2,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        x = DecisionRef(x, 2, 1)
        y1 = DecisionRef(y, 1)
        y2 = DecisionRef(y, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4y2)
        x = sp[1,:x]
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(x, 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x^2 + y1^2 + 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x^2 + y2^2 + 4x + 4y2)
        x = DecisionRef(sp[1,:x])
        @test JuMP.isequal_canonical(JuMP.objective_function(sp), x^2 + 4x)
        @test JuMP.isequal_canonical(x^2 + 4x,
                                     @inferred JuMP.objective_function(sp, DecisionQuadExpr{Float64}))
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x = DecisionRef(x, 2, 1)
        y1 = DecisionRef(y, 1)
        y2 = DecisionRef(y, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x^2 + y1^2 + 4x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x^2 + y2^2 + 4x + 4y2)
        @test JuMP.isequal_canonical(2x^2 + y1^2 + 4x + 4y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(2x^2 + y2^2 + 4x + 4y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        x = sp[1,:x]
        JuMP.set_objective_coefficient(sp, x, 2, 1, 2.0)
        JuMP.set_objective_coefficient(sp, x, 2, 2, 4.0)
        x = DecisionRef(x, 2, 1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), 2x^2 + y1^2 + 2x + 4y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), 2x^2 + y2^2 + 4x + 4y2)
        @test JuMP.isequal_canonical(2x^2 + y1^2 + 2x + 4y1,
                                     @inferred JuMP.objective_function(sp, 2, 1, DecisionQuadExpr{Float64}))
        @test JuMP.isequal_canonical(2x^2 + y2^2 + 4x + 4y2,
                                     @inferred JuMP.objective_function(sp, 2, 2, DecisionQuadExpr{Float64}))
        @test_throws ErrorException JuMP.objective_function(sp)
    end
end

function test_decision_objective_sense_modification(Structure)
    ξ₁ = @scenario a = 2. probability = 0.5
    ξ₂ = @scenario a = 4 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @objective(sp, Min, x)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y)
        @objective(sp, Min, y)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
    # Second-stage
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
    @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
    # Structure specific
    if sp.structure isa StochasticPrograms.DeterministicEquivalent
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2, JuMP.objective_function(sp))
        set_objective_sense(sp, 2, MOI.MAX_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(x - 0.5*y1 - 0.5*y2, JuMP.objective_function(sp))
        set_objective_sense(sp, 1, MOI.MAX_SENSE)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(x + 0.5*y1 + 0.5*y2, JuMP.objective_function(sp))
        set_objective_sense(sp, MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(x - 0.5*y1 - 0.5*y2, JuMP.objective_function(sp))
        set_objective_sense(sp, 2, 1, MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(x + 0.5*y1 - 0.5*y2, JuMP.objective_function(sp))
    end
    if sp.structure isa StochasticPrograms.StageDecompositionStructure
        @test JuMP.objective_function(sp) == x
        @test JuMP.objective_function(sp, 1) == x
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        set_objective_sense(sp, 2, MOI.MAX_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.objective_function(sp) == x
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        set_objective_sense(sp, 1, MOI.MAX_SENSE)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.objective_function(sp) == x
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        set_objective_sense(sp, MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.objective_function(sp) == x
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
        set_objective_sense(sp, 2, 1, MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.objective_function(sp) == x
        @test JuMP.objective_function(sp, 2, 1) == y1
        @test JuMP.objective_function(sp, 2, 2) == y2
    end
    if sp.structure isa StochasticPrograms.ScenarioDecompositionStructure
        x1 = DecisionRef(sp[1,:x], 2, 1)
        x2 = DecisionRef(sp[1,:x], 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 + y2)
        set_objective_sense(sp, 2, MOI.MAX_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 - y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 - y2)
        set_objective_sense(sp, 1, MOI.MAX_SENSE)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 + y2)
        set_objective_sense(sp, MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 - y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 - y2)
        set_objective_sense(sp, 2, 1, MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2)
        @test MOI.MIN_SENSE == @inferred JuMP.objective_sense(sp, 2, 1)
        @test MOI.MAX_SENSE == @inferred JuMP.objective_sense(sp, 2, 2)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 1), x1 + y1)
        @test JuMP.isequal_canonical(JuMP.objective_function(sp, 2, 2), x2 - y2)
    end
end

function runtests()
    @testset "DecisionObjective" begin
        for structure in [(Deterministic(),),
                          (StageDecomposition(),),
                          (ScenarioDecomposition(),),
                          (Deterministic(), () -> MOIU.MockOptimizer(MOIU.Model{Float64}()))]
            name = length(structure) == 1 ? "$(structure[1])" : "$(structure[1]) with decision bridges"
            @testset "$name" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

function run_dtests()
    @testset "DecisionObjective" begin
        for structure in [(DistributedStageDecomposition(),), (DistributedScenarioDecomposition(),)]
            @testset "$(structure[1])" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

end


text:
#!/usr/bin/env julia

# not complete

function main()
    T = parse(Int, readline())
    P = parse(Int, readline())

    for case in 1:T
        corrects = BitMatrix(undef, 100, 10000)
        for i in 1:100
            corrects[i,:] .= map(x -> parse(Bool, x), split(readline()))
        end
        solve(P)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
"""
    sparsity(A)

Calculates the sparsity of the given array.
The sparsity is the number of zero-valued elements over the
total number of elements.
"""
sparsity(A::AbstractArray{<:Number}) = count(iszero, A)/length(A)

"""
    reldiff(A, B)

Relative difference of absolute values of `A` and `B` defined as

``
\\operatorname{reldiff} = 2 \\dfrac{\\operatorname{abs}(A - B)}{\\operatorname{abs}(A+B)}.
``
"""
function reldiff(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
  return 2*abs.(A-B)./abs.(A+B)
end

"""
    effreldiff(A, B, threshold=1e-14)

Same as `reldiff(A,B)` but with all elements set to zero where corresponding element of
`absdiff(A,B)` is smaller than `threshold`. This is useful in avoiding artificially large
relative errors.
"""
function effreldiff(A::AbstractArray{T}, B::AbstractArray{S}, threshold::Float64=1e-14) where T<:Number where S<:Number
  r = reldiff(A,B)
  r[findall(x->abs.(x)<threshold,absdiff(A,B))] .= 0.
  return r
end

"""
    absdiff(A, B)

Difference of absolute values of `A` and `B`.
"""
function absdiff(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
  return abs.(A-B)
end

"""
    compare(A, B)

Compares two matrices `A` and `B`, prints out the maximal absolute and relative differences
and returns a boolean indicating wether `isapprox(A,B)`.
"""
function compare(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
  @printf("max absdiff: %.1e\n", maximum(absdiff(A,B)))
  @printf("mean absdiff: %.1e\n", mean(absdiff(A,B)))
  @printf("max reldiff: %.1e\n", maximum(reldiff(A,B)))
  @printf("mean reldiff: %.1e\n", mean(reldiff(A,B)))

  r = effreldiff(A,B)
  @printf("effective max reldiff: %.1e\n", maximum(r))
  @printf("effective mean reldiff: %.1e\n", mean(r))

  return isapprox(A,B)
end

"""
    compare_full(A, B)

Compares two matrices `A` and `B`, prints out all absolute and relative differences
and returns a boolean indicating wether `isapprox(A,B)`.
"""
function compare_full(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
  compare(A,B)
  println("")
  println("absdiff: ")
  display(absdiff(A,B))
  println("")
  println("reldiff: ")
  display(reldiff(A,B))
  println("")
  return isapprox(A,B)
end

"""
    setrng(rng)
Replaces current `Random.default_rng()` with `rng`.
"""
function setrng(rng::MersenneTwister)
  Random.default_rng().idxI = rng.idxI
  Random.default_rng().idxF = rng.idxF
  Random.default_rng().state = rng.state
  Random.default_rng().vals = rng.vals
  Random.default_rng().seed = rng.seed
  Random.default_rng().ints = rng.ints
  nothing
end

"""
  saverng(filename [, rng::MersenneTwister; group="GLOBAL_RNG"])
  saverng(HDF5.HDF5File [, rng::MersenneTwister; group="GLOBAL_RNG"])
Saves the current state of Julia's random generator (`Random.default_rng()`) to HDF5.
"""
function saverng(f::HDF5.HDF5File, rng::MersenneTwister=Random.default_rng(); group::String="GLOBAL_RNG")
  g = endswith(group, "/") ? group : group * "/"
  try
    if HDF5.exists(f, g)
      HDF5.o_delete(f, g)
    end

    f[g*"idxF"] = rng.idxF
    f[g*"idxI"] = rng.idxI
    f[g*"state_val"] = rng.state.val
    f[g*"vals"] = rng.vals
    f[g*"seed"] = rng.seed
    f[g*"ints"] = Int.(rng.ints)
  catch e
    error("Error while saving RNG state: ", e)
  end
  nothing
end
function saverng(filename::String, rng::MersenneTwister=Random.default_rng(); group::String="GLOBAL_RNG")
  mode = isfile(filename) ? "r+" : "w"
  HDF5.h5open(filename, mode) do f
    saverng(f, rng; group=group)
  end
end

"""
  loadrng(filename [; group="GLOBAL_RNG"]) -> MersenneTwister
  loadrng(f::HDF5.HDF5File [; group="GLOBAL_RNG"]) -> MersenneTwister
Loads a random generator from HDF5.
"""
function loadrng(f::HDF5.HDF5File; group::String="GLOBAL_RNG")::MersenneTwister
  rng = MersenneTwister(0)
  g = endswith(group, "/") ? group : group * "/"
  try
    rng.idxI = read(f[g*"idxI"])
    rng.idxF = read(f[g*"idxF"])
    rng.state = Random.DSFMT.DSFMT_state(read(f[g*"state_val"]))
    rng.vals = read(f[g*"vals"])
    rng.seed = read(f[g*"seed"])
    rng.ints = UInt128.(read(f[g*"ints"]))
  catch e
    error("Error while restoring RNG state: ", e)
  end
  return rng
end
function loadrng(filename::String; group::String="GLOBAL_RNG")
  HDF5.h5open(filename, "r") do f
    loadrng(f; group=group)
  end
end

"""
  restorerng(filename [; group="GLOBAL_RNG"]) -> Void
  restorerng(f::HDF5.HDF5File [; group="GLOBAL_RNG"]) -> Void
Restores a state of Julia's random generator (`Random.GLOBAL_RNG`) from HDF5.
"""
function restorerng(filename::String; group::String="GLOBAL_RNG")
  HDF5.h5open(filename, "r") do f
    restorerng(f; group=group)
  end
  nothing
end
restorerng(f::HDF5.HDF5File; group::String="GLOBAL_RNG") = setrng(loadrng(f; group=group))

"""
    h5repack(src, trg)
Repacks a HDF5 file e.g. to free unused space. If `src == trg`Wrapper to external h5repack
application.
"""
function h5repack(src::String, trg::String)
    if src == trg   h5repack(src)  end
    @static if Sys.iswindows()
        read(`h5repack.exe $src $trg`, String)
    end
    @static if Sys.islinux()
        read(`h5repack $src $trg`, String)
    end
end
function h5repack(filename::String)
    h5repack(filename, "tmp.h5")
    mv("tmp.h5",filename,force=true)
end

##############################################################
# delete this block when https://github.com/JuliaMath/AbstractFFTs.jl/pull/26 is merged
struct Frequencies <: AbstractVector{Float64}
    nreal::Int
    n::Int
    multiplier::Float64
end

unsafe_getindex(x::Frequencies, i::Int) =
    (i-1+ifelse(i <= x.nreal, 0, -x.n))*x.multiplier
function Base.getindex(x::Frequencies, i::Int)
    (i >= 1 && i <= x.n) || throw(BoundsError())
    unsafe_getindex(x, i)
end
if isdefined(Base, :iterate)
    function Base.iterate(x::Frequencies, i::Int=1)
        i > x.n ? nothing : (unsafe_getindex(x,i), i + 1)
    end
else
    Base.start(x::Frequencies) = 1
    Base.next(x::Frequencies, i::Int) = (unsafe_getindex(x, i), i+1)
    Base.done(x::Frequencies, i::Int) = i > x.n
end
Base.size(x::Frequencies) = (x.n,)
Base.step(x::Frequencies) = x.multiplier

"""
    fftfreq(n, fs=1)
Return discrete fourier transform sample frequencies. The returned
Frequencies object is an AbstractVector containing the frequency
bin centers at every sample point. `fs` is the sample rate of the
input signal.
"""
fftfreq(n::Int, fs::Real=1) = Frequencies(((n-1) >> 1)+1, n, fs/n)

"""
    rfftfreq(n, fs=1)
Return discrete fourier transform sample frequencies for use with
`rfft`. The returned Frequencies object is an AbstractVector
containing the frequency bin centers at every sample point. `fs`
is the sample rate of the input signal.
"""
rfftfreq(n::Int, fs::Real=1) = Frequencies((n >> 1)+1, (n >> 1)+1, fs/n)
FFTW.fftshift(x::Frequencies) = (x.nreal-x.n:x.nreal-1)*x.multiplier
##############################################################

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
export IMEXSolverType

"""
# Description
    IMEXSolverType(;
        splitting_type = HEVISplitting(),
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        implicit_solver_adjustable = false,
        solver_method = ARK2GiraldoKellyConstantinescu,
        solver_storage_variant = LowStorageVariant(),
        split_explicit_implicit = false,
        discrete_splitting = true,
    )

This solver type constructs a solver for ODEs with the
additively-partitioned form.  When `split_explicit_implicit == false`
the equation is assumed to be decomposed as

```math
  \\dot{Q} = [l(Q, t)] + [f(Q, t) - l(Q, t)]
```

where `Q` is the state, `f` is the full tendency and `l` is the chosen implicit
operator.

When `split_explicit_implicit == true` the assumed decomposition is

```math
  \\dot{Q} = [l(Q, t)] + [n(Q, t)]
```

where `n` is now only the nonlinear tendency.

# Arguments
- `splitting_type` (DiscreteSplittingType): The type of discrete
    splitting to apply to the right-hand side.
    Default: `HEVISplitting()`
- `implicit_model` (Type): The model describing dynamics to be
    treated implicitly.
    Default: `AtmosAcousticGravityLinearModel`
- `implicit_solver` (Type): A solver for inverting the
    implicit system of equations.
    Default: `ManyColumnLU`
- `implicit_solver_adjustable` (Bool): A flag identifying whether
    or not the `implicit_solver` can be updated as the time-step
    size changes.
    Default: `false`
- `solver_method` (Function): Function defining the particular additive
    Runge-Kutta method to be used for the IMEX method.
    Default: `ARK2GiraldoKellyConstantinescu`
- `solver_storage_variant` (Type): Storage type for the additive
    Runge-Kutta method.
    Default: `LowStorageVariant()`
- `split_explicit_implicit` (Boolean): Whether the tendency is split in explicit
    and implicit parts or not.
- `discrete_splitting` (Boolean): Boolean denoting whether a PDE level or
    discretized level splitting should be used. If `true` then the PDE is
    discretized in such a way that `f_fast + f_slow` is equivalent to
    discretizing the original PDE directly.

### References
 - [Giraldo2013](@cite)
"""
struct IMEXSolverType{DS, ST} <: AbstractSolverType
    # The type of discrete splitting to apply to the right-hand side
    splitting_type::DS
    # The implicit model
    implicit_model::Type
    # Choice of implicit solver
    implicit_solver::Type
    # Can the implicit solver be updated with changing dt?
    implicit_solver_adjustable::Bool
    # Function for the IMEX method
    solver_method::Function
    # Storage type for the ARK scheme
    solver_storage_variant::ST
    # Split tendency or not
    split_explicit_implicit::Bool
    # Whether to use a PDE level or discrete splitting
    discrete_splitting::Bool

    function IMEXSolverType(;
        splitting_type = HEVISplitting(),
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        implicit_solver_adjustable = false,
        solver_method = ARK2GiraldoKellyConstantinescu,
        solver_storage_variant = LowStorageVariant(),
        split_explicit_implicit = false,
        discrete_splitting = true,
    )
        @assert discrete_splitting || split_explicit_implicit

        DS = typeof(splitting_type)
        ST = typeof(solver_storage_variant)

        return new{DS, ST}(
            splitting_type,
            implicit_model,
            implicit_solver,
            implicit_solver_adjustable,
            solver_method,
            solver_storage_variant,
            split_explicit_implicit,
            discrete_splitting,
        )
    end
end

"""
    getdtmodel(ode_solver::IMEXSolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
function getdtmodel(ode_solver::IMEXSolverType, bl)
    # Most restrictive dynamics are treated implicitly
    return ode_solver.implicit_model(bl)
end

"""
# Description
    solversetup(
        ode_solver::IMEXSolverType{HEVISplitting},
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

Creates an ODE solver using a HEVI-type time-integration
scheme. All horizontal acoustic waves are treated explicitly,
while the 1-D vertical problem is treated implicitly. All
other dynamics (advection, diffusion) is treated explicitly
in the additive Runge-Kutta method.

# Comments:
Currently, the only HEVI-type splitting ClimateMachine can currently
do only involves splitting the acoustic processes; it is not currently
possible to perform more fine-grained separation of tendencies
(for example, including vertical advection or diffusion in the 1-D implicit problem)
"""
function solversetup(
    ode_solver::IMEXSolverType{HEVISplitting},
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
)

    # All we need to do is create a DGModel for the
    # vertical acoustic waves (determined from the `implicit_model`)
    vdg = DGModel(
        ode_solver.implicit_model(dg.balance_law),
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        state_gradient_flux = dg.state_gradient_flux,
        states_higher_order = dg.states_higher_order,
        direction = VerticalDirection(),
    )

    if ode_solver.split_explicit_implicit
        if ode_solver.discrete_splitting
            numerical_flux_first_order = (
                dg.numerical_flux_first_order,
                (dg.numerical_flux_first_order,),
            )
        else
            numerical_flux_first_order = dg.numerical_flux_first_order
        end
        rem_dg = remainder_DGModel(
            dg,
            (vdg,);
            numerical_flux_first_order = numerical_flux_first_order,
        )
        solver = ode_solver.solver_method(
            rem_dg,
            vdg,
            LinearBackwardEulerSolver(
                ode_solver.implicit_solver();
                isadjustable = false,
            ),
            Q;
            split_explicit_implicit = true,
            dt = dt,
            t0 = t0,
        )
    else
        solver = ode_solver.solver_method(
            dg,
            vdg,
            LinearBackwardEulerSolver(
                ode_solver.implicit_solver();
                isadjustable = ode_solver.implicit_solver_adjustable,
            ),
            Q;
            dt = dt,
            t0 = t0,
            # NOTE: This needs to be `false` since the ARK method will
            # evaluate the explicit part using the RemainderModel
            # (Difference between full DG model (dg) and the
            # DG model associated with the 1-D implicit problem (vdg))
            split_explicit_implicit = false,
            variant = ode_solver.solver_storage_variant,
        )
    end

    return solver
end

max_forks_repo_name: jm-c/CLIMA