text:
using DiffEqJumpExtensions
using Base.Test

# write your own tests here
tic()
#@time @testset "Functionwrapper tests" begin include("functwrapper_tests.jl") end
toc()

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
"""
# DenseSPNs.jl

Implements dense SPNs for representing images

"""

"""
## Region Graphs

From "Dennis, A. and Ventura, D. Learning the architecture of sum-product networks using clustering on variables.
In Advances in Neural Information Processing Systems 25, 2012":

A *region graph* is a rooted DAG consisting of region nodes and partition nodes.
The root node is a region node.
Partition nodes are restricted to being the children of region nodes and vice versa.
Region and partition nodes have scopes just like nodes in an SPN.
"""

import Base: hash, isequal

export
    RegionGraphNode,
    RegionNode,
    PartitionNode,
    regionGraph,
    buildDenseSPN

abstract type RegionGraphNode end

# Represents rectangular region
# Assume x1,x2,y1,y2 are positive
struct RegionNode  <: RegionGraphNode
    children::Vector{RegionGraphNode}
    x1::Int64
    x2::Int64
    y1::Int64
    y2::Int64
    RegionNode(x1,x2,y1,y2) = new(Vector{PartitionNode}(),x1,x2,y1,y2)
end

# Represents decomposition of regions
struct PartitionNode <: RegionGraphNode
    children::Vector{RegionGraphNode}
    x1::Int64
    x2::Int64
    y1::Int64
    y2::Int64
end

Coordinates = Tuple{Int64,Int64,Int64,Int64}

# The sum and product of two numbers a, b satisfy the equation
# x^2 - (a+b)x + (ab) = 0
# Since the quadratic equation admits only two roots, a and b must be unique
# Hence, if we assume that 0 < a < b, (a+b) + ab == (x+y) + xy iff a=b and x=y
# We can enforce that by using a and b'=a+b (as a > 0)
#Base.hash(r::RegionGraphNode) = 2*(r.x1 + r.x2 + r.x1*r.x2) + r.y1 + r.y2 + r.y1*r.y2 + (r.x1 + r.x2 + r.x1*r.x2)*(r.x1 + r.x2 + r.x1*r.x2+r.y1 + r.y2 + r.y1*r.y2)
#Base.isequal(r1::RegionGraphNode, r2::RegionGraphNode) = (r1.x1 == r2.x1 && r1.x2 == r2.x2 && r1.y1 == r2.y1 && r1.y2 == r2.y2)

# Build region graph on [x1,x2]x[y1,y2], with numparts decompositions per split
function regionGraph(x1::Int64,x2::Int64,y1::Int64,y2::Int64,numparts)
    cache = Dict{Coordinates,RegionNode}()
    r = getOrCreateRegionNode(x1,x2,y1,y2,numparts,cache)
    #println(length(cache), " in cache")
    return r
end

function getOrCreateRegionNode(x1::Int64,x2::Int64,y1::Int64,y2::Int64,numparts::Int64,cache::Dict{Coordinates,RegionNode})
    # root node is region [x1,x2]x[y1,y2] with 3 decompositions
    #spaces = repeat(' ',depth)
    #println("$(hash(r)): $spaces [$(x1),$(x2)]x[$(y1),$(y2)]")
    # if region is pixel, create new node and return it
    if x1 == x2 && y1 == y2
        return RegionNode(x1,x2,y1,y2)
    end
    # if region has already been created
    # return it (and do not decompose it nor duplicated it)
    if haskey(cache,(x1,x2,y1,y2))
        return cache[x1,x2,y1,y2]
    end
    r = RegionNode(x1,x2,y1,y2)
    # create binary decompositions and recur
    # try splitting horizontally if possible
    if y1 < y2
        for i = 1:numparts
            y = rand(y1:y2) # split at y
            if y == y2
                r1 = getOrCreateRegionNode(x1,x2,y1,y-1,numparts,cache)
                r2 = getOrCreateRegionNode(x1,x2,y2,y2,numparts,cache)
            else
                r1 = getOrCreateRegionNode(x1,x2,y1,y,numparts,cache)
                r2 = getOrCreateRegionNode(x1,x2,y+1,y2,numparts,cache)
            end
            part = PartitionNode([r1,r2],x1,x2,y1,y2)
            push!(r.children,part)
        end
    # otherwise split vertically
    elseif x1 < x2
        for i = 1:numparts
            x = rand(x1:x2) # split at
            if x == x2
                r1 = getOrCreateRegionNode(x1,x-1,y1,y2,numparts,cache)
                r2 = getOrCreateRegionNode(x2,x2,y1,y2,numparts,cache)
            else
                r1 = getOrCreateRegionNode(x1,x,y1,y2,numparts,cache)
                r2 = getOrCreateRegionNode(x+1,x2,y1,y2,numparts,cache)
            end
            part = PartitionNode([r1,r2],x1,x2,y1,y2)
            push!(r.children,part)
        end
    end
    cache[x1,x2,y1,y2] = r
    return r
end

function Base.show(io::IO, n::RegionNode)
  show(io, n, 0)
end

function Base.show(io::IO, n::RegionNode, h::Int)
  println(io, " "^h, "(+ [$(n.x1),$(n.x2)]x[$(n.y1),$(n.y2)] $(objectid(n))")
  for ch in n.children
    show(io, ch, h+1)
  end
  println(io, " "^h, ")")
end

function Base.show(io::IO, n::PartitionNode, h::Int)
  println(io, " "^h, "(* ")
  for ch in n.children
    show(io, ch, h+1)
  end
  println(io, " "^h, ")")
end

"""
    Build Discrete Dense SPN from Region Graph.

Takes a region [x1,x2]x[y1,y2], a number numcat of values per pixel, and
a number numpart of decompositions per split.
Use numcat=0 to generate Gaussian leaves.
"""
function buildDenseSPN(x1::Int64,x2::Int64,y1::Int64,y2::Int64,numcat::Int64,numparts=2)
    r = regionGraph(x1::Int64,x2::Int64,y1::Int64,y2::Int64,numparts)
    # obtain topological order of nodes using Kahn's algorithm
    nodes::Vector{Node} = []
    backward::Vector{Vector{Int64}} = Vector() # indexed by [src, dst]
    n2i = Dict{RegionGraphNode,Int64}() # map region graph node into spn node id
    i2n = Dict{Int64,RegionGraphNode}() # inverse map
    l2n = Dict{Int64,Int64}() # map variable indicator to spn node id
    visited = Set{RegionGraphNode}()
    # first compute node in-degree
    indegree = Dict{RegionGraphNode,Int64}(r => 0)
    queue = RegionGraphNode[r]
    while length(queue) > 0
        n = pop!(queue)
        push!(visited,n)
        for ch in n.children
            if haskey(indegree, ch)
                indegree[ch] += 1
            else
                indegree[ch] = 1
            end
        end
        append!(queue, ch for ch in n.children if !in(ch, visited) && !in(ch, queue))
    end
    #@info "#region nodes $(length(visited))"
    tovisit = RegionGraphNode[r]
    # now obtain indices in topo order
    i = 0
    while length(tovisit) > 0
        i += 1
        n = pop!(tovisit)
        if isa(n, RegionNode)
            push!(nodes,SumNode())
        else
            push!(nodes,ProductNode())
        end
        push!(backward, Int64[])
        n2i[n] = i
        i2n[i] = n
        for ch in n.children
            indegree[ch] -= 1
            #@assert indegree[ch] >= 0
            if indegree[ch] == 0
                push!(tovisit, ch)
            end
        end
    end
    @assert length(nodes) == length(keys(indegree)) "Nodes and Indegree size mismatch: $(length(nodes)) $(length(keys(indegree)))"
    @assert length(keys(n2i)) == size(backward,1) "n2i and backward size mismatch"
    # create leaves
    for i = 1:((r.x2-r.x1+1)*(r.y2-r.y1+1))
        l2n[i] = length(nodes)+1
        if numcat > 0
            for k=1:numcat
                values = zeros(Int64, numcat)
                values[k] = 1
                push!(nodes,CategoricalDistribution(i, values))
            end
        else
            push!(nodes,GaussianDistribution(i,0.0,1.0))
        end
    end
    # now obtain DAG
    from = Int64[]
    to = Int64[]
    for i = 1:length(backward)
        n = i2n[i]
        if length(n.children) > 0
            append!( backward[i], n2i[ch] for ch in n.children )
        else
            var = n.x1 + (r.x2-r.x1+1)*(n.y1-1)
            if numcat > 0
                append!( backward[i], [l2n[var]+k-1 for k=1:numcat])
            else
                append!( backward[i], [l2n[var]])
            end
        end
        if isa(nodes[i],SumNode)
            for ch in backward[i]
                push!(to, ch)
                push!(from, i)
            end
        end
    end
    weights = sparse(to, from, ones(Float64, length(to)))
    # start with arbitrary weights
    sumnodes = filter(i -> isa(nodes[i], SumNode), 1:length(nodes))
    for i in sumnodes
        ch = backward[i]
        w = Random.rand(length(ch))
        weights[ch,i] = w/sum(w)
    end
    return SumProductNetwork{Int64}(nodes,backward,weights)
end

ext: jl


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

