text:
const randprims = Set([rand, randn, randexp, randperm, shuffle, sample])
const blacklist = Set([Core.apply_type])
const primnames = Set(Symbol(fn) for fn in randprims)

"""
    Options{recurse::Bool, useslots::Bool, naming::Symbol}

Option type that specifies how Julia methods should be transformed into
generative functions. These options are passed as type parameters so
that they are accessible within `@generated` functions.
"""
struct Options{R, U, S} end
const MinimalOptions = Options{false, false, :static}
const DefaultOptions = Options{true, true, :static}

Options(recurse::Bool, useslots::Bool, naming::Symbol) =
    Options{recurse, useslots, naming}()

"Unpack option type parameters as a tuple."
unpack(::Options{R, U, S}) where {R, U, S} = (R, U, S)

"Transform the IR by wrapping sub-calls in `trace`(@ref)."
function automatic_addressing_transform!(ir::IR, options::Options=DefaultOptions())
    recurse, useslots, naming = unpack(options)
    # Modify arguments
    #optarg = argument!(ir; at=1) # Add argument for options
    rand_addrs = Dict() # Map from rand statements to address statements
    # Iterate over IR
    for (x, stmt) in ir
        !isexpr(stmt.expr, :call) && continue
        fn, args, calltype = unpack_call(stmt.expr)
        (unwrap(fn) == :trace ||
         unwrap(fn) == :learnable ||
         fn isa Variable ||
         !istraced(ir, fn, recurse)) && continue
        # Generate address name from function name and arguments
        addr = genaddr(ir, fn, unpack_args(ir, args, calltype))
        addr = insert!(ir, x, QuoteNode(addr))
        rand_addrs[x] = addr # Remember IRVar for address

        # Rewrite statement by wrapping call within `trace`
        rewrite!(ir, x, calltype, QuoteNode(options), addr, fn, args)
    end
    if (useslots) slotaddrs!(ir, rand_addrs) end # Name addresses using slots
    uniqueaddrs!(ir, rand_addrs) # Ensure uniqueness of random addresses
    loopaddrs!(ir, rand_addrs) # Add loop indices to addresses
    return ir
end

"Unpack calls, special casing `Core._apply` and `Core._apply_iterate`."
function unpack_call(expr::Expr)
    fn, args, calltype = expr.args[1], expr.args[2:end], :call
    if fn == GlobalRef(Core, :_apply)
        fn, args, calltype = args[1], args[2:end], :apply
    elseif fn == GlobalRef(Core, :_apply_iterate)
        fn, args, calltype = args[2], args[3:end], :apply_iterate
    end
    return fn, args, calltype
end

"Unpack tuples in IR."
unpack_tuple(ir, v::Variable) = haskey(ir, v) ?
unpack_tuple(ir, ir[v].expr) : nothing
unpack_tuple(ir, e) = iscall(e, GlobalRef(Core, :tuple)) ?
e.args[2:end] : iscall(e, GlobalRef(Base, :getindex)) ? e : error("Expected tuple, got $e.")

"Unpack arguments, special casing `Core._apply` and `Core._apply_iterate`."
function unpack_args(ir, args, calltype)
    if calltype == :call return args end
    unpacked = [unpack_tuple(ir, a) for a in args]
    filter!(a -> !isnothing(a), unpacked)
    return reduce(vcat, unpacked; init=[])
end

"Determine whether a called function should be traced."
function istraced(ir, fn::GlobalRef, recurse::Bool)
    if !isdefined(fn.mod, fn.name) return error("$fn not defined.") end
    val = getfield(fn.mod, fn.name)
    if val in randprims return true end # Primitives are always traced
    if !recurse return false end # Only trace primitives if not recursing
    for m in (Base, Core, Core.Intrinsics) # Filter out Base, Core, etc.
        if isdefined(m, fn.name) && getfield(m, fn.name) == val return false end
    end
    if val isa Type && val <: Sampleable return false end # Filter distributions
    val isa DataType && return false
    val isa UnionAll && return false
    return true
end
function istraced(ir, fn::Function, recurse::Bool) # Handle injected functions
    !(fn in blacklist) && (fn in randprims || recurse)
end

istraced(ir, fn::Variable, recurse::Bool) = # Handle IR variables
!haskey(ir, fn) || istraced(ir, ir[fn].expr, recurse)

istraced(ir, fn::Expr, recurse::Bool) = # Handle keyword functions
iscall(fn, GlobalRef(Core, :kwfunc)) ?
istraced(ir, fn.args[2], recurse) : true
istraced(ir, fn, recurse::Bool) = # Return true by default, to be safe
true

"Static generation of address names."
function genaddr(ir, fn::Symbol, args)
    if length(args) == 0 || !(fn in primnames) return fn end
    argsym = argaddr(ir, args[1])
    return isnothing(argsym) ? fn : Symbol(fn, :_, argsym)
end
function genaddr(ir, fn::Expr, args)
    if iscall(fn, GlobalRef(Core, :kwfunc))
        return genaddr(ir, fn.args[2], args[3:end])
    elseif isexpr(fn, :call) && fn.args[1] isa GlobalRef
        return genaddr(ir, Symbol(Expr(:call, fn.args[1].name)), args)
    else
        return :unknown
    end
end
genaddr(ir, fn::GlobalRef, args) =
genaddr(ir, fn.name, args)
genaddr(ir, fn::Function, args) =
genaddr(ir, nameof(fn), args)
genaddr(ir, fn::Variable, args) =
genaddr(ir, haskey(ir, fn) ? ir[fn].expr : argname(ir, fn), args)
genaddr(ir, fn, args) =
:unknown

"Generate partial address from argument."
argaddr(ir, arg::Variable) =
haskey(ir, arg) ? argaddr(ir, ir[arg].expr) : argname(ir, arg)
argaddr(ir, arg::Expr) =
isexpr(arg, :call) ? argaddr(ir, arg.args[1]) : nothing
argaddr(ir, arg::GlobalRef) =
arg.name
argaddr(ir, arg) =
Symbol(arg)

"Get argument slot name from IR."
argname(ir::IR, v::Variable) =
argname(ir.meta, v)
argname(meta::IRTools.Meta, v::Variable) =
1 <= v.id <= meta.nargs ? meta.code.slotnames[v.id] : nothing
argname(meta::Any, v::Variable) =
nothing

"Rewrites the statement by wrapping call within `trace`."
function rewrite!(ir, var, calltype, options, addr, fn, args)
    if calltype == :call # Handle basic calls
        if unwrap(fn) == :rand
            ir[var] = xcall(Jaynes, :trace, addr, args...)
        else
            ir[var] = xcall(Jaynes, :trace, addr, fn, args...)
        end
    elseif calltype == :apply # Handle `Core._apply`
        preargs = xcall(Core, :tuple, options, addr, fn)
        preargs = insert!(ir, var, preargs)
        ir[var] = xcall(Core, :_apply,
                        GlobalRef(Jaynes, :trace), preargs, args...)
    elseif calltype == :apply_iterate # Handle `Core._apply_iterate`
        preargs = xcall(Core, :tuple, addr, fn)
        preargs = insert!(ir, var, preargs)
        ir[var] = xcall(Core, :_apply_iterate, GlobalRef(Base, :iterate),
                        GlobalRef(Jaynes, :trace), preargs, args...)
    end
end

"Generate trace addresses from slotnames where possible."
function slotaddrs!(ir::IR, rand_addrs::Dict)
    slotnames = ir.meta.code.slotnames
    for (x, stmt) in ir
        if isexpr(stmt.expr, :(=))
            # Attempt to automatically generate address names from slot names
            slot, v = stmt.expr.args
            if !isa(slot, IRTools.Slot) || !(v in keys(rand_addrs)) continue end
            slot_id = parse(Int, string(slot.id)[2:end])
            addr = slotnames[slot_id] # Look up name in CodeInfo
            ir[rand_addrs[v]] = QuoteNode(addr) # Replace previous address
        end
    end
    return ir
end

"Ensure that all addresses have unique names."
function uniqueaddrs!(ir::IR, rand_addrs::Dict)
    counts = Dict{Symbol,Int}()
    firstuses = Dict{Symbol,Variable}()
    addrvars = sort(collect(values(rand_addrs)), by = v -> ir.defs[v.id])
    for v in addrvars # Number all uses after first occurrence
        addr = ir[v].expr.value
        counts[addr] = get(counts, addr, 0) + 1
        if (counts[addr] == 1) firstuses[addr] = v; continue end
        ir[v] = QuoteNode(Symbol(addr, :_, counts[addr]))
    end
    for (addr, c) in counts # Go back and number first use of name
        if c == 1 continue end
        ir[firstuses[addr]] = QuoteNode(Symbol(addr, :_, 1))
    end
    return ir
end

"Add loop indices to addresses."
function loopaddrs!(ir::IR, rand_addrs::Dict)
    # Add count variables for each loop in IR
    loops, countvars = loopcounts!(ir)
    # Append loop count to addresses in each loop body
    for (loop, count) in zip(loops, countvars)
        for addrvar in values(rand_addrs)
            if !(block(ir, addrvar).id in loop.body) continue end
            if iscall(ir[addrvar].expr, GlobalRef(Base, :Pair))
                head, tail = ir[addrvar].expr.args[2:3]
                tail = insert!(ir, addrvar, xcall(Base, :Pair, tail, count))
                ir[addrvar] = xcall(Base, :Pair, head, tail)
            else
                head = insert!(ir, addrvar, ir[addrvar])
                ir[addrvar] = xcall(Base, :Pair, head, count)
            end
        end
    end
    return ir
end


text:
@testset "Parallel.Stress" begin
    gint = loadgraph(joinpath(testdir, "testdata", "graph-50-500.jgz"), "graph-50-500")
    c = vec(readdlm(joinpath(testdir, "testdata", "graph-50-500-sc.txt"), ','))
    for g in testdigraphs(gint)

        z  = LightGraphs.stress_centrality(g)
        zd = @inferred(Parallel.stress_centrality(g; parallel=:distributed))
        @test z == zd == c
        zt = @inferred(Parallel.stress_centrality(g; parallel=:threads))
        @test z == zt == c

        xd = Parallel.stress_centrality(g, 3; parallel=:distributed)
        @test length(xd) == 50
        xt = Parallel.stress_centrality(g, 3; parallel=:threasd)
        @test length(xt) == 50

        xd2 = Parallel.stress_centrality(g, collect(1:20); parallel=:distributed)
        @test length(xd2) == 50
        xt2 = Parallel.stress_centrality(g, collect(1:20); parallel=:threads)
        @test length(xt2) == 50
    end

    g1 = cycle_graph(4)
    add_vertex!(g1)
    add_edge!(g1, 4, 5)

    for g in testgraphs(g1)
        zd = Parallel.stress_centrality(g; parallel=:distributed)
        @test zd == [4, 2, 4, 10, 0]
        zt = Parallel.stress_centrality(g; parallel=:threads)
        @test zt == [4, 2, 4, 10, 0]
    end
end

ext: jl


text:
let
	lifetime_CH4_2000	= 9.6		# years
	conc_CH4_2000		= 1700.0	# ppb
	lifetime_CH4_strat	= 120.0		# years
	lifetime_CH4_soil	= 160.0		# years
	ppb_per_MtCH4		= 1/2.746	# ppb/MtCH4
	ppb_per_MtN			= 1/4.8		# ppb/MtN
	lifetime_N2O		= 120.0		# years

	global function carbonbalance!(emissions, s::ClimateState, p::ClimateParams)
		@unpack Concentration, NetFluxOcean, NetFluxBiosphere = s
		@unpack timestep = p
		dt_Concentration_CO2 = ppm_per_GtC * (emissions[:CO2] - NetFluxOcean - NetFluxBiosphere)
		Concentration[:CO2] += timestep * dt_Concentration_CO2
		@pack! s = Concentration
	end

	global function concentrations!(emissions, s::ClimateState, p::ClimateParams)
		@unpack Concentration = s
		@unpack timestep = p
		# CH4 loss from tropospheric OH (CH4 feedback on its own lifetime), from IPCC TAR WG1 4.2.1.1 (p250)
		Lifetime_CH4_by_OH = lifetime_CH4_2000 * (Concentration[:CH4]/conc_CH4_2000)^0.28

		# All components of CH4 lifetime, from IPCC TAR WG1 4.2.1.1 (p248)
		Lifetime_CH4 = 1 / (1/lifetime_CH4_strat + 1/lifetime_CH4_soil + 1/Lifetime_CH4_by_OH)

		# Simple 1-box exponential decay
		dt_Concentration_CH4 = emissions[:CH4]*ppb_per_MtCH4 - Concentration[:CH4]/Lifetime_CH4
		Concentration[:CH4] += timestep * dt_Concentration_CH4

		# ditto for N2O
		dt_Concentration_N2O = emissions[:N2O]*ppb_per_MtN - Concentration[:N2O]/lifetime_N2O
		Concentration[:N2O] += timestep * dt_Concentration_N2O
		@pack! s = Concentration
	end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
