text:
@inline alwaystrue(x) = true

# Matcher patterns with Slot and Segment

# matches one term
# syntax:  ~x
struct Slot{P}
    name::Symbol
    predicate::P
end

Slot(s) = Slot(s, alwaystrue)

Base.isequal(s1::Slot, s2::Slot) = s1.name == s2.name

Base.show(io::IO, s::Slot) = (print(io, "~"); print(io, s.name))

# matches zero or more terms
# syntax: ~~x
struct Segment{F}
    name::Symbol
    predicate::F
end

ismatch(s::Segment, t) = s.predicate(t)

Segment(s) = Segment(s, alwaystrue)

Base.show(io::IO, s::Segment) = (print(io, "~~"); print(io, s.name))

makesegment(s::Symbol, keys) = (push!(keys, s); Segment(s))

function makesegment(s::Expr, keys)
    if !(s.head == :(::))
        error("Syntax for specifying a segment is ~~x::\$predicate, where predicate is a boolean function")
    end

    name = s.args[1]

    push!(keys, name)
    :(Segment($(QuoteNode(name)), $(esc(s.args[2]))))
end

makeslot(s::Symbol, keys) = (push!(keys, s); Slot(s))

function makeslot(s::Expr, keys)
    if !(s.head == :(::))
        error("Syntax for specifying a slot is ~x::\$predicate, where predicate is a boolean function")
    end

    name = s.args[1]

    push!(keys, name)
    :(Slot($(QuoteNode(name)), $(esc(s.args[2]))))
end

function makepattern(expr, keys, slots, splat = false)
    res = if expr isa Expr
        if expr.head === :call
            if expr.args[1] === :(~)
                if expr.args[2] isa Expr && expr.args[2].head === :call expr.args[2].args[1] == :(~)
                    # matches ~~x::predicate or ~~x::predicate...
                    return makesegment(expr.args[2].args[2], keys)
                elseif splat
                    # matches ~x::predicate...
                    return makesegment(expr.args[2], keys)
                else
                    # matches ~x::predicate
                    return makeslot(expr.args[2], keys)
                end
            else
                :(term($(map(x->makepattern(x, keys, slots), expr.args)...); type=Any))
            end
        elseif expr.head === :...
            makepattern(expr.args[1], keys, slots, true)
        elseif expr.head == :(::) && expr.args[1] in slots
            return splat ? makesegment(expr, keys) : makeslot(expr, keys)
        elseif expr.head === :ref
            :(term(getindex, $(map(x->makepattern(x, keys, slots), expr.args)...); type=Any))
        elseif expr.head === :$
            esc(expr.args[1])
        else
            Expr(expr.head, makepattern.(expr.args, (keys,), (slots,))...)
        end
    elseif expr in slots
        return splat ? makesegment(expr, keys) : makeslot(expr, keys)
    else
        # treat as a literal
        esc(expr)
    end
    return splat ? Expr(:..., res) : res
end

function makeconsequent(expr)
    if expr isa Expr
        if expr.head === :call
            if expr.args[1] === :(~)
                if expr.args[2] isa Symbol
                    return :(getindex(__MATCHES__, $(QuoteNode(expr.args[2]))))
                elseif expr.args[2] isa Expr && expr.args[2].args[1] == :(~)
                    @assert expr.args[2].args[2] isa Symbol
                    return :(getindex(__MATCHES__, $(QuoteNode(expr.args[2].args[2]))))
                end
            else
                return Expr(:call, map(makeconsequent, expr.args)...)
            end
        else
            return Expr(expr.head, map(makeconsequent, expr.args)...)
        end
    else
        # treat as a literal
        return esc(expr)
    end
end

abstract type AbstractRule end # Currently doesn't really do anything. Can be removed.

#-----------------------------
#### Regular Rewriting Rules

struct Rule{L, M, R} <: AbstractRule
    expr::Expr               # rule pattern stored for pretty printing
    lhs::L                   # the pattern
    matcher::M               # matcher(lhs)
    rhs::R                   # consequent
    depth::Int               # number of levels of expr this rule touches
end

getdepth(r::Rule) = r.depth

function rule_depth(rule, d=0, maxdepth=0)
    if rule isa Term
        maxdepth = reduce(max, (rule_depth(r, d+1, maxdepth) for r in arguments(rule)), init=1)
    elseif rule isa Slot || rule isa Segment
        maxdepth = max(d, maxdepth)
    end
    return maxdepth
end

function Base.show(io::IO, r::Rule)
    Base.print(io, r.expr)
end

const EMPTY_DICT = ImmutableDict{Symbol, Any}(:____, nothing)

function (r::Rule)(term)
    rhs = r.rhs

    try
        # n == 1 means that exactly one term of the input (term,) was matched
        success(bindings, n) = n == 1 ? (@timer "RHS" rhs(bindings)) : nothing
        return r.matcher(success, (term,), EMPTY_DICT)
    catch err
        throw(RuleRewriteError(r, term))
    end
end

"""
    rewrite_rhs(expr::Expr)

Rewrite the `expr` by dealing with `:where` if necessary.
The `:where` is rewritten from, for example, `~x where f(~x)` to `f(~x) ? ~x : nothing`.
"""
function rewrite_rhs(expr::Expr)
    if expr.head == :where
        rhs = expr.args[1]
        predicate = expr.args[2]
        expr = :($predicate ? $rhs : nothing)
    end
    return expr
end
rewrite_rhs(expr) = expr

function addslots(expr, slots)
    if expr isa Expr
        if expr.head === :macrocall && expr.args[1] in [Symbol("@rule"), Symbol("@capture"), Symbol("@slots")]
            Expr(:macrocall, expr.args[1:2]..., slots..., expr.args[3:end]...)
        else
            Expr(expr.head, addslots.(expr.args, (slots,))...)
        end
    else
        expr
    end
end

"""
    @slots [SLOTS...] ex

Declare SLOTS as slot variables for all `@rule` or `@capture` invocations in the expression `ex`.

_Example:_

```julia
julia> @slots x y z a b c SymbolicUtils.Chain([
    (@rule x^2 + 2x*y + y^2 => (x + y)^2),
    (@rule x^a * y^b => (x*y)^a * y^(b-a)),
    (@rule +(x...) => sum(x)),
])
SymbolicUtils.Chain(SymbolicUtils.Rule{SymbolicUtils.Term{Any, Nothing}}[...]))
```

See also: [`@rule`](@ref), [`@capture`](@ref)
"""
macro slots(args...)
    length(args) >= 1 || ArgumentError("@slots requires at least one argument")
    slots = args[1:end-1]
    expr = args[end]

    return esc(addslots(expr, slots))