vstack(plot(sin, 0, 20, Coord.cartesian(xmin=2π, xmax=3π)),
       plot(sin, 0, 20, Coord.cartesian(ymin=-0.5, ymax=0.5)))


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Random
Random.seed!(145)
x, y = 1:2:20, 5*rand(10)
xerr = abs.(randn(10))
lcval = 1:10
@gp x y xerr lcval "w xerr t 'test' lc var" xrange= (0,22) yrange = (0,6)
save(term="pngcairo size 600,400", output="./code/plt_errs18.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# using Quads
using Test
using Random

import Quads.PokerCard
import Quads.PokerDeck
import Quads.Lookup
import Quads.Evaluator

@testset "Quads.jl" begin
    # test the decks
    @test length(PokerDeck.make_deck()) == 52
    @test length(PokerDeck.draw!(5, PokerDeck.make_deck())) == 5

    # some hand and eval tests
    hand_string = ["2h", "3h", "4h", "5h", "Ah"]
    @test length(PokerCard.newCard.(hand_string)) == 5
    hand_ints = PokerCard.newCard.(hand_string)
    @test Evaluator.class_to_string(Evaluator.get_rank_class(Evaluator.evaluate(hand_ints))) == "Straight Flush"

    # seven cards from the deck and eval
    hand_ints = PokerDeck.draw!(7, PokerDeck.make_deck())
    @test typeof(PokerCard.int_to_pretty_str.(hand_ints)) == Array{String,1}
    @test typeof(Evaluator.get_rank_class(Evaluator.evaluate(hand_ints))) == Int

    # some timing tests for 1000 random draws and evals
    t = @elapsed for i in 1:1000
        hand_ints = PokerDeck.draw!(7, PokerDeck.make_deck())
        Evaluator.get_rank_class(Evaluator.evaluate(hand_ints))
    end # elapsed time

    # 1k draws and evals in under 1s
    @assert t < 1

end # tests

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
function OneStepAdjoint!(spt2::SnapShot, spt1::SnapShot, fidMtx::FidMtx)
    spt2.it = spt1.it - 1
    spt2.vz = spt1.vz + (fidMtx.MpzBvz)' * spt1.pz
    spt2.vx = spt1.vx + (fidMtx.MpxBvx)' * spt1.px
    spt2.pz =           (fidMtx.MpzBpz)' * spt1.pz
    spt2.px =           (fidMtx.MpxBpx)' * spt1.px

    spt2.pz = (fidMtx.MvzBp )' * spt2.vz + (fidMtx.MvxBp)' * spt2.vx + spt2.pz
    spt2.px = (fidMtx.MvzBp )' * spt2.vz + (fidMtx.MvxBp)' * spt2.vx + spt2.px
    spt2.vz = (fidMtx.MvzBvz)' * spt2.vz
    spt2.vx = (fidMtx.MvxBvx)' * spt2.vx
    return nothing
end

function MultiStepAdjoint(path::ASCIIString, shot::Shot, fidMtx::FidMtx; interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx   ;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    ot = shot.ot   ; nt    = shot.nt     ; dt = shot.dt;
    if iflag == 1
       zupper = ext
       Nz = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    AddShot2SnapShot!(spt1, shot, Nz, Nx, ext, zupper)
    path_tmp = join([path "_tmp"])
    fid  = writeSnapShot(path_tmp, spt1)
    for it = nt-1: -1: 1
        OneStepAdjoint!(spt2, spt1, fidMtx)
        AddShot2SnapShot!(spt2, shot, Nz, Nx, ext, zupper)
        copy(spt1, spt2)
        writeSnapShot(fid, spt1)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    close(fid)
    reverseSnapShotsOrder(path, path_tmp)
    return nothing
end

function MultiStepAdjoint(ntw::Int64, path::ASCIIString, shot::Shot, fidMtx::FidMtx)
    nz = fidMtx.nz ; nx    = fidMtx.nx   ;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    ot = shot.ot   ; nt    = shot.nt     ; dt = shot.dt;
    if iflag == 1
       zupper = ext
       Nz = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    AddShot2SnapShot!(spt1, shot, Nz, Nx, ext, zupper)
    path_tmp = join([path "_tmp"])
    if ntw > nt
       error("length of src is larger than records")
    elseif ntw == nt
       fid = saveSpt2SrcSpt(path_tmp, spt1, nz, nx, ntw, ext, Nz, Nx, zupper)
       for it = ntw-1 : -1 : 1
           OneStepAdjoint!(spt2, spt1, fidMtx)
           AddShot2SnapShot!(spt2, shot, Nz, Nx, ext, zupper)
           copy(spt1, spt2)
           saveSpt2SrcSpt(fid, spt1, nz, nx, ext, Nz, Nx, zupper)
       end
    elseif ntw < nt
       for it = nt-1: -1: ntw
           OneStepAdjoint!(spt2, spt1, fidMtx)
           AddShot2SnapShot!(spt2, shot, Nz, Nx, ext, zupper)
           copy(spt1, spt2)
       end
       fid = saveSpt2SrcSpt(path_tmp, spt1, nz, nx, ntw, ext, Nz, Nx, zupper)
       for it = ntw-1 : -1 : 1
           OneStepAdjoint!(spt2, spt1, fidMtx)
           AddShot2SnapShot!(spt2, shot, Nz, Nx, ext, zupper)
           copy(spt1, spt2)
           saveSpt2SrcSpt(fid, spt1, nz, nx, ext, Nz, Nx, zupper)
       end
    end
    close(fid)
    reverseSrcSptOrder(path, path_tmp)
    return nothing
end

function MultiStepAdjoint(w::Array{Float64,1}, shot::Shot, fidMtx::FidMtx; interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx   ;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    ot = shot.ot   ; dt    = shot.dt     ;
    nt = shot.nt   ; ntsrc = length(w)   ;
    if iflag == 1
       zupper = ext
       Nz = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    dis= zeros(nz*nx)
    spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    AddShot2SnapShot!(spt1, shot, Nz, Nx, ext, zupper)
    for it = nt-1: -1: 1
        OneStepAdjoint!(spt2, spt1, fidMtx)
        AddShot2SnapShot!(spt2, shot, Nz, Nx, ext, zupper)
        copy(spt1, spt2)
        if it <= ntsrc
           dis = dis + spt2dis(w, spt1)
        end
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return dis
end

function MultiStepAdjoint(ntw::Int64, dis::Array{Float64,1}, shot::Shot, fidMtx::FidMtx; interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx   ;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    ot = shot.ot   ; dt    = shot.dt     ;
    nt = shot.nt   ;
    if iflag == 1
       zupper = ext
       Nz = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    w  = zeros(ntw)
    spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, nt)
    AddShot2SnapShot!(spt1, shot, Nz, Nx, ext, zupper)
    for it = nt-1: -1: 1
        OneStepAdjoint!(spt2, spt1, fidMtx)
        AddShot2SnapShot!(spt2, shot, Nz, Nx, ext, zupper)
        copy(spt1, spt2)
        if it <= ntw
           spt2wlet(w, dis, spt1)
        end
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return w
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
include("MapTransform/MapTransform.jl")

using JSON
using .MapTransform
using NCDatasets
using Formatting
using ArgParse
using LaTeXStrings

missing2nan  = (x,) -> nomissing(x, NaN)
missing2zero = (x,) -> nomissing(x, 0.0)

ρcp = 1026.0 * 3996.0

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--input-file"
            help = "Ocean history file."
            arg_type = String
            required = true

        "--output-img"
            help = "Output image name."
            arg_type = String
            default = ""

        "--no-display"
            help = ""
            action = :store_true

    end

    return parse_args(s)
end

parsed = parse_commandline()

JSON.print(parsed, 4)

Dataset(parsed["input-file"], "r") do ds

    global Nx = ds.dim["Nx"]
    global Ny = ds.dim["Ny"]
    global Nz = ds.dim["Nz"]
    global Nt = ds.dim["time"]

    global dz_cT   = reshape(nomissing(ds["dz_cT"][:], 0), 1, 1, Nz)
    global dσ_sT   = reshape(nomissing(ds["area_sT"][:], 0), Nx, Ny)
    global ϕ_sT    = reshape(nomissing(ds["lat_sT"][:], 0), Nx, Ny)
    global mask_sT = reshape(nomissing(ds["mask_sT"][:], 0), Nx, Ny)

    global ADVT   = nomissing(ds["ADVT"][:], 0.0)
    global QFLXT = nomissing(ds["QFLXT"][:], 0.0)

end

lat_v = range(-90.0, 90.0, length=181) |> collect
lat_t = (lat_v[1:end-1] + lat_v[2:end]) / 2.0
r = MapTransform.Relation(;
    lat     = ϕ_sT,
    area    = dσ_sT,
    mask    = mask_sT,
    lat_bnd = lat_v,
)

using PyCall
import PyPlot as plt

if parsed["no-display"]
    pygui(false)
else
    pygui(true)
end

fig, ax = plt.subplots(2, 1, sharex=true, constrained_layout=true)

fig.suptitle(parsed["input-file"])

for t=1:Nt

    _ADVT   = view(ADVT, :, :, :, t)
    HFC_ADV = sum(_ADVT .* dz_cT, dims=3)[:, :, 1]  * ρcp # HFC = Heat Flux Convergence

    HFC_ADV_t =   MapTransform.transform(r, HFC_ADV)
    OHT_ADV_v = - MapTransform.∫∂a(r, HFC_ADV)

    _QFLXT = view(QFLXT, :, :, :, t)
    HFC_QFLX = sum(_QFLXT .* dz_cT, dims=3)[:, :, 1]  * ρcp # HFC = Heat Flux Convergence

    HFC_QFLX_t =   MapTransform.transform(r, HFC_QFLX)
    OHT_QFLX_v = - MapTransform.∫∂a(r, HFC_QFLX)

#    ax[1].scatter(lat_t, HFC_ADV_t)

    HFC_ADV_t[isnan.(HFC_ADV_t)] .= 0

    ax[1].plot(lat_t, HFC_ADV_t,   "r--", label="ADV")
    ax[1].plot(lat_t, HFC_QFLX_t, "b--", label="QFLX")
    ax[1].plot(lat_t, HFC_ADV_t + HFC_QFLX_t, "k-", label="SUM")

    ax[2].plot(lat_v, OHT_ADV_v / 1e15, "r--")
    ax[2].plot(lat_v, OHT_QFLX_v / 1e15, "b--")
    ax[2].plot(lat_v, (OHT_ADV_v + OHT_QFLX_v) / 1e15, "k-")

end
ax[1].legend()
ax[1].grid(true)
ax[2].grid(true)
ax[2].set_ylim([-6, 8])

ax[1].set_ylabel(L"[$ \mathrm{W} / \mathrm{m}^2 $]")
ax[2].set_ylabel(L"[$ \mathrm{PW} $]")

if parsed["output-img"] != ""
    println("Output image: ", parsed["output-img"])
    fig.savefig(parsed["output-img"], dpi=150)
end

if ! parsed["no-display"]
    println("Showing image...")
    plt.show()
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module testSep_Model
using Test
using JuMP

srcDir = pwd()*"\\src\\"
include(srcDir*"Examples.jl")
include(srcDir*"sep_Model.jl")
using .Examples
using .sep_Model

@testset "Model test" begin
    ρ_dict = Examples.get_examples()
    ρ_sep  = ρ_dict["sep"]
    ρ = ρ_sep["sep3d3r3"]
    t = 3
    con_list = "S₁S₂S₃sG"

    sep_mod = Modelξₜˢᵉᵖ(ρ,t,con_list)

    @test typeof(sep_mod) == JuMP.Model
    @test size(sep_mod[:Lx])[1] == 924
# end

# @testset "Model export test" begin
    export_model(sep_mod,pwd()*"\\Defualt_sep_mod.dat-s")
#end

#@testset "Model export test" begin
    sep_mod_read = read_model(pwd()*"\\Defualt_sep_mod.dat-s")
    # @test sep_mod_read  == sep_mod
end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
function inverse(side::HorizonSide)
    return HorizonSide(mod(Int(side)+2, 4))
end

function main(r::Robot)
    for side in sides
        putmarkers!(r,side)
        followMarkers(r,inverse(side))
    end
    putmarker!(r)
end

function putmarkers!(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

function followMarkers(r::Robot,side::HorizonSide)
    while ismarker(r)==true
        move!(r,side)
    end
end

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
# ---
# title: 424. Longest Repeating Character Replacement
# id: problem424
# author: Indigo
# date: 2020-02-02
# difficulty: Medium
# categories: Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/longest-repeating-character-replacement/description/>
# hidden: true
# ---
#
# Given a string `s` that consists of only uppercase English letters, you can
# perform at most `k` operations on that string.
#
# In one operation, you can choose **any** character of the string and change it
# to any other uppercase English character.
#
# Find the length of the longest sub-string containing all repeating letters you
# can get after performing the above operations.
#
# **Note:**
# Both the string's length and _k_ will not exceed 104.
#
# **Example 1:**
#
#
#
#     Input:
#     s = "ABAB", k = 2
#
#     Output:
#     4
#
#     Explanation:
#     Replace the two 'A's with two 'B's or vice versa.
#
#
#
#
# **Example 2:**
#
#
#
#     Input:
#     s = "AABABBA", k = 1
#
#     Output:
#     4
#
#     Explanation:
#     Replace the one 'A' in the middle with 'B' and form "AABBBBA".
#     The substring "BBBB" has the longest repeating letters, which is 4.
#
#
#
#
#
## @lc code=start
using LeetCode

function character_replacement(s::String, k::Int)
    i, j, len, maxn = 1, 1, length(s), 0
    cntr = fill(0, 26)
    while j ≤ len
        maxn = max(maxn, cntr[s[j] - 'A' + 1] += 1)
        if maxn + k < j - i + 1
            cntr[s[i] - 'A' + 1] -= 1
            i += 1
        end
        j += 1
    end
    j - i
end
## @lc code=end


text:
module SparseVectors

using Compat

import Base:
    length, size, nnz, countnz, nonzeros, getindex, setindex!,
    showarray, show, writemime,
    convert, full, vec, copy, reinterpret, float, complex,
    scale, scale!, dot, vecnorm, hcat, vcat,
    sum, sumabs, sumabs2, maximum, minimum, maxabs, minabs,
    sparse, sprand, sprandn,
    A_mul_B!, At_mul_B, At_mul_B!, Ac_mul_B, Ac_mul_B!

import Base:
    +, .+, -, .-, *, .*, max, min,
    abs, abs2, conj, real, imag,
    floor, ceil, trunc, round,
    exp, exp2, exp10, log, log2, log10,
    log1p, expm1,
    sin, cos, tan, csc, cot, sec, sinpi, cospi,
    sind, cosd, tand, cscd, cotd, secd,
    asin, acos, atan, acot,
    asind, acosd, atand, acotd,
    sinh, cosh, tanh, csch, coth, sech,
    asinh, atanh, acsch, asech

import Base.LinAlg

import ArrayViews: view
import Base.LinAlg: axpy!

export
    # reexport
    view, axpy!,

    # types
    SparseVector,
    SparseVectorView,

    # functions
    nonzeroinds,
    sparsevector,
    getcol,
    unsafe_colrange,
    densemv

# sources
include("common.jl")

include("sparsevec.jl")
include("sparsevecview.jl")
include("sparsematview.jl")

include("generics.jl")
include("math.jl")
include("linalg.jl")

end # module


text:
# scale factors are standard deviations
# square them for (co)variance
struct PercentileMatch{T} <: AbstractVector{T}
    scale_factors::Vector{T}
end
Base.size(pm::PercentileMatch) = size(pm.scale_factors)
Base.length(pm::PercentileMatch) = length(pm.scale_factors)
Base.getindex(pm::PercentileMatch, i) = pm.scale_factors[i]
Base.setindex!(pm::PercentileMatch, v, i) = pm.scale_factors[i] = v
Base.resize!(pm::PercentileMatch, N) = resize!(pm.scale_factors, N)

PercentileMatch{T}(undef, N) where T = PercentileMatch(Vector{T}(undef, N))
function PercentileMatch(s::MahalanobisDistances)
    match_percentiles!(PercentileMatch(similar(s.mahal)), s)
end
function PercentileMatch!(pm::PercentileMatch, s::MahalanobisDistances)
    resize!(pm, length(s))
    match_percentiles!(pm, s)
end

function match_percentiles!(pm::PercentileMatch{T}, s::MahalanobisDistances{T}) where T
    sort!(s)
    N = length(s)
    denom = 1 / (N + 1)
    @inbounds for i ∈ eachindex(s)
        pm[i] = s[i] / quantile(Chi(3), i * denom)
    end
    pm
end

function sample_distances!(rng::AbstractRNG, d::UniformSamples, res::PercentileMatch)
    @fastmath @inbounds for n ∈ eachindex(d)
        d[n] = sqrt( (abs2(randn(rng)) + abs2(randn(rng)) + abs2(randn(rng))) ) * rand(rng, res)
    end
    d
end
sample_distances!(d::UniformSamples, res::PercentileMatch) = sample_distances!(Random.GLOBAL_RNG, d, res)

function sample_Pc!(rng::AbstractRNG, pc_array::AbstractArray, res::AbstractVector, conj, n::Int)
    c1 = SymmetricM3(conj[n,73],conj[n,74],conj[n,79],conj[n,75],conj[n,80],conj[n,84]) / 1e6
    c2 = SymmetricM3(conj[n,133],conj[n,124],conj[n,139],conj[n,135],conj[n,140],conj[n,144]) / 1e6

    r1 = SVector(ntuple(i -> conj[n,i+171], Val(3)))
    v1 = SVector(ntuple(i -> conj[n,i+174], Val(3)))
    r2 = SVector(ntuple(i -> conj[n,i+177], Val(3)))
    v2 = SVector(ntuple(i -> conj[n,i+180], Val(3)))

    HBR = conj[n,157] / 1e3
    sample_Pc!(rng, pc_array, res, r1, v1, c1, r2, v2, c2, HBR)
end

function sample_Pc!(rng::AbstractRNG, pc_array::AbstractArray, res::AbstractVector, r1, v1, c1, r2, v2, c2, HBR)
    resize!(pc_array, length(res))
    @inbounds for i ∈ eachindex(pc_array,res)
        pc_array[i] = pc2dfoster_RIC(r1, v1, c1, r2, v2, c2 * (res[i])^2, HBR)
    end
    pc_array
end

function Distributions.logpdf(res::PercentileMatch{T}, x::SVector{3}) where {T}
    N = length(res)
    m = T(-0.5)*(x' * x)
    ms = zero(T)
    @fastmath @simd for s ∈ res.scale_factors
        ms += SLEEFwrap.exp( m / s^2 - T(3)*SLEEFwrap.log(s) )
    end
    # log(ms) - T(log(N)) - T(1.5log(2π))
    log(ms) - log(N) - 1.5log(2π)
end


text:
if basename(pwd()) != "9"
    cd("2021/9")
end

function loadheightmap(filename::AbstractString)
    len = length(readline(filename))
    transpose(reshape([parse(Int, d) for line in eachline(filename) for d in line], len, :))
end

function part1(heights::AbstractMatrix{Int})
    sum(
        heights[i] + 1
        for i in CartesianIndices(heights)
        if all(>(heights[i]), neighbours(heights, i))
    )
end

function neighbourcoords(A, coord)
    (
        coord + offset
        for offset in (
            CartesianIndex(0, 1), CartesianIndex(0, -1),
            CartesianIndex(1, 0), CartesianIndex(-1, 0)
        )
        if coord + offset in CartesianIndices(A)
    )
end

neighbours(A, coord) = (A[c] for c in neighbourcoords(A, coord))

function part2(heights::AbstractMatrix{Int})
    prod(sort!(vec([floodnines(heights, i) for i in CartesianIndices(heights)]), rev=true)[1:3])
end

"""
Flood fills a basin with 9s and returns the number of locations flooded.
"""
function floodnines(heights::AbstractMatrix{Int}, i::CartesianIndex)
    if heights[i] != 9
        heights[i] = 9
        1 + sum(floodnines(heights, n) for n in neighbourcoords(heights, i))
    else
        0
    end
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase
using Test

# Calculus
@parameters t σ ρ β
@variables x y z

eqs = [σ*(y-x),
       x*(ρ-z)-y,
       x*y - β*z]

simpexpr = [
   :(σ * (y - x))
   :(x * (ρ - z) - y)
   :(x * y - β * z)
   ]

for i in 1:3
   @test ModelingToolkit.simplified_expr.(eqs)[i] == simpexpr[i]
   @test ModelingToolkit.simplified_expr.(eqs)[i] == simpexpr[i]
end

∂ = ModelingToolkit.jacobian(eqs,[x,y,z])
for i in 1:3
    ∇ = ModelingToolkit.gradient(eqs[i],[x,y,z])
    @test isequal(∂[i,:],∇)
end
@test all(isequal.(ModelingToolkit.gradient(eqs[1],[x,y,z]),[σ * -1,σ,0]))
@test all(isequal.(ModelingToolkit.hessian(eqs[1],[x,y,z]),0))

# Function building

@parameters σ() ρ() β()
@variables x y z
eqs = [σ*(y-x),
       x*(ρ-z)-y,
       x*y - β*z]
f1,f2 = ModelingToolkit.build_function(eqs,[x,y,z],[σ,ρ,β])
f = eval(f1)
out = [1.0,2,3]
o1 = f([1.0,2,3],[1.0,2,3])
f = eval(f2)
f(out,[1.0,2,3],[1.0,2,3])
@test all(o1 .== out)

function test_worldage()
   @parameters σ() ρ() β()
   @variables x y z
   eqs = [σ*(y-x),
          x*(ρ-z)-y,
          x*y - β*z]
   f, f_iip = ModelingToolkit.build_function(eqs,[x,y,z],[σ,ρ,β],(),ModelingToolkit.simplified_expr,Val{false})
   out = [1.0,2,3]
   o1 = f([1.0,2,3],[1.0,2,3])
   f_iip(out,[1.0,2,3],[1.0,2,3])
end
test_worldage()

## No parameters
@variables x y z
eqs = [(y-x)^2,
       x*(x-z)-y,
       x*y - y*z]
f1,f2 = ModelingToolkit.build_function(eqs,[x,y,z])
f = eval(f1)
out = zeros(3)
o1 = f([1.0,2,3])
f = eval(f2)
f(out,[1.0,2,3])
@test all(out .== o1)

function test_worldage()
   @variables x y z
   eqs = [(y-x)^2,
          x*(x-z)-y,
          x*y - y*z]
   f, f_iip = ModelingToolkit.build_function(eqs,[x,y,z],(),(),ModelingToolkit.simplified_expr,Val{false})
   out = zeros(3)
   o1 = f([1.0,2,3])
   f_iip(out,[1.0,2,3])
end
test_worldage()


text:
# Here we calculate the effective wavenumber and effective wave amplitudes, without any restriction to the volume fraction of particles and incident wave frequency.
using NLsolve

kernelN(n,x,y) = x*diffhankelh1(n,x)*besselj(n,y) - y*hankelh1(n,x)*diffbesselj(n,y)

# include depricated function to find a single effective wavenumber, when in fact there are many. The code is still used in tests and gives many correct results
include("wavenumber_single.jl")

" Returns all the transmitted effective wavenumbers"
wavenumbers(ω::T, medium::Medium{T}, specie::Specie{T}; kws...) where T<:Number = wavenumbers(ω, medium, [specie]; kws...)

function wavenumbers(ω::T, medium::Medium{T}, species::Vector{Specie{T}}; tol::T = 1e-6,
        hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
        mesh_points::Int = 7, mesh_size::T = 0.65,
        max_Imk::T = 0.0, max_Rek::T = 0.0,
        time_limit::T = 1.0,
        iterations::Int = 200,
        radius_multiplier::T = 1.005,
        t_vecs = t_vectors(ω, medium, species; hankel_order = hankel_order),
        kws...) where T<:Number

    k = ω/medium.c
    S = length(species)
    ho = hankel_order

    # Z_l_n = Zn_matrix(ω, medium, species; hankel_order = ho)
    # r = maximum(s.r for s in species)
    # φ = sum(volume_fraction.(species))

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    function M_component(keff,j,l,m,n)
        - (n==m ? 1.0:0.0)*(j==l ? 1.0:0.0) + 2.0pi*species[l].num_density*t_vecs[l][n+ho+1]*
            kernelN(n-m,k*as[j,l],keff*as[j,l])/(k^2.0-keff^2.0)
    end

    # this matrix is needed to calculate the eigenvectors
    MM(keff::Complex{T}) = reshape(
        [M_component(keff,j,l,m,n) for m in -ho:ho, j = 1:S, n in -ho:ho, l = 1:S]
    , ((2ho+1)*S, (2ho+1)*S))

    # the constraint uses keff_vec[2] < -tol to better specify solutions where imag(k_effs)<0
    constraint(keff_vec::Array{T}) = ( (keff_vec[2] < zero(T)) ? one(T) : zero(T) )*(-1 + exp(-T(100.0)*keff_vec[2]))
    # detMM2(keff_vec::Array{T}) =  constraint(keff_vec) + map(x -> real(x*conj(x)), det(MM(keff_vec[1]+im*keff_vec[2])))

    function detM!(F, x)
        cons = constraint(x)
        detM = det(MM(x[1]+im*x[2]))
        F[1] = real(detM) + cons
        F[2] = imag(detM) + cons
    end

    kφ = wavenumber_low_volfrac(ω, medium, species; verbose = false)
    eff_medium = effective_medium(medium, species)
    k0 = ω/eff_medium.c
    if isnan(k0) k0 = kφ end

    if max_Rek == 0.0
        dk_x = max(real(k0),abs(real(kφ))) * mesh_size
        max_Rek = mesh_points * dk_x
    else dk_x = max_Rek/mesh_points
    end
    if max_Imk == 0.0
        dk_y = abs(imag(kφ)) * mesh_size
        max_Imk = mesh_points * dk_y
    else dk_y = max_Imk/mesh_points
    end

    kx = -max_Rek:dk_x:max_Rek
    ky = 0.0:dk_y:max_Imk
    kins = [[x,y] for x in kx, y in ky][:]
    # add both the low volfrac and frequency as initial guesses
    push!(kins, [real(kφ),imag(kφ)], [real(k0),imag(k0)])

    # Find all wavenumbers
    # results = map(kins) do kin
    #    optimize(detMM2, kin; time_limit = time_limit)
    # end

    k_vecs = map(kins) do kin
        # println(kin)
        # try
            res = mcpsolve(detM!, [-T(10)*max_Rek, zero(T)], [T(10)*max_Rek, T(10)*max_Imk],
                       iterations = iterations, kin, xtol = sqrt(tol), factor=0.1)
            # res = nlsolve(detM!,kin; iterations = iterations, xtol = sqrt(tol))
            res.zero
        # catch
        #     [zero(T),-one(T)]
        # end
        # if res.residual_norm > sqrt(tol)
        #     [zero(T),-one(T)] # discard
        # else res.zero
        # end
    end

    # Take the best length(kx) candidates, plus those that are smaller than tolerance T(1e-4).
    # As tol does not affect the above results, it is best to use an imperical first pass tolerance T(1e-4).
    # sort!(results, by = r->r.minimum)
    # len = max(length(kx), findfirst(r -> sqrt(r.minimum) > T(1e-4), results))
    # k_vecs = [r.minimizer for r in results[1:len]]

    # remove unphysical wavenumbers
    # deleteat!(k_vecs, find(k_vec[2] < -tol for k_vec in k_vecs))

    function reduce_vecs(vecs::Vector{Vector{T}},tol::T)
        all_inds = collect(eachindex(vecs))
        vecs = map(vecs) do k_vec
            ind_ins = find(norm(ks - k_vec) < sqrt(tol) for ks in vecs[all_inds])
            inds = all_inds[ind_ins]
            deleteat!(all_inds,ind_ins)
            isempty(inds) ? [zero(T),-one(T)] :  mean(vecs[inds])
        end
        vecs = deleteat!(vecs, find(k_vec[2] < -tol for k_vec in vecs))
    end

    # group together wavenumbers which are closer than sqrt(tol)
    k_vecs = reduce_vecs(k_vecs,sqrt(tol))

    # Here we refine the effective wavenumbers
    # k_vecs = map(k_vecs) do k_vec    # Here we refine the effective wavenumbers
    #     res = optimize(detMM2, k_vec;  g_tol = tol^2.0, f_tol = tol^4.0)
    #     res.minimizer
    # end

    k_vecs = map(k_vecs) do kin
        res = mcpsolve(detM!, [-T(10)*max_Rek, zero(T)], [T(10)*max_Rek, T(10)*max_Imk],
                   iterations = 10*iterations, kin, xtol = tol, factor=0.1)
        # res = nlsolve(detM!,kin; iterations = 10*iterations, xtol = tol, method = :newton)
        if res.residual_norm > sqrt(tol)  # not sure about this cristeria
            [zero(T),-one(T)] # discard
        else res.zero
        end
    end

    # group together wavenumbers which are closer than tol
    k_vecs = reduce_vecs(k_vecs,tol)
    # nl_kvecs = reduce_vecs(nl_kvecs,tol)

    # Finally delete unphysical waves, including waves travelling backwards with almost no attenuation. This only is important in the limit of very low frequency or very weak scatterers.
    deleteat!(k_vecs, find(k_vec[2] < -tol for k_vec in k_vecs))
    deleteat!(k_vecs, find(k_vec[2] < tol && k_vec[1] < tol for k_vec in k_vecs))

    k_effs = [ k_vec[1] + k_vec[2]*im for k_vec in k_vecs]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end

# Alternative grouping
# all_inds = collect(eachindex(k_vecs))
# k_vecs = map(k_vecs) do k_vec
#     nk = norm(k_vec)*sqrt(tol)
#     ind_ins = find(norm(ks - k_vec) < nk for ks in k_vecs[all_inds])
#     inds = all_inds[ind_ins]
#     deleteat!(all_inds,ind_ins)
#     isempty(inds) ? [zero(T),-one(T)] :  mean(k_vecs[inds])
# end

# alternative grouping methods

    # digs = 3 - Int(round(log(10,sqrt(tol)))) # number of digits to round to check if equal
    # k_vecs = map(groupby(k_vec -> round.(k_vec,digs), k_vecs)) do g
    #     mean(g)
    # end

# detMM2(keff_vec::Array{T}) =  map(x -> real(x*conj(x)), det(MM(keff_vec[1]+im*keff_vec[2])))

# function detMM!(F,x)
#     F[1] = abs(det(MM(x[1]+im*x[2])))
# end

# Alternative solvers
# res = nlsolve(detMM!,initial_k_eff; iterations = 10000, factor=2.0)
# k_eff_nl = res.zero[1] + im*res.zero[2]
# lower = [0.,0.]; upper = [T(2)*k0,k0]
# result = optimize(detMM2, initial_k0, lower, upper; g_tol = tol^2.0, f_tol = tol^4.0)


text:
"""
    to_matrix(operator)

Converts to a dense matrix the KLocal Operator
"""
function to_matrix(op::AbsLinearOperator)
    hilb = basis(op)
    @assert indexable(hilb)

    N = spacedimension(hilb)
    mat = zeros(ComplexF64, N, N)

    for (i, σ) = enumerate(states(hilb))
        fun = (mel, cngs, σ) -> begin
            σp = apply(σ, cngs)
            j = index(hilb, σp)
            mat[i, j] += mel
        end

        map_connections(fun, op, σ)
    end
    return mat
end

Base.Matrix(op::AbsLinearOperator) = to_matrix(op)

"""
    to_sparse(operator)

Converts to a sparse matrix the KLocal Operator.
"""
function to_sparse(op::AbsLinearOperator)
    hilb = basis(op)
    @assert indexable(hilb)

    i_vals   = Vector{Int}()
    j_vals   = Vector{Int}()
    mel_vals = Vector{ComplexF64}()

    for (i, σ) = enumerate(states(hilb))
        conns = row_valdiff(op, σ)

        for (mel, cngs) = conns
            σp = apply(σ, cngs)
            j = index(hilb, σp)

            push!(i_vals, i)
            push!(j_vals, j)
            push!(mel_vals, mel)
        end
    end

    N = spacedimension(hilb)
    return sparse(i_vals, j_vals, mel_vals, N, N)
end

SparseArrays.sparse(op::AbsLinearOperator) = to_sparse(op)

function to_map(op::AbsLinearOperator)
    hilb = basis(op)
    N = spacedimension(hilb)

    σ = state(hilb)
    function op_v_prod!(y::AbstractVector, x::AbstractVector)
        y .= zero(eltype(y))

        for (i, x_val) = enumerate(x)
            x_val == 0 && continue

            state!(σ, hilb, i)
            conns = row_valdiff(op, σ)

            for (mel, cngs) = conns
                mel == 0 && continue

                σp = apply(σ, cngs)
                j = index(hilb, σp)
                println("$i, $j  => $mel -- $σ - $σp")
                y[j] += x_val*mel
            end
        end
        return y
    end

    return LinearMap{ComplexF64}(op_v_prod!, N, N; ismutating=true)
end

QuantumOpticsBase.DenseOperator(op::AbsLinearOperator) =
    DenseOperator(convert(CompositeBasis, basis(op)), to_matrix(op))

QuantumOpticsBase.SparseOperator(op::AbsLinearOperator) =
    SparseOperator(convert(CompositeBasis, basis(op)), to_sparse(op))

QuantumOpticsBase.DenseSuperOperator(op::KLocalLiouvillian) =
    DenseSuperOperator(nq_to_qo_basis(basis(op)), to_matrix(op))

QuantumOpticsBase.SparseSuperOperator(op::KLocalLiouvillian) =
    SparseSuperOperator(nq_to_qo_basis(basis(op)), to_sparse(op))


text:
@testset "W" begin
    @testset "wncard/wninsd" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 3
    end
    @testset "wncomd" begin
        window1 = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        for d in darray
            wninsd!(window1, d...)
        end
        @test wncard(window1) == 3
        window2 = wncomd(window1, 2.0, 20.0)
        @test wncard(window2) == 2
        @test wnfetd(window2, 1) == (3.0, 7.0)
        @test wnfetd(window2, 2) == (11.0, 20.0)
    end
    @testset "wncond!" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 3
        window = wncond!(window, 2.0, 1.0)
        @test wncard(window) == 2
        @test wnfetd(window, 1) == (9.0, 10.0)
        @test wnfetd(window, 2) == (25.0, 26.0)
    end
    @testset "wndifd" begin
        window1 = SpiceDoubleCell(8)
        window2 = SpiceDoubleCell(8)
        darray1 = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        darray2 = [[2.0, 6.0], [8.0, 10.0], [16.0, 18.0]]
        for d in darray1
            wninsd!(window1, d...)
        end
        @test wncard(window1) == 3
        for d in darray2
            wninsd!(window2, d...)
        end
        @test wncard(window2) == 3
        window3 = wndifd(window1, window2)
        @test wncard(window3) == 4
        @test wnfetd(window3, 1) == (1.0, 2.0)
        @test wnfetd(window3, 2) == (7.0, 8.0)
        @test wnfetd(window3, 3) == (10.0, 11.0)
        @test wnfetd(window3, 4) == (23.0, 27.0)
    end
    @testset "wnelmd" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 3
        array = [0.0, 1.0, 9.0, 13.0, 29.0]
        expected = [false, true, true, false, false]
        @test wnelmd.(Ref(window), array) == expected
    end
    @testset "wnexpd" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0], [29.0, 29.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 4
        wnexpd!(window, 2.0, 1.0)
        @test wncard(window) == 3
        @test wnfetd(window, 1) == (-1.0, 4.0)
        @test wnfetd(window, 2) == (5.0, 12.0)
        @test wnfetd(window, 3) == (21.0, 30.0)
    end
    @testset "wnextd" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0], [29.0, 29.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 4
        wnextd!(window, :L)
        @test wncard(window) == 4
        @test wnfetd(window, 1) == (1.0, 1.0)
        @test wnfetd(window, 2) == (7.0, 7.0)
        @test wnfetd(window, 3) == (23.0, 23.0)
        @test wnfetd(window, 4) == (29.0, 29.0)
        @test_throws SpiceError wnextd!(window, :test)
    end
    @testset "wnfild" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0], [29.0, 29.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 4
        wnfild!(window, 3.0)
        @test wncard(window) == 3
        @test wnfetd(window, 1) == (1.0, 3.0)
        @test wnfetd(window, 2) == (7.0, 11.0)
        @test wnfetd(window, 3) == (23.0, 29.0)
    end
    @testset "wnfltd" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0], [29.0, 29.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 4
        wnfltd!(window, 3.0)
        @test wncard(window) == 2
        @test wnfetd(window, 1) == (7.0, 11.0)
        @test wnfetd(window, 2) == (23.0, 27.0)
    end
    @testset "wnincd" begin
        window = SpiceDoubleCell(8)
        darray = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        for d in darray
            wninsd!(window, d...)
        end
        @test wncard(window) == 3
        array = [
            [1.0, 3.0],
            [9.0, 10.0],
            [0.0, 2.0],
            [13.0, 15.0],
            [29.0, 30.0],
        ]
        expected = [true, true, false, false, false]
        @testset for i in eachindex(array, expected)
            @test wnincd(window, array[i]...) == expected[i]
        end
    end
    @testset "wnintd" begin
        window1 = SpiceDoubleCell(8)
        window2 = SpiceDoubleCell(8)
        darray1 = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        darray2 = [[2.0, 6.0], [8.0, 10.0], [16.0, 18.0]]
        for d in darray1
            wninsd!(window1, d...)
        end
        @test wncard(window1) == 3
        for d in darray2
            wninsd!(window2, d...)
        end
        @test wncard(window2) == 3
        window3 = wnintd(window1, window2)
        @test wncard(window3) == 2
        @test wnfetd(window3, 1) == (2.0, 3.0)
        @test wnfetd(window3, 2) == (8.0, 10.0)
    end
    @testset "wnreld" begin
        window1 = SpiceDoubleCell(8)
        window2 = SpiceDoubleCell(8)
        darray1 = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        darray2 = [[1.0, 2.0], [9.0, 9.0], [24.0, 27.0]]
        for d in darray1
            wninsd!(window1, d...)
        end
        @test wncard(window1) == 3
        for d in darray2
            wninsd!(window2, d...)
        end
        @test wncard(window2) == 3
        ops = [:(==), :!=, :<=, :<, :>=, :>]
        expected = [false, true, false, false, true, true]
        for i in eachindex(ops, expected)
            @test wnreld(window1, ops[i], window2) == expected[i]
        end
        @test !(window1 == window2)
        @test window1 != window2
        @test !(window1 <= window2)
        @test !(window1 ⊆ window2)
        @test !(window1 < window2)
        @test !(window1 ⊊ window2)
        @test window1 >= window2
        @test window1 ⊇ window2
        @test window1 > window2
        @test window1 ⊋ window2
    end
    @testset "wnsumd" begin
        window = SpiceDoubleCell(12)
        darray = [
            [1.0, 3.0],
            [7.0, 11.0],
            [18.0, 18.0],
            [23.0, 27.0],
            [30.0, 69.0],
            [72.0, 80.0],
        ]
        for d in darray
            wninsd!(window, d...)
        end
        meas, avg, stddev, shortest, longest = wnsumd(window)
        @test meas == 57.0
        @test avg == 9.5
        @test stddev ≈ 13.413302
        @test shortest == 5
        @test longest == 9
    end
    @testset "wnunid" begin
        window1 = SpiceDoubleCell(8)
        window2 = SpiceDoubleCell(8)
        darray1 = [[1.0, 3.0], [7.0, 11.0], [23.0, 27.0]]
        darray2 = [[2.0, 6.0], [8.0, 10.0], [16.0, 18.0]]
        for d in darray1
            wninsd!(window1, d...)
        end
        @test wncard(window1) == 3
        for d in darray2
            wninsd!(window2, d...)
        end
        @test wncard(window2) == 3
        window3 = wnunid(window1, window2)
        @test wncard(window3) == 4
        @test wnfetd(window3, 1) == (1.0, 6.0)
        @test wnfetd(window3, 2) == (7.0, 11.0)
        @test wnfetd(window3, 3) == (16.0, 18.0)
        @test wnfetd(window3, 4) == (23.0, 27.0)
    end
    @testset "wnvald" begin
        window = SpiceDoubleCell(30)
        darray = [[0.0, 0.0], [10.0, 12.0], [2.0, 7.0],
                [13.0, 15.0], [1.0, 5.0], [23.0, 29.0],
                [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]
        for d in darray
            wninsd!(window, d...)
        end
        result = wnvald!(window)
        @test wncard(result) == 5
        @test wnfetd(result, 1) == (0.0, 0.0)
        @test wnfetd(result, 2) == (1.0, 7.0)
        @test wnfetd(result, 3) == (10.0, 12.0)
        @test wnfetd(result, 4) == (13.0, 15.0)
        @test wnfetd(result, 5) == (23.0, 29.0)
    end
end


text:
function pipeline(df::Table2Process)
    cols, kw = df.table, copy(df.kw)
    kw[:axis_type] = get(kw, :axis_type, :auto)
    kw[:compute_axis] = get(kw, :compute_axis, :auto)
    kw[:compute_error] = get(kw, :compute_error, :none)
    if isa(kw[:compute_error], Number)
        kw[:summarize] = get(kw, :summarize, (mean,std))
    else
        kw[:summarize] = get(kw, :summarize, (mean,sem))
    end
    kw[:acrossall] = get(kw, :acrossall, false)
    kw[:fkwargs] = get(kw, :fkwargs, [])
    kw[:xreduce] = get(kw, :xreduce, false)
    kw[:yreduce] = get(kw, :yreduce, false)
    kw[:plot_kwargs] = get(kw, :plot_kwargs, [])
    t = process_axis_type!(cols, kw)
    (kw[:axis_type] == :pointbypoint) || process_function!(t)
    return ProcessedTable(_group_apply(t), t.kw)
end

function process_axis_type!(cols, kw)
    n_grp = length(cols)- 3 - kw[:compare]
    pkey_vec = vcat([Symbol(:s, i) for i in 1:n_grp], kw[:compare] ? [:compare, :across] : [:across])
    col_names = vcat(pkey_vec, [:x, :y])
    pkey_tup = tuple(pkey_vec...)
    x, y = cols[end-1:end]
    kw[:acrossall] && (cols[end-2] .= collect(1.:Float64(length(cols[end-2]))))
    at = kw[:axis_type]
    (at == :pointbypoint) && return Table2Process(table(cols[1:end-2]..., x, y, names = col_names, pkey = pkey_tup), kw)

    if !(eltype(x)<:Real)
        (kw[:axis_type] in [:discrete, :auto]) || warn("Changing to discrete axis, x values are not real numbers!")
        kw[:axis_type] = :discrete
    end
    bin_width = 1.0
    if kw[:axis_type] == :binned
        nbins = get(kw, :nbins, 30)
        start, stop = extrema(x)
        edges = range(start, stop = stop, length = nbins+1)
        middles = (edges[2:end] .+ edges[1:end-1])./2
        indices = [searchsortedfirst(edges[2:end], s) for s in x]
        x_binned = middles[indices]
        bin_width = step(edges)
        kw[:axis_type] = :discrete
        cols = tuple(cols[1:end-2]..., x_binned, cols[end])
    end
    (kw[:axis_type] == :auto) && (kw[:axis_type] = :continuous)
    kw[:axis_type] in [:discrete, :continuous] ||
        error("Axis type $(kw[:axis_type]) is not supported")
    if all(t -> (t isa Real) && isnan(t), y)
        if kw[:axis_type] == :discrete
            y .= bin_width
        else
            y .= 0.0
        end
    end
    Table2Process(table(cols..., names = col_names, pkey = pkey_tup), kw)
end

function process_function!(t::Table2Process)
    if isa(t.kw[:f], Symbol)
        t.kw[:compute_axis] = get(builtin_axis, t.kw[:f], :auto)
        t.kw[:f] = builtin_funcs[t.kw[:f]]
    end
    t.kw[:fclosure] = (args...) -> t.kw[:f](Val{t.kw[:axis_type]}(), args...; t.kw[:fkwargs]...)
end

function get_grouped_error(trend, variation, f, xaxis, split_table, compute_error)
    if !isa(compute_error, Integer)
        subject_table = groupby(tt -> f(xaxis, table(tt)), split_table, flatten = true, select = IndexedTables.valuenames(split_table))
    else
        ns = compute_error
        large_table = table(repeat(collect(1:ns), inner = length(xaxis)),
            repeat(xaxis, outer = ns), zeros(ns*length(xaxis)), names = [:across, :x, :y],
            pkey = :across, presorted = true)
        subject_table = groupby(large_table, flatten = true, select = IndexedTables.valuenames(large_table)) do tt
            nd = length(split_table)
            perm = rand(1:nd, nd)
            permuted_data = table(columns(split_table, :x)[perm], columns(split_table, :x)[perm],
                names = [:x, :y])
            f(xaxis, permuted_data)
        end
    end
    nanfree = filter(_isfinite, subject_table)
    if compute_error == :none
        return select(nanfree, (:x, :y))
    else
        return groupby((:y => trend, :err => variation), nanfree, :x, select = :y)
    end
end

_isfinite(t) = true
_isfinite(t::Number) = isfinite(t)
_isfinite(t::Union{Tuple, NamedTuple}) = all(_isfinite, t)

function _group_apply(t::Table2Process)
    n = length(columns(t.table))-3
    if t.kw[:axis_type] == :pointbypoint
        if (t.kw[:xreduce] == false) || (t.kw[:yreduce] == false)
            t.kw[:compare] && error("can't compare without xreduce and yreduce")
            g = t.table
        elseif !t.kw[:compare]
            g = groupby((:x => :x => t.kw[:xreduce], :y => :y => t.kw[:yreduce]), t.table)
        else
            single_w = groupby((:xy => t.kw[:xreduce], ), t.table, select = :x)
            a, b = unique(column(single_w, :compare))
            xtable = renamecol(filter(i -> i.compare == a, single_w), :xy, :x)
            ytable = renamecol(filter(i -> i.compare == b, single_w), :xy, :y)
            pkey_tup = (listsplits(t)..., :across)
            g = join(xtable, ytable, lkey = pkey_tup, rkey = pkey_tup,
                lselect = :x, rselect = :y)
        end
    else
        g = groupby(t.table, (listsplits(t)...,), select = (:across, :x, :y), flatten = true) do tt
            split_table = table(tt, pkey = :across, presorted = true)
            xaxis = get_axis(columns(split_table, :x), t.kw[:axis_type], t.kw[:compute_axis])
            get_grouped_error(t.kw[:summarize]..., t.kw[:fclosure], xaxis, split_table, t.kw[:compute_error])
        end
    end
    return filter(_isfinite, g)
end


text:
using Decimals
using Test

@testset "Decimals" begin
    include("decimal.jl")
    include("arithmetic.jl")
end


text:
#=
Non-dominated Sorting Genetic Algorithm (NSGA-II) for Multi-objective Optimization

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `crossoverRate`: The fraction of the population at the next generation, that is created by the crossover function
- `mutationRate`: Probability of chromosome to be mutated
- `selection`: [Selection](@ref) function (default: `tournament`)
- `crossover`: [Crossover](@ref) function (default: `SBX`)
- `mutation`: [Mutation](@ref) function (default: `PLM`)
- `metrics` is a collection of convergence metrics.
=#
Base.@kwdef struct NSGAII <: AbstractOptimizer
    populationSize::Int = 50
    crossoverRate::Float64 = 0.9
    mutationRate::Float64 = 0.2
    selection = tournament(2, select = twowaycomp)
    crossover = SBX()
    mutation = PLM()
    metrics::ConvergenceMetrics = ConvergenceMetric[GD(), GD(true)]
end

population_size(method::NSGAII) = method.populationSize
default_options(method::NSGAII) = (iterations = 1000,)
summary(m::NSGAII) =
    "NSGA-II[P=$(m.populationSize),x=$(m.crossoverRate),μ=$(m.mutationRate)]"
show(io::IO, m::NSGAII) = print(io, summary(m))

mutable struct NSGAIIIndividual <: AbstractIndividual
    variables::AbstractVector
    objectives::Union{AbstractVector, Number, Nothing}
    constraints::Union{AbstractVector, Nothing}
    rank
    distance

    NSGAIIIndividual(var) = new(var, nothing, nothing)
    NSGAIIIndividual(var, obj) = new(var, obj, nothing)
    NSGAIIIndividual(var, obj, cst) = new(var, obj, cst)
    NSGAIIIndividual(var, obj, cst, rank, distance) = new(var, obj, cst, rank, distance)
end
variables(i::NSGAIIIndividual) = i.variables
objectives(i::NSGAIIIndividual) = i.objectives
constraints(i::NSGAIIIndividual) = i.constraints
rank(i::NSGAIIIndividual) = i.rank
distance(i::NSGAIIIndividual) = i.distance
copy(i::NSGAIIIndividual) = NSGAIIIndividual(copy(i.variables), copy(i.objectives), copy(i.constraints), copy(i.rank), copy(i.distance))
NSGAIIIndividual(i::Individual) = NSGAIIIndividual(copy(i.variables), copy(i.objectives), copy(i.constraints), nothing, nothing)

mutable struct NSGAIIState <: AbstractOptimizerState
    iteration
    start_time
    stop_time
    metrics::ConvergenceMetrics # collection of convergence metrics

    population                  # population
    pfront                      # individuals of the first Pareto front
end
pfront(s::NSGAIIState) = s.pfront
value(s::NSGAIIState) = objectives(pfront(s))
minimizer(s::NSGAIIState) = variables(pfront(s))
copy(s::NSGAIIState) = NSGAIIState(s.iteration, s.start_time, s.stop_time, copy(s.metrics),
                                 copy(s.population),
                                 copy(s.pfront))

function ranks(P::Population)
    map(rank, P)
end
function distances(P::Population)
    map(distance, P)
end
function maxF(P::Population)
    maxF, _ = findmax(map(rank, P))
    maxF
end
function calcF(P::Population)
    F = [Vector{Int64}() for _ in 1:maxF(P)]
    for i = 1:length(P)
        push!(F[rank(P[i])], i)
    end
    F
end

function initial_state(method::NSGAII, objective, population, options)
    population = map(NSGAIIIndividual, population)
    value!(objective, population)
    nondominatedsort!(population)
    crowding_distance!(population)
    pfront = filter(i -> rank(i) == 1, population)
    return NSGAIIState(0, 0, 0, copy(method.metrics), population, pfront)
end

function update_state!(
    method::NSGAII,
    state,
    objective,
    constraints,
    options
)
    parents = state.population

    populationSize = method.populationSize
    rng = options.rng

    F = calcF(parents)
    for f = F
        d = sortperm(distances(parents[f]))
        f = shuffle(f)
    end
    selected = collect(Iterators.flatten(F))

    # select offspring
    # selected = method.selection(hcat(ranks(parents), distances(parents))', populationSize; rng = rng)

    offspring = copy(parents)

    # perform mating
    recombine!(offspring, parents, selected, method)

    # perform mutation
    mutate!(offspring, method, constraints, rng = rng)

    # handle constraint
    map(x -> apply!(constraints, variables(x)), offspring)

    evaluate!(objective, offspring, constraints)

    population = vcat(parents, offspring)

    # calculate ranks & crowding for population
    nondominatedsort!(population)
    crowding_distance!(population)

    ### println(filter(d -> d == -Inf, map(distance, filter(i -> rank(i) == 1, population))))

    # select best individuals
    fitidx = Int[]
    F = calcF(population)
    for f in F
        if length(fitidx) + length(f) > populationSize
            idxs = sortperm(map(rank, population[f]))
            append!(fitidx, idxs[1:(populationSize-length(fitidx))])
            break
        else
            append!(fitidx, f)
        end
    end

    state.pfront = filter(i -> rank(i) == 1, population)

    # construct new parent population
    parents .= population[fitidx]

    return false
end


text:
## --- Load external packages
    if VERSION>=v"0.7"
        using Statistics
        using DelimitedFiles
        using SpecialFunctions
    end

    using Plots; gr(); default(fmt=:svg);
    using StatGeochem

## --- Plot and compare different estimates of continental coverage

    # Area of the continents
    area_cont_km2 = 1.4894E8
    # Excluding antarctica
    area_cc_km2 = 1.4894E8 - 1.4E7

    # Digitized from scan of Egyed (1956) Nature with PlotDigitizer.
    # Paleozoic, Mesozoic, and "Caniozoic" rescaled to their modern boundaries (541,252,66,0 Myr)
    egyed_termier_myr = [529.076, 507.928, 486.289, 458.775, 440.531, 417.673, 394.023, 375.125, 361.855, 345.028, 330.946, 313.166, 306.193, 296.568, 282.748, 275.076, 258.413, 248.355, 227.428, 212.373, 198.061, 187.34, 175.453, 151.163, 139.336, 110.362, 85.9226, 73.1421, 64.2713, 52.789, 41.7133, 20.8795, 9.55939, 3.51455]
    egyed_termier_km2 = 1E6 .* [34.7499, 31.6007, 33.1326, 28.0706, 45.8783, 34.0062, 52.383, 55.2615, 28.7052, 42.7889, 32.8643, 32.09, 42.7662, 17.5917, 31.5325, 38.1309, 20.8666, 19.1131, 22.9933, 16.2493, 25.0294, 26.9874, 27.3472, 15.349, 29.2083, 36.0758, 34.9838, 30.7257, 15.583, 18.1485, 20.9982, 10.9498, 3.95325, 0.508505]

    egyed_strahov_myr = [484.608, 399.113, 329.952, 309.219, 279.509, 256.341, 225.505, 177.268, 151.933, 92.0874, 51.7846, 8.63077]
    egyed_strahov_km2 = 1E6 .* [65.7566, 72.6922, 58.6154, 53.227, 49.878, 19.2318, 26.9204, 39.7824, 24.1789, 54.748, 33.151, 1.23662]

    # From Ronov, 1992
    ronov_myr =  [3.9605, 14.1815, 28.465, 44.95, 61.0, 83.25, 122.75, 154.25, 168.8, 187.7, 219.15, 242.1, 249.685, 262.235, 285.6, 311.05, 341.05, 370.8, 388.0, 406.25, 423.3, 435.6, 451.1, 464.2, 477.7, 493.2, 506.0, 526.0]
    ronov_coverage = 1/100 .* [5, 10, 9, 18, 15, 29, 24, 24, 22, 18, 18, 18, 17, 18, 22, 25, 28, 33, 35, 28, 28, 38, 34, 39, 37, 33, 32, 30]

    # Macrostrat
    data = readdlm("../macrostrat/macrostrat.csv",',')
    macrostrat = elementify(data)

    # Plot results
    h = plot(xflip=true, xlims=(0,1600), ylims=(0, 0.5), fg_color_legend=:white, legend=:topleft, framestyle=:box)
    plot!(h, egyed_strahov_myr,egyed_strahov_km2/area_cont_km2,color=:lightblue,label="Egyed/Strahov (global)")
    plot!(h, egyed_strahov_myr,egyed_strahov_km2/area_cont_km2,color=:lightblue,label="", seriestype=:scatter, markersize=2, markerstrokecolor=:auto)

    plot!(h, egyed_termier_myr,egyed_termier_km2/area_cc_km2,color=:darkblue,label="Egyed/Termier (global)")
    plot!(h, egyed_termier_myr,egyed_termier_km2/area_cc_km2,color=:darkblue,label="", seriestype=:scatter, markersize=2, markerstrokecolor=:auto)

    plot!(h, ronov_myr, ronov_coverage, color=:black, label = "Ronov (global)")
    plot!(h, ronov_myr, ronov_coverage, color=:black, label = "", seriestype=:scatter,  markersize=2, markerstrokecolor=:auto)

    plot!(h, macrostrat["Time"], macrostrat["Coverage"], color=:darkred, label = "Macrostrat (N. Am.)")
    plot!(h, xlabel="Age (Ma)", ylabel="Fraction covered")
    savefig(h, "CoverageFraction.pdf")
    display(h);

max_stars_repo_name: brenhinkeller/GreatUnconformity

max_issues_repo_name: brenhinkeller/GreatUnconformity

max_forks_repo_name: brenhinkeller/GreatUnconformity


text:
immutable BED <: Bio.IO.FileFormat end
export BED

immutable BigWig <: Bio.IO.FileFormat end
immutable BigBed <: Bio.IO.FileFormat end
export BigWig, BigBed


text:
module UserSitgnature
import MLStyle.AbstractPatterns: PComp
export and, or, literal, and, wildcard, decons,
       guard, effect, self

PatternUse = NamedTuple{
    (
        :and,
        :or,
        :literal,
        :wildcard,
        :decons,
        :guard,
        :effect
    )
}

@nospecialize
_empty_ntuple = NamedTuple()
and(ps::Vector, config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    xs = [p(impl) for p in ps]
    impl.and(xs, config)
end

or(ps::Vector, config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    xs = [p(impl) for p in ps]
    impl.or(xs, config)
end
literal(val, config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    impl.literal(val, config)
end
wildcard(config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    impl.wildcard(config)
end

decons(comp::PComp, extract::Function, ps, config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    xs = [p(impl) for p in ps]
    impl.decons(comp, extract, xs, config)
end

guard(pred, config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    impl.guard(pred, config)
end

effect(ctx_perf, config::NamedTuple=_empty_ntuple) = function apply(impl::PatternUse)
    impl.effect(ctx_perf, config)
end

const self = (
    and = and,
    or = or,
    literal = literal,
    wildcard = wildcard,
    decons = decons,
    guard = guard,
    effect = effect
)

end
@specialize


text:
module TestCopulas

using Copulas
using Base.Test

## checkSameLength
##----------------

@test Copulas.checkSameLength([0.3, 0.4], [0.8, 0.9])
@test_throws Exception Copulas.checkSameLength([0.3, 0.4], [0.8])

## params
##-------

cop = Copulas.GaussianPC_Cpp([-0.5])
rho = Copulas.params(cop)
@test rho == [-0.5]

## getFamAndParams
##----------------

famId, rho = Copulas.getFamAndParams(cop)
@test famId == 10.0
@test rho == [-0.5]

## triangular indexing
##--------------------

inds = Copulas.triang2sub(4, 1)
@test (1, 2) == inds

inds = Copulas.triang2sub(4, 2)
@test (1, 3) == inds

inds = Copulas.triang2sub(4, 3)
@test (1, 4) == inds

inds = Copulas.triang2sub(6, 5)
@test (1, 6) == inds

inds = Copulas.triang2sub(4, 4)
@test (2, 3) == inds

inds = Copulas.triang2sub(8, 24)
@test (5, 7) == inds

inds = Copulas.triang2sub(7, 12)
@test (3, 4) == inds

## reverse direction
ind = Copulas.sub2triang(4, 1, 3)
@test ind == 2

ind = Copulas.sub2triang(4, 3, 4)
@test ind == 6

ind = Copulas.sub2triang(6, 3, 5)
@test ind == 11

@test_throws Exception Copulas.sub2triang(5, 3, 3)
@test_throws Exception Copulas.sub2triang(5, 3, 2)
@test_throws Exception Copulas.sub2triang(5, 3, 6)

end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

equations = IdealGlmMhdEquations3D(5/3)

initial_condition = initial_condition_convergence_test

volume_flux = flux_central
solver = DGSEM(polydeg=5, surface_flux=flux_hll,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

##############################################################################
# Create the mesh

# Note, we use the domain [-1, 1]^3 for the Alfvén wave convergence test case so the
# warped mapping simplifies (quite a bit)

# Mapping as described in https://arxiv.org/abs/2012.12040
function mapping(xi, eta, zeta)
  y = eta + 0.125 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta))

  x = xi + 0.125 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta))

  z = zeta + 0.125 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta))

  return SVector(x, y, z)
end

cells_per_dimension = (4, 4, 4)
mesh = CurvedMesh(cells_per_dimension, mapping)

###############################################################################
# create the semi discretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

cfl = 1.2
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
import MixedSubdivisions
const MS = MixedSubdivisions

@testset "Binomial System Test" begin

    if Sys.iswindows()
        @testset "GMP set Int64 on 32-bit systems" begin
            x = BigInt()
            a = round(UInt64, rand() * typemax(Int32))
            for y in [a, a + typemax(Int32), a + UInt64(4) * typemax(UInt32)]
                HC.set_ui64!(x, y)
                @test x == y
                HC.set_si64!(x, Int64(y))
                @test x == Int64(y)
                HC.set_si64!(x, -Int64(y))
                @test x == -Int64(y)
            end
            HC.set_si64!(x, Int64(typemin(Int32)))
            @test x == typemin(Int32)
            HC.set_si64!(x, typemin(Int64))
            @test x == typemin(Int64)
        end
    end

    @testset "Hermite Normal Form" begin
        A = [0 3 1 0; -2 2 -1 2; 1 -1 2 3; -3 3 3 2]
        H, U = HC.hnf(A)
        @test A * U == H

        D = zeros(Int, 1, 1)
        D[1, 1] = -2
        H, U = HC.hnf(D)
        @test H[1, 1] == 2
        @test U[1, 1] == -1

        # overflow A
        A = [
            2 4 4 2 -4 -3 -4 3 0 3
            -4 -2 -5 -5 1 0 -3 0 -1 -2
            -3 5 -3 1 4 3 -1 -4 -1 0
            1 3 3 4 -2 -3 -2 -5 -4 -3
            2 2 0 2 -4 4 3 -4 -2 -4
            -3 -3 1 1 -4 4 0 4 4 -4
            -3 -1 -1 -1 2 -4 -3 -4 4 -4
            0 1 -2 4 5 4 3 1 -5 2
            -5 -4 -5 3 1 5 0 -3 -3 -1
            2 -5 -3 -1 -1 5 -4 5 -3 1
        ]
        @test_throws OverflowError HC.hnf(A)
        H, U = HC.hnf(A, BigInt)
        @test A * U == H
    end

    @testset "Solving (unit b)" begin
        A = [0 1 1 0 -1; 0 0 0 1 -1; -1 0 -1 0 -1; 0 -1 0 -1 -1; 1 0 0 0 -1]
        b = [
            0.9053223983046926 + 0.4247250347316951im,
            0.7000429487508004 - 0.7141007421255663im,
            0.018811552539476483 - 0.9998230470893611im,
            -0.9983533473373446 + 0.057363698105327716im,
            0.9999409355072137 - 0.010868555421874721im,
        ]
        X = HC.solve!(HC.BinomialSystemSolver(A, b))
        @test maximum(eachcol(X)) do x
            maximum(abs.([prod(x .^ a) for a in eachcol(A)] - b) ./ abs.(b))
        end < 1e-12

        A = [
            3 2 3 4 2 5
            6 -3 8 -3 8 7
            -2 -5 7 3 6 5
            1 2 3 4 5 6
            1 0 2 0 5 0
            1 2 0 -3 0 5
        ]
        b = randn(ComplexF64, 6)
        b ./= abs.(b)
        X = HC.solve!(HC.BinomialSystemSolver(A, b))
        @test maximum(eachcol(X)) do x
            maximum(abs.([prod(x .^ a) for a in eachcol(A)] - b) ./ abs.(b))
        end < 1e-12
    end

    @testset "General solving" begin
        A = [0 1 1 0 -1; 0 0 0 1 -1; -1 0 -1 0 -1; 0 -1 0 -1 -1; 1 0 0 0 -1]
        b = randn(ComplexF64, size(A, 1))
        X = HC.solve!(HC.BinomialSystemSolver(A, b))
        @test maximum(eachcol(X)) do x
            maximum(abs.([prod(x .^ a) for a in eachcol(A)] - b) ./ abs.(b))
        end < 1e-12

        A = [
            3 2 3 4 2 5
            6 -3 8 -3 8 7
            -2 -5 7 3 6 5
            1 2 3 4 5 6
            1 0 2 0 5 0
            1 2 0 -3 0 5
        ]
        b = randn(ComplexF64, 6)
        X = HC.solve!(HC.BinomialSystemSolver(A, b))
        @test maximum(eachcol(X)) do x
            maximum(abs.([prod(x .^ a) for a in eachcol(A)] - b) ./ abs.(b))
        end < 1e-12
    end

    @testset "MixedCell" begin
        f = cyclic(5)
        supp = first.(ModelKit.exponents_coefficients.(f, Ref(variables(f))))

        cells, lift = MS.fine_mixed_cells(supp)
        coeffs = map(l -> randn(ComplexF64, length(l)), lift)

        cell = first(cells)
        BSS = HC.BinomialSystemSolver(5)
        X = HC.solve(BSS, supp, coeffs, cell)
        @test maximum(eachcol(X)) do x
            maximum(abs.([prod(x .^ a) for a in eachcol(BSS.A)] - BSS.b) ./ abs.(BSS.b))
        end < 1e-12
    end
end


text:
import Stripeline
const Sl = Stripeline

##########################################################################
# Simple ADC, with gain 1, no offset and 20 bit of dynamic range

adc = Sl.ADC()

@test Sl.adc_response(adc, 0.0) == 0
@test Sl.adc_response(adc, 1.0) == 1
@test Sl.adc_response(adc, -1.0) == -1
@test Sl.adc_response(adc, 2^21) == 2^19

@test Sl.adc_inv_response(adc, 0) ≈ 0.0
@test Sl.adc_inv_response(adc, 1) ≈ 1.0
@test Sl.adc_inv_response(adc, -1) ≈ -1.0

##########################################################################
# The same as above, but with a 2× gain and an offset of 5 K

adc = Sl.ADC(offset_k = 5.0, gain_k_over_adu = 2.0)

@test Sl.adc_response(adc, 0.0) == -10
@test Sl.adc_response(adc, 1.0) == -8
@test Sl.adc_response(adc, -1.0) == -12
@test Sl.adc_response(adc, 2^21) == 2^19

@test Sl.adc_inv_response(adc, 0) ≈ 5.0
@test Sl.adc_inv_response(adc, 1) ≈ 5.5
@test Sl.adc_inv_response(adc, -1) ≈ 4.5

##########################################################################
# Non-linearity checks

adc = Sl.ADC(min_output_adu = -10,
    max_output_adu = 10,
    non_linearities_x_adu = [-10, 0, 10],
    non_linearities_y_adu = [0, 10, 0],
)

@test Sl.adc_response.(Ref(adc), [-10., -5., 0., 5., 10.]) == [-10, 0, 10, 10, 10]

##########################################################################
# Additional checks

adc = Sl.ADC()

@test Sl.adc_filter(adc, 1.1) ≈ 1.0


text:
@add_problem mip function mip_lp_fallback_interface(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable()
    p = minimize(x, x >= 4.3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4.3 atol = atol rtol = rtol
    end

    x = Variable(2)
    p = minimize(norm(x, 1), x[1] >= 4.3; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4.3 atol = atol rtol = rtol
    end
end

@add_problem mip function mip_integer_variables(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(:Int)
    p = minimize(x, x >= 4.3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol = atol rtol = rtol
    end

    x = Variable(2, :Int)
    p = minimize(sum(x), x >= 4.3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 10 atol = atol rtol = rtol
    end

    x = Variable(:Int)
    y = Variable()
    p = minimize(sum(x + y), x >= 4.3, y >= 7; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 12 atol = atol rtol = rtol
    end

    x = Variable(2, :Int)
    p = minimize(norm(x, 1), x[1] >= 4.3; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol = atol rtol = rtol
    end

    x = Variable(2, :Int)
    p = minimize(sum(x), x[1] >= 4.3, x >= 0; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol = atol rtol = rtol
    end

    x = Variable(2, :Int)
    p = minimize(sum(x), x >= 0.5; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
    end
end

@add_problem mip function mip_binary_variables(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, :Bin)
    p = minimize(sum(x), x >= 0.5; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
    end

    x = Variable(2, :Bin)
    p = minimize(sum(x), x[1] >= 0.5, x >= 0; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
    end
end


text:
"""
Returns a signal, which is true as long as the window is open.
returns `Signal{Bool}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#gaade9264e79fae52bdb78e2df11ee8d6a)
"""
function window_open(window, s::Signal{Bool}=Signal(true))
    GLFW.SetWindowCloseCallback(window, (window,) -> begin
        push!(s, false)
    end)
    s
end
"""
Size of the window. Must not correlate to the real pixel size.
This is why there is also framebuffer_size.
returns `Signal{Vec{2,Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#gaaca1c2715759d03da9834eac19323d4a)
"""
function window_size(window, s::Signal{Vec{2,Int}}=Signal(Vec{2,Int}( GLFW.GetWindowSize(window))))
    GLFW.SetWindowSizeCallback(window, (window, w::Cint, h::Cint,) -> begin
        push!(s, Vec{2,Int}(w, h))
    end)
    s
end
"""
Size of window in pixel.
returns `Signal{Vec{2,Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#ga311bb32e578aa240b6464af494debffc)
"""
function framebuffer_size(window, s::Signal{Vec{2, Int}}=Signal(Vec{2, Int}(GLFW.GetFramebufferSize(window))))
    GLFW.SetFramebufferSizeCallback(window, (window, w::Cint, h::Cint) -> begin
        push!(s, Vec(Int(w), Int(h)))
    end)
    s
end
"""
Position of the window in screen coordinates.
returns `Signal{Vec{2,Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#ga1c36e52549efd47790eb3f324da71924)
"""
function window_position(window, s::Signal{Vec{2,Int}}=Signal(Vec(0,0)))
    GLFW.SetWindowPosCallback(window, (window, x::Cint, y::Cint,) -> begin
        push!(s, Vec(Int(x), Int(y)))
    end)
    s
end
"""
Registers a callback for the mouse buttons + modifiers
returns `Signal{NTuple{4, Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function keyboard_buttons(window, s::Signal{NTuple{4, Int}}=Signal((0,0,0,0)))
    keydict = Dict{Int, Bool}()
    GLFW.SetKeyCallback(window, (window, button::Cint, scancode::Cint, action::Cint, mods::Cint) -> begin
        push!(s, (Int(button), Int(scancode), Int(action), Int(mods)))
    end)
    s
end
"""
Registers a callback for the mouse buttons + modifiers
returns an `Signal{NTuple{3, Int}}`,
containing the pressed button the action and modifiers.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function mouse_buttons(window, s::Signal{NTuple{3, Int}}=Signal((0,0,0)))
    GLFW.SetMouseButtonCallback(window, (window, button::Cint, action::Cint, mods::Cint) -> begin
        push!(s, (Int(button), Int(action), Int(mods)))
    end)
    s
end
"""
Registers a callback for drag and drop of files.
returns `Signal{Vector{Compat.UTF8String}}`, which are absolute file paths
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#gacc95e259ad21d4f666faa6280d4018fd)
"""
function dropped_files(window, s::Signal{Vector{Compat.UTF8String}}=Signal(Compat.UTF8String[]))
    GLFW.SetDropCallback(window, (window, files) -> begin
        push!(s, map(Compat.String, files))
    end)
    s
end

"""
Registers a callback for keyboard unicode input.
returns an `Signal{Vector{Char}}`,
containing the pressed char. Is empty, if no key is pressed.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function unicode_input(window, s::Signal{Vector{Char}}=Signal(Char[]))
    GLFW.SetCharCallback(window, (window, c::Char) -> begin
        push!(s, Char[c])
        push!(s, Char[])
    end)
    s
end
"""
Registers a callback for the mouse cursor position.
returns an `Signal{Vec{2, Float64}}`,
which is not in screen coordinates, with the upper left window corner being 0
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function cursor_position(window, s::Signal{Vec{2, Float64}}=Signal(Vec(0.,0.)))
    GLFW.SetCursorPosCallback(window, (window, x::Cdouble, y::Cdouble) -> begin
        push!(s, Vec{2, Float64}(x, y))
    end)
    s
end
"""
Registers a callback for the mouse scroll.
returns an `Signal{Vec{2, Float64}}`,
which is an x and y offset.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#gacc95e259ad21d4f666faa6280d4018fd)
"""
function scroll(window, s::Signal{Vec{2, Float64}}=Signal(Vec(0.,0.)))
    GLFW.SetScrollCallback(window, (window, xoffset::Cdouble, yoffset::Cdouble) -> begin
        push!(s, Vec{2, Float64}(xoffset, yoffset))
        push!(s, Vec{2, Float64}(0))
    end)
    s
end
"""
Registers a callback for the focus of a window.
returns an `Signal{Bool}`,
which is true whenever the window has focus.
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#ga6b5f973531ea91663ad707ba4f2ac104)
"""
function hasfocus(window, s::Signal{Bool}=Signal(false))
    GLFW.SetWindowFocusCallback(window, (window, focus::Bool) -> begin
        push!(s, focus)
    end)
    s
end
"""
Registers a callback for if the mouse has entered the window.
returns an `Signal{Bool}`,
which is true whenever the cursor enters the window.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga762d898d9b0241d7e3e3b767c6cf318f)
"""
function entered_window(window, s::Signal{Bool}=Signal(false))
    GLFW.SetCursorEnterCallback(window, (window, entered::Bool) -> begin
        push!(s, entered)
    end)
    s
end

"""
Takes a screen and registers a list of callback functions.
Returns a dict{Symbol, Signal}(name_of_callback => signal)
"""
function register_callbacks(window::GLFW.Window, callbacks::Vector{Function})
    tmp = map(callbacks) do f
        (Symbol(last(split(string(f),"."))), f(window))
    end
    Dict{Symbol, Any}(tmp)
end
function register_callbacks(window::Screen, callbacks::Vector{Function})
    register_callbacks(window.nativewindow, callbacks)
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds Base
# The underlying abstraction layer for DynamicBounds.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/problem.jl
# Defines attributes and functions used by relaxation problems.
#############################################################################

"""
$(TYPEDEF)

Abstract supertype for all relaxation problems.
"""
abstract type AbstractDERelaxProblem  end

"""
$(TYPEDEF)

Abstract supertype for all problem attributes.
"""
abstract type AbstractRelaxProblemAttribute end

#This allows to use attributes in broadcast calls without the need to
# embed it in a `Ref`
Base.broadcastable(attribute::AbstractRelaxProblemAttribute) = Ref(attribute)

"""
$(TYPEDEF)

A problem attribute used to check whether state bounds are defined.
"""
struct HasStateBounds <: AbstractRelaxProblemAttribute end

"""
$(TYPEDEF)

A problem attribute used to check whether state bounds are defined as constant.
"""
struct HasConstantStateBounds <: AbstractRelaxProblemAttribute end

"""
$(TYPEDEF)

A problem attribute used to store constant state bounds.

$(TYPEDFIELDS)
"""
struct ConstantStateBounds <: AbstractRelaxProblemAttribute
    "Lower state variable bounds"
    xL::Vector{Float64}
    "Upper state variable bounds"
    xU::Vector{Float64}
    function ConstantStateBounds(xL::Vector{Float64}, xU::Vector{Float64})
        @assert length(xL) == length(xU)
        new(xL, xU)
    end
end
ConstantStateBounds() = ConstantStateBounds(Float64[], Float64[])

"""
$(TYPEDEF)

A problem attribute used to check whether variable state bounds are defined.
"""
struct HasVariableStateBounds <: AbstractRelaxProblemAttribute end

"""
$(TYPEDEF)

A problem attribute used to store time-varying state bounds.

$(TYPEDFIELDS)
"""
struct VariableStateBounds{F1 <: Function, F2 <: Function} <: AbstractRelaxProblemAttribute
    "Lower state variable bounds"
    xL::F1
    "Upper state variable bounds"
    xU::F2
end
VariableStateBounds() = VariableStateBounds(Base.isempty, Base.isempty)

"""
$(TYPEDEF)

A problem attribute used to check whether a user-defined jacobian was specified.
"""
struct HasUserJacobian <: AbstractRelaxProblemAttribute end

"""
$(TYPEDEF)

An abstract supertype for all manners of path constraints.
"""
abstract type AbstractPathConstraint <: AbstractRelaxProblemAttribute end

"""
$(TYPEDEF)

An object used to specify that the invariant `Ax(t,p) <= b` is valid fo all `t`.
"""
struct PolyhedralConstraint <: AbstractPathConstraint
    A::Array{Float64,2}
    b::Vector{Float64}
    function PolyhedralConstraint(A::Array{Float64,2}, b::Vector{Float64})
        @assert size(A,1) == length(b)
        return new(A, b)
    end
end
PolyhedralConstraint() = PolyhedralConstraint(zeros(Float64,1,1), zeros(Float64,1))


text:
# Represents a graph which with different lincomb functionality:
# It can add more than two matrices. Only used in code generation.
struct MultiLincombCompgraph{T}
    operations::Dict{Symbol,Symbol}
    parents::Dict{Symbol,NTuple{<:Any,Symbol}}
    coeffs::Dict{Symbol,NTuple{<:Any,T}}
    outputs::Vector{Symbol}
end

function MultiLincombCompgraph(g::Compgraph)
    T = eltype(g)
    extgraph = MultiLincombCompgraph(
        Dict{Symbol,Symbol}(),
        Dict{Symbol,NTuple{<:Any,Symbol}}(),
        Dict{Symbol,NTuple{<:Any,T}}(),
        Vector{Symbol}(),
    )
    Z = extract_sums(g)
    for k in keys(g.operations)
        if (g.operations[k] == :mult)
            add_mult!(extgraph, k, g.parents[k][1], g.parents[k][2])
        end
        if (g.operations[k] == :ldiv)
            add_ldiv!(extgraph, k, g.parents[k][1], g.parents[k][2])
        end
    end
    for s in Z
        coeff_list = s[1]
        symbol_list = s[2]
        key = s[3][end]
        p = size(coeff_list, 1)
        extgraph.operations[key] = :lincomb
        extgraph.coeffs[key] = NTuple{p,T}(coeff_list)
        extgraph.parents[key] = NTuple{p,Symbol}(symbol_list)
    end
    for k in g.outputs
        push!(extgraph.outputs, k)
    end
    return extgraph
end

"""
    sums=extract_sums(graph::Compgraph{T})

`sums::Vector{Tuple{Vector{T},Vector{Symbol},Vector{Symbol}}}`

Returns a representation of sums in `graph` which may potentially be merged, in a dot-fusion.
The vector `sums` contains a tuple for each of these sums. The three entries of
the tuple are:

  - a vector of `T` values that represent the coefficients of the summands;
  - a vector of `Symbol`s that correspond to the summands; and
  - a vector of intermediate `Symbol`s (i.e, nodes) that can be merged.

The first two vectors have the same number of entries, one for each element that
can be merged in the sum.
"""
function extract_sums(graph)
    coeff, nodes, merged, sums =
        find_mergeable_sums(graph, graph.outputs[1], [])
    return sums
end

# Return true if `node` has multiple parents. Such nodes cannot be freed.
function has_multiple_parents(graph, node)
    return sum(map(x -> any(x .== node), values(graph.parents))) > 1
end

function find_mergeable_sums(graph, node, processed, curr_coeff = 1)
    # Extract sums in the subgraph of `graph` with root `node`.
    #
    # The function accepts four parameters:
    #     * `graph`: current graph.
    #     * `node`: current node.
    #     * `processed`: set of nodes the algorithm has already processed..
    #     * `curr_coeff`: coefficient `node` if parent is lincomb, 1 otherwise.
    #
    # The functions returns four vectors:
    #     * `pcoeffs`: coefficients of sum currently being constructed.
    #     * `pnodes`: corresponding nodes of sum currently being constructed.
    #     * `pmerged`: nodes merged in the current sum (these will disappear).
    #     * `sums`: completely extracted sums.
    #
    # The algorithm starts from the first output node, which is seen as the root
    # of a spanning tree with edges defined in `graph.parents`. The graph may
    # have cycles, but the vector `processed` ensures that each node is
    # processed only once, the first time it is visited.
    #
    # If `node` is an input node, that is, a node without parents, then the
    # algorithm returns four empty vectors, as 1) the subtree rooted at `node`,
    # being empty, does not have extracted sums, and 2) no sum is being
    # constructed.
    #
    # If the node is not a leaf, the function is called recursively on the
    # two parents, and three cases are possible:
    #
    # 1) If `node` is not a `:lincomb`, then the function returns the union of
    # the sums extracted in the two subgraphs rooted at the parents. If either
    # parent is a `:lincomb` the sum that parent was constructing is added to
    # the vector of extracted sums. The three other output vectors are empty.
    #
    # 2) If `node` is a `:lincomb`, is parent to only one node, and is not an
    # output node, then the function merges the two sums being constructed by
    # the parents, if any, adds `node` to it, and returns the data accordingly.
    # The union of the vectors of extracted sums is also returned.
    #
    # 3) Otherwise, `node` is added to the union of the (possibly empty) sums
    # being constructed by the parents. In particular, the algorithm will add
    # `node` to the list of nodes to be merged, will updated the coefficients of
    # the constructed sum accordingly, and will add the current sum to the
    # vector of extracted sums, which will also include the union of the sum
    # extracted in the subgraph rooted at the parents.

    if !(node in keys(graph.operations)) || (node in processed)
        # Nothing to do for leaf nodes and nodes already processed.
        return Float64[], Symbol[], Symbol[], []
    else
        # Call function recursively on both parents.
        push!(processed, node)
        (parent1, parent2) = graph.parents[node]
        curr_lincomb = graph.operations[node] == :lincomb
        (coeff1, coeff2) = curr_lincomb ? graph.coeffs[node] : (1, 1)
        pcoeffs1, pnodes1, pmerged1, sums1 =
            find_mergeable_sums(graph, parent1, processed, coeff1)
        pcoeffs2, pnodes2, pmerged2, sums2 =
            find_mergeable_sums(graph, parent2, processed, coeff2)
        if curr_lincomb
            # Grow the sum by adjoining terms coming from parents, if any.
            has_multiple_parents(graph, node) && (curr_coeff = 1)
            new_coeffs = vcat(curr_coeff * pcoeffs1, curr_coeff * pcoeffs2)
            new_nodes = vcat(pnodes1, pnodes2)
            new_merged = vcat(pmerged1, pmerged2)
            # Lincomb parents are added to the mergeable nodes.
            # Non-lincomb parents are added to the sum and put, and their
            # coefficients are added to the vector of coefficients.
            if haskey(graph.operations, parent1) &&
               graph.operations[parent1] == :lincomb &&
               !has_multiple_parents(graph, parent1)
                new_merged = vcat(new_merged, parent1)
            else
                new_coeffs = vcat(new_coeffs, curr_coeff * coeff1)
                new_nodes = vcat(new_nodes, parent1)
            end
            if haskey(graph.operations, parent2) &&
               graph.operations[parent2] == :lincomb &&
               !has_multiple_parents(graph, parent2)
                new_merged = vcat(new_merged, parent2)
            else
                new_coeffs = vcat(new_coeffs, curr_coeff * coeff2)
                new_nodes = vcat(new_nodes, parent2)
            end
            if node in graph.outputs || has_multiple_parents(graph, node)
                return Float64[],
                Symbol[],
                Symbol[],
                vcat(
                    sums1,
                    sums2,
                    (new_coeffs, new_nodes, vcat(new_merged, node)),
                )
            else
                return new_coeffs, new_nodes, new_merged, vcat(sums1, sums2)
            end
        else
            # Current node is not a lincomb node.
            sums = vcat(sums1, sums2)
            # Add sum of lincomb parents, if any, to sums.
            sums =
                !isempty(pcoeffs1) ?
                vcat(sums, (pcoeffs1, pnodes1, vcat(pmerged1, parent1))) : sums
            sums =
                !isempty(pcoeffs2) ?
                vcat(sums, (pcoeffs2, pnodes2, vcat(pmerged2, parent2))) : sums
            return Float64[], Symbol[], Symbol[], sums
        end
    end
end


text:
module TestFlattenedArray

include("preamble.jl")

function makedata(V::Type, M::Integer, N::Integer)
    dims = testdims(M + N)
    inner_size, outersize = tuple_split(dims, Val(M))

    nested = Array{Array{V,M},N}(undef, outersize...)
    for I in eachindex(nested)
        nested[I] = rand!(Array{V,M}(undef, inner_size...))
    end

    return (
        nested = nested,
        flat = reshape(reduce(hcat, nested), dims),
        dims = dims,
        inner_size = inner_size,
        outersize = outersize,
        inner_axes = axes(first(nested)),
        outeraxes = axes(nested),
    )
end

@testset "M = $M, N = $N, V = $V" for M in (0, 2), N = 1:2, V in (Float64,)
    @testset "constructors" begin
        data = makedata(V, M, N)
        Expected = FlattenedArray{V,M+N,M,N,typeof(data.nested),typeof(data.inner_axes)}
        @test flatten(data.nested) isa Expected
        @test_inferred flatten(data.nested)

        @test flatten(data.flat) === data.flat
        @test flatten(data.nested) == data.flat
    end

    @testset "basic" begin
        data = makedata(V, M, N)
        F = flatten(data.nested)
        F[:] .= 1:length(F)
        len = prod(data.inner_size)
        for i = 1:length(data.nested)
            a = data.nested[i]
            offs = (i - 1) * len + 1
            @test F[offs:offs+len-1] == vec(data.nested[i])
        end
    end

    let V = V, M = M, N = N
        test_array_AB() do
            data = makedata(V, M, N)
            A = FlattenedArray(data.nested)
            B = data.flat
            return A, B
        end
    end
end

@testset "aliasing" begin
    x1 = rand(2, 3)
    x2 = rand(2, 3)
    x3 = rand(2, 3)
    F = flatten([x1, x2])

    @test dataids(F) === (dataids(F.parent)..., Iterators.flatten(dataids.(F.parent))...)

    @test mightalias(F, x1)
    @test !mightalias(F, rand(2, 3))

    @test mightalias(x1, F)
    @test !mightalias(rand(2, 3), F)

    @test mightalias(F, flatten([x1, x3]))
    @test !mightalias(F, flatten([x3, x3]))

    let F2 = unalias(x1, F)
        @test !mightalias(x1, F2)
        @test !mightalias(F2, x1)
        @test F2.parent == F.parent
        @test F2.parent !== F.parent
        @test all(zip(F2.parent, F.parent)) do (f2, f)
            f2 !== f
        end
    end
end

end # module


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+b*csc(c+d*x)^2)^4, x) == :(x*a^4+-1//7*b^4*d^-1*cot(c+d*x)^7+-1//3*b^2*d^-1*cot(c+d*x)^3*(3*b^2+6*a^2+8*a*b)+-1//5*b^3*d^-1*cot(c+d*x)^5*(3b+4a)+-1*b*d^-1*(b+2a)*(b^2+2*a^2+2*a*b)*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^3, x) == :(x*a^3+-1//5*b^3*d^-1*cot(c+d*x)^5+-1*b*d^-1*(b^2+3*a^2+3*a*b)*cot(c+d*x)+-1//3*b^2*d^-1*cot(c+d*x)^3*(2b+3a))
@test integrate((a+b*csc(c+d*x)^2)^2, x) == :(x*a^2+-1//3*b^2*d^-1*cot(c+d*x)^3+-1*b*d^-1*(b+2a)*cot(c+d*x))
@test integrate(a+b*csc(c+d*x)^2, x) == :(a*x+-1*b*d^-1*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^-1, x) == :(x*a^-1+-1*a^-1*b^(1/2)*d^-1*(a+b)^-1//2*arctan(b^-1//2*(a+b)^(1/2)*tan(c+d*x)))
@test integrate((a+b*csc(c+d*x)^2)^-2, x) == :(x*a^-2+(1/2)*b*a^-1*d^-1*(a+b)^-1*(a+b+b*cot(c+d*x)^2)^-1*cot(c+d*x)+(1/2)*a^-2*b^(1/2)*d^-1*(a+b)^-3//2*(2b+3a)*arctan(b^(1/2)*(a+b)^-1//2*cot(c+d*x)))
@test integrate((a+b*csc(c+d*x)^2)^-3, x) == :(x*a^-3+1//4*b*a^-1*d^-1*(a+b)^-1*(a+b+b*cot(c+d*x)^2)^-2*cot(c+d*x)+1//8*a^-3*b^(1/2)*d^-1*(a+b)^-5//2*(8*b^2+15*a^2+20*a*b)*arctan(b^(1/2)*(a+b)^-1//2*cot(c+d*x))+1//8*b*a^-2*d^-1*(a+b)^-2*(a+b+b*cot(c+d*x)^2)^-1*(4b+7a)*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^-4, x) == :(x*a^-4+1//6*b*a^-1*d^-1*(a+b)^-1*(a+b+b*cot(c+d*x)^2)^-3*cot(c+d*x)+1//16*a^-4*b^(1/2)*d^-1*(a+b)^-7//2*(16*b^3+35*a^3+56*a*b^2+70*b*a^2)*arctan(b^(1/2)*(a+b)^-1//2*cot(c+d*x))+1//16*b*a^-3*d^-1*(a+b)^-3*(a+b+b*cot(c+d*x)^2)^-1*(8*b^2+19*a^2+22*a*b)*cot(c+d*x)+1//24*b*a^-2*d^-1*(a+b)^-2*(a+b+b*cot(c+d*x)^2)^-2*(6b+11a)*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^5//2, x) == :(-1*a^5//2*d^-1*arctan(a^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1//4*b*d^-1*(a+b+b*cot(c+d*x)^2)^3//2*cot(c+d*x)+-1//8*b^(1/2)*d^-1*(3*b^2+15*a^2+10*a*b)*arctanh(b^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1//8*b*d^-1*(a+b+b*cot(c+d*x)^2)^(1/2)*(3b+7a)*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^3//2, x) == :(-1*a^3//2*d^-1*arctan(a^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1//2*b*d^-1*(a+b+b*cot(c+d*x)^2)^(1/2)*cot(c+d*x)+-1//2*b^(1/2)*d^-1*(b+3a)*arctanh(b^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x)))
@test integrate((a+b*csc(c+d*x)^2)^(1/2), x) == :(-1*a^(1/2)*d^-1*arctan(a^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1*b^(1/2)*d^-1*arctanh(b^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x)))
@test integrate((a+b*csc(c+d*x)^2)^-3//2, x) == :(-1*a^-3//2*d^-1*arctan(a^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+b*a^-1*d^-1*(a+b)^-1*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^-5//2, x) == :(-1*a^-5//2*d^-1*arctan(a^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+1//3*b*a^-1*d^-1*(a+b)^-1*(a+b+b*cot(c+d*x)^2)^-3//2*cot(c+d*x)+1//3*b*a^-2*d^-1*(a+b)^-2*(a+b+b*cot(c+d*x)^2)^-1//2*(3b+5a)*cot(c+d*x))
@test integrate((a+b*csc(c+d*x)^2)^-7//2, x) == :(-1*a^-7//2*d^-1*arctan(a^(1/2)*(a+b+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+1//5*b*a^-1*d^-1*(a+b)^-1*(a+b+b*cot(c+d*x)^2)^-5//2*cot(c+d*x)+1//15*b*a^-3*d^-1*(a+b)^-3*(a+b+b*cot(c+d*x)^2)^-1//2*(15*b^2+33*a^2+40*a*b)*cot(c+d*x)+1//15*b*a^-2*d^-1*(a+b)^-2*(a+b+b*cot(c+d*x)^2)^-3//2*(5b+9a)*cot(c+d*x))
@test integrate((1+csc(x)^2)^3//2, x) == :(-1*arctan((2+cot(x)^2)^-1//2*cot(x))+-2*arcsinh((1/2)*2^(1/2)*cot(x))+-1//2*(2+cot(x)^2)^(1/2)*cot(x))
@test integrate((1+csc(x)^2)^(1/2), x) == :(-1*arcsinh((1/2)*2^(1/2)*cot(x))+-1*arctan((2+cot(x)^2)^-1//2*cot(x)))
@test integrate((1+csc(x)^2)^-1//2, x) == :(-1*arctan((2+cot(x)^2)^-1//2*cot(x)))
@test integrate((1+-1*csc(x)^2)^3//2, x) == :((1/2)*(-1*cot(x)^2)^(1/2)*cot(x)+(-1*cot(x)^2)^(1/2)*log(sin(x))*tan(x))
@test integrate((1+-1*csc(x)^2)^(1/2), x) == :((-1*cot(x)^2)^(1/2)*log(sin(x))*tan(x))
@test integrate((1+-1*csc(x)^2)^-1//2, x) == :(-1*(-1*cot(x)^2)^-1//2*cot(x)*log(cos(x)))
@test integrate((-1+csc(x)^2)^3//2, x) == :(-1//2*(cot(x)^2)^3//2*tan(x)+-1*(cot(x)^2)^(1/2)*log(sin(x))*tan(x))
@test integrate((-1+csc(x)^2)^(1/2), x) == :((cot(x)^2)^(1/2)*log(sin(x))*tan(x))
@test integrate((-1+csc(x)^2)^-1//2, x) == :(-1*(cot(x)^2)^-1//2*cot(x)*log(cos(x)))
@test integrate((-1+-1*csc(x)^2)^3//2, x) == :(-1*arctanh((-2+-1*cot(x)^2)^-1//2*cot(x))+-2*arctan((-2+-1*cot(x)^2)^-1//2*cot(x))+(1/2)*(-2+-1*cot(x)^2)^(1/2)*cot(x))
@test integrate((-1+-1*csc(x)^2)^(1/2), x) == :(arctan((-2+-1*cot(x)^2)^-1//2*cot(x))+arctanh((-2+-1*cot(x)^2)^-1//2*cot(x)))
@test integrate((-1+-1*csc(x)^2)^-1//2, x) == :(-1*arctanh((-2+-1*cot(x)^2)^-1//2*cot(x)))

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.7 (d trig)^m (a+b (c csc)^n)^p.jl


text:
# using Revise
# using CairoMakie, AbstractPlotting, ZigZagBoomerang, SparseArrays, LinearAlgebra

# Essence of the model section procedure.
# a = Dict((1,0)=>1, (0,0)=>2)
# push!(a, (2,0) => 3)
# (1,0) ∈ keys(a)
# key1 = (1,0)
# a[key1] = a[key1] + 1

# function ϕ(x, i, μ)
#     x[i] - μ[i]
# end
# κ = 1.0
# n = 100
# μ = zeros(n)
# x0 = randn(n)
# θ0 = rand((-1.0,1.0), n)
# T = 1000.0
# c = 10ones(n)
#@time trace0, _ = ZigZagBoomerang.spdmp(ϕ, 0.0, x0, θ0, T, c, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)
# @time trace0, _ = ZigZagBoomerang.sspdmp(ϕ, 0.0, x0, θ0, T, c, κ, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)

is_stuck(xi,vi) = xi == 0 && vi == 0 ? 1 : 0

"""
    model_trace(x0, v0, events)
given an initial position x0, v0 and a trace of events, assign to each
event, a vector k corresponding to the open sets as in Appendix 1.
`m0` is the initial model, `mm` is the trace of models at each event
"""
function model_trace(x0, v0, events)
    m0 = is_stuck.(x0, v0)
    m = copy(m0)
    mm = Vector{Int64}[]
    k = 1
    for event in events
        t, i, x, v = event
        m[i] = is_stuck(x, v)
        push!(mm, copy(m))
        k += 1
    end
    m0, mm
end

"""
    model_probs(m0, mm, events)
Given the initial model and the trace of models and trace of events,
computes the time spent by the sampler in each model explored. returns
a Ditcionary `probs` containing as `key` the model and as value its total time
"""
function model_probs(m0, mm, events)
    probs = Dict(m0=> events[1][1],)
    t = 0.0
    for i in eachindex(events)[2:end]
        tnew, i, x, v = events[i]
        δt = tnew - t
        t = tnew
        key = mm[i]
        if key ∈ keys(probs)
            probs[key] += δt
        else
            push!(probs, key => δt)
        end
    end
    probs
end

# m0, mm = model_trace(x0, θ0, trace0.events)
# model_probs(m0, mm, trace0.events)


text:
@testset "Auto-generated source files" begin
    @testset "R4" begin
        @testset "Patient" begin
            p = FHIRClient.R4Types.Patient()
            @test p.resourceType isa String
            @test p.resourceType == "Patient"
        end
    end
end