export iauBp00
"""
Frame bias and precession, IAU 2000.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
    date1,date2  double         TT as a 2-part Julian Date (Note 1)

Returned:
    rb           double[3][3]   frame bias matrix (Note 2)
    rp           double[3][3]   precession matrix (Note 3)
    rbp          double[3][3]   bias-precession matrix (Note 4)

Notes:

    1. The TT date date1+date2 is a Julian Date, apportioned in any
    convenient way between the two arguments.  For example,
    JD(TT)=2450123.7 could be expressed in any of these ways,
    among others:

            date1         date2

        2450123.7           0.0       (JD method)
        2451545.0       -1421.3       (J2000 method)
        2400000.5       50123.2       (MJD method)
        2450123.5           0.2       (date & time method)

    The JD method is the most natural and convenient to use in
    cases where the loss of several decimal digits of resolution
    is acceptable.  The J2000 method is best matched to the way
    the argument is handled internally and will deliver the
    optimum resolution.  The MJD method and the date & time methods
    are both good compromises between resolution and convenience.

    2. The matrix rb transforms vectors from GCRS to mean J2000.0 by
    applying frame bias.

    3. The matrix rp transforms vectors from J2000.0 mean equator and
    equinox to mean equator and equinox of date by applying
    precession.

    4. The matrix rbp transforms vectors from GCRS to mean equator and
    equinox of date by applying frame bias then precession.  It is
    the product rp x rb.

    5. It is permissible to re-use the same array in the returned
    arguments.  The arrays are filled in the order given.

Called:
    iauBi00      frame bias components, IAU 2000
    iauPr00      IAU 2000 precession adjustments
    iauIr        initialize r-matrix to identity
    iauRx        rotate around X-axis
    iauRy        rotate around Y-axis
    iauRz        rotate around Z-axis
    iauCr        copy r-matrix
    iauRxr       product of two r-matrices

Reference:
    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

    n.b. The celestial ephemeris origin (CEO) was renamed "celestial
        intermediate origin" (CIO) by IAU 2006 Resolution 2.

This revision:  2013 August 21

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauBp00(double date1, double date2,
#              double rb[3][3], double rp[3][3], double rbp[3][3])

function iauBp00(date1::Real, date2::Real)
    # Allocate return value
    rb  = zeros(Float64, 3, 3)
    rp  = zeros(Float64, 3, 3)
    rbp = zeros(Float64, 3, 3)

    ccall((:iauBp00, libsofa_c), Cvoid,
        (Cdouble, Cdouble,
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
        convert(Float64, date1), convert(Float64, date2),
        rb, rp, rbp)

    return SMatrix{3,3}(rb'), SMatrix{3,3}(rp'), SMatrix{3,3}(rbp')
end


text:
using Test
using MiniKalman, GaussianDistributions, StaticArrays, FillArrays, Statistics, Random,
    Unitful, Parameters

@with_kw struct RandomWalk <: MiniKalman.Model
    step_size
    obs_noise
end
MiniKalman.observation_noise(rw::RandomWalk, _, _) = white_noise1(rw.obs_noise)
MiniKalman.transition_noise(rw::RandomWalk, _, _) = white_noise1(rw.step_size)
MiniKalman.initial_state(rw::RandomWalk) = white_noise1(0.5) # centered around 0
MiniKalman.nobservations(rw::RandomWalk, _) = 1

# Check that sampling from the true model, then inferring the parameters of RandomWalk
# through maximum likelihood applied to that sample, yields the true model.
true_model = RandomWalk(step_size=2.0, obs_noise=0.5)
sr = sample_and_recover(true_model, nothing, 100000,
                        start_model=RandomWalk(step_size=1.0, obs_noise=1.0),
                        rng=MersenneTwister(1))

@test sr.estimated_model.step_size ≈ sr.true_model.step_size rtol=0.02
@test sr.estimated_model.obs_noise ≈ sr.true_model.obs_noise rtol=0.02

# Needs revision:
# N = 100
# rng = MersenneTwister(1)

# # Need to generate it as regular vectors/matrices first to use `\` below.
# xx = rand(rng, 1, 1,N)
# XX = hcat(ones(1,1,N), xx)
# XX_mat = [convert(SMatrix{1,2}, XX[:,:,i]) for i in 1:N]
# YY = 7*xx[:] + randn(rng, N) .* sqrt(1) .+ 15
# YY_vec = map(SVector, YY)

# xxf2, ll2 = kalman_filter(Gaussian([20., 20.0], [1000.0 0; 0 1000.0]), YY_vec,
#                           Fill(MiniKalman.white_noise1(10.0), N),
#                           observation_mats=XX_mat)

# xf = mean(xxf2[end])
# @test xf ≈ XX[1, :, :]' \ YY rtol=0.01
# @assert N==100 # otherwise the following test is not exact
# @test xf ≈ [15.0589, 6.73069] rtol=0.001

# @test sum(ll2) ≈ -328.39288090529027
# rand(white_noise2(1.0u"m^2"))


text:
export Kernel
module Kernel

using CUDA
using GemmKernels
using GemmKernels.Tiling
using GPUifyLoops: @unroll
using StaticArrays

function matmul_singlestage(a, b, c, d,
                          transf_gl2sh_a, transf_gl2sh_b, transf_gl2sh_c, transf_sh2gl_d,
                          transf_sh2rf_a, transf_sh2rf_b, transf_sh2rf_c, transf_rf2sh_d,
                          epilogue,
                          ::Type{conf}) where {conf <: GemmKernels.Config}
    # Calculate the number of fragments needed to fully cover a warp tile
    num_fragments_m = conf.compute_warp.M ÷ conf.compute_op_shape.M
    num_fragments_n = conf.compute_warp.N ÷ conf.compute_op_shape.N

    # Constants
    block_i = (blockIdx().x - 1) * conf.block_shape.M
    block_j = (blockIdx().y - 1) * conf.block_shape.N

    warpId = (threadIdx().x - 1) ÷ 32 + 1
    laneId = (threadIdx().x - 1) % 32 + 1

    gemm_sz = Tile(conf.matmul_shape)
    block_tile = Tile(conf.block_shape)

    # (1) Cooperatively load a block_shape.M x block_shape.N tile of C from global to shared memory within one threadblock
    shmem_c = @cuDynamicSharedMem(Layout.eltype(conf.shared_c_layout), Layout.physical_size(conf.shared_c_layout, block_tile.MN.size))

    @unroll for warp_tile = parallellise(block_tile.MN, Tile(conf.mem_cd_warp), warpId, conf.warps_per_block)
        @unroll for thread_tile = parallellise(warp_tile, Tile(conf.mem_cd_thread), laneId, 32)
            x = Layout.load(conf.global_c_layout, c, translate_base(thread_tile, (M = block_i, N = block_j)))
            x = transf_gl2sh_c(x, thread_tile)
            Layout.store!(conf.shared_c_layout, shmem_c, x, thread_tile)
        end
    end

    sync_threads()

    # (2) Load a compute_warp.M x compute_warp.N tile of C from shared memory into registers
    warp_tile = subdivide(block_tile.MN, Tile(conf.compute_warp).MN, warpId, conf.warps_per_block)

    c_frags = MArray{Tuple{num_fragments_m, num_fragments_n}, Operator.fragtype_accum(conf.operator, conf.shared_c_layout)}(undef)

    @unroll for i = 1 : num_fragments_m
        @unroll for j = 1 : num_fragments_n
            tile = translate_offset(warp_tile, (M = (i-1)*conf.compute_op_shape.M, N = (j-1)*conf.compute_op_shape.N))
            @inbounds c_frags[i, j] = transf_sh2rf_c(Operator.load_c(conf.operator, conf.shared_c_layout, shmem_c, tile), tile)
        end
    end

    sync_threads()

    # (3) Compute a block_shape.M x block_shape.N x block_shape.K matrix product within one threadblock
    shmem_a = @cuDynamicSharedMem(Layout.eltype(conf.shared_a_layout), Layout.physical_size(conf.shared_a_layout, block_tile.MK.size))
    shmem_b = @cuDynamicSharedMem(Layout.eltype(conf.shared_b_layout), Layout.physical_size(conf.shared_b_layout, block_tile.KN.size),
                                    length(shmem_a) * sizeof(Layout.eltype(conf.shared_a_layout)))

    @unroll for block_k = 0 : block_tile.size.K : gemm_sz.size.K - 1
        if Layout.threadblock_condition(conf.global_a_layout, conf.global_b_layout, block_i, block_j, block_k, block_tile)
            # (3.1) Cooperatively load a block_shape.M x block_shape.K tile of A from global to shared memory within one threadblock
            @unroll for warp_tile = parallellise(block_tile.MK, Tile(conf.mem_a_warp), warpId, conf.warps_per_block, conf.is_a_col_major)
                @unroll for thread_tile = parallellise(warp_tile, Tile(conf.mem_a_thread), laneId, 32, conf.is_a_col_major)
                    x = Layout.load(conf.global_a_layout, a, translate_base(thread_tile, (M = block_i, K = block_k)))
                    x = transf_gl2sh_a(x, thread_tile)
                    Layout.store!(conf.shared_a_layout, shmem_a, x, thread_tile)
                end
            end

            # (3.2) Cooperatively load a block_shape.K x block_shape.N tile of B from global to shared memory within one threadblock
            @unroll for warp_tile = parallellise(block_tile.KN, Tile(conf.mem_b_warp), warpId, conf.warps_per_block, conf.is_b_col_major)
                @unroll for thread_tile = parallellise(warp_tile, Tile(conf.mem_b_thread), laneId, 32, conf.is_b_col_major)
                    x = Layout.load(conf.global_b_layout, b, translate_base(thread_tile, (K = block_k, N = block_j)))
                    x = transf_gl2sh_b(x, thread_tile)
                    Layout.store!(conf.shared_b_layout, shmem_b, x, thread_tile)
                end
            end

            sync_threads()

            # (3.3) Calculate a compute_warp.M x compute_warp.N tile of D, using a compute_warp.M x compute_warp.N x compute_warp.K operation
            @unroll for warp_tile = parallellise(block_tile, Tile(conf.compute_warp), warpId, conf.warps_per_block)
                # (3.3.1) Load a compute_warp.M x compute_warp.K tile of A from shared memory into registers
                a_frags = MArray{Tuple{num_fragments_m}, Operator.fragtype_a(conf.operator, conf.shared_a_layout)}(undef)

                @unroll for i = 1 : num_fragments_m
                    a_tile = translate_offset(warp_tile.MK, (M = (i-1)*conf.compute_op_shape.M, K = 0))
                    @inbounds a_frags[i] = transf_sh2rf_a(Operator.load_a(conf.operator, conf.shared_a_layout, shmem_a, a_tile), a_tile)
                end

                # (3.3.2) Load a compute_warp.K x compute_warp.N tile of B from shared memory into registers
                b_frags = MArray{Tuple{num_fragments_n}, Operator.fragtype_b(conf.operator, conf.shared_b_layout)}(undef)

                @unroll for j = 1 : num_fragments_n
                    b_tile = translate_offset(warp_tile.KN, (K = 0, N = (j-1)*conf.compute_op_shape.N))
                    @inbounds b_frags[j] = transf_sh2rf_b(Operator.load_b(conf.operator, conf.shared_b_layout, shmem_b, b_tile), b_tile)
                end

                # (3.3.3) Compute a compute_warp.M x compute_warp.N x compute_warp.K matrix product within one warp
                @unroll for i = 1 : num_fragments_m
                    @unroll for j = 1 : num_fragments_n
                        @inbounds c_frags[i, j] = Operator.mma(conf.operator, a_frags[i], b_frags[j], c_frags[i, j])
                    end
                end
            end
        end

        sync_threads()
    end

    # (4) Store the compute_warp.M x compute_warp.N tile of D from registers to shared memory
    shmem_d = @cuDynamicSharedMem(Layout.eltype(conf.shared_d_layout), Layout.physical_size(conf.shared_d_layout, block_tile.MN.size))

    warp_tile = subdivide(block_tile.MN, Tile(conf.compute_warp).MN, warpId, conf.warps_per_block)

    @unroll for i = 1 : num_fragments_m
        @unroll for j = 1 : num_fragments_n
            tile = translate_offset(warp_tile, (M = (i-1)*conf.compute_op_shape.M, N = (j-1)*conf.compute_op_shape.N))
            Operator.store_d(conf.operator, conf.shared_d_layout, shmem_d, transf_rf2sh_d(c_frags[i, j], tile), tile)
        end
    end

    sync_threads()

    # (5) Run the epilogue
    epilogue(d, shmem_d, transf_sh2gl_d, conf)

    return
end

function matmul_pipelined(a, b, c, d,
                          transf_gl2sh_a, transf_gl2sh_b, transf_gl2sh_c, transf_sh2gl_d,
                          transf_sh2rf_a, transf_sh2rf_b, transf_sh2rf_c, transf_rf2sh_d,
                          epilogue,
                          ::Type{conf}) where {conf <: GemmKernels.Config}
    # Calculate the number of fragments needed to fully cover a warp tile
    num_fragments_m = conf.compute_warp.M ÷ conf.compute_op_shape.M
    num_fragments_n = conf.compute_warp.N ÷ conf.compute_op_shape.N

    # Constants
    block_i = (blockIdx().x - 1) * conf.block_shape.M
    block_j = (blockIdx().y - 1) * conf.block_shape.N

    warpId = (threadIdx().x - 1) ÷ 32 + 1
    laneId = (threadIdx().x - 1) % 32 + 1

    gemm_sz = Tile(conf.matmul_shape)
    block_tile = Tile(conf.block_shape)

    # (1) Cooperatively load a block_shape.M x block_shape.N tile of C from global to shared memory within one threadblock
    shmem_c = @cuDynamicSharedMem(Layout.eltype(conf.shared_c_layout), Layout.physical_size(conf.shared_c_layout, block_tile.MN.size))

    @unroll for warp_tile = parallellise(block_tile.MN, Tile(conf.mem_cd_warp), warpId, conf.warps_per_block)
        @unroll for thread_tile = parallellise(warp_tile, Tile(conf.mem_cd_thread), laneId, 32)
            x = Layout.load(conf.global_c_layout, c, translate_base(thread_tile, (M = block_i, N = block_j)))
            x = transf_gl2sh_c(x, thread_tile)
            Layout.store!(conf.shared_c_layout, shmem_c, x, thread_tile)
        end
    end

    sync_threads()

    # (2) Load a compute_warp.M x compute_warp.N tile of C from shared memory into registers
    warp_tile = subdivide(block_tile.MN, Tile(conf.compute_warp).MN, warpId, conf.warps_per_block)

    c_frags = MArray{Tuple{num_fragments_m, num_fragments_n}, Operator.fragtype_accum(conf.operator, conf.shared_c_layout)}(undef)

    @unroll for i = 1 : num_fragments_m
        @unroll for j = 1 : num_fragments_n
            tile = translate_offset(warp_tile, (M = (i-1)*conf.compute_op_shape.M, N = (j-1)*conf.compute_op_shape.N))
            @inbounds c_frags[i, j] = transf_sh2rf_c(Operator.load_c(conf.operator, conf.shared_c_layout, shmem_c, tile), tile)
        end
    end

    sync_threads()

    # (3) Compute a block_shape.M x block_shape.N x block_shape.K matrix product within one threadblock
    shmem_a = @cuDynamicSharedMem(Layout.eltype(conf.shared_a_layout), Layout.physical_size(conf.shared_a_layout, block_tile.MK.size))
    shmem_b = @cuDynamicSharedMem(Layout.eltype(conf.shared_b_layout), Layout.physical_size(conf.shared_b_layout, block_tile.KN.size),
                                    length(shmem_a) * sizeof(Layout.eltype(conf.shared_a_layout)))

    # Sizes of a_fragment and b_fragment
    a_frag_i = (block_tile.size.M * block_tile.size.K) ÷ (conf.mem_a_warp.M * conf.mem_a_warp.K * conf.warps_per_block)
    a_frag_j = (conf.mem_a_warp.M * conf.mem_a_warp.K) ÷ (conf.mem_a_thread.M * conf.mem_a_thread.K * 32)
    b_frag_i = (block_tile.size.K * block_tile.size.N) ÷ (conf.mem_b_warp.K * conf.mem_b_warp.N * conf.warps_per_block)
    b_frag_j = (conf.mem_b_warp.K * conf.mem_b_warp.N) ÷ (conf.mem_b_thread.K * conf.mem_b_thread.N * 32)

    a_fragment = MArray{Tuple{a_frag_i, a_frag_j}, Layout.fragtype(conf.global_a_layout, conf.mem_a_thread)}(undef)
    b_fragment = MArray{Tuple{b_frag_i, b_frag_j}, Layout.fragtype(conf.global_b_layout, conf.mem_b_thread)}(undef)

    a_frags = MArray{Tuple{2, num_fragments_m}, Operator.fragtype_a(conf.operator, conf.shared_a_layout)}(undef)
    b_frags = MArray{Tuple{2, num_fragments_n}, Operator.fragtype_b(conf.operator, conf.shared_b_layout)}(undef)

    warp_tile_mn = subdivide(block_tile, Tile(conf.compute_warp), warpId, conf.warps_per_block)

    # ld.global(0 : block_shape.K)
    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.MK, Tile(conf.mem_a_warp), warpId, conf.warps_per_block, conf.is_a_col_major))
        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_a_thread), laneId, 32, conf.is_a_col_major))
            @inbounds a_fragment[i, j] = Layout.load(conf.global_a_layout, a, translate_base(thread_tile, (M = block_i, K = 0)))
        end
    end

    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.KN, Tile(conf.mem_b_warp), warpId, conf.warps_per_block, conf.is_b_col_major))
        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_b_thread), laneId, 32, conf.is_b_col_major))
            @inbounds b_fragment[i, j] = Layout.load(conf.global_b_layout, b, translate_base(thread_tile, (K = 0, N = block_j)))
        end
    end

    # st.shared()
    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.MK, Tile(conf.mem_a_warp), warpId, conf.warps_per_block, conf.is_a_col_major))
        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_a_thread), laneId, 32, conf.is_a_col_major))
            @inbounds x = transf_gl2sh_a(a_fragment[i, j], thread_tile)
            Layout.store!(conf.shared_a_layout, shmem_a, x, thread_tile)
        end
    end

    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.KN, Tile(conf.mem_b_warp), warpId, conf.warps_per_block, conf.is_b_col_major))
        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_b_thread), laneId, 32, conf.is_b_col_major))
            @inbounds x = transf_gl2sh_b(b_fragment[i, j], thread_tile)
            Layout.store!(conf.shared_b_layout, shmem_b, x, thread_tile)
        end
    end

    sync_threads()

    # ld.shared(0 : compute_op_shape.K, stage = 1)
    warp_tile = translate_offset(warp_tile_mn, (M = 0, N = 0, K = 0))

    @unroll for i = 1 : num_fragments_m
        a_tile = translate_offset(warp_tile.MK, (M = (i-1)*conf.compute_op_shape.M, K = 0))
        @inbounds a_frags[1, i] = transf_sh2rf_a(Operator.load_a(conf.operator, conf.shared_a_layout, shmem_a, a_tile), a_tile)
    end

    @unroll for j = 1 : num_fragments_n
        b_tile = translate_offset(warp_tile.KN, (K = 0, N = (j-1)*conf.compute_op_shape.N))
        @inbounds b_frags[1, j] = transf_sh2rf_b(Operator.load_b(conf.operator, conf.shared_b_layout, shmem_b, b_tile), b_tile)
    end

    # ld.global(block_shape.K : 2 * block_shape.K)
    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.MK, Tile(conf.mem_a_warp), warpId, conf.warps_per_block, conf.is_a_col_major))
        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_a_thread), laneId, 32, conf.is_a_col_major))
            @inbounds a_fragment[i, j] = Layout.load(conf.global_a_layout, a, translate_base(thread_tile, (M = block_i, K = block_tile.size.K)))
        end
    end

    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.KN, Tile(conf.mem_b_warp), warpId, conf.warps_per_block, conf.is_b_col_major))
        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_b_thread), laneId, 32, conf.is_b_col_major))
            @inbounds b_fragment[i, j] = Layout.load(conf.global_b_layout, b, translate_base(thread_tile, (K = block_tile.size.K, N = block_j)))
        end
    end

    @unroll 1 for block_k = 0 : block_tile.size.K : gemm_sz.size.K - 1
        @unroll for (i, warp_k) = enumerate(0 : conf.compute_op_shape.K : block_tile.size.K - 1)
            cur_stage = mod1(i, 2)
            nxt_stage = mod1(i + 1, 2)

            if i == block_tile.size.K ÷ conf.compute_op_shape.K # last iteration
                sync_threads()

                # st.shared()
                @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.MK, Tile(conf.mem_a_warp), warpId, conf.warps_per_block, conf.is_a_col_major))
                    @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_a_thread), laneId, 32, conf.is_a_col_major))
                        @inbounds x = transf_gl2sh_a(a_fragment[i, j], thread_tile)
                        Layout.store!(conf.shared_a_layout, shmem_a, x, thread_tile)
                    end
                end

                @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.KN, Tile(conf.mem_b_warp), warpId, conf.warps_per_block, conf.is_b_col_major))
                    @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_b_thread), laneId, 32, conf.is_b_col_major))
                        @inbounds x = transf_gl2sh_b(b_fragment[i, j], thread_tile)
                        Layout.store!(conf.shared_b_layout, shmem_b, x, thread_tile)
                    end
                end

                sync_threads()

                # avoid out of bounds access for global memory
                if block_k < (gemm_sz.size.K - 2 * block_tile.size.K)
                    # ld.global(block_k + 2 * block_shape.K : block_k + 3 * block_shape.K)
                    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.MK, Tile(conf.mem_a_warp), warpId, conf.warps_per_block, conf.is_a_col_major))
                        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_a_thread), laneId, 32, conf.is_a_col_major))
                            @inbounds a_fragment[i, j] = Layout.load(conf.global_a_layout, a, translate_base(thread_tile, (M = block_i, K = block_k + 2 * block_tile.size.K)))
                        end
                    end

                    @unroll for (i, warp_tile) = enumerate(parallellise(block_tile.KN, Tile(conf.mem_b_warp), warpId, conf.warps_per_block, conf.is_b_col_major))
                        @unroll for (j, thread_tile) = enumerate(parallellise(warp_tile, Tile(conf.mem_b_thread), laneId, 32, conf.is_b_col_major))
                            @inbounds b_fragment[i, j] = Layout.load(conf.global_b_layout, b, translate_base(thread_tile, (K = block_k + 2 * block_tile.size.K, N = block_j)))
                        end
                    end
                end
            end

            # ld.shared((warp_k + compute_op_shape.K) % block_shape.K, stage = nxt_stage)
            warp_tile = translate_offset(warp_tile_mn, (M = 0, N = 0, K = (warp_k + conf.compute_op_shape.K) % block_tile.size.K))

            @unroll for i = 1 : num_fragments_m
                a_tile = translate_offset(warp_tile.MK, (M = (i-1)*conf.compute_op_shape.M, K = 0))
                @inbounds a_frags[nxt_stage, i] = transf_sh2rf_a(Operator.load_a(conf.operator, conf.shared_a_layout, shmem_a, a_tile), a_tile)
            end

            @unroll for j = 1 : num_fragments_n
                b_tile = translate_offset(warp_tile.KN, (K = 0, N = (j-1)*conf.compute_op_shape.N))
                @inbounds b_frags[nxt_stage, j] = transf_sh2rf_b(Operator.load_b(conf.operator, conf.shared_b_layout, shmem_b, b_tile), b_tile)
            end

            # mma(cur_stage)
            @unroll for i = 1 : num_fragments_m
                @unroll for j = 1 : num_fragments_n
                    @inbounds c_frags[i, j] = Operator.mma(conf.operator, a_frags[cur_stage, i], b_frags[cur_stage, j], c_frags[i, j])
                end
            end
        end

        sync_threads()
    end

    # (4) Store the compute_warp.M x compute_warp.N tile of D from registers to shared memory
    shmem_d = @cuDynamicSharedMem(Layout.eltype(conf.shared_d_layout), Layout.physical_size(conf.shared_d_layout, block_tile.MN.size))

    warp_tile = subdivide(block_tile.MN, Tile(conf.compute_warp).MN, warpId, conf.warps_per_block)

    @unroll for i = 1 : num_fragments_m
        @unroll for j = 1 : num_fragments_n
            tile = translate_offset(warp_tile, (M = (i-1)*conf.compute_op_shape.M, N = (j-1)*conf.compute_op_shape.N))
            Operator.store_d(conf.operator, conf.shared_d_layout, shmem_d, transf_rf2sh_d(c_frags[i, j], tile), tile)
        end
    end

    sync_threads()

    # (5) Run the epilogue
    epilogue(d, shmem_d, transf_sh2gl_d, conf)

    return
end

end


text:
struct PrimalLieTrotter end
struct DualLieTrotter end
struct Strang end

struct ProjectorSplitting_Params{sType, lType, kType}
    S_rhs # rhs of S step (core projected rhs)
    L_rhs # rhs of L step (range projected rhs)
    K_rhs # rhs of K step (corange projected rhs)
    S_kwargs
    L_kwargs
    K_kwargs
    S_alg::sType
    L_alg::lType
    K_alg::kType
end

struct ProjectorSplitting_Cache{uType,SIntegratorType,LIntegratorType,KIntegratorType,yType} <: AbstractDLRAlgorithm_Cache
    US::Matrix{uType}
    VS::Matrix{uType}
    QRK::LinearAlgebra.QRCompactWY{uType, Matrix{uType}}
    QRL::LinearAlgebra.QRCompactWY{uType, Matrix{uType}}
    SIntegrator::SIntegratorType
    LIntegrator::LIntegratorType
    KIntegrator::KIntegratorType
    y
    ycurr::yType
    yprev::yType
    Δy::yType
end

struct ProjectorSplitting{oType, sType, lType, kType} <: AbstractDLRAlgorithm
    order::oType
    alg_params::ProjectorSplitting_Params{sType, lType, kType}
end

function ProjectorSplitting(order = PrimalLieTrotter();S_rhs = nothing, L_rhs = nothing, K_rhs = nothing,
                             S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(),
                             S_alg=Tsit5(), L_alg = Tsit5(), K_alg = Tsit5())
    params = ProjectorSplitting_Params(S_rhs,L_rhs,K_rhs,S_kwargs,L_kwargs,K_kwargs,S_alg,L_alg,K_alg)
    return ProjectorSplitting(order, params)
end

function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1])
    # allocate memory for frequently accessed arrays
    tspan = (t0,t0+dt)

    US = u.U*u.S
    VS = u.V*u.S'
    QRK = qr(US)
    QRL = qr(VS)

    if isnothing(alg.alg_params.K_rhs)
        K_rhs = function (US, V, t)
                    return Matrix(prob.f(TwoFactorRepresentation(US,V),t)*V)
                end
    else
        K_rhs = alg.alg_params.K_rhs
    end
    KProblem = ODEProblem(K_rhs, US, tspan, u.V)
    KIntegrator = init(KProblem, alg.alg_params.K_alg, save_everystep=false, alg.alg_params.K_kwargs...)

    if isnothing(alg.alg_params.S_rhs)
        S_rhs = function (S, (U,V), t)
                    return Matrix(-U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V)
                end
    else
        S_rhs = alg.alg_params.S_rhs
    end
    SProblem = ODEProblem(S_rhs, QRK.R, tspan, (u.U, u.V))
    SIntegrator = init(SProblem, alg.alg_params.S_alg, save_everystep=false, alg.alg_params.S_kwargs...)

    if isnothing(alg.alg_params.L_rhs)
        L_rhs = function (VS, U, t)
                    return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U)
                end
    else
        L_rhs = alg.alg_params.L_rhs
    end
    LProblem = ODEProblem(L_rhs, VS, tspan, u.U)
    LIntegrator = init(LProblem, alg.alg_params.L_alg, save_everystep=false, alg.alg_params.L_kwargs...)

    return ProjectorSplitting_Cache(US, VS, QRK, QRL,
                                    SIntegrator, LIntegrator, KIntegrator,
                                    nothing, nothing, nothing, nothing)
end

function alg_cache(prob::MatrixDataProblem, alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1])
    # creates caches for frequently used arrays by performing the first time step
    @unpack y = prob

    yprev = y isa AbstractArray ? deepcopy(y[1]) : y(t0)
    ycurr = deepcopy(yprev)
    Δy = similar(yprev)
    VS = u.V*u.S'
    US = u.U*u.S
    QRK = qr(US)
    QRL = qr(VS)
    KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1)
    SIntegrator = MatrixDataIntegrator(Δy, QRK.R, u.U, u.V, -1)
    LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1)

    return ProjectorSplitting_Cache(US, VS, QRK, QRL,
                                    SIntegrator, LIntegrator, KIntegrator,
                                    prob.y, ycurr, yprev, Δy)
end

function init(prob::AbstractDLRProblem, alg::ProjectorSplitting, dt)
    t0, tf = prob.tspan
    @assert tf > t0 "Integration in reverse time direction is not supported"
    u = deepcopy(prob.u0)
    sol = init_sol(dt, t0, tf, prob.u0)
    cache = alg_cache(prob, alg, u, dt, t0 = t0)
    sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object
    return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0)
end

function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem})
    @unpack y, ycurr, yprev, Δy = cache
    update_data!(ycurr, y, t, dt)
    Δy .= ycurr - yprev
    yprev .= ycurr
    primal_LT_step!(u, cache, t, dt)
end

function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem})
    primal_LT_step!(u, cache, t, dt)
end

function primal_LT_step!(u, cache, t, dt)
    @unpack US, VS, QRK, QRL, KIntegrator, SIntegrator, LIntegrator = cache

    # K step
    mul!(US,u.U,u.S)
    set_u!(KIntegrator, US)
    step!(KIntegrator, dt, true)
    US .= KIntegrator.u
    QRL = qr!(US)
    u.U .= Matrix(QRL.Q)

    # S step
    set_u!(SIntegrator, QRL.R)
    step!(SIntegrator, dt, true)

    # L step
    mul!(VS,u.V,SIntegrator.u')
    set_u!(LIntegrator, VS)
    step!(LIntegrator, dt, true)
    VS .= LIntegrator.u
    QRL = qr!(VS)
    u.V .= Matrix(QRL.Q)
    u.S .= QRL.R'
end

function dual_LT_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem})
    @unpack y, ycurr, yprev, Δy = cache
    update_data!(ycurr, y, t, dt)
    Δy .= ycurr - yprev
    yprev .= ycurr
    dual_LT_step!(u, cache, t, dt)
end

function dual_LT_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem})
    dual_LT_step!(u, cache, t, dt)
end

function dual_LT_step!(u, cache, t, dt)
    @unpack US, VS, QRK, QRL, KIntegrator, SIntegrator, LIntegrator = cache

    # L step
    mul!(VS,u.V,u.S')
    set_u!(LIntegrator, VS)
    step!(LIntegrator, dt, true)
    VS .= LIntegrator.u
    QRL = qr!(VS)
    u.V .= Matrix(QRL.Q)

    # S step
    set_u!(SIntegrator, Matrix(QRL.R'))
    step!(SIntegrator, dt, true)

    # K step
    mul!(US,u.U,SIntegrator.u)
    set_u!(KIntegrator, US)
    step!(KIntegrator, dt, true)
    US .= KIntegrator.u
    QRK = qr!(US)
    u.U .= Matrix(QRK.Q)
    u.S .= QRK.R
end

function step!(integrator::DLRIntegrator, ::ProjectorSplitting{PrimalLieTrotter,S,L,K}, dt) where {S,L,K}
    @unpack u, t, iter, cache, probType = integrator
    primal_LT_step!(u, cache, t, dt, probType)
    integrator.t += dt
    integrator.iter += 1
end

function step!(integrator::DLRIntegrator, ::ProjectorSplitting{DualLieTrotter,S,L,K}, dt) where {S,L,K}
    @unpack u, t, iter, cache, probType = integrator
    dual_LT_step!(u, cache, t, dt, probType)
    integrator.t += dt
    integrator.iter += 1
end

function step!(integrator::DLRIntegrator, ::ProjectorSplitting{Strang,S,L,K}, dt) where {S,L,K}
    @unpack u, t, iter, cache, probType = integrator
    primal_LT_step!(u, cache, t, dt/2, probType)
    dual_LT_step!(u, cache, t + dt/2, dt/2, probType)
    integrator.t += dt
    integrator.iter += 1
end


text:
module AIMASamples

include("VacuumEnvironment.jl")
include("Romania.jl")
include("SlideBlock.jl")
include("Queens.jl")

end # module


text:
#=
% Example: Find the ocean volume that has originated from each    %
%            surface box.                                           %
%                                                                   %
% This is equivalent to solving a sensitivity problem:              %
% The total volume is V = v^T c , where v is the volume of each box %
% and c is the fraction of volume from a given source which         %
% satisfies the equation A c = d.                                   %
% Next, dV/d(d) = A^(-T) v, and dV/d(d) is exactly the volume       %
% originating from each source.
%
% See Section 3 and Supplementary Section 4, Gebbie & Huybers 2011.
=#
using Revise
using TMI, PyPlot, PyCall

TMIversion = "modern_90x45x33_GH10_GH12"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)

volume = volumefilled(TMIversion,Alu,γ)

# view the surface
cntrs = 1:0.25:6
tlabel="Volume filled by surface patch [log10(m)]"
planviewplot(volume,cntrs,γ,titlelabel=tlabel)


text:
module Sessions
using DataFrames, NaturalSort, DataStructures, HDF5

export PROCESSED_DATA_DIR, REGIONS, REGION_LABELS, Session, ClusterMetadata, sessionnames, readsession, patient, siteregion

const RAW_DATA_DIR = joinpath(dirname(@__FILE__), "..", "raw")
const PROCESSED_DATA_DIR = joinpath(dirname(@__FILE__), "..", "processed")
const REGIONS = OrderedDict(
    :H => r"^..H$",
    :EC => r"^.EC",
    :A => r"^.A$",
    :PHG => r"^.(?:PHG|PHC|PG.?)$",
    :MTL => r"^.(?:.H|EC|A|PHG|PHC|PG.?)$",
    :All => r""
)
const REGION_LABELS = OrderedDict(
    :H => "Hippocampus",
    :EC => "Entorhinal Cortex",
    :A => "Amygdala",
    :PHG => "Parahippocampal Cortex",
    :MTL => "Medial Temporal Lobe",
    :All => "All"
)

type Session
    folder::String
    trials::DataFrame
    spikes::Vector{Vector{Vector{Float64}}}
    unit_type::Vector{Vector{Symbol}}
    sites::Vector{String}
    stimuli::Vector{String}
end

sessionnames() = sort!([x[1:rsearch(x, '.')-1] for x in readdir(RAW_DATA_DIR) if endswith(x, ".h5")], lt=natural)
splitcols(x) = [x[:, i] for i = 1:size(x, 2)]
function readsession(sname)::Session
    h5open(joinpath(RAW_DATA_DIR, "$sname.h5")) do f
        ftrials = f["trials"]
        stim = read(ftrials, "stimuli")
        nstim = maximum(stim)
        trials = DataFrame(correct=read(ftrials, "correct") .== 1,
                          isi=read(ftrials, "isi"),
                          probe1=read(ftrials, "probe1"),
                          probe2=read(ftrials, "probe2"),
                          stimuli=splitcols(stim),
                          onsets=splitcols(read(ftrials, "onsets")),
                          offsets=splitcols(read(ftrials, "offsets")),
                          maint=read(ftrials, "maint"),
                          probes=read(ftrials, "probes"))
        close(ftrials)
        fchannels = f["channels"]

        sites = read(f, "sites")
        spikes = Vector{Vector{Vector{Float64}}}(length(sites))
        unit_type = Vector{Vector{Symbol}}(length(sites))
        for i = 1:length(sites)
            st = "ch$(i)_spike_times"
            if exists(fchannels, st)
                spikes[i] = read(fchannels, st)
                unit_type[i] = map(Symbol, read(fchannels, "ch$(i)_unit_types"))
            end
        end

        # Not included in data files for reasons of confidentiality
        stimuli = ["Stimulus $i" for i in 1:nstim]

        Session(sname, trials, spikes, unit_type, sites, stimuli)
    end
end

immutable ClusterMetadata
    session::String
    site::String
    channel::Int
    cluster::Int
    unit_type::Symbol
end

ClusterMetadata(s::Session, ch::Int, cluster::Int) =
    ClusterMetadata(s.folder, s.sites[ch], ch, cluster, s.unit_type[ch][cluster])

Base.:(==)(x::ClusterMetadata, y::ClusterMetadata) =
    x.cluster == y.cluster && x.channel == y.channel && x.session == y.session
Base.hash(x::ClusterMetadata) = hash(x.session, hash(x.site, hash(x.channel, hash(x.cluster))))
Base.isless(x::ClusterMetadata, y::ClusterMetadata) =
    x.session != y.session ? natural(x.session, y.session) : x.channel != y.channel ? isless(x.channel, y.channel) :
    isless(x.cluster, y.cluster)
Base.show(io::IO, m::ClusterMetadata) =
    print(io, m.session, ' ', m.site, " ch", m.channel, " clu", m.cluster, ' ', m.unit_type)

function siteregion(m::ClusterMetadata)
    for (k, v) in REGIONS
        (k == :All || k == :MTL) && continue
        ismatch(v, m.site) && return k
    end
    ismatch(REGIONS[:MTL], m.site) && return :MTL
    :All
end
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
module TestPredicatedSets
using Test
using Catlab.CategoricalAlgebra.PredicatedSets
using Catlab.Theories

odds = PredicatedSet(Int, x -> x % 2 == 1)
evens = PredicatedSet(Int, x -> x % 2 == 0)

plus_one_to_odd = PredicatedHom(odds,evens,x -> x+1)
plus_one_to_even = PredicatedHom(evens,odds,x -> x+1)

@test plus_one_to_odd(1) == 2
@test_throws AssertionError plus_one_to_odd(2) == 3
@test plus_one_to_even(2) == 3
@test_throws AssertionError plus_one_to_even(3) == 4

@test compose(plus_one_to_odd, plus_one_to_even)(3) == 5

plus_two = PredicatedHom(odds,evens,x -> x+2)

@test_throws AssertionError plus_two(3)
end


text:
using Optimize
using NLPModels

using Plots
# pyplot()
plotly()

n_min = 1
n_max = 1

# Math Prog Base collection
include("MPBProblems-1D.jl")
test_probs = mpb_probs_1D

using ScalarSolvers
solvers = [bissect, bissect_Cub, bissect_nwt, bissect_sec, bissect_secA, ARC_Cub, ARC_Nwt, ARC_Sec, ARC_SecA, TR_Cub, TR_Nwt, TR_Sec, TR_SecA, zoom_Cub, zoom_Nwt, zoom_Sec, zoom_SecA]
labels = ["Biss"; "Biss_Cub"; "Biss_Nwt"; "Biss_Sec"; "Biss_SecA"; "ARC_Cub"; "ARC_Nwt"; "ARC_Sec"; "ARC_SecA"; "TR_Cub";"TR_Nwt"; "TR_Sec"; "TR_SecA"; "zoom_Cub"; "zoom_Nwt"; "zoom_Sec"; "zoom_SecA"]

options = [Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           Dict{Symbol,Any}(:verbose=>false, :tol => 0.1)
           ]

include("compare_solvers.jl")
# s1, P1, t1, Pt1 = compare_solvers_with_options2(solvers, options, labels, test_probs, n_min, n_max, printskip = false)
P = compare_solvers_with_options2(solvers, options, labels, test_probs, n_min, n_max, printskip = false)[2]


text:
# This file was generated, do not modify it. # hide
mtm = machine(tm, Xhp, y)
fit!(mtm)
rep = report(mtm)

res = rep.plotting

@show round.(res.measurements.^2, digits=2)
@show argmin(res.measurements)


text:
module JFVMvis

using JFVM, PyPlot

export visualizeCells, visualizeCellVectors

include("VisualizePDE.jl")

end # module


text:
# Tue 13 Nov 2012 04:13:36 PM EST
#
# NeHe Tut 10 - Move around in a 3D world
#
# NOTE: This example does not run slowly.  Rather, the key repeat rate limits
# the speed at which the "player's" position is updated, giving the illusion of
# a slow example.  Try lowering the key repeat interval to see it run at
# "normal" speed.
#
#
# Q - quit
# B - turn texture alpha-blending on/off
# L - turn lights on/off
# F - change texture filter (linear, nearest, mipmap)
# PageUp/Down - look up/down
# Up/Down - move forward/backward
# Left/Right - turn left/right

# load necessary GL/SDL routines

using OpenGL
@OpenGL.version "1.0"
@OpenGL.load
using SDL

# initialize variables

bpp                = 16
wintitle           = "NeHe Tut 10"
icontitle          = "NeHe Tut 10"
width              = 640
height             = 480

walkbias           = 0.0
walkbiasangle      = 0.0

lookupdown         = 0.0

xpos               = 0.0
zpos               = 0.0

yrot               = 0.0

LightAmbient       = [0.5f0, 0.5f0, 0.5f0, 1.0f0]
LightDiffuse       = [1.0f0, 1.0f0, 1.0f0, 1.0f0]
LightPosition      = [0.0f0, 0.0f0, 2.0f0, 1.0f0]

filter             = 3
light              = true
blend              = false

keystate_checked   = false
lastkeycheckTime   = 0
key_repeatinterval = 75 #ms

# open SDL window with an OpenGL context

SDL_Init(SDL_INIT_VIDEO)
#videoInfo = SDL_GetVideoInfo()
videoFlags = (SDL_OPENGL | SDL_GL_DOUBLEBUFFER | SDL_HWPALETTE | SDL_RESIZABLE)
#if videoInfo.hw_available
    videoFlags = (videoFlags | SDL_HWSURFACE)
#else
    #videoFlags = (videoFlags | SDL_SWSURFACE)
#end
#if videoInfo.blit_hw
    videoFlags = (videoFlags | SDL_HWACCEL)
#end
SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1)
SDL_SetVideoMode(width, height, bpp, videoFlags)
SDL_WM_SetCaption(wintitle, icontitle)

glViewport(0, 0, width, height)
glClearColor(0.0, 0.0, 0.0, 0.0)
glClearDepth(1.0)
glDepthFunc(GL_LESS)
glEnable(GL_DEPTH_TEST)
glShadeModel(GL_SMOOTH)

glMatrixMode(GL_PROJECTION)
glLoadIdentity()

gluPerspective(45.0,width/height,0.1,100.0)

glMatrixMode(GL_MODELVIEW)

### auxiliary functions

# initialize global variables

global numtriangles = 0

function SetupWorld(world_map::String)

    global numtriangles

    filein       = open(world_map)
    world_data   = readlines(filein)

    numtriangles = parseint(chomp(split(world_data[1],' ')[2]))

    sector       = zeros(numtriangles,3,5)

    loop = 1
    vert = 1
    line = 1

    while line <= length(world_data)-2
        if world_data[2+line][1] != '/' && world_data[2+line][1] != '\n'
            while vert <= 3
                (x, y, z, u, v)      = split(chomp(world_data[2+line]),' ')
                x                    = parsefloat(x)
                y                    = parsefloat(y)
                z                    = parsefloat(z)
                u                    = parsefloat(u)
                v                    = parsefloat(v)
                sector[loop,vert,:]  = [x,y,z,u,v]
                vert                 += 1
                line                 += 1
            end
            vert = 1
            loop += 1
        else
            line += 1
        end
    end

    return sector

end

### end of auxiliary functions

# initialize sector1 with SetupWorld

sector1 = SetupWorld(expanduser("~/.julia/SDL/Examples/NeHe/tut10/world.txt"))

# load textures from images

tex   = Array(Uint32,3) # generating 3 textures

img, w, h = glimread(expanduser("~/.julia/SDL/Examples/NeHe/tut10/mud.bmp"))

glGenTextures(3,tex)
glBindTexture(GL_TEXTURE_2D,tex[1])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)

glBindTexture(GL_TEXTURE_2D,tex[2])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)

glBindTexture(GL_TEXTURE_2D,tex[3])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)

gluBuild2DMipmaps(GL_TEXTURE_2D, 3, w, h, GL_RGB, GL_UNSIGNED_BYTE, img)

# initialize lights

glLightfv(GL_LIGHT1, GL_AMBIENT, LightAmbient)
glLightfv(GL_LIGHT1, GL_DIFFUSE, LightDiffuse)
glLightfv(GL_LIGHT1, GL_POSITION, LightPosition)

glEnable(GL_LIGHT1)
glEnable(GL_LIGHTING)

# enable texture mapping and setup alpha blending

glEnable(GL_TEXTURE_2D)
glBlendFunc(GL_SRC_ALPHA, GL_ONE)

# main drawing loop

while true
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    xtrans = -xpos
    ztrans = -zpos
    ytrans = -walkbias-0.25
    sceneroty = 360.0-yrot

    glRotate(lookupdown, 1.0, 0.0, 0.0)
    glRotate(sceneroty, 0.0, 1.0, 0.0)
    glTranslate(xtrans, ytrans, ztrans)

    glBindTexture(GL_TEXTURE_2D,tex[filter])

    for face = 1:numtriangles
        glBegin(GL_TRIANGLES)
            glNormal(0.0, 0.0, 1.0)
            x_m = sector1[face,1,1]
            y_m = sector1[face,1,2]
            z_m = sector1[face,1,3]
            u_m = sector1[face,1,4]
            v_m = sector1[face,1,5]
            glTexCoord(u_m,v_m)
            glVertex(x_m,y_m,z_m)

            x_m = sector1[face,2,1]
            y_m = sector1[face,2,2]
            z_m = sector1[face,2,3]
            u_m = sector1[face,2,4]
            v_m = sector1[face,2,5]
            glTexCoord(u_m,v_m)
            glVertex(x_m,y_m,z_m)

            x_m = sector1[face,3,1]
            y_m = sector1[face,3,2]
            z_m = sector1[face,3,3]
            u_m = sector1[face,3,4]
            v_m = sector1[face,3,5]
            glTexCoord(u_m,v_m)
            glVertex(x_m,y_m,z_m)
        glEnd()
    end

    SDL_GL_SwapBuffers()

    SDL_PumpEvents()
    if SDL_GetTicks() - lastkeycheckTime >= key_repeatinterval
        keystate         = SDL_GetKeyState()
        keystate_checked = true
        lastkeycheckTime = SDL_GetTicks()
    end

    # Sampling rates for event processing are so fast that a single key press
    # lasts through several iterations of this loop.  This means that one press
    # can be seen as 50 or more presses by the SDL event system.  To correct
    # for this, we only check keypresses every 75ms.

    if keystate_checked == true
        if keystate[SDLK_q] == true
            break
        elseif keystate[SDLK_b] == true
            println("Blend was: $blend")
            blend = (blend ? false : true)
            if blend
                glEnable(GL_BLEND)
                glDisable(GL_DEPTH_TEST)
            else
                glDisable(GL_BLEND)
                glEnable(GL_DEPTH_TEST)
            end
            println("Blend is now: $blend")
        elseif keystate[SDLK_l] == true
            println("Light was: $light")
            light = (light ? false : true)
            println("Light is now: $light")
            if light
                glEnable(GL_LIGHTING)
            else
                glDisable(GL_LIGHTING)
            end
        elseif keystate[SDLK_f] == true
            println("Filter was: $filter")
            filter +=1
            if filter > 3
                filter = 1
            end
            println("Filter is now: $filter")
        elseif keystate[SDLK_PAGEUP] == true
            lookupdown -=0.2
        elseif keystate[SDLK_PAGEDOWN] == true
            lookupdown +=1.0
        elseif keystate[SDLK_UP] == true
            xpos -=sin(degrees2radians(yrot))*0.05
            zpos -=cos(degrees2radians(yrot))*0.05
            walkbias +=10
            if walkbiasangle <= 359.0
                walkbiasangle =0.0
            else
                walkbiasangle +=10
            end
            walkbias = sin(degrees2radians(walkbiasangle))/20.0
        elseif keystate[SDLK_DOWN] == true
            xpos +=sin(degrees2radians(yrot))*0.05
            zpos +=cos(degrees2radians(yrot))*0.05
            walkbias -=10
            if walkbiasangle <= 1.0
                walkbiasangle =359.0
            else
                walkbiasangle -=10
            end
            walkbias = sin(degrees2radians(walkbiasangle))/20.0
        elseif keystate[SDLK_LEFT] == true
            yrot +=1.5
        elseif keystate[SDLK_RIGHT] == true
            yrot -=1.5
        end
        keystate_checked = false
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: marketplace_entitlement_service
using AWS.Compat
using AWS.UUIDs

"""
    get_entitlements(product_code)
    get_entitlements(product_code, params::Dict{String,<:Any})