end

"""
    @rule [SLOTS...] LHS => RHS

Creates a `Rule` object. A rule object is callable, and takes an expression and
rewrites it if it matches the LHS pattern to the RHS expression, returns
`nothing` otherwise. The rule language is described below.

LHS can be any possibly nested function call expression where any of the arguments can
optionally be a Slot (`~x`) or a Segment (`~x...`) (described below).

SLOTS is an optional list of symbols to be interpted as slots or segments
directly (without using `~`).  To declare slots for several rules at once, see
the `@slots` macro.

If an expression matches LHS entirely, then it is rewritten to the result of the
expression RHS, whose local scope includes the slot matches as variables.

**Slot**:

A Slot variable is written as `~x` and matches a single expression. `x` is the name of the variable. If a slot appears more than once in an LHS expression then expression matched at every such location must be equal (as shown by `isequal`).

_Example:_

Simple rule to turn any `sin` into `cos`:

```julia
julia> @syms a b c
(a, b, c)

julia> r = @rule sin(~x) => cos(x)
sin(~x) => cos(x)

julia> r(sin(1+a))
cos((1 + a))
```

A rule with 2 segment variables

```julia
julia> r = @rule sin(~x + ~y) => sin(x)*cos(y) + cos(x)*sin(y)
sin(~x + ~y) => sin(x) * cos(y) + cos(x) * sin(y)

julia> r(sin(a + b))
cos(a)*sin(b) + sin(a)*cos(b)
```

A rule that matches two of the same expressions:

```julia
julia> r = @rule sin(~x)^2 + cos(~x)^2 => 1
sin(~x) ^ 2 + cos(~x) ^ 2 => 1

julia> r(sin(2a)^2 + cos(2a)^2)
1

julia> r(sin(2a)^2 + cos(a)^2)
# nothing
```

A rule without `~`

```julia
julia> r = @slots x y z @rule x(y + z) => xy + xz
x(y + z) => xy + xz
```

**Segment**:

A Segment variable matches zero or more expressions in the function call.
Segments may be written by splatting slot variables (`~x...`).

_Example:_

This implements the distributive property of multiplication: `+(~ys...)` matches expressions
like `a + b`, `a+b+c` and so on. On the RHS `ys` presents as any old julia array.

```julia
julia> r = @rule ~x * +((~ys...)) => sum(map(y-> x * y, ys));

julia> r(2 * (a+b+c))
((2 * a) + (2 * b) + (2 * c))
```

A segment without `~`.

```julia
julia> r = @slots xs @rule min(xs...) => foldl(min, xs, Inf);

julia> r(min(a, b, c))
min(min(a, b), c)

**Predicates**:

There are two kinds of predicates, namely over slot variables and over the whole rule.
For the former, predicates can be used on both `~x` and `~x...` like `~x::f` or `~x::f...`.
Here `f` can be any julia function. In the case of a slot the function gets a single
matched subexpression, in the case of segment, it gets an array of matched expressions.

The predicate should return `true` if the current match is acceptable, and `false`
otherwise.

```julia
julia> two_πs(x::Number) = abs(round(x/(2π)) - x/(2π)) < 10^-9
two_πs (generic function with 1 method)

julia> two_πs(x) = false
two_πs (generic function with 2 methods)

julia> r = @slots x y z @rule sin(+(x..., y::two_πs, z...)) => sin(+(x..., z...))
sin(+(x..., y::two_πs, z...)) => sin(+(x..., z...))

julia> r(sin(a+3π))

julia> r(sin(a+6π))
sin(a)

julia> r(sin(a+6π+c))
sin((a + c))
```

For the predicate over the whole rule, use `@rule <LHS> => <RHS> where <predicate>`:

```
julia> @syms a b;

julia> predicate(x) = x === a;

julia> r = @rule ~x => x where f(x);

julia> r(a)
a

julia> r(b) === nothing
true
```

Note that this is syntactic sugar and that it is the same as something like
`@rule ~x => f(x) ? x : nothing`.

**Context**:

_In predicates_: Contextual predicates are functions wrapped in the `Contextual` type.
The function is called with 2 arguments: the expression and a context object
passed during a call to the Rule object (maybe done by passing a context to `simplify` or
a `RuleSet` object).

The function can use the inputs however it wants, and must return a boolean indicating
whether the predicate holds or not.

_In the consequent pattern_: Use `(@ctx)` to access the context object on the right hand side
of an expression.

**Compatibility**:
Segment variables may still be written as (`~~x`), and slot (`~x`) and segment (`~x...` or `~~x`) syntaxes on the RHS will still substitute the result of the matches.

See also: [`@capture`](@ref), [`@slots`](@ref)
"""
macro rule(args...)
    length(args) >= 1 || ArgumentError("@rule requires at least one argument")
    slots = args[1:end-1]
    expr = args[end]

    @assert expr.head == :call && expr.args[1] == :(=>)
    lhs = expr.args[2]
    rhs = rewrite_rhs(expr.args[3])
    keys = Symbol[]
    lhs_term = makepattern(lhs, keys, slots)
    unique!(keys)
    bind = map(key-> :($(esc(key)) = getindex(__MATCHES__, $(QuoteNode(key)))), keys)
    quote
        $(__source__)
        lhs_pattern = $(lhs_term)
        Rule($(QuoteNode(expr)),
             lhs_pattern,
             matcher(lhs_pattern),
             __MATCHES__ -> let $(bind...); $(makeconsequent(rhs)) end,
             rule_depth($lhs_term))
    end
end

"""
    @capture [SLOTS...] ex pattern

Uses a `Rule` object to capture an expression if it matches the `pattern`. Returns `true` and injects
slot variable match results into the calling scope when the `pattern` matches, otherwise returns false. The
rule language for specifying the `pattern` is the same in @capture as it is in `@rule`. Contextual matching
is not yet supported.

```julia
julia> @syms a; ex = a^a;

julia> if @capture ex (~x)^(~x)
           @show x
       elseif @capture ex 2(~y)
           @show y
       end;
x = a
```

See also: [`@rule`](@ref), [`@slots`](@ref)
"""
macro capture(args...)
    length(args) >= 2 || ArgumentError("@capture requires at least two arguments")
    slots = args[1:end-2]
    ex = args[end-1]
    lhs = args[end]

    keys = Symbol[]
    lhs_term = makepattern(lhs, keys, slots)
    unique!(keys)
    bind = Expr(:block, map(key-> :($(esc(key)) = getindex(__MATCHES__, $(QuoteNode(key)))), keys)...)
    quote
        $(__source__)
        lhs_pattern = $(lhs_term)
        __MATCHES__ = Rule($(QuoteNode(lhs)),
             lhs_pattern,
             matcher(lhs_pattern),
             identity,
             rule_depth($lhs_term))($(esc(ex)))
        if __MATCHES__ !== nothing
            $bind
            true
        else
            false
        end
    end