GetEntitlements retrieves entitlement values for a given product. The results can be
filtered based on customer identifier or product dimensions.

# Arguments
- `product_code`: Product code is used to uniquely identify a product in AWS Marketplace.
  The product code will be provided by AWS Marketplace when the product listing is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filter"`: Filter is used to return entitlements for a specific customer or for a
  specific dimension. Filters are described as keys mapped to a lists of values. Filtered
  requests are unioned for each value in the value list, and then intersected for each filter
  key.
- `"MaxResults"`: The maximum number of items to retrieve from the GetEntitlements
  operation. For pagination, use the NextToken field in subsequent calls to GetEntitlements.
- `"NextToken"`: For paginated calls to GetEntitlements, pass the NextToken from the
  previous GetEntitlementsResult.
"""
get_entitlements(ProductCode; aws_config::AbstractAWSConfig=global_aws_config()) = marketplace_entitlement_service("GetEntitlements", Dict{String, Any}("ProductCode"=>ProductCode); aws_config=aws_config)
get_entitlements(ProductCode, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = marketplace_entitlement_service("GetEntitlements", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ProductCode"=>ProductCode), params)); aws_config=aws_config)


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	configure_scip(solver_settings_path::String)

Reads user-specified solver settings from scip\_settings.yml in the directory specified by the string solver\_settings\_path.

Returns a MathOptInterface OptimizerWithAttributes SCIP optimizer instance to be used in the GenX.generate_model() method.

The SCIP optimizer instance is configured with the following default parameters if a user-specified parameter for each respective field is not provided:

 - Dispverblevel = 0
 - limitsgap = 0.05

"""
function configure_scip(solver_settings_path::String)

    solver_settings = YAML.load(open(solver_settings_path))

    # Optional solver parameters ############################################
    display_verblevel = 0
    if (haskey(solver_settings, "display_verblevel"))
        display_verblevel = solver_settings["display_verblevel"]
    end

    limits_gap = 0.05
    if (haskey(solver_settings, "limits_gap"))
        limits_gap = solver_settings["limits_gap"]
    end
    ########################################################################

    OPTIMIZER = optimizer_with_attributes(
        SCIP.Optimizer,
        "display_verblevel" => display_verblevel,
        "limits_gap" => limits_gap
    )

    return OPTIMIZER
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
push!(LOAD_PATH, ".")

module Parallel

using Random
using Future

const GLOBAL_THREAD_RNG = [MersenneTwister(rand(1:(1<<30))) for i=1:Threads.nthreads()];

function threadRand(args...)
    return rand(GLOBAL_THREAD_RNG[Threads.threadid()], args...)
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
using Clustering
using Test
using CodecZlib
using Distances

@testset "hclust() (hierarchical clustering)" begin

@testset "param checks" begin
    D = rand(5, 5)
    Dsym = D + D'
    Dnan = copy(Dsym)
    Dnan[1, 3] = Dnan[3, 1] = NaN
    @testset "hclust()" begin
        @test_throws ArgumentError hclust(Dsym, linkage=:typo)
        @test_throws ArgumentError hclust(D, linkage=:single)
        @test_throws ArgumentError hclust(Dnan, linkage=:single)
    end
    hclu = @inferred(hclust(Dsym, linkage=:single))
    @test hclu isa Clustering.Hclust
    @testset "cutree()" begin
        @test_throws ArgumentError cutree(hclu)
        @test_throws ArgumentError cutree(hclu, k=-1)
        @test_throws ArgumentError cutree(hclu, k=0)
        @test cutree(hclust(Dsym), k=10) isa Vector{Int}
        @test cutree(hclust(fill(0.0, 0, 0)), k=0) == Int[]
        @test cutree(hclust(fill(0.0, 0, 0)), k=1) == Int[]
    end