end

#-----------------------------
#### Associative Commutative Rules

struct ACRule{F,R} <: AbstractRule
    sets::F
    rule::R
    arity::Int
end

Rule(acr::ACRule)   = acr.rule
getdepth(r::ACRule) = getdepth(r.rule)

macro acrule(expr)
    arity = length(expr.args[2].args[2:end])
    quote
        ACRule(permutations, $(esc(:(@rule($(expr))))), $arity)
    end
end

macro ordered_acrule(expr)
    arity = length(expr.args[2].args[2:end])
    quote
        ACRule(combinations, $(esc(:(@rule($(expr))))), $arity)
    end
end

Base.show(io::IO, acr::ACRule) = print(io, "ACRule(", acr.rule, ")")

function (acr::ACRule)(term)
    r = Rule(acr)
    if !istree(term)
        r(term)
    else
        f =  operation(term)
        # Assume that the matcher was formed by closing over a term
        if f != operation(r.lhs) # Maybe offer a fallback if m.term errors.
            return nothing
        end

        T = symtype(term)
        args = unsorted_arguments(term)

        itr = acr.sets(eachindex(args), acr.arity)

        for inds in itr
            result = r(Term{T}(f, @views args[inds]))
            if result !== nothing
                # Assumption: inds are unique
                length(args) == length(inds) && return result
                return similarterm(term, f, [result, (args[i] for i in eachindex(args) if i ∉ inds)...], symtype(term))
            end
        end
    end
end

struct RuleRewriteError
    rule
    expr
end

getdepth(::Any) = typemax(Int)

@noinline function Base.showerror(io::IO, err::RuleRewriteError)
    msg = "Failed to apply rule $(err.rule) on expression "
    msg *= sprint(io->showraw(io, err.expr))
    print(io, msg)
end

function timerewrite(f)
    if !TIMER_OUTPUTS
        error("timerewrite must be called after enabling " *
              "TIMER_OUTPUTS in the main file of this package")
    end
    reset_timer!()
    being_timed[] = true
    x = f()
    being_timed[] = false
    print_timer()
    println()
    x
end

"""
    @timerewrite expr

If `expr` calls `simplify` or a `RuleSet` object, track the amount of time
it spent on applying each rule and pretty print the timing.

This uses [TimerOutputs.jl](https://github.com/KristofferC/TimerOutputs.jl).

## Example:

```julia

julia> expr = foldr(*, rand([a,b,c,d], 100))
(a ^ 26) * (b ^ 30) * (c ^ 16) * (d ^ 28)

julia> @timerewrite simplify(expr)
 ────────────────────────────────────────────────────────────────────────────────────────────────
                                                         Time                   Allocations
                                                 ──────────────────────   ───────────────────────
                Tot / % measured:                     340ms / 15.3%           92.2MiB / 10.8%

 Section                                 ncalls     time   %tot     avg     alloc   %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────────────────
 ACRule((~y) ^ ~n * ~y => (~y) ^ (~n ...    667   11.1ms  21.3%  16.7μs   2.66MiB  26.8%  4.08KiB
   RHS                                       92    277μs  0.53%  3.01μs   14.4KiB  0.14%     160B
 ACRule((~x) ^ ~n * (~x) ^ ~m => (~x)...    575   7.63ms  14.6%  13.3μs   1.83MiB  18.4%  3.26KiB
 (*)(~(~(x::!issortedₑ))) => sort_arg...    831   6.31ms  12.1%  7.59μs    738KiB  7.26%     910B
   RHS                                      164   3.03ms  5.81%  18.5μs    250KiB  2.46%  1.52KiB
   ...
   ...
 ────────────────────────────────────────────────────────────────────────────────────────────────
(a ^ 26) * (b ^ 30) * (c ^ 16) * (d ^ 28)
```
"""
macro timerewrite(expr)
    :(timerewrite(()->$(esc(expr))))
end

Base.@deprecate RuleSet(x) Postwalk(Chain(x))


text:
# compare two jld4 files

using HDF5,JLD,KUparser
a=load(ARGS[1])
b=load(ARGS[2])
@show a["deprel"]==b["deprel"]
@show a["postag"]==b["postag"]
aform=map(x->x.form, a["corpus"]); bform=map(x->x.form, b["corpus"])
@show aform==bform
apostag=map(x->x.postag, a["corpus"]); bpostag=map(x->x.postag, b["corpus"])
@show apostag==bpostag
ahead=map(x->x.head, a["corpus"]); bhead=map(x->x.head, b["corpus"])
@show ahead==bhead
adeprel=map(x->x.deprel, a["corpus"]); bdeprel=map(x->x.deprel, b["corpus"])
@show adeprel==bdeprel
awvec=(isdefined(a["corpus"][1],:wvec) ? map(x->x.wvec, a["corpus"]) : Any[])
bwvec=(isdefined(b["corpus"][1],:wvec) ? map(x->x.wvec, b["corpus"]) : Any[])
@show awvec==bwvec

ext: jl


text:
@testset "Wavefront Operations" begin
    function reduce_kernel(op,X,Y)
        idx = AMDGPU.workitemIdx().x
        Y[1] = AMDGPU.wfred(op,X[idx])
        nothing
    end
    function scan_kernel(op,X,Y)
        idx = AMDGPU.workitemIdx().x
        Y[idx] = AMDGPU.wfscan(op,X[idx],true)
        nothing
    end
    function bool_kernel(X,Y)
        idx = AMDGPU.workitemIdx().x
        Y[1] = AMDGPU.wfany(X[idx])
        Y[2] = AMDGPU.wfall(X[idx])
        Y[3] = AMDGPU.wfsame(X[idx])
        nothing
    end

    for T in (Cint, Clong, Cuint, Culong)
        X = rand(T(1):T(100), 64)
        for op in (Base.:+, max, min, Base.:&, Base.:|, Base.:⊻)
            RX, RY = ROCArray(X), ROCArray(zeros(T,1))
            wait(@roc groupsize=64 reduce_kernel(op,RX,RY))
            @test Array(RY)[1] == reduce(op,X)

            RX, RY = ROCArray(X), ROCArray(zeros(T,64))
            wait(@roc groupsize=64 scan_kernel(op,RX,RY))
            @test Array(RY) == accumulate(op,X)
        end
    end
    for T in (Float32, Float64)
        X = rand(T, 64)
        for op in (Base.:+, max, min)
            RX, RY = ROCArray(X), ROCArray(zeros(T,1))
            wait(@roc groupsize=64 reduce_kernel(op,RX,RY))
            @test Array(RY)[1] ≈ reduce(op,X)

            RX, RY = ROCArray(X), ROCArray(zeros(T,64))
            wait(@roc groupsize=64 scan_kernel(op,RX,RY))
            @test Array(RY) ≈ accumulate(op,X)
        end
    end
    for X in (rand(Cint(0):Cint(1), 64),
              zeros(Cint, 64),
              ones(Cint, 64),
              )
        RX, RY = ROCArray(X), ROCArray(zeros(Bool,3))
        wait(@roc groupsize=64 bool_kernel(RX,RY))
        Y = Array(RY)
        @test_skip Y[1] == all(x->x==1,X)
        @test_skip Y[2] == any(x->x==1,X)
        @test_skip Y[3] == (length(unique(X)) == 1)
    end