end

@testset "R hclust() generated examples" begin
# load the examples
hclu_examples_filename = joinpath(@__DIR__, "data", "hclust_generated_examples.jl.gz")
Base.include_string(@__MODULE__,
                    open(io -> read(io, String), GzipDecompressorStream, hclu_examples_filename),
                    hclu_examples_filename)

# test to make sure many random examples match R's implementation
@testset "example #$i (linkage=:$(example["linkage"]), n=$(size(example["D"], 1)))" for
        (i, example) in enumerate(examples)

    linkage = example["linkage"]
    hclu = @inferred(hclust(example["D"], linkage=linkage))
    @test hclu isa Clustering.Hclust
    @test Clustering.nnodes(hclu) == size(example["D"], 1)
    @test Clustering.nmerges(hclu) == Clustering.nnodes(hclu)-1
    @test Clustering.height(hclu) ≈ maximum(example["height"]) atol=1e-5
    @test hclu.merges == example["merge"]
    @test hclu.heights ≈ example["height"] atol=1e-5
    @test hclu.order == example["order"]

    # compare hclust_nn_lw() (the default) and hclust_nn() (slower) methods
    if linkage ∈ [:complete, :average]
        @testset "hclust_nn()" begin
            linkage_fun = linkage == :complete ? Clustering.slicemaximum : Clustering.slicemean
            hclu2 = Hclust(Clustering.orderbranches_r!(Clustering.hclust_nn(example["D"], linkage_fun)),
                linkage)
            @test hclu2.merges == hclu.merges
            @test hclu2.heights ≈ hclu.heights atol=1e-5
            @test hclu2.order == hclu.order
        end
    end

    local cut_k = example["cut_k"]
    local cut_h = example["cut_h"]
    if cut_h !== nothing
        # due to small arithmetic differences between R and Julia heights might be slightly different
        # find the matching height
        cut_h_r = cut_h
        cut_h_ix = findmin(abs.(hclu.heights .- cut_h_r))[2]
        cut_h = hclu.heights[cut_h_ix]
        @assert isapprox(cut_h, cut_h_r, atol=1e-6) "h=$cut_h ≈ h_R=$cut_h_r"
    end
    @testset "cutree(hclu, k=$(repr(cut_k)), h=$(repr(cut_h)))" begin
        cutt = @inferred(cutree(hclu, k=cut_k, h=cut_h))
        @test cutt isa Vector{Int}
        @test length(cutt) == Clustering.nnodes(hclu)
        @test cutt == example["cutree"]
    end
end
end

@testset "hclust_n3()" begin
    # no thorough testing (it's O(N³)), just test one example
    example_n3 = examples[10]
    hclu_n3 = @inferred(Clustering.orderbranches_r!(Clustering.hclust_n3(example_n3["D"], Clustering.slicemaximum)))
    @test hclu_n3.mleft == example_n3["merge"][:, 1]
    @test hclu_n3.mright == example_n3["merge"][:, 2]
    @test hclu_n3.heights ≈ example_n3["height"] atol=1e-5
end

local hclust_linkages = [:single, :average, :complete]

@testset "hclust(0×0 matrix, linkage=$linkage)" for linkage in hclust_linkages
    hclu = hclust(fill(0.0, 0, 0), linkage=linkage)
    @test Clustering.nnodes(hclu) == 0
    @test Clustering.nmerges(hclu) == 0
    @test Clustering.height(hclu) == -Inf
    cut1 = @inferred(cutree(hclu, h=1))
    @test cut1 == Int[]
end

@testset "hclust([$dist] 1×1 matrix, linkage=$linkage)" for
    linkage in hclust_linkages, dist in [-Inf, 0, 1, 2, Inf]
    hclu = hclust(fill(dist, 1, 1), linkage=linkage)
    @test Clustering.nnodes(hclu) == 1
    @test Clustering.nmerges(hclu) == 0
    @test Clustering.height(hclu) == -Inf
    cut1 = @inferred(cutree(hclu, h=1))
    @test cut1 == [1]
end

@testset "hclust(linkage=$linkage) when data contains an isolated point (#109)" for linkage in hclust_linkages
    # point #2 is isolated: distances to all the other points are Inf
    mdist = [
        0.0 Inf 0.1  Inf  Inf  Inf  Inf  Inf  Inf  Inf;
        Inf 0.0 Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf;
        0.1 Inf 0.0  0.11 Inf  Inf  Inf  Inf  Inf  Inf;
        Inf Inf 0.11 0.0  0.86 0.86 Inf  Inf  Inf  Inf;
        Inf Inf Inf  0.86 0.0  0.67 0.72 1.93 Inf  Inf;
        Inf Inf Inf  0.86 0.67 0.0  0.72 Inf  Inf  Inf;
        Inf Inf Inf  Inf  0.72 0.72 0.0  Inf  0.72 Inf;
        Inf Inf Inf  Inf  1.93 Inf  Inf  0.0  3.91 Inf;
        Inf Inf Inf  Inf  Inf  Inf  0.72 3.91 0.0  0.52;
        Inf Inf Inf  Inf  Inf  Inf  Inf  Inf  0.52 0.0]

    hclu = hclust(mdist, linkage=linkage)
    @test Clustering.nnodes(hclu) == 10
    @test Clustering.nmerges(hclu) == 9
    @test Clustering.height(hclu) == Inf
    cut1 = @inferred(cutree(hclu, h=1))
    @test cut1 isa Vector{Int}
    @test length(cut1) == Clustering.nnodes(hclu)
    if linkage == :single
        @test cut1 == [1, 2, 1, 1, 1, 1, 1, 3, 1, 1]
    end
end

@testset "cutree(hclust, h=h) when the height of all subtrees greater than h (#141)" begin
    A = [0.0 0.7; 0.7 0.0]
    hA = hclust(A, linkage=:average)
    @test cutree(hA, h=0.5) == [1, 2]
    @test cutree(hA, h=0.7) == [1, 1]
    @test cutree(hA, h=0.9) == [1, 1]
end

@testset "Leaf ordering methods" begin
    n = 10
    mat = zeros(Int, n, n)

    for i in 1:n
        last = minimum([i+Int(floor(n/5)), n])
        for j in i:last
            mat[i,j] = 1
        end
    end

    dm = pairwise(Euclidean(), mat, dims=2)

    hcl_r = hclust(dm, linkage=:average)
    hcl_barjoseph = hclust(dm, linkage=:average, branchorder=:barjoseph)
    hcl_optimal = hclust(dm, linkage=:average, branchorder=:optimal)

    @test hcl_r.order == [3, 1, 2, 4, 5, 9, 10, 6, 7, 8]
    @test hcl_r.merges == [-1 -2; -3 1; -4 -5; -9 -10; -7 -8; -6 5; 2 3; 4 6; 7 8]

    @test hcl_barjoseph.order == collect(1:10)
    @test hcl_barjoseph.merges == [-1 -2; 1 -3; -4 -5; -9 -10; -7 -8; -6 5; 2 3; 6 4; 7 8]

    @test hcl_barjoseph.merges == hcl_optimal.merges
    @test hcl_barjoseph.order == hcl_optimal.order

    @test_throws ArgumentError hclust(dm, linkage=:average, branchorder=:wrong)

    hcl_zero = hclust(fill(0.0, 0, 0), linkage=:average, branchorder=:barjoseph)
    @test Clustering.nnodes(hcl_zero) == 0

    hcl_one = hclust(fill(0.0, 1, 1), linkage=:average, branchorder=:barjoseph)
    @test Clustering.nnodes(hcl_one) == 1

    # Larger matrix to make sure all swaps are tested
    Random.seed!(1)
    D = rand(50,50)
    Dm = D + D'
    hcl_rand = hclust(Dm, linkage=:average, branchorder=:optimal)

    @test hcl_rand.merges == [-29 -1; -32 -24; -46 -44; -10 -41; -17 -12; -40 5;
                              -8 -28; -13 -35; -19 -20; -43 -42; -34 -18; -15 -39;
                              -30 -49; -26 -22; -36 -31; -38 -4; -5 -9; 1 16; -6 13;
                              17 10; -33 -45; -7 -2; -23 -21; -48 -27; -37 -16; -14 8;
                              2 -50; 19 20; -47 -11; 9 -3; 6 15; 14 22; 31 18; 23 25;
                              11 3; 32 29; 33 30; 7 12; 35 38; -25 4; 21 28; 26 34;
                              41 27; 24 40; 39 42; 37 43; 45 46; 47 36; 48 44]

    @test hcl_rand.order == [34, 18, 46, 44, 8, 28, 15, 39, 14, 13, 35, 23, 21, 37,
                             16, 40, 17, 12, 36, 31, 29, 1, 38, 4, 19, 20, 3, 33, 45,
                             6, 30, 49, 5, 9, 43, 42, 32, 24, 50, 26, 22, 7, 2, 47,
                             11, 48, 27, 25, 10, 41]

    @test Clustering.nnodes(hcl_rand) == 50