end

@testset "Wavefront Information" begin
    function kernel(X)
        X[1] = wavefrontsize()
        nothing
    end
    RX = ROCArray(zeros(UInt32, 1))
    wait(@roc kernel(RX))
    @allowscalar @test RX[1] == 64
end


text:
#############################################################################
#############################################################################
#
# This file contains units tests for polynomial operations
#
#############################################################################
#############################################################################
"""
Test Polynomial Addition and Subtraction
"""

function mod_sum_test_poly(mod::Int ; num_sums::Int = 100, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:num_sums
        p1 = rand(PolynomialModP,mod)
        p2 = rand(PolynomialModP,mod)
        sum = p1+p2
        @assert iszero((p2+p1) - sum)
    end
    println("Addition Test - PASSED")
end
"""
Test product of polynomials.
"""
function mod_prod_test_poly(m::Int ;N::Int = 100, N_prods::Int = 20, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(PolynomialModP,m)
        p2 = rand(PolynomialModP,m)
        prod = p1*p2
        prodtest= mod(leading(p1)*leading(p2),m)
        if !iszero(prodtest - leading(prod))
            println("There is a multiplication problem for $p1 & $p2")
            @assert iszero(prodtest - leading(prod))
        end
    end

    for _ in 1:N
        p_base = PolynomialModP(Term(1,0),m)
        for _ in 1:N_prods
            p = rand(PolynomialModP,m)
            prod = p_base*p
            test = leading(p_base)*leading(p)
            test.coeff = (test.coeff)%m
            if !iszero(test - leading(prod))
                println("There is a multiplication problem for $p1 & $p2")
                @assert iszero(test - leading(prod))
            end
            p_base = prod
        end
    end
    println("Multiplication Test - PASSED")
end

function new_prod_test(p1,p2)
    println("")
    a = p1*p2 - p1*̄p1
    if a != 0
        println("The multiplication does not work for:")
        println(p1*p2)
        println(p1*̄p2)
        return
    end
    return "Multiplication Works"
end

"""
Test derivative of polynomials (as well as product).
"""
function mod_prod_derivative_test_poly(mod::Int;N::Int = 10^2,  seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(PolynomialModP,mod)
        p2 = rand(PolynomialModP,mod)
        p1d = derivative(p1)
        p2d = derivative(p2)
        @assert iszero((p1d*p2) + (p1*p2d)-derivative(p1*p2))
    end
    println("Derivative Test - PASSED")
end

"""
Test division of polynomials modulo p.
"""
function mod_division_test_poly(mod::Int; N::Int = 10^4, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(PolynomialModP,mod)
        p2 = rand(PolynomialModP,mod)
        p_prod = p1*p2
        q, r = PolynomialModP(mod), PolynomialModP(mod)
        try
            q, r = divide(p_prod, p2)
            if (q, r) == (nothing,nothing)
                println("Unlucky prime: $p1 is reduced to $(p1 % prime) modulo $prime")
                continue
            end
        catch e
            if typeof(e) == DivideError
                @assert p2, prime == 0
            else
                throw(e)
            end
        end
        @assert iszero( q*p2+r - p_prod )
    end
    println("Division Test - PASSED")
end

"""
Test the extended euclid algorithm for polynomials modulo p.
"""
function mod_ext_euclid_test_poly(mod::Int; N::Int = 100, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(PolynomialModP,mod)
        p2 = rand(PolynomialModP,mod)
        g, s, t = extended_euclid_alg(p1, p2)
        @assert s*p1 + t*p2 - g == 0
    end
    println("Extended Euclid Algorithm Test - PASSED")
end


text:
#=

Code for solving Laplacian and Diagnally Dominant Systems
by augmented spanning tree preconditioners.

Started by Dan Spielman

=#

#=========================================
  Augmented Spanning Tree Preconditioner
=========================================#

mutable struct AugTreeParams
    treeAlg::Function
    opt::Bool
    nnzL_fac::Float64
    flops::Float64
    verbose::Bool
end

AugTreeParams() = AugTreeParams(akpw, true, 4.0, 200.0, false)
AugTreeParamsOld() = AugTreeParams(akpw, false, 0, 0, false)

"""
    B = augmentTree{Tv,Ti}(tree, A, k)

Takes as input a tree and an adjacency matrix of a graph.
It then computes the stretch of every edge of the graph wrt
the tree.  It then adds back the k edges of highest stretch,
and k edges sampled according to stretch.

This is the old alg.  We now recommend using augmentTreeOpt.
"""
function augmentTree(tree::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}, k::Ti) where {Tv,Ti}

    st = comp_stretches(tree, A)

    # just to be safe, remove the tree from this
    #=
    (ti,tj) = findnz(tree)
    for i in 1:length(ti)   # most of the time
        st[ti[i],tj[i]] = 0
        st[tj[i],ti[i]] = 0
    end
    =#

    (ai,aj,av) = findnz(triu(st))

    ord = sortperm(av, rev=true)

    edgeinds = zeros(Bool,length(av))
    for i in 1:k
        edgeinds[ord[i]] = true
    end

    s = sum(av[ord[(k+1):end]])
    probs = av * k / s
    probs[ord[1:k]] .= 0
    edgeinds[rand(length(av)) .< probs] .= true

    augi = ai[edgeinds]
    augj = aj[edgeinds]
    augm = length(augi)
    augv = zeros(Tv, augm)
    for i in 1:augm,
        augv[i] = A[augi[i],augj[i]]
    end

    n = size(A)[1]
    aug = sparse(augi, augj, augv, n, n)
    aug = aug + aug'

    return tree + aug

end

"""
    B = augmentTreeOpt{Tv,Ti}(tree, A, params)

Takes as input a tree and an adjacency matrix of a graph.
It then computes the stretch of every edge of the graph wrt
the tree.  It uses cholmod to decide how many edge to add back,
shooting for nnzL_fac times n entries in the factored augmented tree,
with a number of flops to factor equal to nnz(a)*flops_fac.
The edges to add back are then choen at random.
"""
function augmentTreeOpt(tree::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}; params=AugTreeParams()) where {Tv,Ti}

    Aminus = A - tree

    ai,aj,av = findnz(triu(Aminus))

    n = A.n
    m = length(ai)

    st = comp_stretches(tree, Aminus)
    _,_,sv = findnz(triu(st))

    r = -log.(rand(m)) ./ sv
    ord = sortperm(r)

    nnzLTooBig(nnzL) = (nnzL-2*(n-1)) > n*params.nnzL_fac
    nnzLTooSmall(nnzL) = (nnzL-2*(n-1)) < n*params.nnzL_fac/2
    flopsTooBig(flops) = flops > (n+m)*params.flops
    flopsTooSmall(flops) = flops < (n+m)*params.flops/4

    k = Int(round(2*sqrt(n)))

    first = true
    direction = 0
    done = false

    while ~done

        edgeinds = ord[1:min(k,m)]
        augi = ai[edgeinds]
        augj = aj[edgeinds]
        augv = av[edgeinds]

        n = size(A,1)
        aug = sparse(augi, augj, augv, n, n)
        aug = aug + aug'

        augTree = tree+aug

        nnzL, flops = ask_cholmod(lap(augTree))

        if first
            first = false
            if ~(nnzLTooBig(nnzL) || nnzLTooSmall(nnzL) ||
                 flopsTooBig(flops) || flopsTooSmall(flops))
                done = true
            elseif nnzLTooBig(nnzL) || flopsTooBig(flops)
                k = div(k,2)
                direction = -1
            else
                k = k * 2
                direction = 1
                if k >= m
                    done = true
                end

            end
        else
            if direction == -1
                if nnzLTooBig(nnzL) || flopsTooBig(flops)
                    k = div(k,2)
                else
                    done = true
                end

            else # direction == 1

                if nnzLTooSmall(nnzL) && flopsTooSmall(flops)
                    k = k * 2
                    if k >= m
                        done = true
                    end
                else
                    done = true
                end
            end
        end

        if done
            if params.verbose
                println("Increased k by $(k/round(2*sqrt(n))) to $(k).")
            end

            return augTree
        end
    end
end

"""
    pre = augTreePrecon{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; params=AugTreeParams())

This is an augmented spanning tree preconditioner for diagonally dominant
linear systems.  It takes as optional input a tree growing algorithm.
It adds back 2sqrt(n) edges via augmentTree: the sqrt(n) of highest stretch
and another sqrt(n) sampled according to stretch.
For most purposes, one should directly call `augTreeSolver`.
"""
function augTreePrecon(ddmat::SparseMatrixCSC{Tv,Ti};  params=AugTreeParams()) where {Tv,Ti}

  adjmat = -triu(ddmat,1)
  adjmat = adjmat + adjmat'

  tree = params.treeAlg(adjmat)

  n = size(ddmat)[1]

  if params.opt
      augtree = augmentTreeOpt(tree,adjmat,params=params)
  else
      augtree = augmentTree(tree,adjmat,convert(Int,round(sqrt(n))))
  end

  Dx = sparse(Diagonal(ddmat*ones(n)))

  augDD = Dx + sparse(Diagonal(augtree*ones(n))) - augtree

  F = cholesky(augDD)

  return x -> (F\x)

end

"""
    solver = augTreeSddm(sddm; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[],  params=AugTreeParams())

An "augmented spanning tree" solver for positive definite diagonally dominant matrices.  It works by adding edges to a low stretch spanning tree.  It calls `augTreePrecon` to form the preconditioner.  `params` has entries

* `params.treeAlg` default to `akpw`
* `params.opt` if true, it interacts with cholmod to choose a good number of edges to add back.  If false, it adds back 2*sqrt(n).
"""
function augTreeSddm(sddm::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[],  params=AugTreeParams()) where {Tv,Ti}

    t1 = time()

    F = augTreePrecon(sddm; params=params)
    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_,pcgIts=pcgIts_) = pcg(sddm, b, F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)

    if verbose
        println("Solver build time: ", round((time() - t1),digits=3), " seconds.")
    end

    return f

end

"""
    pre = augTreeLapPrecon{Tv,Ti}(A; params=AugTreeParams())

This is an augmented spanning tree preconditioner for Laplacians.
It takes as optional input a tree growing algorithm.
It adds back 2sqrt(n) edges via `augmentTree`: the sqrt(n) of highest stretch
and another sqrt(n) sampled according to stretch.
For most purposes, one should directly call `augTreeLapSolver`."""
function augTreeLapPrecon(a::SparseMatrixCSC{Tv,Ti}; params=AugTreeParams()) where {Tv,Ti}

    tree = params.treeAlg(a)

    if params.opt
        augtree = augmentTreeOpt(tree,a, params=params)
    else
        augtree = augmentTree(tree,a,convert(Int,round(sqrt(a.n))))
    end

    F = (h -> h)

    try
        F = chol_lap(augtree)
    catch
        println("cholfact failed in augTreeLapPrecon.  Going to backup routine")
        F = augTreeFactor(augtree, tree)
    end

    return F

end

"""
    solver = augTreeLap(A; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params=AugTreeParams())

An "augmented spanning tree" solver for Laplacians.  It works by adding edges to a low stretch spanning tree.  It calls `augTreePrecon` to form the preconditioner.  `params` has entries

* `params.treeAlg` default to `akpw`
* `params.opt` if true, it interacts with cholmod to choose a good number of edges to add back.  If false, it adds back 2*sqrt(n).
"""
function augTreeLap(a::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params=AugTreeParams()) where {Tv,Ti}

    return lapWrapComponents(augTreeLap1, a, verbose=verbose, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, params=params)

end

function augTreeLap1(a::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params=AugTreeParams()) where {Tv,Ti}

  F = augTreeLapPrecon(a, params=params)

  tol_ =tol
  maxits_ =maxits
  maxtime_ =maxtime
  verbose_ =verbose
  pcgIts_ =pcgIts

  la = lap(a)

  f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) = pcg(la, b .- mean(b), F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose)

  return f

end

function augTreeFactor(a, tree; verbose=false,
                      tol::Real=1e-6, maxits::Integer=1000, maxtime=Inf, pcgIts=Int[])
    n = size(a,1);

    ord::Array{Int64,1} = Laplacians.dfsOrder(tree)

    # these lines could be MUCH faster
    aord = symPermuteCSC(a,ord)
    tord = symPermuteCSC(tree,ord)

    laord = lap(aord)

    Aminus = aord - tord

    (aminusi, aminusj, aminusv) = findnz(Aminus)

    marked = ones(Int64,n)
    marked[aminusi] .= 0
    marked[aminusj] .= 0

    if (nnz(a) == 2*(n - 1))
        if verbose
            println("The graph is a tree.")
        end

        # mark something so this code stops somewhere
        i = 1
        j = nbri(aord, i, 1)
        marked[i] = 0
        marked[j] = 0
    end

    elims1, elims2, ind::Array{Int64,1}, subtree = elimDeg12(tord, marked)

    map = zeros(Int64,n)

    n1 = length(ind)

    map[ind] = collect(1:n1)
    mapi = map[aminusi]
    mapj = map[aminusj]

    rest = sparse(mapi, mapj, aminusv, n1,n1)
    a1 = rest + subtree
    la1 = lap(a1)

    fsub = chol_lap(a1)

    f1 = function(b::Array{Float64,1})

        subMean!(b) # b = b - mean(b)

        y = forwardSolve(b, elims1, elims2)
        ys = y[ind]

        xs = fsub(ys)

        x = zeros(Float64,n)
        x[ind] = xs

        backSolve(x, y, elims1, elims2)
        subMean!(x) # x = x - mean(x)

        return x
    end

    f = function(b)

        bord = b[ord]

        xord = f1(bord)

        #xord = pcg(laord, bord, f1, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)

        x = zeros(Float64,n)
        x[ord] = xord .- mean(xord)

        return x
    end

    return f

end


text:
module Day1
export day1, doublePointerSearch

function readInput(inputfile="inputs/input01.txt")
    open(inputfile) do f
        lines = Int[];
        while (l=readline(f)) !== ""
            push!(lines, parse(Int, l))
        end
        sort!(lines)
    end
end

function doublePointerSearch(values, targetSum)
    front = 1
    back = length(values)
    sum = values[front] + values[back]
    while (sum !== targetSum) && (front < back)
        if sum < targetSum
            front += 1
        else
            back -= 1
        end
        sum = values[front] + values[back]
    end
    return (front, back) # caller has to handle the fail case where front == back
end

function day1(inputfile="inputs/input01.txt")
    values = readInput(inputfile)
    (day1Part1(values), day1Part2(values))
end

function day1Part1(values)
    front, back = doublePointerSearch(values, 2020)
    @assert front < back "No 2 terms sum to 2020!"
    values[front]*values[back]
end

function day1Part2(values)
    for index3 in length(values):-1:3
        index1, index2 = doublePointerSearch(view(values, 1:index3), 2020-values[index3])
        if index1 < index2 # search successful
            return values[index1]*values[index2]*values[index3]
        end
    end
    error("No triplet summing to 2020 found!")
end

end #module


text:
jtype(tag::BSONType)::DataType =
  tag == null ? Nothing :
  tag == boolean ? Bool :
  tag == int32 ? Int32 :
  tag == int64 ? Int64 :
  tag == double ? Float64 :
  error("Unsupported tag $tag")

parse_cstr(io::IOT) where {IOT <: IO} =
  Base.readuntil_string(io, 0x0, false)

parse_cstr_unsafe(io::IO)::Vector{UInt8} = readuntil(io, 0x0, keep=false)

# Not really unsafe, but we do access internal fields directly. Avoids
# allocating a temporary string/vector.
function parse_cstr_unsafe(io::IOBuffer)::SubArray{UInt8, 1}
  st = io.ptr

  while read(io, UInt8) ≠ 0x0
  end

  view(io.data, st:(io.ptr - 2))
end

function parse_tag(io::IOT, tag::BSONType) where {IOT <: IO}
  if tag == null
    nothing
  elseif tag == document
    parse_doc(io)
  elseif tag == array
    parse_array(io)
  elseif tag == string
    len = read(io, Int32)-1
    s = String(read(io, len))
    eof = read(io, 1)
    s
  elseif tag == binary
    len = read(io, Int32)
    subtype = read(io, 1)
    read(io, len)
  else
    read(io, jtype(tag))
  end
end

function parse_array(io::IOT)::BSONArray where {IOT <: IO}
  len = read(io, Int32)
  ps = BSONArray()

  while (tag = read(io, BSONType)) ≠ eof
    # Note that arrays are dicts with the index as the key
    while read(io, UInt8) != 0x00
      nothing
    end
    push!(ps, parse_tag(io::IOT, tag))
  end

  ps
end

const SEEN_REF = 1
const SEEN_DATA = 1 << 1
const SEEN_TYPE = 1 << 2
const SEEN_TYPENAME = 1 << 3
const SEEN_TAG = 1 << 4
const SEEN_NAME = 1 << 5
const SEEN_PARAMS = 1 << 6
const SEEN_PATH = 1 << 7
const SEEN_OTHER = 1 << 8
const SEEN_SIZE = 1 << 9
const SEEN_VAR = 1 << 10
const SEEN_BODY = 1 << 11

const SEEN_TAG_STRUCT = 1 << 12
const SEEN_TAG_BACKREF = 1 << 13
const SEEN_TAG_DATATYPE = 1 << 14
const SEEN_TAG_SYMBOL = 1 << 15
const SEEN_TAG_TUPLE = 1 << 16
const SEEN_TAG_SVEC = 1 << 17
const SEEN_TAG_UNION = 1 << 18
const SEEN_TAG_ANON = 1 << 19
const SEEN_TAG_REF = 1 << 20
const SEEN_TAG_ARRAY = 1 << 21
const SEEN_TAG_UNIONALL = 1 << 22

function classify_doc_tag(tag::AbstractVector{UInt8})::Union{Int64, String}
  if tag == b"backref"
    SEEN_TAG_BACKREF
  elseif tag == b"struct"
    SEEN_TAG_STRUCT
  elseif tag == b"datatype"
    SEEN_TAG_DATATYPE
  elseif tag == b"symbol"
    SEEN_TAG_SYMBOL
  elseif tag == b"tuple"
    SEEN_TAG_TUPLE
  elseif tag == b"svec"
    SEEN_TAG_SVEC
  elseif tag == b"jl_bottom_type"
    SEEN_TAG_UNION
  elseif tag == b"jl_anonymous"
    SEEN_TAG_ANON
  elseif tag == b"ref"
    SEEN_TAG_REF
  elseif tag == b"array"
    SEEN_TAG_ARRAY
  elseif tag == b"unionall"
    SEEN_TAG_UNIONALL
  else
    String(tag)
  end
end

function parse_doc_tag(io::IO)::Union{Int64, String}
  len = read(io, Int32) - 1
  tag = read(io, len)
  eof = read(io, 1)

  classify_doc_tag(tag)
end

function parse_doc_tag(io::IOBuffer)::Union{Int64, String}
  len = read(io, Int32) - 1
  spos = position(io)
  tag = parse_cstr_unsafe(io)

  if length(tag) == len
    classify_doc_tag(tag)
  else
    seek(io, spos)
    s = String(read(io, len))
    eof = read(io, 1)
    s
  end
end

function parse_doc(io::IOT) where {IOT <: IO}
  #@debug "parse_doc"
  len = read(io, Int32)

  seen::Int64 = 0
  see(it::Int64) = seen = seen | it
  saw(it::Int64)::Bool = seen & it != 0
  only_saw(it::Int64)::Bool = seen == it

  # First try to parse this document as a Tagged* intermediate type. Note that
  # both nothing and missing are valid data values.
  local tref, tdata, ttype, ttypename, ttag, tname, tparams, tpath, tsize, tvar, tbody
  local other
  local k::AbstractVector{UInt8}
  backrefs = nothing

  for _ in 1:6
    if (tag = read(io, BSONType)) == eof
      break
    end
    k = parse_cstr_unsafe(io)
    #@debug "Read key" String(k)

    if k == b"ref"
      see(SEEN_REF)
      tref = parse_tag(io, tag)
    elseif k == b"data"
      see(SEEN_DATA)
      tdata = parse_tag(io, tag)
      #@debug "Read" tdata
    elseif k == b"type"
      see(SEEN_TYPE)
      ttype = parse_tag(io, tag)
      #@debug "Read" ttype
    elseif k == b"typename"
      see(SEEN_TYPENAME)
      ttypename = parse_tag(io, tag)
      #@debug "Read" ttype
    elseif k == b"tag"
      see(SEEN_TAG)
      if tag == string
        if (dtag = parse_doc_tag(io)) isa Int64
          see(dtag)
        else
          ttag = dtag
          #@debug "Read" dtag
        end
      else
        ttag = parse_tag(io, tag)
        #@debug "Read" ttag
      end
    elseif k == b"name"
      see(SEEN_NAME)
      tname = parse_tag(io, tag)
      #@debug "Read" tname
    elseif k == b"params"
      see(SEEN_PARAMS)
      tparams = parse_tag(io, tag)
      #@debug "Read" tparams
    elseif k == b"path"
      see(SEEN_PATH)
      tpath = parse_tag(io, tag)
    elseif k == b"size"
      see(SEEN_SIZE)
      tsize = parse_tag(io, tag)
    elseif k == b"var"
      see(SEEN_VAR)
      tvar = parse_tag(io, tag)
    elseif k == b"body"
      see(SEEN_BODY)
      tbody = parse_tag(io, tag)
    elseif k == b"_backrefs"
      backrefs = parse_tag(io, tag)
    else
      see(SEEN_OTHER)
      other = parse_tag(io, tag)
      #@debug "Read" other
      break
    end
  end

  ret = if only_saw(SEEN_TAG | SEEN_REF | SEEN_TAG_BACKREF)
    TaggedBackref(tref)
  elseif only_saw(SEEN_TAG | SEEN_TYPE | SEEN_DATA | SEEN_TAG_STRUCT)
    TaggedStruct(ttype, tdata)
  elseif only_saw(SEEN_TAG | SEEN_NAME | SEEN_PARAMS | SEEN_TAG_DATATYPE)
    TaggedType(tname, tparams)
  elseif only_saw(SEEN_TAG | SEEN_NAME | SEEN_TAG_SYMBOL)
    Symbol(tname)
  elseif only_saw(SEEN_TAG | SEEN_DATA | SEEN_TAG_TUPLE)
    TaggedTuple(tdata)
  elseif only_saw(SEEN_TAG | SEEN_DATA | SEEN_TAG_SVEC)
    TaggedSvec(tdata)
  elseif only_saw(SEEN_TAG | SEEN_TAG_UNION)
    Union{}
  elseif only_saw(SEEN_TAG | SEEN_TYPENAME | SEEN_PARAMS | SEEN_TAG_ANON)
    TaggedAnonymous(ttypename, tparams)
  elseif only_saw(SEEN_TAG | SEEN_PATH | SEEN_TAG_REF)
    TaggedRef(tpath)
  elseif only_saw(SEEN_TAG | SEEN_TYPE | SEEN_SIZE | SEEN_DATA | SEEN_TAG_ARRAY)
    TaggedArray(ttype, tsize, tdata)
  elseif only_saw(SEEN_TAG | SEEN_VAR | SEEN_BODY | SEEN_TAG_UNIONALL)
    TaggedUnionall(tvar, tbody)
  else
    # It doesn't look like a Tagged*, so just allocate a Dict
    dic = BSONDict()
    saw(SEEN_REF) && (dic[:ref] = tref)
    saw(SEEN_DATA) && (dic[:data] = tdata)
    saw(SEEN_TYPE) && (dic[:type] = ttype)
    saw(SEEN_TYPENAME) && (dic[:typename] = ttypename)
    saw(SEEN_TAG)  && (dic[:tag] = ttag)
    saw(SEEN_NAME) && (dic[:name] = tname)
    saw(SEEN_PARAMS) && (dic[:params] = tparams)
    saw(SEEN_PATH) && (dic[:path] = tpath)
    saw(SEEN_SIZE) && (dic[:size] = tsize)
    saw(SEEN_VAR) && (dic[:var] = tvar)
    saw(SEEN_BODY) && (dic[:body] = tbody)
    saw(SEEN_OTHER) && (dic[Symbol(String(k))] = other)

    if tag != eof
      while (tag = read(io, BSONType)) ≠ eof
        local k = Symbol(parse_cstr(io))
        #@debug "Read key" k
        dic[k] = parse_tag(io::IOT, tag)
      end
    end

    if haskey(dic, :_backrefs)
      backrefs = dic[:_backrefs]
      delete!(dic, :_backrefs)
    end

    dic
  end

  if backrefs != nothing
    BackRefsWrapper(ret, backrefs)
  else
    ret
  end
end

backrefs!(x, refs) = applychildren!(x -> backrefs!(x, refs), x)

backrefs!(bref::TaggedBackref, refs) = refs[bref.ref]

backrefs!(x) = x
function backrefs!(dict::BackRefsWrapper)
  backrefs!(dict.refs, dict.refs)
  backrefs!(dict.root, dict.refs)
end

parse(io::IOT) where {IOT <: IO} = backrefs!(parse_doc(io))
parse(path::String) = open(parse, path)

function load_compat(x)
  @warn "load_compat() is deprecated and may be removed in future versions!"
  raise_recursive(parse(x), IdDict{Any, Any}())
end

function roundtrip(x, flags...)
  buf = IOBuffer()
  bson(buf, Dict(:stuff => x))
  seek(buf, 0)
  if :set_root_type in flags
    doc = Document(buf, typeof(x))
    doc[:stuff]
  else
    (:compat in flags ? load_compat(buf) : load(buf))[:stuff]
  end
end

function halftrip(x)
  buf = IOBuffer()
  bson(buf, Dict(:stuff => x))
  parse_doc(seek(buf, 0))
end


text:
function conv(signal1::Vector{Complex}, signal2::Vector{Complex})
    n = length(signal1) + length(signal2)
    out = Vector{Complex}(n)
    sum = 0

    for i = 0:n
        for j = 0:i
            if(j < length(signal1))
                sum += signal1[j] * signal2[i-j]
            end
        end
        out[i] = sum
        sum = 0
    end

    return out
end

function conv_fft(signal1::Vector{Complex}, signal2::Vector{Complex})
    return ifft(fft(signal1).*fft(signal2))
end

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
# test_utils.jl : utilities for testing SimpleMultigrid

# log resnorm
function log(mg,d)
    println(string("+","-"^30,"+"))
	println(string("|"," "^11,@sprintf("n = %4i",(size(mg.grids[1].A,1))^(1/d)+1)," "^11,"|"))
    println(string("+","-"^9,"+","-"^20,"+"))
    println(string("| V-cycle |  |rʰ|"," "^7,"ratio","  |"))
    println(string("+","-"^9,"+","-"^20,"+"))
    println(string(@sprintf("|   %4i  |",0),@sprintf("  %6.3e",mg.resnorm[1])," "^8," |"))
    for i in 2:length(mg.resnorm)
        println(string(@sprintf("|   %4i  |",i-1),@sprintf("  %6.3e",mg.resnorm[i]),@sprintf("  %3.2f",mg.resnorm[i]/mg.resnorm[i-1]),"   |"))
    end
    println(string("+","-"^9,"+","-"^20,"+"))
end

# plot
@require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
    function mplot(x...)
        n = maximum(length.(x))+1
        figure(figsize=(5,5))
        for i in 1:length(x)
            factor = round(Int, n/(length(x[i])+1))
            plot(factor:factor:n-1,x[i],"-o")
        end
    end

    function msurf(s...)
        n = maximum(size.(s,1))+1
        m = maximum(size.(s,2))+1
        figure(figsize=(5,5))
        alphas = linspace(0.2,1,length(s))
        for i in 1:length(s)
            nfactor = round(Int, n/(size(s[i],1)+1))
            mfactor = round(Int, m/(size(s[i],2)+1))
            x = nfactor:nfactor:n-1
            y = mfactor:mfactor:m-1
            xgrid = repmat(x,1,length(y))
            ygrid = repmat(y',length(x),1)
            surf(xgrid,ygrid,s[i],alpha=alphas[i])
        end
    end
end


text:
"""
Resorption. Parameters for reabsorbtion of nutrients from structures when metabolic rates fall.
"""
abstract type AbstractResorption end

K_resorption(p::AbstractResorption) = p.K_resorption

resorption!(o, u) = resorption!(resorption_pars(o), o, u)
resorption!(p::Nothing, o, u) = nothing

"""
    StructuralLossResorption(K_resorption)

Structure is lost while reserves are retained.

$(FIELDDOCTABLE)
"""
@columns struct StructuralLossResorption{R} <: AbstractResorption
    K_resorption::R | 0.000001 | _    | (1e-8, 1e-3) | true | "Half saturation metabolic rate for resorption of tissues."
end

resorption!(p::StructuralLossResorption, o, u) = begin
    o.J[:V,:res] -= resorption(p, o, u)
    nothing
end

"""
    LosslessResorption(K_resorption)

Structure is distributed back to C an N reserves without loss.

$(FIELDDOCTABLE)
"""
@columns struct LosslessResorption{R} <: AbstractResorption
    K_resorption::R | 0.000001 | _    | (1e-8, 1e-3) | true | "Half saturation metabolic rate for resorption of tissues."
end

resorption!(p::LosslessResorption, o, u) = begin
    aph = resorption(p, o, u)
    o.J[:C,:res] += aph
    o.J[:N,:res] += aph * n_N_V(o)
    o.J[:V,:res] -= aph
    nothing
end

"""
    DissipativeResorption(r_EN_V, r_EN_V, K_resorption)

Some structure is distributed back to C an N reserves with a proportion lost to the environment.

This model has parameters for controlling differential fractions of C an N resorption.

$(FIELDDOCTABLE)
"""
@columns struct DissipativeResorption{R,P} <: AbstractResorption
    # Field         | Default  | Unit       | Bounds       | Log  | Description
    K_resorption::R | 0.000001 | _          | (1e-8, 1e-3) | true | "Half saturation metabolic rate for resorption of tissues."
    r_EC_V::P       | 0.0      | mol*mol^-1 | (0.0, 1.0)   | _    | "Proportion of C recovered from structure"
    r_EN_V::P       | 0.5      | mol*mol^-1 | (0.0, 1.0)   | _    | "Proportion of N recovered from structure"
end

resorption!(p::DissipativeResorption, o, u) = begin
    aph = resorption(m, o, u)
    o.J[:C,:res] += aph * p.r_EC_V
    o.J[:N,:res] += aph * p.r_EN_V * n_N_V(o)
    o.J[:V,:res] -= aph
    nothing
end

"""
    resorption(p::AbstractResorption, o::Organ, u)

Resoption of structure `V` with shape-adjusted metabolic rate
"""
@inline resorption(p::AbstractResorption, o, u) =
    resorption(u[:V], rate(o) / scaling(o), K_resorption(p))

@inline resorption(v::Number, r::Number, a::Number) =
    v * (oneunit(r) - 1 / (oneunit(1 / r) + a / r))