end

end # testset "hclust()"


text:
function mkcgent2dgent(dgnodes,tol)

dig = Int(log10(1/tol));
dgnodes = round.(dgnodes,digits=dig);

# CG element-to-node connectivities
cgnodes,cgelcon = mkelconcg(dgnodes);

# CG node-to-element connectivities
rowent2elem,colent2elem,~,~,nent = mkent2elem(cgelcon);

cgent2dgent = colent2elem;
npe = size(dgnodes,1);
dim = size(dgnodes,2);
#nent = maximum(cgelcon[:]);

# display(size(dgnodes))
# display(size(cgnodes))
# display(nent)
# display(size(rowent2elem))
for i = 1:nent # for CG node i
    nelem = rowent2elem[i+1]-rowent2elem[i]; # number of elements connected to CG node i
    elem = colent2elem[(rowent2elem[i]+1):rowent2elem[i+1]]; # elements connected to CG node i
    xcg = reshape(cgnodes[:,i],1,dim); # coordinates of CG node i
    for j = 1:nelem # for element j connected to CG node i
        xdg = dgnodes[:,:,elem[j]]; # coordinates of DG nodes on element j
        in = xiny(xcg,xdg,1); # match CG node i to one of the DG nodes on element j
        cgent2dgent[rowent2elem[i]+j] = (elem[j]-1)*npe+in[1]; # index of the matched DG node
    end
end

return cgelcon,rowent2elem,colent2elem,cgent2dgent,cgnodes

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using FFTW

function periodogram(y,f_range)
    # The periodogram estimator computes the estimate of the spectrum of a
    # signal based on multiple samples in time. For line spectral estimation,
    # the estimates of the frequencies of the sinusoids can be found by
    # searching for local maxima in the spectrum.
    #
    # Inputs:
    # y: measurements
    # f_range: sample point of the spectrum \in [0,1)
    #
    # Output:
    # spec: estimated spectrum

    y=y[:,1]
    N = length(y)
    spec(f) = 1/N*abs( sum(y[i]*exp(-1im*2*pi*f*i) for i = 1:N ))^2
    spec = spec.(f_range)/maximum(spec.(f_range))
end


text:
@eval using Zygote: unbroadcast, Numeric, @adjoint

@eval begin

# lenseflow

@adjoint function *(Lϕ::LenseFlowOp, f::Field{B}) where {B}
    f̃ = Lϕ * f
    function back(Δ)
        δf, δϕ = δf̃ϕ_δfϕ(Lϕ,f̃,f)' * FΦTuple(Δ,zero(getϕ(Lϕ)))
        LenseFlow(δϕ), B(δf)
    end
    f̃, back
end

@adjoint function \(Lϕ::LenseFlowOp, f̃::Field{B}) where {B}
    f = Lϕ \ f̃
    function back(Δ)
        δf, δϕ = δfϕ_δf̃ϕ(Lϕ,f,f̃)' * FΦTuple(Δ,zero(getϕ(Lϕ)))
        LenseFlow(δϕ), B(δf)
    end
    f, back
end

@adjoint (::Type{L})(ϕ) where {L<:LenseFlowOp} = L(ϕ), Δ -> (Δ.ϕ,)
@adjoint cache(Lϕ::LenseFlowOp, f) = cache(Lϕ,f), Δ->(Δ,nothing)
@adjoint cache!(cL::CachedLenseFlow, ϕ) where {L<:LenseFlowOp} = cache!(cL,ϕ), Δ -> (nothing,Δ.ϕ)

# algebra

@adjoint (::Type{B})(f::Field{B′}) where {B<:Basis, B′} = B(f), Δ -> (B′(Δ),)
@adjoint (*)(L::LinOp, f::Field{B}) where {B} = L*f, Δ -> (nothing, B(L'*Δ))
@adjoint (\)(L::LinOp, f::Field{B}) where {B} = L\f, Δ -> (nothing, B(L'\Δ))
@adjoint (*)(a::Adjoint{<:Any,<:Field{B1}}, b::Field{B2}) where {B1,B2} = a*b, Δ -> (B1(Δ*b)',  B2(Δ*a'))
@adjoint function (*)(x::Adjoint{<:Any,<:Field{B1}}, D::Union{Diagonal,OuterProdOp}, y::Field{B2}) where {B1,B2}
    z = x*D*y
    back = if parent(x)===y
        function (Δ)
            g = B1(Δ*(D*y))
            (g', Δ*z, g)
        end
    else
        function (Δ)
            (B1(Δ*(D*y))', Δ*z, B2(Δ*(D*x')))
        end
    end
    z, back
end

# eventually we need to implement these to allow gradients w.r.t. θ:
@adjoint pinv(D::LinOp) = pinv(D), Δ->nothing
@adjoint logdet(L::LinOp, θ) = logdet(L,θ), Δ->nothing
@adjoint (ds::DataSet)(args...; kwargs...) = ds(args...; kwargs...), Δ -> nothing

# some stuff which arguably belongs in Zygote or ChainRules
# see also: https://github.com/FluxML/Zygote.jl/issues/316
@adjoint *(D::Diagonal, v::AbstractVector) = D.diag .* v,
    Δ -> (Diagonal(unbroadcast(D.diag, Δ .* conj.(v))), unbroadcast(v, Δ .* conj.(D.diag)))

@adjoint broadcasted(::typeof(-), x ::Numeric, y::Numeric) =
    broadcast(-, x, y), Δ -> (nothing, unbroadcast(x, Δ), unbroadcast(y, -Δ))

@adjoint broadcasted(::typeof(\), x ::Numeric, y::Numeric) =
    broadcast(\, x, y), Δ -> (nothing, unbroadcast(x, @. -Δ*y/x^2), unbroadcast(y, @. Δ/x))

end


text:
push!(LOAD_PATH,"../src/")
using TopologicalDataAnalysis
using Documenter

makedocs(
         sitename = "TopologicalDataAnalysis.jl",
         modules  = [TopologicalDataAnalysis],
         pages=[
                "Home" => "index.md"
               ])

deploydocs(;
    repo="github.com/davibarreira/TopologicalDataAnalysis.jl",
)


text:
function plotIRF(estimate, H)
	irfs = Psi(estimate, H)
	n, n, H = size(irfs)
	fig = figure("pyplot_subplot_mixed",figsize=(10,10))
	for j = 1:n
		for i = 1:n
			subplot(n,n,(i-1)*n+j)
			grid("on")
			plot(vec(irfs[i,j,:]))
		end
	end
	fig[:canvas][:draw]()
	fig.close()
end


text:
import PAPI

function slow_code()
    tmp = 1.1
    for i=1:200_000_000
        tmp = (tmp+100)/i
    end
    return tmp
end

function main()
    precompile(slow_code, ())

    flpins, mflips, real_time, proc_time = @PAPI.flips slow_code()

    @printf("Real_time: %f Proc_time: %f Total flpins: %lld MFLIPS: %f\n",
             real_time, proc_time, flpins, mflips)
end

main()


text:
using Seismic, PyCall
@pyimport matplotlib.pyplot as plt

nt = 256;
nx = convert(Int64,128);
h = Array(Header,nx);
for ix = 1:nx
  h[ix] = Seismic.InitSeisHeader();
  h[ix].tracenum = ix;
  h[ix].d1 = 0.004;
  h[ix].n1 = nt;
  h[ix].imx = ix;
  h[ix].imy = 0;
  h[ix].ihx = 0;
  h[ix].ihy = 0;
end

dt=0.004;
dx=25;
f_niquest=Int(floor(1/dt))
param=["dt"=>dt, "dx"=>dx]

data=rand(nt,nx)

f1 = plt.figure(1)
plt.imshow(data[1:nt,:],cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
println("1")
plt.title("data before fk filter")
println("1")
plt.xlabel("X (m)")
println("1")
plt.ylabel("T (s)")
println("1")
plt.show()
println("1")

m=fft(data,1)/sqrt(size(data,1));
m=fft(m,2)/sqrt(size(data,2))
m=fftshift(m)

m=m[Int(floor(f_niquest/2)):f_niquest,:]

println("2")

f11 = plt.figure(11)
plt.imshow(abs(m),cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("fk spectrum before filter")
plt.xlabel("K (1/m)")
plt.ylabel("W (1/s)")
plt.show()

println("3")

# apply FK filter
data_fk,fk_spec,h = SeisFKFilter(data,h,param)

println("4")

f2 = plt.figure(2)
plt.imshow(data_fk[1:nt,:],cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("data after fk filter")
plt.xlabel("X (m)")
plt.ylabel("T (s)")
plt.show()

println("5")

f21 = plt.figure(21)
plt.imshow(fk_spec,cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("fk spectrum after filter")
plt.xlabel("K (1/m)")
plt.ylabel("W (1/s)")
plt.show()


text:
###############################################################################
#
# We already illustrated the use of Stopping for optimization algorithm,
# however, in the case where one algorithm/solver is not Stopping-compatible,
# a buffer solver is required to unify the formalism.
# We illustrate this situation here with the Ipopt solver.
#
# Remark in the buffer function: in case the solver stops with success
# but the stopping condition is not satisfied, one option is to iterate
# and reduce the various tolerances.
#
# Documentation for Ipopt options can be found here:
# https://coin-or.github.io/Ipopt/OPTIONS.html#OPTIONS_REF
###############################################################################

using Ipopt, NLPModels, NLPModelsIpopt, Stopping

include("../test-stopping/rosenbrock.jl")
x0  = 1.5 * ones(6)
nlp = ADNLPModel(rosenbrock,  x0)

#The traditional way to solve an optimization problem using NLPModelsIpopt
#https://github.com/JuliaSmoothOptimizers/NLPModelsIpopt.jl
printstyled("Oth scenario:\n")

stats = ipopt(nlp, print_level = 0, x0 = x0)
#Use y0 (general), zL (lower bound), zU (upper bound)
#for initial guess of Lagrange multipliers.
@show stats.solution, stats.status

#Using Stopping, the idea is to create a buffer function
function solveIpopt(stp :: NLPStopping)

 #xk = solveIpopt(stop.pb, stop.current_state.x)
 stats = ipopt(nlp, print_level     = 0,
                    tol             = stp.meta.rtol,
                    x0              = stp.current_state.x,
                    max_iter        = stp.meta.max_iter,
                    max_cpu_time    = stp.meta.max_time,
                    dual_inf_tol    = stp.meta.atol,
                    constr_viol_tol = stp.meta.atol,
                    compl_inf_tol   = stp.meta.atol)

 #Update the meta boolean with the output message
 if stats.status == :first_order stp.meta.suboptimal      = true end
 if stats.status == :acceptable  stp.meta.suboptimal      = true end
 if stats.status == :infeasible  stp.meta.infeasible      = true end
 if stats.status == :small_step  stp.meta.stalled         = true end
 if stats.status == :max_iter    stp.meta.iteration_limit = true end
 if stats.status == :max_time    stp.meta.tired           = true end

 stp.meta.nb_of_stop = stats.iter
 #stats.elapsed_time

 x = stats.solution

 #Not mandatory, but in case some entries of the State are used to stop
 fill_in!(stp, x)

 stop!(stp)

 return stp
end

nlp_at_x = NLPAtX(x0)
stop = NLPStopping(nlp, unconstrained_check, nlp_at_x)

#1st scenario, we solve again the problem with the buffer solver
printstyled("1st scenario:\n")
solveIpopt(stop)
@show stop.current_state.x, status(stop)
nbiter = stop.meta.nb_of_stop

#2nd scenario: we check that we control the maximum iterations.
printstyled("2nd scenario:\n")
#rstate is set as true to allow reinit! modifying the State
reinit!(stop, rstate = true, x = x0)
stop.meta.max_iter = max(nbiter-4,1)

solveIpopt(stop)
#Final status is :IterationLimit
@show stop.current_state.x, status(stop)

printstyled("The End.\n")


text:
using AbstractGPs

using Distributions: Poisson, LogNormal
using SampleChainsDynamicHMC
using Soss: Soss, For
using Test
using Turing: Turing

@testset "PPLs" begin
    include("turing.jl")
    println(" ")
    @info "Ran Turing tests"

    include("soss.jl")
    println(" ")
    @info "Ran Soss tests"
end


text:
# resample(uvals::Vector{AbstractUncertainValue})
@test resample(uvals) isa Vector
@test length(resample(uvals)) == length(uvals)

# resample(uvals::Vector{AbstractUncertainValue}, n::Int)
n = 2
@test resample(uvals, n) isa Vector{Vector{T}} where T <: Real
@test length(resample(uvals, n)) == n

# resample_elwise(uvals::Vector{AbstractUncertainValue}, n::Int)
n = 2
@test resample_elwise(uvals, n) isa Vector{Vector}
@test length(resample_elwise(uvals, n)) == length(uvals)
@test length(resample_elwise(uvals, n)[1]) == n


text:
"""
Convert from a hypergraph to a dyadic graph by turning hyperedges into cliques.
Convert between hypergraph and bipartite representations.
"""

using MatrixNetworks
using SparseArrays
using LinearAlgebra

"""
`to_adjacency_matrix`
=====================

Convert a hypergraph to an adjacency matrix

Arguments
----------
    - `H::Hypergraphs`: The hypergraph that will be converted into an adjacency matrix
    - `simple::Bool (=false)`: Whether to exclude multi-edges and self-loops from the matrix (they are included by default)

Examples
----------
~~~~
to_adjacency_matrix(Hypergraphs([[1,2,3,4], [3,4,1], [1,4,5]], ["One", "Two", "Three", "Four", "Five"], 5, 3); simple=true)
~~~~
"""
function to_adjacency_matrix(H::Hypergraph; simple::Bool=false)
    A = zeros(Int64, H.n, H.n) # Create adjacency matrix of size n x n

    # Loop through each hyperedge and create edges in the new adjacency matrix
    for i = 1:H.m
        for j = 1:length(H.edges[i])-1
            for k = j+1:length(H.edges[i])
                index1 = H.edges[i][j]
                index2 = H.edges[i][k]
                A[index1,index2] += 1
                A[index2,index1] += 1
            end
        end
    end

    if simple
        A -= Diagonal(A) # Remove self-loops
        A = min.(A, 1) # Remove multi-edges
    end

    return A
end

"""
`hypergraph_to_matrixnetworks`
=====================

Convert a hypergraph to a MatrixNetwork

Arguments
----------
    - `H::Hypergraphs`: The hypergraph that will be converted into a MatrixNetwork
    - `simple::Bool (=false)`: Whether to exclude multi-edges and self-loops from the MatrixNetwork (they are included by default)

Examples
----------
~~~~
hypergraph_to_matrixnetwork(Hypergraphs([[1,2,3,4], [3,4,1], [1,4,5]], ["One", "Two", "Three", "Four", "Five"], 5, 3); simple=true)
~~~~
"""
function hypergraph_to_matrixnetwork(H::Hypergraph; simple=false)
    return MatrixNetwork(sparse(to_adjacency_matrix(H; simple=simple)))
end

"""
`hyper_to_bipartite`
====================

Converts a hypergraph into its bipartite representation. For all vertices v and
edges e in the hypergraph, an edge (v, e) exists in the bipartite if and only if
v ∈ e in the hypergraph.
"""
function hyper_to_bipartite(H::Hypergraph)
   di_edges::Vector{Vector{Int64}} = []

   for i = 1:h.m
      for j in h.edges[i]
         push!(di_edges, [i + h.n, j])
      end
   end

   return MatrixNetwork(di_edges, h.n + h.m)
end

"""
`is_bipartite`
==============

Checks whether a graph is bipartite by checking whether it is 2-colorable.
"""
function is_bipartite(A::MatrixNetwork)
   group1::Set{Int64} = Set()
   group2::Set{Int64} = Set()

   function color_neighbors(v::Int64, is_group1::Bool)
      for n in A.ci[A.rp[v]:A.rp[v+1]-1]
         if is_group1 # v is in group 2, so all neighbors should be in group 1
            if n in group2
               return false
            end
            push!(group1, n)
         else
            if n in group1
               return false
            end
            push!(group2, n)
         end
      end
      return true
   end

   is_group1 = true

   for i = 1:A.n
      is_group1 = !is_group1
      if !color_neighbors(i, is_group1)
         return false
      end
   end
   return (group1, group2)
end

"""
`bipartite_to_hyper`
====================

Returns the left and right projections of a bipartite graph into a hypergraph.
In the left projection, all vertices in group two that are adjacent to the same
vertex in group one are put together in a hyperedge. The right projection is the
same except groups one and two are swapped in their roles.
"""
function bipartite_to_hyper(A::MatrixNetwork)
   b = is_bipartite(A)
   if typeof(b) <: Bool
      error("ArgumentError: The network must be bipartite to be projected to a hypergraph.")
   end

   g1, g2 = collect.(b)
   left_edges = [map(x -> findfirst(y -> x == y, g2), A.ci[A.rp[i]:A.rp[i+1]-1]) for i in g1]
   right_edges = [map(x -> findfirst(y -> x == y, g1), A.ci[A.rp[j]:A.rp[j+1]-1]) for j in g2]

   return Hypergraphs(left_edges, A.vals[g2], length(g2), length(left_edges)),
          Hypergraphs(right_edges, A.vals[g1], length(g1), length(right_edges))
end


text:
abstract type Summarizer end

"""
    update!(s::Summarizer, nt::NamdedTuple)

Update the summary statistics with a named tuple of one row of data.
"""
function update!(::Summarizer, nt) end

"""
    get(s::Summarizer)

Extract the summary statistics from a summarizer.  Returns a stat-value `Pair`.
"""
function Base.get(::Summarizer) end

"""
    name(s::Summarizer)

The name of the variable(s) summarized.  Defaults to the `name` field.
"""
name(s::Summarizer) = s.name


text:
# function discrete_value_estimate(π, px, s)
#         xs = support(px)
#         pdfs = exp.(logpdf(px, s, xs))
#         sum(value(π, s) .* pdfs, dims=1)
# end
#
# function CERL_DQN_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
#         px = 𝒫[:px]
#         term = 𝒟[:done] .| 𝒟[:fail]
#         term .* 𝒟[:fail] .+ (1.f0 .- term) .* (discrete_value_estimate(π, px, 𝒟[:sp]) ./ discrete_value_estimate(π, px, 𝒟[:s]))
# end
#
# function IS_estimate_log_pfail(π, px, s, Nsamples=10)
#         xs_and_logpdfs = [exploration(π, s) for _ in 1:Nsamples]
#         xs = [e[1] for e in xs_and_logpdfs]
#         νlogpdfs = vcat([e[2] for e in xs_and_logpdfs]...)
#         logpfails = vcat([value(π, s, x) for x in xs]...)
#         logpxs = Float32.(vcat([logpdf.(px, x) for x in xs]...))
#
#         -Float32(Base.log(Nsamples)) .+ logsumexp(νlogpdfs .+ logpfails .- logpxs, dims=1)
# end
#
# function estimate_log_pfail(π, px, s, Nsamples=10)
#         xs = [reshape(rand(px, size(s,2)), :, size(s,2)) for _ in 1:Nsamples]
#         logpfails = vcat([value(π, s, x) for x in xs]...)
#
#         -Float32(Base.log(Nsamples)) .+ logsumexp(logpfails, dims=1)
# end
#
# function estimate_pfail(π, px, s, Nsamples=10)
#         xs = [reshape(rand(px, size(s,2)), :, size(s,2)) for _ in 1:Nsamples]
#         pfails = vcat([value(π, s, x) for x in xs]...)
#
#         mean(pfails, dims=1)
# end
#
# function IS_Continuous_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
#         term = 𝒟[:done] .| 𝒟[:fail]
#         # term .* (𝒟[:fail] .== false) .* -100f0 .+ (1.f0 .- term) .* estimate_log_pfail(π, 𝒫[:px], 𝒟[:sp], 𝒫[:N_IS_Samples])
#         norm = estimate_pfail(π, 𝒫[:px], 𝒟[:s], 𝒫[:N_IS_Samples])
#         term .* 𝒟[:fail] .+ (1.f0 .- term) .* (estimate_pfail(π, 𝒫[:px], 𝒟[:sp], 𝒫[:N_IS_Samples])) ./ norm
# end
#
# function IS_L_KL_log(π, 𝒫, 𝒟; kwargs...)
#         x, logνx = exploration(π, 𝒟[:s])
#         logpfail_s = Zygote.ignore() do
#                 estimate_log_pfail(π, 𝒫[:px], 𝒟[:s], 𝒫[:N_IS_Samples])
#         end
#         logpx = logpdf(𝒫[:px], x)
#         logpfail_sx = value(π, 𝒟[:s], x)
#
#
#         -mean(exp.( logpx .+ logpfail_sx .- logνx .- logpfail_s) .* logνx)
# end
#
# function IS_L_KL(π, 𝒫, 𝒟; kwargs...)
#         x, logνx = exploration(π, 𝒟[:s])
#         νx = exp.(logνx)
#         pfail_s = Zygote.ignore() do
#                 estimate_pfail(π, 𝒫[:px], 𝒟[:s], 𝒫[:N_IS_Samples])
#         end
#         px = exp.(logpdf(𝒫[:px], x))
#         pfail_sx = value(π, 𝒟[:s], x)
#
#
#         -mean(logνx .* px .* pfail_sx ./ (νx .* pfail_s))
# end
#
#
# function compute_IS_weight(𝒟, 𝒫; info=Dict())
#         𝒟[:weight] .= exp.(sum(𝒫[:xlogprobs] .* 𝒟[:x], dims=1) .- 𝒟[:xlogprob])
#         info[:mean_is_weight] = mean(𝒟[:weight])
#         info[:std_is_weight] = std(𝒟[:weight])
#         info[:min_is_weight] = minimum(𝒟[:weight])
#         info[:max_is_weight] = maximum(𝒟[:weight])
# end
#
# function compute_IS_weight_continuous(𝒟, 𝒫; info=Dict())
#         if 𝒫[:px] isa UnivariateDistribution
#                 𝒟[:weight] .= clamp.(exp.(reshape(logpdf.(𝒫[:px], 𝒟[:x]), 1, :) .- 𝒟[:xlogprob]), 0f0, 5f0)
#         else
#                 𝒟[:weight] .= clamp.(exp.(reshape(logpdf(𝒫[:px], 𝒟[:x]), 1, :) .- 𝒟[:xlogprob]), 0f0, 5f0)
#         end
#         info[:mean_is_weight] = mean(𝒟[:weight])
#         info[:std_is_weight] = std(𝒟[:weight])
#         info[:min_is_weight] = minimum(𝒟[:weight])
#         info[:max_is_weight] = maximum(𝒟[:weight])
# end
#
# ISARL_DQN(;kwargs...) = DQN(;c_loss=td_loss(name=:x_Qavg, a_key=:x), prefix="x_", target_fn=IS_DQN_target)
# ISARL_DDPG(;kwargs...) = DQN(;a_loss=IS_L_KL, c_loss=td_loss(name=:x_Qavg, a_key=:x), prefix="x_", target_fn=IS_Continuous_target, )
#
# function CERL_DQN(;π::DiscreteNetwork,
#               N::Int,
#               ΔN=4,
#               π_explore=ϵGreedyPolicy(LinearDecaySchedule(1., 0.1, floor(Int, N/2)), π.outputs),
#               c_opt::NamedTuple=(;),
#               log::NamedTuple=(;),
#               c_loss=td_loss(),
#               target_fn=DQN_target,
#               prefix="",
#               kwargs...)
#
#
#                     OffPolicySolver(;agent=PolicyParams(π=π, π_explore=π_explore, π⁻=deepcopy(π)),
#                                      log=LoggerParams(;dir="log/cerl_dqn", log...),
#                                      N=N,
#                                      ΔN=ΔN,
#                                      c_opt = TrainingParams(;loss=c_loss, name=string(prefix, "critic_"), epochs=ΔN, c_opt...),
#                                      target_fn=target_fn,
#                                      kwargs...)
# end