text:
# ./src/lif_conductance.jl

############

#=
coductance-based models, generalized so that they can have
different sypatic kernels (Exp or Expdiff)
and different spike generating functions (None, or Exp)
=#
# abstract type SynapticKernel end moved to main
struct SKExp <: SynapticKernel
  τ::Float64
end

# Input update for conductance based IF neurons, assuming the kernel takes
# at most two traces arguments (as in diff of exponentials)
@inline function postsynaptic_kernel_update!(input::Vector{R},post_state::Vector{R},
    post_trace1::Vector{R},post_trace2::Vector{R},pre_synkernel::SynapticKernel,
    pre_v_reversal::R,idx::Integer) where R
  ker_term = pre_synkernel(post_trace1[idx],post_trace2[idx])
	input[idx] += ker_term*(pre_v_reversal - post_state[idx])
  return nothing
end

# # this function is not needed
# @inline function (sk::SKExp)(x::Real)
#   return x/sk.τ
# end
@inline function (sk::SKExp)(x::Real,::Real)
  return x/sk.τ
end
@inline function trace_spike_update!(w::R,trace1::Vector{R},::Vector{R},
    ::SKExp,idx::Integer) where R
  trace1[idx] += w
end

# difference of exponentials, needs two traces
struct SKExpDiff <: SynapticKernel
  τ_plus::Float64
  τ_minus::Float64
end
@inline function (sk::SKExpDiff)(x_plus::Real,x_minus::Real)
  return (x_plus-x_minus) / (sk.τ_plus-sk.τ_minus)
end
@inline function trace_spike_update!(w::R,trace1::Vector{R},trace2::Vector{R},
    ::SKExpDiff,idx::Integer) where R
  trace1[idx] += w
  trace2[idx] += w
end

# abstract type SpikeGenFunction end # moved to main file
struct SpikeGenNone <: SpikeGenFunction end
@inline function (sk::SpikeGenNone)(::R) where R
  return zero(R)
end

struct SpikeGenEIF <: SpikeGenFunction
  vth_exp::Float64
  s::Float64
end
@inline function (sk::SpikeGenEIF)(v::R) where R<:Real
  return sk.s*exp((v-sk.vth_exp)/sk.s)
end

# abstract type NTConductance <:NeuronType end moved to main file
# ALL conductance neurons must have spikes (neuron.isfiring::BitArray{1})
# a synaptic kernel (neuron.synaptic_kernel::SynapticKernel)
# and a reversal potential (neuron.v_reversal::Float64)

struct NTLIFConductance{SK<:SynapticKernel} <: NTConductance
  synaptic_kernel::SK
  spikegen::SpikeGenFunction
  τ_post::Float64 # time constant for postsynaptic update
  capacitance::Float64 # membrane capacitance
  v_threshold::Float64 # max voltage before reset
  v_reset::Float64 # reset after spike
  v_leak::Float64 # reversal potential for leak term
  τ_refractory::Float64 # refractory time
  v_reversal::Float64 # reversal potential that affects postsynaptic neurons
end

struct PSLIFConductance{NT} <: PSSpikingType{NT}
  neurontype::NT
  n::Int64 # pop size
  state_now::Vector{Float64}
  input::Vector{Float64}
	alloc_dv::Vector{Float64}
	last_fired::Vector{Float64}
	isfiring::BitArray{1}
	isrefractory::BitArray{1}
	pre_reverse_potentials::Vector{Float64}
	pre_conductances_now::Vector{Float64}
end
function PSLIFConductance(p::NTLIFConductance,n)
  zz = _ -> zeros(Float64,n)
  ff = _ ->  falses(n)
  PSLIFConductance(p,n,ntuple(zz,4)...,ntuple(ff,2)...,ntuple(zz,2)...)
end

function reset!(ps::PSLIFConductance)
  fill!(ps.state_now,0.0)
  fill!(ps.input,0.0)
  fill!(ps.last_fired,-1E6)
  fill!(ps.isfiring,false)
  fill!(ps.isrefractory,false)
  fill!(ps.pre_reverse_potentials,0.0)
  fill!(ps.pre_conductances_now,0.0)
  return nothing
end

# connection
# use ConnGeneralIF2

@inline function trace_spike_update!(conn::ConnGeneralIF2,
    w::Real,pre_sker::SynapticKernel,post_idx::Integer)
  # not normalized by taus
  return trace_spike_update!(w,conn.post_trace1,conn.post_trace2,pre_sker,post_idx)
end

@inline function postsynaptic_kernel_update!(pspost::PopulationState,
     conn::ConnGeneralIF2,pre_synkernel::SynapticKernel,pre_v_reversal::Float64,idx::Integer)
  return postsynaptic_kernel_update!(pspost.input,pspost.state_now,
    conn.post_trace1,conn.post_trace2,pre_synkernel,
    pre_v_reversal,idx)
end

@inline function trace_decay!(dt::Real,conn::ConnGeneralIF2,synker::SynapticKernel)
  return trace_decay!(dt,conn.post_trace1,conn.post_trace2,synker)
end
@inline function trace_decay!(dt::R,tr1::Vector{R},::Vector{R},sk::SKExp) where R
  tr1 .*= exp(-dt/sk.τ)
  return nothing
end
@inline function trace_decay!(dt::R,tr1::Vector{R},tr2::Vector{R},sk::SKExpDiff) where R
  tr1 .*= exp(-dt/sk.τ_plus)
  tr2 .*= exp(-dt/sk.τ_minus)
  return nothing
end

## define two main functions here
function local_update!(t_now::Float64,dt::Float64,ps::PSLIFConductance)
	# computes the update to internal voltage, given the total input
  # dv =  (v_leak - v ) dt / τ + I dt / Cap
  dttau =  dt / ps.neurontype.τ_post
  dtCap = dt / ps.neurontype.capacitance
  @inbounds for i in eachindex(ps.state_now)
    state_now = ps.state_now[i]
    state_now += (  (ps.neurontype.v_leak - state_now)*dttau +
      ps.input[i]*dtCap +
      ps.neurontype.spikegen(state_now)*dttau )
    ps.state_now[i] = state_now
  end
	# update spikes and refractoriness, and end
  return _spiking_state_update!(ps.state_now,ps.isfiring,ps.isrefractory,ps.last_fired,
    t_now,ps.neurontype.τ_refractory,ps.neurontype.v_threshold,ps.neurontype.v_reset)
end

function forward_signal!(t_now::Real,dt::Real,
      pspost::PSLIFConductance,conn::ConnGeneralIF2,pspre::PopulationState)
	post_idxs = rowvals(conn.weights) # postsynaptic neurons
	weightsnz = nonzeros(conn.weights) # direct access to weights
  pre_synker = pspre.neurontype.synaptic_kernel
  pre_v_reversal = pspre.neurontype.v_reversal
	for _pre in findall(pspre.isfiring)
		_posts_nz = nzrange(conn.weights,_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
			post_idx = post_idxs[_pnz]
      trace_spike_update!(conn,weightsnz[_pnz],pre_synker,post_idx)
		end
	end
	# add the currents to postsynaptic input
	# ONLY non-refractory neurons
	@inbounds @simd for i in eachindex(pspost.input)
		if ! pspost.isrefractory[i]
      postsynaptic_kernel_update!(pspost,conn,pre_synker,pre_v_reversal,i)
		end
	end
  # finally, all postsynaptic conductances decay in time
  trace_decay!(dt,conn,pre_synker)
  return nothing
end


text:
ENV["GKSwstype"] = "mp4"  # or "webm", "mov"

using GR

x = LinRange(0, 2π, 200)
for i = 0:100
    plot(x, sin.(x .+ i / 10.0))
end

ext: jl


text:
@testset "synchronous.jl" begin
    x, A = demodata()
    Ã = A .- mean(A, dims=1)

    Φ = synchronous(Ã)
    @test Φ ≈ cov(Ã)

    @plottest begin
        heatmap(x, x, Φ)
        plot!(xflip=true, yflip=true)
        xlabel!("Peak position")
        ylabel!("Peak position")
        title!("Synchronous")
    end "synchronous.png"

    @plottest begin
        plot(x, diag(Φ), label=:none)
        xlabel!("x")
        ylabel!("I")
    end "autocorrelation.png"
end


text:
"""
Base type for operations on a tape
"""
abstract type AbstractOp end

########################################################################
#                             VARIABLE                                 #
########################################################################

"""
Variable represents a reference to an operation on a tape.
Variables can be used to index tape or keep reference to
a specific operation on the tape.

Variables (also aliesed as `V`) can be:

* free, created as `V(id)` - used for indexing into tape
* bound, created as `V(op)`` - used to keep a robust reference
  to an operation on the tape
"""
mutable struct Variable
    _id::Union{<:Integer,Nothing}
    _op::Union{AbstractOp,Nothing}
end

Variable(id::Integer) = Variable(id, nothing)
Variable(op::AbstractOp) = Variable(nothing, op)

Base.show(io::IO, v::Variable) = print(io, "%$(v.id)")

function Base.getproperty(v::Variable, p::Symbol)
    if p == :id
        if v._op !== nothing
            # variable bound to a specific operation on a tapea
            return v._op.id
        else
            # free variable with only ID
            return v._id
        end
    else
        return getfield(v, p)
    end
end

function Base.setproperty!(v::Variable, p::Symbol, x)
    if p == :id
        if v._op !== nothing
            # variable bound to a specific operation on a tapea
            v._op.id = x
        else
            # free variable with only ID
            v.id = x
        end
    else
        return setfield!(v, p, x)
    end
end

function Base.:(==)(v1::Variable, v2::Variable)
    # variables are equal if:
    # * both are bound to the same operation, or
    # * both are unbound and their IDs are equal
    return v1._op === v2._op && v1.id == v2.id
end

Base.hash(v::Variable, h::UInt) = hash(v.id, hash(v._op, h))

const V = Variable

########################################################################
#                            OPERATIONS                                #
########################################################################

function Base.getproperty(op::AbstractOp, f::Symbol)
    if f == :typ
        return typeof(op.val)
    elseif f == :var
        return Variable(nothing, op)
    else
        getfield(op, f)
    end
end

## Input

"Operation representing input data of a tape"
mutable struct Input <: AbstractOp
    id::Int
    val::Any
end

Input(val::Any) = Input(0, val)

Base.show(io::IO, op::Input) = print(io, "inp %$(op.id)::$(op.typ)")

## Constant

"Operation representing a constant value on a tape"
mutable struct Constant <: AbstractOp
    id::Int
    typ::Type
    val
end

Constant(id::Int, val) = Constant(id, typeof(val), val)
Constant(val) = Constant(0, typeof(val), val)
Base.show(io::IO, op::Constant) = print(io, "const %$(op.id) = $(op.val)::$(op.typ)")

## Call

"""
Operation represening function call on tape. Typically, calls
are constructed using [`mkcall`](@ref) function.

Important fields of a Call{T}:

* `fn::T` - function or object to be called
* `args::Vector` - vector of variables or values used as arguments
* `val::Any` - the result of the function call
"""
mutable struct Call{T} <: AbstractOp
    id::Int
    val::Any
    fn::T
    args::Vector{Any}   # vector of Variables or const values
end
Call(id, val, fn::T, args) where T = Call{T}(id, val, fn, args)

pretty_type_name(T) = string(T)
pretty_type_name(T::Type{<:Broadcast.Broadcasted}) = "Broadcasted{}"

function Base.show(io::IO, op::Call)
    arg_str = join(["$v" for v in op.args], ", ")
    typ_str = pretty_type_name(op.typ)
    print(io, "%$(op.id) = $(op.fn)($arg_str)::$typ_str")
end

"""
Helper function to map a function only to Variable arguments of a Call
leaving constant values as is
"""
function map_vars(fn::Function, args::Union{Vector,Tuple})
    return map(v -> v isa Variable ? fn(v) : v, args)
end

"""
    mkcall(fn, args...; val=missing)

Convenient constructor for Call operation. If val is `missing` (default)
and call value can be calculated from (bound) variables and constants,
they are calculated. To prevent this behavior, set val to some neutral value.
"""
function mkcall(fn, args...; val=missing)
    fargs = (fn, args...)
    calculable = all(
        a -> !isa(a, Variable) ||                      # not variable
        (a._op !== nothing && a._op.val !== missing),  # bound variable
        fargs
    )
    if val === missing && calculable
        fargs_ = map_vars(v -> v._op.val, fargs)
        fn_, args_ = fargs_[1], fargs_[2:end]
        val_ = fn_(args_...)
    else
        val_ = val
    end
    return Call(0, val_, fn, [args...])
end

########################################################################
#                                 TAPE                                 #
########################################################################

"""
Linearized representation of a function execution.

Fields
======

* `ops` - vector of operations on the tape
* `result` - variable pointing to the operation to be used as the result
* `parent` - parent tape if any
* `meta` - internal metadata
* `c` - application-specific context
"""
mutable struct Tape{C}
    # linearized execution graph
    ops::Vector{<:AbstractOp}
    # result variable
    result::Variable
    # for subtapes - parent tape
    parent::Union{Tape,Nothing}
    # tape metadata (depends on the context)
    meta::Dict
    # application-specific context
    c::C
end

Tape(c::C) where C = Tape(AbstractOp[], Variable(0), nothing, Dict(), c)
# by default context is just a Dict{Any, Any}
Tape() = Tape(Dict{Any,Any}())

function Base.show(io::IO, tape::Tape{C}) where C
    println(io, "Tape{$C}")
    for op in tape.ops
        println(io, "  $op")
    end
end

function Base.getproperty(tape::Tape, p::Symbol)
    if p == :retval
        return tape[tape.result].val
    else
        return getfield(tape, p)
    end
end

"Get list of a tape input variables"
inputs(tape::Tape) = [V(op) for op in tape.ops if op isa Input]

"Set values of a tape inputs"
function inputs!(tape::Tape, vals...)
    @assert(isempty(tape) || length(inputs(tape)) == length(vals),
            "This tape contains $(length(inputs(tape))) inputs, but " *
            "$(length(vals)) value(s) were provided")
    if isempty(tape)
        # initialize inputs
        for val in vals
            push!(tape, Input(val))
        end
    else
        # rewrite input values
        for (i, val) in enumerate(vals)
            tape[V(i)].val = val
        end
    end
    return [V(op) for op in tape.ops[1:length(vals)]]
end

Base.getindex(tape::Tape, v::Variable) = tape.ops[v.id]

function Base.setindex!(tape::Tape, op::AbstractOp, v::Variable)
    op.id = v.id
    tape.ops[v.id] = op
    v._op = op   # bind to op, overriding v.id
end

Base.lastindex(tape::Tape) = lastindex(tape.ops)
Base.length(tape::Tape) = length(tape.ops)
Base.iterate(tape::Tape) = iterate(tape.ops)       # exclude inputs?
Base.iterate(tape::Tape, s) = iterate(tape.ops, s)

"""
    push!(tape::Tape, op::AbstractOp)

Push a new operation to the end of the tape.
"""
function Base.push!(tape::Tape, op::AbstractOp)
    new_id = length(tape) + 1
    op.id = new_id
    push!(tape.ops, op)
    return V(op)
end

"""
    insert!(tape::Tape, idx::Integer, ops::AbstractOp...)

Insert new operations into tape starting from position idx.
"""
function Base.insert!(tape::Tape, idx::Integer, ops::AbstractOp...)
    num_new_ops = length(ops)
    old_ops = tape.ops
    new_ops = Vector{AbstractOp}(undef, length(tape) + num_new_ops)
    # copy old ops before insertion point
    for i = 1:idx - 1
        new_ops[i] = old_ops[i]
    end
    # insert target ops, assign ids
    for i = 1:num_new_ops
        id = idx + i - 1
        new_ops[id] = ops[i]
        new_ops[id].id = id
    end
    # insert the rest of old ops
    for i = idx:length(old_ops)
        id = i + num_new_ops
        new_ops[id] = old_ops[i]
        new_ops[id].id = id
    end
    tape.ops = new_ops
    return [V(op) for op in ops]
end

"""
    replace!(tape, op  => new_ops; rebind_to=length(new_ops), old_new=Dict())

Replace specified operation with 1 or more other operations,
rebind variables in the reminder of the tape to ops[rebind_to].

Operation can be specified directly, by a variable or by ID.
"""
function Base.replace!(tape::Tape, idx_ops::Pair{<:Integer,<:Union{Tuple,Vector}};
                       rebind_to=length(idx_ops[2]), old_new=Dict{Int,Int}())
    idx, ops = idx_ops
    tape[V(idx)] = ops[1]
    if idx < length(tape)
        insert!(tape, idx + 1, ops[2:end]...)
    else
        for op in ops[2:end]
            push!(tape, op)
        end
    end

    st = merge(old_new, Dict(idx => ops[rebind_to].id))
    rebind!(tape, st; from=idx + length(ops))
    return V(ops[rebind_to])
end

Base.replace!(
    tape::Tape,
    idx_ops::Pair{Variable, <:Union{Tuple,Vector}};
    kwargs...) = replace!(tape, idx_ops[1].id => idx_ops[2]; kwargs...)

Base.replace!(
        tape::Tape,
        idx_ops::Pair{<:AbstractOp, <:Union{Tuple,Vector}};
        kwargs...) = replace!(tape, idx_ops[1].id => idx_ops[2]; kwargs...)

"""
    deleteat!(tape::Tape, idx; rebind_to = nothing)

Remove `tape[V(idx)]` from the `tape`.
If `rebind_to` is not `nothing`, then
replace all references to `V(idx)` with `V(rebind_to)`.

`idx` may be an index or `Variable`/`AbstractOp` directly.
"""
function Base.deleteat!(tape::Tape, idx::Integer; rebind_to = nothing)
    # delete and rebind
    deleteat!(tape.ops, idx)
    isnothing(rebind_to) || rebind!(tape, Dict(idx => rebind_to))

    # shift indices for outputs up
    for i in idx:length(tape)
        tape.ops[i].id -= 1
    end

    return tape
end
Base.deleteat!(tape::Tape, idx::Variable; kwargs...) =
    deleteat!(tape, idx.id; kwargs...)
Base.deleteat!(tape::Tape, idx::AbstractOp; kwargs...) =
    deleteat!(tape, idx.id; kwargs...)

########################################################################
#                       SPECIAL OPERATIONS                             #
########################################################################

## Loop

"""
Operation representing a loop in an computational graph.
See the online documentation for details.
"""
mutable struct Loop <: AbstractOp
    id::Int
    parent_inputs::Vector{Variable}
    condition::Variable
    cont_vars::Vector{Variable}
    exit_vars::Vector{Variable}
    subtape::Tape
    val::Any
end

function Base.show(io::IO, loop::Loop)
    input_str = join(map(string, loop.parent_inputs), ", ")
    print(io, "%$(loop.id) = Loop($input_str)")
end

###############################################################################
#                                 REBIND                                      #
###############################################################################

"""Returned version of the var bound to the tape op"""
bound(tape::Tape, v::Variable) = Variable(tape[v])

"""
    rebind!(tape::Tape, op, st::Dict)
    rebind!(tape::Tape, st::Dict; from, to)

Rebind all variables according to substitution table. Example:

    tape = Tape()
    v1, v2 = inputs!(tape, nothing, 3.0, 5.0)
    v3 = push!(tape, mkcall(*, v1, 2))
    st = Dict(v1.id => v2.id)
    rebind!(tape, st)
    @assert tape[v3].args[1].id == v2.id

See also: rebind_context!()
"""
function rebind!(tape::Tape, v::Variable, st::Dict)
    if haskey(st, v.id)
        # rebind to a new op
        v._op = tape[V(st[v.id])]
end
end

rebind!(::Tape, ::Input, ::Dict) = ()
rebind!(::Tape, ::Constant, ::Dict) = ()

function rebind!(tape::Tape, op::Call, st::Dict)
    for v in op.args
        if v isa Variable
            rebind!(tape, v, st)
        end
    end
    return op
end

"""
    rebind_context!(tape::Tape, st::Dict)

Rebind variables in the tape's context according to substitution table.
By default does nothing, but can be overwitten for specific Tape{C}
"""
rebind_context!(tape::Tape, st::Dict) = ()

function rebind!(tape::Tape, st::Dict; from=1, to=length(tape))
    for id = from:to
        rebind!(tape, tape[V(id)], st)
    end
    rebind!(tape, tape.result, st)
    rebind_context!(tape, st)
    return tape
end

########################################################################
#                              EXECUTION                               #
########################################################################

exec!(::Tape, ::Input) = ()
exec!(::Tape, ::Constant) = ()

function exec!(tape::Tape, op::Call)
    fn = op.fn isa V ? tape[op.fn].val : op.fn
    arg_vals = map_vars(v -> tape[v].val, op.args)
    op.val = fn(arg_vals...)
end

"""
Collect variables which will be used at loop exit if it happens
at this point on tape.
"""
function loop_exit_vars_at_point(op::Loop, id::Int)
    input_vars = inputs(op.subtape)
    exit_idxs = findall(v -> v in op.exit_vars, op.cont_vars)
    vars = Vector{Variable}(undef, length(exit_idxs))
    for (i, idx) in enumerate(exit_idxs)
        if id > op.cont_vars[idx].id
            # if condition is checked after this continue var is changed,
            # use continue var
            vars[i] = op.cont_vars[idx]
        else
            # otherwise use input var
            vars[i] = input_vars[idx]
        end
    end
    return vars
end

function exec!(tape::Tape, op::Loop)
    subtape = op.subtape
    # initialize inputs
    inputs!(subtape, [tape[v].val for v in op.parent_inputs]...)
    # run the loop strictly while continue condition is true
    # note that subtape execution may finish before the full
    # iteration is done
    cond_var = op.condition
    vi0 = length(op.parent_inputs) + 1
    vi = vi0
    while true
        # @show vi
        # @show subtape[V(1)].val
        # @show subtape[V(2)].val
        # @show subtape[V(7)].val
        # sleep(1)
        exec!(subtape, subtape[V(vi)])
        if vi == cond_var.id && subtape[V(vi)].val == false
            actual_exit_vars = loop_exit_vars_at_point(op, vi)
            op.val = ([v._op.val for v in actual_exit_vars]...,)
            break
        end
        vi += 1
        if vi > length(subtape)
            vi = vi0
            inputs!(subtape, [subtape[v].val for v in op.cont_vars]...)
        end
    end
    # # exit_var is special - it's a tuple combining all the exit variables
    # # since it doesn't exist in the original code, it may be not executed
    # # by loop logic at the last iteration; hence, we execute it manually
    # exec!(subtape, subtape[op.exit_var])
    # op.val = subtape[op.exit_var].val
end

"""
    play!(tape::Tape, args...; debug=false)

Execute operations on the tape one by one.
If `debug=true`, print each operation before execution.
"""
function play!(tape::Tape, args...; debug=false)
    for (i, val) in enumerate(args)
        @assert(tape[V(i)] isa Input, "More arguments than the original function had")
        tape[V(i)].val = val
    end
    for op in tape
        if debug
            println(op)
        end
        exec!(tape, op)
    end
    return tape[tape.result].val
end

########################################################################
#                                 UTILS                                #
########################################################################

"""
    call_signature(fn, args...)
    call_signature(tape::Tape, op::Call)

Get a signature of a function call. The obtain signature is suitable
for `is_primitive(sig)`.
"""
function call_signature(tape::Tape, op::Call)
    farg_vals = map_vars(v -> tape[v].val, [op.fn, op.args...])
    return Tuple{map(typeof, farg_vals)...}
end

function call_signature(fn, args...)
    return Tuple{map(typeof, (fn, args...))...}
end

"""
    primitivize!(tape::Tape; is_primitive=is_primitive)

Trace non-primitive function calls on a tape and decompose them
into a list of corresponding primitive calls.

# Example

    f(x) = 2x - 1
    g(x) = f(x) + 5

    tape = Tape()
    _, x = inputs!(tape, g, 3.0)
    y = push!(tape, mkcall(f, x))
    z = push!(tape, mkcall(+, y, 5))
    tape.result = z

    primitivize!(tape)

    # output

    Tape{Dict{Any, Any}}
      inp %1::typeof(g)
      inp %2::Float64
      %3 = *(2, %2)::Float64
      %4 = -(%3, 1)::Float64
      %5 = +(%4, 5)::Float64
"""
function primitivize!(tape::Tape, op::AbstractOp)
    id = op.id
    fn = op.fn isa V ? tape[op.fn].val : op.fn
    args = map_vars(a -> tape[a].val, op.args)
    _, sub = trace(fn, args...)

    new_ops = sub.ops[length(inputs(sub))+1:end]
    old_new = Dict{Int, Int}()
    for (i, v) in enumerate((op.fn, op.args...))
        if v isa V
            old_new[i] = v.id
        end
    end
    replace!(tape, id => new_ops, old_new=old_new)
    # note: not touching the context since replacement
    # may be ambiguous for it
end

function primitivize!(tape::Tape; is_primitive=is_primitive)
    # note: referencing concrete operations on the original tape
    # they will stay the same even when we modify the tape
    vars = [V(op) for op in tape]
    for v in vars
        op = tape[v]
        if op isa Call && !is_primitive(call_signature(tape, op))
            primitivize!(tape, op)
        end
    end
end


text:
using Quadhps
using ForwardDiff, QuadGK, Test

@testset "Quadhps tests" begin

  @testset "Basic" begin
    a = Quadhps.quadhp(x->1.0, 0.0, 2.0, rtol=2*eps())
    @test a[1] ≈ 2.0 atol=10 * eps()
    @test a[2] ≈ 0.0 atol=10 * eps()
    b = Quadhps.quadhp(x->x, 1.0, 2.0, rtol=2*eps())
    @test b[1] ≈ 1.5 atol=10 * eps()
    @test b[2] ≈ 0.0 atol=10 * eps()
    c = Quadhps.quadhp(x->3 .* x.^2 .+ 1.0, 2.0, 5.0, rtol=120*eps())
    @test c[1] ≈ 120.0 atol=1000 * eps()
    @test c[2] ≈ 0.0 atol=1000 * eps()
  end

  @testset "Versus QuadGK" begin
    f(x) = 1 ./ (1 .- x.^2)
    lower = -0.999
    upper = 0.999
    a = Quadhps.quadhp(x->f(x), lower, upper, rtol=2*eps())[1]
    b = QuadGK.quadgk(x->f(x), lower, upper, rtol=2*eps())[1]
    @test a ≈ b atol=1000 * eps()
  end

  @testset "Types" begin
    f(x) = x
    lower = BigFloat(0.1)
    upper = 10.0
    a = Quadhps.quadhp(x->f(x), lower, upper)[1]
    b = QuadGK.quadgk(x->f(x), lower, upper)[1]
    @test a ≈ b rtol=eps()
  end

  @testset "Vectorised" begin
    a = Quadhps.quadhp(x->[1, 2], 0.0, 1.0)[1]
    @test a[1] ≈ 1 atol=2 * eps()
    @test a[2] ≈ 2 atol=2 * eps()
  end

  @testset "Intervals" begin
    a = Quadhps.quadhp(x->3 .* x.^2 .+ 1.0, 0.0, 5.0, rtol=120*eps())
    b = Quadhps.quadhp(x->3 .* x.^2 .+ 1.0, 0.0, 2.0, 5.0, rtol=120*eps())
    c = Quadhps.quadhp(x->3 .* x.^2 .+ 1.0, 0.0, 2.0, 4.0, 5.0, rtol=120*eps())
    @test a[1] ≈ b[1]
    @test a[1] ≈ c[1]
  end

  @testset "ForwadDiff" begin
    integrand = sin
    differentiand(x) = Quadhps.quadhp(y->integrand(y[1]), 0.0, x[1])[1]
    b = π/2
    @test integrand(b) ≈ ForwardDiff.gradient(differentiand, [b])[1]
    a = ForwardDiff.gradient(x->
      Quadhps.quadhp(y->sin(y[1]*x[1]), 0.0, 1.0)[1],
          [0.5])[1]
    @test a ≈ 0.469181324769897
  end
end


text:
##############################################################################
##
## Iterate on terms
##
##############################################################################

eachterm(@nospecialize(x::AbstractTerm)) = (x,)
eachterm(@nospecialize(x::NTuple{N, AbstractTerm})) where {N} = x
TermOrTerms = Union{AbstractTerm, NTuple{N, AbstractTerm} where N}

##############################################################################
##
## Parse FixedEffect
##
##############################################################################
fesymbol(t::FixedEffectModels.FixedEffectTerm) = t.x
fesymbol(t::FunctionTerm{typeof(fe)}) = Symbol(t.args_parsed[1])


text:
"Validates the given plot_attributes according to their type"
function _validate_plot_attributes!(plot_attributes::Dict{Symbol, Any})
    for attr in keys(plot_attributes)
      if !haskey(default_plot_attributes, attr)
        Memento.warn(_LOGGER, "Ignoring unexpected attribute $(repr(attr))")
      end
    end

    # validate color attributes
    for attr in _color_attributes
        color = plot_attributes[attr]
        if !(typeof(color) <: Union{String, Symbol, AbstractVector})
          Memento.warn(_LOGGER, "Color value for $(repr(attr)) should be given as symbol or string")
        else
          try
            if typeof(color) <: AbstractVector
                parse.(Colors.Colorant, color) # parses all colors as CSS color
            else
                parse(Colors.Colorant, color) # try to parse the color as a CSS color
                plot_attributes[attr] = [color] # package color into an array
            end
          catch e
            Memento.warn(_LOGGER, "Invalid color $(repr(color)) given for $(repr(attr))")
          end
        end
    end

    # validate numeric attributes
    for attr in _numeric_attributes
      value = plot_attributes[attr]
      if !(typeof(value) <: Union{Number, String})
        Memento.warn(_LOGGER, "Value for $(repr(attr)) should be given as a number or numeric String")
      elseif typeof(value) <: String
        try
            parse(Float64, value)
        catch e
            Memento.warn(_LOGGER, "Invalid number $(repr(value)) given for $(repr(attr))")
        end
      end
    end

    # validate data label attributes
    for attr in _label_attributes
      value = plot_attributes[attr]
      if !(typeof(value) <: Union{String, Symbol})
        Memento.warn(_LOGGER, "Value for $(repr(attr)) should be given as a String or Symbol")
      end
    end
end

"Checks that the given column plot_attributes[data_attr] exists in the data"
function _validate_data(data::DataFrames.DataFrame, data_column::Any, data_name::String)
    if !(typeof(data_column) <: Union{String, Symbol})
        return
    end
    if !(data_column in names(data) || data_column in propertynames(data))
        Memento.warn(_LOGGER, "Data column $(repr(data_column)) does not exist for $(data_name)")
    end
end

"Checks that the given data type attribute is a valid VegaLite data type"
function _validate_data_type(plot_attributes::Dict{Symbol, Any}, attr::Symbol)
    valid_types = Set([:quantitative, :temporal, :ordinal, :nominal, :geojson])
    data_type = plot_attributes[attr]
    if !(Symbol(data_type) in valid_types)
        Memento.warn(_LOGGER, "Data type $(repr(data_type)) not a valid VegaLite data type")
    end
end


text:
immutable ElError <: Exception
end

immutable ElRuntimeError <: Exception
end

immutable ElLogicError <: Exception
end

immutable ElArgumentError <: Exception
end

immutable ElBoundsError <: Exception
end

immutable ElAllocError <: Exception
end

function ElError(code::Integer)
    if code == 1
        return OutOfMemoryError()
    elseif code == 2
        return ArgumentError()
    elseif code == 3
        return ElLogicError()
    elseif code == 4
        return ElRuntimeError()
    elseif code == -1
        # catchall error code
        return ElError()
    else
        throw(ArgumentError("Unknown Elemental Error Code: $code"))
    end
end

ErrorString(code::Integer) =
    btyestring(ccall((:ElErrorString, libEl), Ptr{Cchar}, (Cint,), code))


text:
# estimators for ratio of means of two time series (including blocking)

"""
    RatioBlockingResult(ratio, f, σ_f, δ_y, k, success)
Result of [`ratio_of_means()`](@ref).

### Fields:
- `ratio::P`: ratio with uncertainties propagated by MonteCarloMeasurements
- `f::T`: ratio of means
- `σ_f::T`: std from linear propagation
- `δ_y::T`: coefficient of variation for denominator (≤ 0.1 for normal approx)
- `k::Int`: k-1 blocking steps were used to uncorrelate time series
- `blocks::Int`: number of data values after blocking
- `success::Bool`: false if any of the blocking steps failed

Has methods for [`NamedTuple`](@ref), [`val_and_errs`](@ref), [`val`](@ref), [`errs`](@ref).

Note: to compute statistics on the `RatioBlockingResult`, use functions `pmedian`,
`pquantile`, `pmiddle`, `piterate`, `pextrema`, `pminimum`, `pmaximum`, `pmean`, and `pcov`.
"""
struct RatioBlockingResult{T,P}
    ratio::P    # ratio with uncertainties propagated by MonteCarloMeasurements
    f::T        # ratio of means
    σ_f::T      # std from linear propagation
    δ_y::T      # coefficient of variation for denominator (≤ 0.1 for normal approx)
    k::Int      # k-1 blocking steps were used to uncorrelate time series
    blocks::Int # number of data values after blocking
    success::Bool # false if any of the blocking steps failed
end
# make it behave like Particles
# Base.iterate(r::RatioBlockingResult, args...) = iterate(r.ratio, args...)
# Base.length(r::RatioBlockingResult) = length(r.ratio)
# Base.eltype(r::RatioBlockingResult) = eltype(r.ratio)
MacroTools.@forward RatioBlockingResult.ratio MonteCarloMeasurements.pmedian, MonteCarloMeasurements.pquantile
MacroTools.@forward RatioBlockingResult.ratio MonteCarloMeasurements.pmiddle, MonteCarloMeasurements.piterate, MonteCarloMeasurements.pextrema
MacroTools.@forward RatioBlockingResult.ratio MonteCarloMeasurements.pminimum, MonteCarloMeasurements.pmaximum
MacroTools.@forward RatioBlockingResult.ratio MonteCarloMeasurements.pmean, MonteCarloMeasurements.pcov

import Statistics: median, quantile, mean, cov
import Base: iterate, extrema, minimum, maximum
@deprecate median(r::RatioBlockingResult) pmedian(r)
@deprecate quantile(r::RatioBlockingResult, args...) pquantile(r, args...)
@deprecate iterate(r::RatioBlockingResult, args...) piterate(r, args...)
@deprecate extrema(r::RatioBlockingResult) pextrema(r)
@deprecate minimum(r::RatioBlockingResult) pminimum(r)
@deprecate maximum(r::RatioBlockingResult) pmaximum(r)
@deprecate mean(r::RatioBlockingResult) pmean(r)
@deprecate cov(r::RatioBlockingResult) pcov(r)

function MonteCarloMeasurements.pmedian(r::RatioBlockingResult{<:Complex})
    complex(pmedian(real(r.ratio)), pmedian(imag(r.ratio)))
end
function MonteCarloMeasurements.pquantile(r::RatioBlockingResult{<:Complex}, args...)
    throw(ErrorException("""
    `pquantile(r, args...)` called with `Complex` data type.
    Try `pquantile(real(r.ratio), args...)` and `pquantile(imag(r.ratio), args...)`!
    """
    ))
end

function Base.show(io::IO, r::RatioBlockingResult{T,P}) where {T<:Real, P}
    q = pquantile(r.ratio, [0.16,0.5,0.84])
    qr95 = pquantile(r.ratio, [0.025,0.975])
    println(io, "RatioBlockingResult{$T,$P}")
    println(io, f"  ratio = \%g(q[2]) ± (\%g(q[3]-q[2]), \%g(q[2]-q[1])) (MC)")
    println(io, f"  95% confidence interval: [\%g(qr95[1]), \%g(qr95[2])]) (MC)")
    println(io, f"  linear error propagation: \%g(r.f) ± \%g(r.σ_f)")
    println(io, f"  |δ_y| = |\%g(r.δ_y)| (≤ 0.1 for normal approx)")
    if r.success
        println(io, "  Blocking successful with $(r.blocks) blocks after $(r.k-1) transformations (k = $(r.k)).")
    else
        println(io, "  Blocking unsuccessful.")
    end
end

function Base.show(io::IO, r::RatioBlockingResult{T,P}) where {T<:Complex, P}
    qr = pquantile(real(r.ratio), [0.16,0.5,0.84])
    qi = pquantile(imag(r.ratio), [0.16,0.5,0.84])
    println(io, "RatioBlockingResult{$T,$P}")
    println(io, f"  ratio = \%g(qr[2]) ± (\%g(qr[3]-qr[2]), \%g(qr[2]-qr[1])) + [\%g(qi[2]) ± (\%g(qi[3]-qi[2]), \%g(qi[2]-qi[1]))]*im (MC)")
    qr95 = pquantile(real(r.ratio), [0.025,0.975])
    qi95 = pquantile(imag(r.ratio), [0.025,0.975])
    println(io, f"  95% confidence interval real: [\%g(qr95[1]), \%g(qr95[2])] (MC)")
    println(io, f"  95% confidence interval imag: [\%g(qi95[1]), \%g(qi95[2])] (MC)")
    println(io, "  linear error propagation: ($(r.f)) ± ($(r.σ_f))")
    println(io, f"  |δ_y| = \%g(abs(r.δ_y)) (≤ 0.1 for normal approx)")
    if r.success
        println(io, "  Blocking successful with $(r.blocks) blocks after $(r.k-1) transformations (k = $(r.k)).")
    else
        println(io, "  Blocking unsuccessful.")
    end
end

"""
    ratio_of_means(num, denom; α=0.01, corrected=true, mc_samples=nothing, skip=0) -> r
Estimate the ratio of `mean(num)/mean(denom)` assuming that `num` and `denom` are possibly
correlated time series, skipping the first `skip` elements. A blocking analysis with
m-test is used to uncorrelate the time series, see [`blocking_analysis()`](@ref). The
remaining standard error and correlation of the means is propagated using
`MonteCarloMeasurements`. The results are reported as a [`RatioBlockingResult`](@ref).

Robust estimates for the ratio are obtained from `pmedian(r)` and confidence intervals from
`pquantile()`, e.g. `pquantile(r, [0.025, 0.975])` for the 95% confidence interval.

Estimates from linear uncertainty propagation are returned as `r.f` and `r.σ_f` using
[`x_by_y_linear()`](@ref).

The keyword `mc_samples` controls the number of samples used for error propagation
by `MonteCarloMeasurements`. Use `nothing` for the default and `Val(1000)` to set the number
to 1000 samples in a type-consistent way.

Note: to compute statistics on the [`RatioBlockingResult`](@ref), use functions `pmedian`,
`pquantile`, `pmiddle`, `piterate`, `pextrema`, `pminimum`, `pmaximum`, `pmean`, and `pcov`.
"""
function ratio_of_means(num, denom; α=0.01, corrected=true, mc_samples=nothing, skip=0)
    num = @view num[skip+1:end]
    denom = @view denom[skip+1:end]
    # determine how many blocking steps are needed to uncorrelate data
    bt_num = blocking_analysis(num; α)
    bt_den = blocking_analysis(denom; α)
    # choose largst k from mtest on blocking analyses on numerator and denominator
    ks = (bt_num.k, bt_den.k)
    success = any(k->k<0, ks) ? false : true
    k = max(ks...)

    # MC and linear error propagation
    r, f, σ_f, δ_y, blocks = ratio_estimators(num, denom, k; corrected, mc_samples)

    # more accurately before blocking
    f = bt_num.mean/bt_den.mean
    return RatioBlockingResult(r, f, σ_f, δ_y, k, blocks, success)
end

"""
    x_by_y_linear(μ_x,μ_y,σ_x,σ_y,ρ) -> f, σ_f
Linear error propagation for ratio `f = x/y` assuming `x` and `y` are correlated normal
random variables and assuming the ratio can be approximated as a normal distribution.
See [wikipedia](https://en.wikipedia.org/wiki/Propagation_of_uncertainty) and
[Díaz-Francés, Rubio (2013)](http://link.springer.com/10.1007/s00362-012-0429-2).
```math
σ_f = \\sqrt{\\frac{σ_x}{μ_y}^2 + \\frac{μ_x σ_y}{μ_y^2}^2 - \\frac{2 ρ μ_x}{μ_y^3}}
```
"""
function x_by_y_linear(μ_x, μ_y, σ_x, σ_y, ρ)
    f = μ_x/μ_y
    # σ_f = abs(f)*sqrt((σ_x/μ_x)^2 + (σ_y/μ_y)^2 - 2*ρ/(μ_x*μ_y))
    σ_f = sqrt((σ_x/μ_y)^2 + (μ_x*σ_y/μ_y^2)^2 - 2*ρ*μ_x/μ_y^3)
    return f, σ_f
end

"""
    particles(samples, d)
    particles(::Nothing, d)
    particles(::Val{T}, d) where T
Return `Particles` object from  `MonteCarloMeasurements` using  a type-stable constructor
if possible. Pass `nothing` for the default number of particles or `Val(1_000)` for using
1000 particles in a type-stable manner. If `d` is a `Particles` object it is passed through
without re-sampling.
"""
particles(samples, d::Distribution) = particles(Val(samples),d)
particles(::Nothing, d::Distribution) = Particles(d)
particles(::Val{T}, d::Distribution) where T = Particles{eltype(d),T}(Random.GLOBAL_RNG,d)
function particles(samples, m::Measurements.Measurement)
    particles(samples, Measurements.value(m), Measurements.uncertainty(m))
end
particles(_, p::Particles) = p # don't re-sample if it is already a Particles object
"""
    particles(samples, μ, σ)
    particles(samples, μ::AbstractVector, Σ::AbstractMatrix)
Return `Particles` object from `MonteCarloMeasurements` with single- or multivariate
normal distribution. Zero variance parameters are supported.
"""
particles(samples, μ, σ) = particles(samples, Normal(μ, σ))
function particles(samples, μ::AbstractVector, Σ::AbstractMatrix)
    singular_dim = map(x->abs(x)<√eps(x), diag(Σ))
    s = sum(singular_dim) # number of near-zero variances
    if iszero(s) # assume Σ is positive definite
        return particles(samples, MvNormal(μ, Σ))
    else
        l = length(singular_dim)
        Σ_reg = [Σ[i,j] for i in 1:l, j in 1:l if !singular_dim[i] && !singular_dim[j]]
        Σ_reg = reshape(Σ_reg, (l-s,l-s)) # make sure it's a matrix!
        μ_reg = [μ[i] for i in 1:l if !singular_dim[i]]
        p_reg = particles(samples, MvNormal(μ_reg, Σ_reg)) # Vector{Particles}
        p = similar(p_reg, 0) # has the right type but is empty
        count_reg = 0
        for i in 1:l
            if singular_dim[i]
                push!(p, particles(samples, μ[i], √Σ[i,i]))
            else
                count_reg += 1
                push!(p, p_reg[count_reg])
            end
        end
        return p
    end
end

"""
    ratio_estimators(x, y, [k]; corrected=true, mc_samples=10_000) -> (; r, f, σ_f, δ_y, n)
Estimators for the ratio of means `mean(x)/mean(y)`.
If `k` is given, `k-1` blocking steps are performed to remove internal correlations in
the time series `x` and `y`. Otherwise these are assumed to be free of internal
correlations. Correlations between `x` and `y` may be present and are taken into
account.

### Return values:
- `r::Particles` is the Monte Carlo sampled ratio estimator, see [`Particles`](@ref)
- `f = mean(x)/mean(y)`
- `σ_f` standard deviation of `f` from linear error propagation (normal approximation)
- `δ_y = std(y)/mean(y)` coefficient of variation; < 0.1 for normal approximation to work
- `n`: number of uncorrelated data used for uncertainty estimation
"""
function ratio_estimators(
    x::AbstractVector{<:Real}, y::AbstractVector{<:Real};
    corrected = true, mc_samples = nothing,
)
    n = length(x)
    @assert n == length(y)
    μ_x = mean(x)
    var_x = var(x; corrected)/n # variance of mean
    μ_y = mean(y)
    var_y = var(y; corrected)/n # variance of mean
    ρ = cov(x, y; corrected)/n # estimated correlation of sample means μ_x and μ_y

    # Monte Carlo sampling of correlated normal distribution of sample means for x and y
    # x_y_ps = particles(mc_samples, (μ_x, μ_y, var_x, var_y, ρ))
    x_y_ps = particles(mc_samples, [μ_x,μ_y],[var_x ρ; ρ var_y])
    # Note: type instability creeps in here through `Particles`
    r = x_y_ps[1]/x_y_ps[2] # MC sampled ratio of means

    # linear error propagation
    f, σ_f = x_by_y_linear(μ_x, μ_y, √var_x, √var_y, ρ)

    # coefficient of variation, should be <0.1 for normal approximation
    # [Kuethe(2000), Diaz-Frances & Rubio (2013)]
    δ_y = √var_y/μ_y
    return (; r, f, σ_f, δ_y, n)
end

function ratio_estimators(num, denom, k; corrected = true, mc_samples = nothing)
    for i in 1:(k-1) # decorrelate time series by `k-1` blocking steps
        num = blocker(num)
        denom = blocker(denom)
    end
    return ratio_estimators(num, denom; corrected, mc_samples)
end

# x or y could be complex
function ratio_estimators(x, y; corrected = true, mc_samples = nothing)
    n = length(x)
    @assert n == length(y)
    μ_x = mean(x)
    var_x = var(x; corrected)/n # variance of mean
    μ_y = mean(y)
    var_y = var(y; corrected)/n # variance of mean
    ρ = cov(x, y; corrected)/n # estimated correlation of sample means μ_x and μ_y

    Σ = [
        var(real(x); corrected)          cov(real(x), imag(x); corrected) cov(real(x), real(y); corrected) cov(real(x), imag(y); corrected)
        cov(imag(x), real(x); corrected) var(imag(x); corrected)          cov(imag(x), real(y); corrected) cov(imag(x), imag(y); corrected)
        cov(real(y), real(x); corrected) cov(real(y), imag(x); corrected) var(real(y); corrected)          cov(real(y), imag(y); corrected)
        cov(imag(y), real(x); corrected) cov(imag(y), imag(x); corrected) cov(imag(y), real(y); corrected) var(imag(y); corrected)
    ]/n
    # Monte Carlo sampling of correlated normal distribution of sample means for x and y
    x_y_ps = particles(mc_samples, [real(μ_x), imag(μ_x), real(μ_y), imag(μ_y)], Σ)
    # MC sampled ratio of means
    r = Base.FastMath.div_fast(x_y_ps[1] + im*x_y_ps[2], x_y_ps[3] + im*x_y_ps[4])
    # r = (x_y_ps[1] + im*x_y_ps[2]) / (x_y_ps[3] + im*x_y_ps[4])

    # linear error propagation
    f, σ_f = x_by_y_linear(μ_x, μ_y, √var_x, √var_y, ρ)

    # coefficient of variation, should be <0.1 for normal approximation
    # [Kuethe(2000), Diaz-Frances & Rubio (2013)]
    T = promote_type(typeof(μ_x),typeof(μ_y))
    δ_y = T(√var_y/μ_y)
    return (; r, f, σ_f, δ_y, n)
end


text:
@testset "consistency_bars.jl" begin
    # default values
    bars = ConsistencyBars()
    @test bars.samples == 1_000
    @test bars.quantiles[1] ≈ 0.025
    @test bars.quantiles[2] ≈ 0.975

    # exceptions
    @test_throws ErrorException ConsistencyBars(; samples=0)
    @test_throws ErrorException ConsistencyBars(; samples=-1)
    @test_throws ErrorException ConsistencyBars(; coverage=-0.5)
    @test_throws ErrorException ConsistencyBars(; coverage=1.4)

    # random values
    samples = ceil(Int, abs(randn()))
    coverage = rand()
    bars = ConsistencyBars(; samples=samples, coverage=coverage)
    @test bars.samples == samples
    @test bars.quantiles == promote((1 - coverage) / 2, (1 + coverage) / 2)
end


text:
using Optimisers
using ChainRulesCore, Functors, StaticArrays, Zygote
using LinearAlgebra, Statistics, Test, Random
using Optimisers: @.., @lazy

Random.seed!(1)

struct Foo; x; y; end
Functors.@functor Foo
Optimisers.trainable(x::Foo) = (x.y, x.x)

struct TwoThirds a; b; c; end
Functors.@functor TwoThirds (a, c)
Optimisers.trainable(x::TwoThirds) = (a = x.a,)

@testset verbose=true "Optimisers.jl" begin
  @testset verbose=true "Features" begin

    @testset "very basics" begin
      m = ([1.0, 2.0],)
      mid = objectid(m[1])
      g = ([25, 33],)
      o = Descent(0.1)
      s = Optimisers.setup(o, m)

      s2, m2 = Optimisers.update(s, m, g)
      @test m[1] == 1:2  # not mutated
      @test Optimisers.iswriteable(m[1])
      @test m2[1] ≈ [1,2] .- 0.1 .* [25, 33]

      s3, m3 = Optimisers.update!(s, m, g)
      @test objectid(m3[1]) == mid
      @test m3[1] ≈ [1,2] .- 0.1 .* [25, 33]

      g4 = Tangent{typeof(m)}(g...)
      s4, m4 = Optimisers.update!(s, ([1.0, 2.0],), g4)
      @test m4[1] ≈ [1,2] .- 0.1 .* [25, 33]
    end

    @testset "gradient clipping" begin
      m = (α = ([0.0], sin), γ = rand(3))
      s1 = Optimisers.setup(ClipGrad(13), m)
      _, m1 = Optimisers.update(s1, m, (α = nothing, γ = [1,10,100],))
      @test m.γ .- m1.γ ≈ [1, 10, 13]

      s2 = Optimisers.setup(ClipNorm(10), m)
      _, m2 = Optimisers.update(s2, m, (α = ([0.1], nothing), γ = [1,10,100],))
      @test only(m.α[1] .- m2.α[1]) ≈ 0.1
      @test norm(m.γ .- m2.γ) ≈ 10
      @test_throws DomainError Optimisers.update(s2, m, (α = [0.1], γ = [1,10,NaN],))

      s3 = Optimisers.setup(ClipNorm(5, 1; throw=false), m)
      _, m3 = Optimisers.update(s3, m, (α = ([0.1], nothing), γ = [1,10,100],))
      @test only(m.α[1] .- m3.α[1]) ≈ 0.1
      @test norm(m.γ .- m3.γ, 1) ≈ 5
      _, m3n = Optimisers.update!(s3, m, (α = nothing, γ = [1,10,Inf],))
      @test isnan(m3n.γ[3])
    end

    @testset "OptimiserChain" begin
      x = [1, 10, 100.0]; dx = [1, 2, 3.0];
      @test Optimisers.update(Optimisers.setup(WeightDecay(0.1), x), x, dx)[2] ≈ [1-0.1-1, 10-1-2, 100-10-3]
      @test Optimisers.update(Optimisers.setup(ClipGrad(2), x), x, dx)[2] ≈ [1-1, 10-2, 100-2]

      o2 = OptimiserChain(ClipGrad(2), WeightDecay(0.1))
      @test Optimisers.update(Optimisers.setup(o2, x), x, dx)[2] ≈ [1-0.1-1, 10-1-2, 100-10-2]

      o2n = OptimiserChain(OptimiserChain(ClipGrad(2), WeightDecay(0.1)))  # nested
      @test Optimisers.update(Optimisers.setup(o2n, x), x, dx)[2] ≈ [1-0.1-1, 10-1-2, 100-10-2]

      o2r = OptimiserChain(WeightDecay(0.1), ClipGrad(2))  # reversed
      @test Optimisers.update(Optimisers.setup(o2r, x), x, dx)[2] != [1-0.1-1, 10-2, 100-2]

      # Trivial cases
      o1 = OptimiserChain(Descent(0.1))
      @test Optimisers.update(Optimisers.setup(o1, x), x, dx)[2] ≈ [0.9, 9.8, 99.7]

      o0 = OptimiserChain()
      @test Optimisers.update(Optimisers.setup(o0, x), x, dx)[2] ≈ [1-1,10-2,100-3]
    end

    @testset "trainable subset" begin
      @info "ignore these warnings about trainable, testing the old path"
      # Foo has an old-style tuple trainable, both elements
      mf = Foo([1.0, 2.0], (a = sin, b = [3.0, 4.0], c = 5))
      sf = Optimisers.setup(Descent(0.1), mf)
      gf = (x = nothing, y = (a = nothing, b = [1,1], c = 1))
      _, mf2 = Optimisers.update(sf, mf, gf)
      @test mf2.x == [1,2]
      @test mf2.y == (a = sin, b = [2.9, 3.9], c = 5)

      gf3 = Tangent{typeof(mf)}(; x = NoTangent(), y = Tangent{typeof(mf.y)}(; a = NoTangent(), b = [1,1], c = 1))
      _, mf3 = Optimisers.update(sf, mf, gf3)  # the same, but with ChainRules types
      @test mf3.x == [1,2]
      @test mf3.y == (a = sin, b = [2.9, 3.9], c = 5)

      # TwoThirds has functor a,c only, and trainable a only
      mt = TwoThirds(Float32[1,2], Float32[3,4], Float32[5,6])
      mt10 = fmap(x -> 10x, mt)
      @test mt10.a == [10, 20]
      @test mt10.b == [3, 4]
      @test mt10.c == [50, 60]
      st = Optimisers.setup(Momentum(0.1, 0.9), mt)
      gt = gradient(m -> sum(abs2, m.a) + 100sum(abs2, m.b), mt)
      _, mtup = Optimisers.update(st, mt, gt...)
      @test mtup.a ≈ [0.8, 1.6]
      @test mtup.b == [3, 4]
      @test mtup.c == [5, 6]

      # Various kinds of missing branches together:
      m = Foo(
          TwoThirds(Foo(1.0, Float32[2,3,4]), 5.0, Float32[6,7]),
          TwoThirds((p = Float32[1,2,3],), sin, (q = 4.0, r = cos,)),
          )
      s = Optimisers.setup(Momentum(0.1, 0.9), m)
      g = gradient(m -> sum(abs2, m.x.a.y) + m.x.b^2 + log(m.y.c.q), m)
      @test Optimisers.update!(s, m, g...)[2] isa Foo
    end

    @testset "eltype preservation" begin
      m = (Float16[1,2], Float32[3,4])
      s1 = Optimisers.setup(Descent(0.1), m)
      s2, m2 = Optimisers.update(s1, m, m)
      @test eltype(m2[1]) == Float16  # because update copies & calls update!
      @test eltype(m2[2]) == Float32

      staticm = (SA{Float16}[1,2], SA{Float32}[3,4])
      s3 = Optimisers.setup(Descent(0.1), staticm)
      s4, m4 = Optimisers.update(s3, staticm, staticm)
      @test eltype(m4[1]) == Float16  # because of explicit broadcast in subtract!
      @test eltype(m4[2]) == Float32
    end

    @testset "forgotten gradient" begin
      x = [1.0, 2.0]
      sx = Optimisers.setup(Descent(), x)
      @test_throws MethodError Optimisers.update(sx, x)

      m = (x = x, y = sin)
      sm = Optimisers.setup(Descent(), m)
      @test_throws MethodError Optimisers.update(sm, m)
    end

    @testset "broadcasting macros" begin
      x = [1.0, 2.0]; y = [3,4]; z = [5,6]
      @test (@lazy x + y * z) isa Broadcast.Broadcasted
      bc = @lazy x + y * z
      @test (y .+ 2 .* bc) == [35,56]

      @test (@.. x = y * z) isa Array
      @test x == y .* z  # mutated
      r = 1.0:2.0  # immutable
      @test (@.. r = y * z) isa Array
      @test (@.. r = y * z) == y .* z
      @.. r = y * z
      @test r == y .* z  # attaches name r to result
    end

    @testset "tied weights" begin
      ok = (1.0:3.0, sin, "abc", :abc)
      m = (α = ok, β = rand(3), γ = ok)
      m1 = (rand(3), m, rand(3))
      @test Optimisers.setup(ADAMW(), m1) isa Tuple
      m2 = (rand(3), m, rand(3), m, rand(3))  # illegal
      @test_throws ArgumentError Optimisers.setup(ADAMW(), m2)
    end

  end
  @testset verbose=true "Destructure" begin
    include("destructure.jl")
  end
  @testset verbose=true "Optimisation Rules" begin
    include("rules.jl")
  end
end


text:
"""
测试重新切分
"""

using Test
using MARY_fRG.Refine, MARY_fRG.Basics, MARY_fRG.Triangulated
using MARY_fRG.Drawers

@testset "三角形的切分" begin
    rtri = RtTriangle(Point2D(0., 0.), Point2D(1., 1.), Point2D(1., -1.))
    newtris = split_triangle(rtri)
    println(newtris)
    plt = draw_polygon(newtris)
    @savepng plt "newtris"
    eqtri = EqTriangle(Point2D(0., 0.), Point2D(1., sqrt(3)), Point2D(2., 0.))
    newtris = split_triangle(eqtri)
    println(newtris)
    plt = draw_polygon(newtris)
    @savepng plt "newtris2"
end

@testset "正方形的细化" begin
    squ1 = Square(Point2D(0., 0.), 1)
    ltris1, ladjs1 = split_square(squ1, 3)
    plt = draw_polygon(ltris1)
    @savepng plt "refine1"
    ltris2 = refine_triangles(ltris1)
    plt = draw_polygon(ltris2)
    @savepng plt "refine2"
end

@testset "六边型的细化" begin
    hexa1 = EqHexagon(Point2D(0., 0.), 1)
    ltris1, ladjs1 = split_hexagon(hexa1, 4)
    plt = draw_polygon(ltris1)
    @savepng plt "refine3"
    ltris2 = refine_triangles(ltris1)
    plt = draw_polygon(ltris2)
    @savepng plt "refine4"
end

@testset "寻找最近邻" begin
    #正方格子
    squ1 = Square(Point2D(0., 0.), 1.)
    ltris, ladjs = split_square(squ1, 11)
    ladjs2 = find_adjs_by_adjoint(ltris)
    save_triangulated(squ1, ltris, ladjs, "squ1")
    save_triangulated(squ1, ltris, ladjs2, "squ2")
    fl1 = open("squ1.adj", "r")
    fl2 = open("squ2.adj", "r")
    @test read(fl1) == read(fl2)
    #六角格子
    hexa1 = EqHexagon(Point2D(0., 0.), 1.)
    ltris, ladjs = split_hexagon(hexa1, 13)
    ladjs2 = find_adjs_by_adjoint(ltris)
    save_triangulated(squ1, ltris, ladjs, "hexa1")
    save_triangulated(squ1, ltris, ladjs2, "hexa2")
    fl1 = open("hexa1.adj", "r")
    fl2 = open("hexa2.adj", "r")
    @test read(fl1) == read(fl2)
    rm("hexa1.adj")
    rm("hexa2.adj")
    rm("hexa1.ply")
    rm("hexa2.ply")
    rm("hexa1.tri")
    rm("hexa2.tri")
    rm("squ1.adj")
    rm("squ2.adj")
    rm("squ1.ply")
    rm("squ2.ply")
    rm("squ1.tri")
    rm("squ2.tri")
end


text:
# ---
# title: 895. Maximum Frequency Stack
# id: problem895
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Stack
# link: <https://leetcode.com/problems/maximum-frequency-stack/description/>
# hidden: true
# ---
#
# Implement `FreqStack`, a class which simulates the operation of a stack-like
# data structure.
#
# `FreqStack` has two functions:
#
#   * `push(int x)`, which pushes an integer `x` onto the stack.
#   * `pop()`, which **removes** and returns the most frequent element in the stack.
#     * If there is a tie for most frequent element, the element closest to the top of the stack is removed and returned.
#
#
#
# **Example 1:**
#
#
#
#     Input:
#     ["FreqStack","push","push","push","push","push","push","pop","pop","pop","pop"],
#     [[],[5],[7],[5],[7],[4],[5],[],[],[],[]]
#     Output: [null,null,null,null,null,null,null,5,7,5,4]
#     **Explanation** :
#     After making six .push operations, the stack is [5,7,5,7,4,5] from bottom to top.  Then:
#
#     pop() -> returns 5, as 5 is the most frequent.
#     The stack becomes [5,7,5,7,4].
#
#     pop() -> returns 7, as 5 and 7 is the most frequent, but 7 is closest to the top.
#     The stack becomes [5,7,5,4].
#
#     pop() -> returns 5.
#     The stack becomes [5,7,4].
#
#     pop() -> returns 4.
#     The stack becomes [5,7].
#
#
#
#
# **Note:**
#
#   * Calls to `FreqStack.push(int x)` will be such that `0 <= x <= 10^9`.
#   * It is guaranteed that `FreqStack.pop()` won't be called if the stack has zero elements.
#   * The total number of `FreqStack.push` calls will not exceed `10000` in a single test case.
#   * The total number of `FreqStack.pop` calls will not exceed `10000` in a single test case.
#   * The total number of `FreqStack.push` and `FreqStack.pop` calls will not exceed `150000` across all test cases.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function nfkbOde!(delta, x, v, t)

    # nf-kb part
    IkBa           = x[1]
    IkBan          = x[2]
    IkBaNFkB       = x[3]
    IkBaNFkBn      = x[4]
    IkBat          = x[5]
    IkBb           = x[6]
    IkBbn          = x[7]
    IkBbNFkB       = x[8]
    IkBbNFkBn      = x[9]
    IkBbt          = x[10]
    IkBe           = x[11]
    IkBen          = x[12]
    IkBeNFkB       = x[13]
    IkBeNFkBn      = x[14]
    IkBet          = x[15]
    NFkB           = x[16]
    NFkBn          = x[17]

    # TLR4&IKK part
    LPSo           = x[18]; # input LPS
    LPS            = x[19];
    LPSen          = x[20];
    TLR4           = x[21];
    TLR4en         = x[22];
    TLR4LPS        = x[23]
    TLR4LPSen      = x[24];
    MyD88          = x[25]
    MyD88s         = x[26]
    TRIF           = x[27];
    TRIFs          = x[28];
    TRAF6          = x[29]
    TRAF6s         = x[30]
    IKKK_off       = x[31];
    IKKK           = x[32];
    IKK_off        = x[33];
    IKK            = x[34];
    IKK_i          = x[35];

    # TNF part
    TNFnas        = x[36]
    TNFmRNA       = x[37]
    TNFpro        = x[38]
    TNF           = x[39]

    # TNFR part
    tnfrm          = x[40]
    TNFR           = x[41]
    TNFRtnf        = x[42]
    C1             = x[43]
    C1_off         = x[44]
    C1tnf          = x[45]
    C1tnf_off      = x[46]
    TTR            = x[47]
    # A20
    a20            = x[48]
    a20t           = x[49]

    ## set IKK flux
    IKK_flux = IKK / (v.IKK_TOTAL); # ikk modification

    # TNF mRNA stabilization
    fa = (TRIFs + v.TP[10]) / (TRIFs + v.TP[11])
    fi = v.TP[12] / (TRIFs + v.TP[12])

    kdegm = v.TP[5] * fi
    kpr = v.TP[6] * fa
    ktl = v.TP[7] * fa
    ksec = v.TP[9] * fa

    # golgi plug effect
    kdeg_tnf = v.TP[8]
    ktl_a = v.NP[16]
    ktl_b = v.NP[17]
    ktl_e = v.NP[18]

    # TODO:  is t > 0 nessessary?
    if (v.flag_noTnfFeedback) & (t > 0.0)
        golgi = t / (t + v.TP[14])
        ksec = ksec * (1 - golgi)
        ktnfr = v.AP[2] * (1 - golgi)
        ktlr4 = v.IP[9] * (1 - golgi)
        kshut = v.IP[2] * (1 - golgi);
        kshut_com = v.IP[14] * (1 - golgi)
    else
        ktnfr = v.AP[2]
        ktlr4 = v.IP[9]
        kshut = v.IP[2];
        kshut_com = v.IP[14]
    end

    ## Calculate Iteration Fluxes

    # -- NFkB Activation Module
    f_rsu_a         = v.NP[1]
    f_rsu_b         = v.NP[2]
    f_rsu_e         = v.NP[3]
    f_rsu_f         = v.TP[1] ; # tnf basal transcription [NEW]

    f_rsr_an        = v.NP[4]   * (NFkBn^v.NP[7])
    f_rsr_bn        = v.NP[5]   * (NFkBn^v.NP[8])
    f_rsr_en        = v.NP[6]   * (NFkBn^v.NP[9])

    f_rsr_fn        = v.TP[2]  * (NFkBn^v.TP[3]) / (NFkBn^v.TP[3] + v.TP[4]^v.TP[3]); # tnf transcription
    f_rp_f          = kpr * TNFnas; # tnf RNA processing

    f_rd_a          = v.NP[13]  * IkBat
    f_rd_b          = v.NP[14]  * IkBbt
    f_rd_e          = v.NP[15]  * IkBet
    f_rd_f          = kdegm * TNFmRNA; # tnf mRNA degradation

    f_ps_c_a        = ktl_a	* IkBat
    f_ps_c_b        = ktl_b	* IkBbt
    f_ps_c_e        = ktl_e	* IkBet
    f_ps_c_f        = ktl	* TNFmRNA;  # tnf protein synthsis
    f_pd_f          = kdeg_tnf    * TNFpro;  # tnf protein deg. (inside the cell)
    f_sc_c_f        = ksec   * TNFpro;  # tnf protein secrection

    f_in_a          = v.NP[19]	* IkBa
    f_in_b          = v.NP[20]	* IkBb
    f_in_e          = v.NP[21]	* IkBe
    f_in_n          = v.NP[22]	* NFkB
    f_ex_a          = v.NP[23]	* IkBan
    f_ex_b          = v.NP[24]	* IkBbn
    f_ex_e          = v.NP[25]	* IkBen
    f_ex_n          = v.NP[26]	* NFkBn
    f_in_2an        = v.NP[27]	* IkBaNFkB
    f_in_2bn        = v.NP[28]	* IkBbNFkB
    f_in_2en        = v.NP[29]	* IkBeNFkB
    f_ex_2an        = v.NP[30]	* IkBaNFkBn
    f_ex_2bn        = v.NP[31]	* IkBbNFkBn
    f_ex_2en        = v.NP[32]	* IkBeNFkBn

    f_pd_c_a        = v.NP[33]	* IkBa
    f_pd_c_b        = v.NP[34]	* IkBb
    f_pd_c_e        = v.NP[35]	* IkBe
    f_pd_n_a        = v.NP[36]	* IkBan
    f_pd_n_b        = v.NP[37]	* IkBbn
    f_pd_n_e        = v.NP[38]	* IkBen
    f_pd_c_2an      = v.NP[39]	* IkBaNFkB
    f_pd_c_2bn      = v.NP[40]	* IkBbNFkB
    f_pd_c_2en      = v.NP[41]	* IkBeNFkB
    f_pd_n_2an      = v.NP[42]	* IkBaNFkBn
    f_pd_n_2bn      = v.NP[43]	* IkBbNFkBn
    f_pd_n_2en      = v.NP[44]	* IkBeNFkBn

    f_a_c_an        = v.NP[45]  * IkBa	* NFkB
    f_a_c_bn        = v.NP[46]  * IkBb	* NFkB
    f_a_c_en        = v.NP[47]  * IkBe	* NFkB
    f_a_n_an        = v.NP[48]	* IkBan	* NFkBn
    f_a_n_bn        = v.NP[49]	* IkBbn	* NFkBn
    f_a_n_en        = v.NP[50]  * IkBen	* NFkBn
    f_d_c_an        = v.NP[51]	* IkBaNFkB
    f_d_c_bn        = v.NP[52]	* IkBbNFkB
    f_d_c_en        = v.NP[53]	* IkBeNFkB
    f_d_n_an        = v.NP[54]	* IkBaNFkBn
    f_d_n_bn        = v.NP[55]	* IkBbNFkBn
    f_d_n_en        = v.NP[56]	* IkBeNFkBn

    # IKK Mediated IkB Degradation [free & bound]
    f_ph_c_a        = v.NP[57]  * IkBa * IKK_flux
    f_ph_c_b        = v.NP[58]	* IkBb * IKK_flux
    f_ph_c_e        = v.NP[59]	* IkBe * IKK_flux

    f_ph_c_an       = v.NP[60]  * IkBaNFkB * IKK_flux
    f_ph_c_bn       = v.NP[61]	* IkBbNFkB * IKK_flux
    f_ph_c_en       = v.NP[62]	* IkBeNFkB * IKK_flux

    # A20 Fluxes
    f_rsu_a20       = v.NP[63]
    f_rsr_a20       = v.NP[64]  * (NFkBn^v.NP[65])
    f_rd_a20        = v.NP[66]  * a20t
    f_ps_c_a20      = v.NP[67]  * NFkBn
    f_pd_c_a20      = v.NP[68]  * a20

    if t > v.NP[69]  # disables inducible A20 txn to match exp data
        f_rsr_a20 *= 0.0
    end

    # -- IKK Activation Module & TLR4 module
    # LPS binding to the surface()
    f_b_LPS         = v.IP[1]   * LPSo

    # ligand binding
    f_b_TLR4LPS     = v.IP[5]    * TLR4 * LPS
    f_ub_TLR4LPS    = v.IP[6]    * TLR4LPS
    f_b_TLR4LPSen   = v.IP[7]    * TLR4en * LPSen
    f_ub_TLR4LPSen  = v.IP[8]    * TLR4LPSen

    # Activation modified on 0426
    f_a_MyD88       = v.IP[18]    * (TLR4LPS)^v.IP[19] / ((TLR4LPS)^v.IP[19] + (v.IP[20])^v.IP[19])  * MyD88
    f_i_MyD88       = v.IP[21]    * MyD88s
    f_a_TRIF        = v.IP[22]    * TLR4LPSen * TRIF
    f_i_TRIF        = v.IP[23]   * TRIFs

    # shuttling; modified on 0426
    f_in_LPS        = kshut   * LPS
    f_out_LPS       = v.IP[3]   * LPSen
    f_in_TLR4       = v.IP[12]   * TLR4
    f_out_TLR4      = v.IP[13]   * TLR4en
    f_in_TLR4LPS    = kshut_com   * TLR4LPS
    f_out_TLR4LPS   = v.IP[15]   * TLR4LPSen

    # generation & degradation
    f_g_TLR4        = ktlr4
    f_d_TLR4        = v.IP[10]   * TLR4
    f_d_TLR4en      = v.IP[11]   * TLR4en
    f_d_LPSen       = v.IP[4]   * LPSen
    f_d_TLR4LPSen   = v.IP[17]   * TLR4LPSen
    # modified 0515
    f_d_TLR4LPS     = v.IP[16]   * TLR4LPS
    # -----------------------IKK module
    # TRAF6
    f_a_TRAF6_MyD88s = v.IP[24]   * TRAF6 * MyD88s
    f_a_TRAF6_TRIFs = v.IP[25]   * TRAF6 * TRIFs
    f_i_TRAF6       = v.IP[26]   * TRAF6s

    # IKKK

    f_IKKK_on       = v.IP[28]   * IKKK_off
    f_IKKK_on_TRAF6s = v.IP[27]   * IKKK_off * TRAF6s
    f_IKKK_off      = v.IP[29]   * IKKK

    # IKK
    f_IKK_on        = v.IP[31]   * IKK_off
    f_IKK_on_IKKK   = v.IP[30]   * IKK_off * IKKK
    f_IKK_off       = v.IP[32]   * IKK
    f_IKK_off_i     = v.IP[33]   * IKK
    f_IKKi          = v.IP[34]   * IKK_i

    ## TNFR part
    f_pd_m_tnf      = v.AP[1]    * TNF ; # pd_m_tnf 45' half life of exogenous TNF
    f_syn_tnfrm     = ktnfr;  # tnfrm synthesis [txn, tsl, localization]
    f_pd_tnfrm      = v.AP[3]    * tnfrm;# tnfrm --> deg  -- 120' halflife
    f_a_tnfrm       = v.AP[4]    * tnfrm;# 3tnfrm --> TNFR
    f_d_TNFR        = v.AP[5]    * TNFR;# TNFR   --> 3tnfrm
    f_i_TNFR        = v.AP[6]    * TNFR;# TNFR internalization -- 30' halflife
    f_a_C1_off      = v.AP[7]    * TNFR * TTR;# TNFR + TTR --> C1_off
    f_d_C1_off      = v.AP[8]    * C1_off;# C1_off --> TNFR + TTR
    f_a_C1          = v.AP[9]    * C1_off;# C1_off --> C1
    f_C1_off        = v.AP[10]   * C1;# C1     --> C1_off
    f_C1_A20        = v.AP[11]   * C1 * a20;# C1     --> C1_off [A20 Mediated]
    f_d_C1          = v.AP[12]   * C1; # UPDATE # C1     --> TNFR + TTR
    f_i_C1_off      = v.AP[13]   * C1_off; # C1_off internalization
    f_i_C1          = v.AP[14]   * C1; # C1 internalization
    f_a_tnfrmtnf    = v.AP[15]   * tnfrm * TNF; # 3tnfrm + tnf --> TNFRtnf
    f_a_TNFRtnf     = v.AP[16]   * TNFR * TNF; # TNFR + tnf --> TNFRtnf
    f_d_TNFRtnf     = v.AP[17]   * TNFRtnf; # TNFRtnf   --> TNFR + tnf
    f_i_TNFRtnf     = v.AP[18]   * TNFRtnf;  # TNFRtnf internalization -- 30' halflife

    f_a_C1tnf_off   = v.AP[19]   * TNFRtnf * TTR; # TNFRtnf + TTR --> C1tnf_off
    f_d_C1tnf_off   = v.AP[20]   * C1tnf_off; # C1tnf_off --> TNFRtnf + TTR
    f_a_C1tnf       = v.AP[21]   * C1tnf_off; # C1tnf_off --> C1tnf
    f_C1tnf_off     = v.AP[22]   * C1tnf; # C1tnf     --> C1tnf_off
    f_C1tnf_A20     = v.AP[23]   * C1tnf * a20;  # C1tnf     --> C1tnf_off [A20 Mediated]
    f_d_C1tnf       = v.AP[24]   * C1tnf; # updated  # C1tnf     --> TNFRtnf + TTR
    f_i_C1tnf       = v.AP[25]   * C1tnf; # C1tnf internalization
    f_i_C1tnf_off   = v.AP[26]   * C1tnf_off;  # C1tnf_off internalization
    f_d_tnf_C1_off  = v.AP[27]   * C1tnf_off; # C1tnf_off --> C1_off + tnf
    f_a_tnf_C1_off  = v.AP[28]   * C1_off * TNF; # C1_off + tnf --> C1tnf_off
    f_d_tnf_C1      = v.AP[29]   * C1tnf; # C1tnf    --> C1 + tnf
    f_a_tnf_C1      = v.AP[30]   * C1 * TNF; # C1 + tnf --> C1tnf
    f_IKKK_on_C1    = v.AP[31]   * IKKK_off * C1; # IKKK_off --> IKKK [C1 mediated]?500
    f_IKKK_on_C1tnf = v.AP[32]   * IKKK_off * C1tnf;# IKKK_off --> IKKK [C1tnf mediated]

    ## Add Fluxes to component concentrations; then Save

    # IkB Transcrv.IPtion
    d_IkBat     = d_IkBat + f_rsu_a
    d_IkBbt     = d_IkBbt + f_rsu_b
    d_IkBet     = d_IkBet + f_rsu_e
    d_IkBat     = d_IkBat + f_rsr_an
    d_IkBbt     = d_IkBbt + f_rsr_bn
    d_IkBet     = d_IkBet + f_rsr_en
    d_TNFnas    = d_TNFnas + f_rsu_f; # NEW
    d_TNFnas    = d_TNFnas + f_rsr_fn; # NEW
    d_TNFnas    = d_TNFnas - f_rp_f; # NEW
    d_TNFmRNA   = d_TNFmRNA + f_rp_f; # NEW

    # IkB Transcrv.IPt Degradation
    d_IkBbt     = d_IkBbt - f_rd_b
    d_IkBet     = d_IkBet - f_rd_e
    d_IkBat     = d_IkBat - f_rd_a
    d_TNFmRNA   = d_TNFmRNA  - f_rd_f;   # tnft degradation

    # IkB Translation
    d_IkBat     = d_IkBat - f_ps_c_a
    d_IkBat     = d_IkBat + f_ps_c_a
    d_IkBa      = + f_ps_c_a

    d_IkBbt     = d_IkBbt - f_ps_c_b
    d_IkBbt     = d_IkBbt + f_ps_c_b
    d_IkBb      = d_IkBb  + f_ps_c_b

    d_IkBet     = d_IkBet - f_ps_c_e
    d_IkBet     = d_IkBet + f_ps_c_e
    d_IkBe      = d_IkBe  + f_ps_c_e

    d_TNFpro    = d_TNFpro + f_ps_c_f;

    # TNF secretion + protein degradation
    d_TNFpro     = d_TNFpro - f_sc_c_f;
    d_TNF        = d_TNF    + f_sc_c_f;
    d_TNFpro     = d_TNFpro - f_pd_f;

    # IkB:NFkB Shuttling [Free & Bound]
    d_IkBa      = d_IkBa  - f_in_a
    d_IkBan     = + f_in_a

    d_IkBb      = d_IkBb  - f_in_b
    d_IkBbn     = d_IkBbn + f_in_b

    d_IkBe      = d_IkBe  - f_in_e
    d_IkBen     = d_IkBen + f_in_e

    d_NFkB      = d_NFkB  - f_in_n
    d_NFkBn     = d_NFkBn + f_in_n

    d_IkBan     = d_IkBan - f_ex_a
    d_IkBa      = d_IkBa  + f_ex_a

    d_IkBbn     = d_IkBbn - f_ex_b
    d_IkBb      = d_IkBb  + f_ex_b

    d_IkBen     = d_IkBen - f_ex_e
    d_IkBe      = d_IkBe  + f_ex_e

    d_NFkBn     = d_NFkBn - f_ex_n
    d_NFkB      = d_NFkB  + f_ex_n

    d_IkBaNFkB  = - f_in_2an
    d_IkBaNFkBn = d_IkBaNFkBn + f_in_2an

    d_IkBbNFkB  = d_IkBbNFkB  - f_in_2bn
    d_IkBbNFkBn = d_IkBbNFkBn + f_in_2bn

    d_IkBeNFkB  = d_IkBeNFkB  - f_in_2en
    d_IkBeNFkBn = d_IkBeNFkBn + f_in_2en

    d_IkBaNFkBn = d_IkBaNFkBn - f_ex_2an
    d_IkBaNFkB  = d_IkBaNFkB  + f_ex_2an

    d_IkBbNFkBn = d_IkBbNFkBn - f_ex_2bn
    d_IkBbNFkB  = d_IkBbNFkB  + f_ex_2bn

    d_IkBeNFkBn = d_IkBeNFkBn - f_ex_2en
    d_IkBeNFkB  = d_IkBeNFkB  + f_ex_2en

    # IkB:NFkB Association [Cytoplasm & Nucleus]
    d_IkBa      = d_IkBa - f_a_c_an
    d_NFkB      = d_NFkB - f_a_c_an
    d_IkBaNFkB  = d_IkBaNFkB + f_a_c_an

    d_IkBb      = d_IkBb - f_a_c_bn
    d_NFkB      = d_NFkB - f_a_c_bn
    d_IkBbNFkB  = d_IkBbNFkB + f_a_c_bn

    d_IkBe      = d_IkBe - f_a_c_en
    d_NFkB      = d_NFkB - f_a_c_en
    d_IkBeNFkB  = d_IkBeNFkB + f_a_c_en

    d_IkBan     = d_IkBan - f_a_n_an
    d_NFkBn     = d_NFkBn - f_a_n_an
    d_IkBaNFkBn = d_IkBaNFkBn + f_a_n_an

    d_IkBbn     = d_IkBbn - f_a_n_bn
    d_NFkBn     = d_NFkBn - f_a_n_bn
    d_IkBbNFkBn = d_IkBbNFkBn + f_a_n_bn

    d_IkBen     = d_IkBen - f_a_n_en
    d_NFkBn     = d_NFkBn - f_a_n_en
    d_IkBeNFkBn = d_IkBeNFkBn + f_a_n_en

    # IkB:NFkB Dissociation [Cytoplasm & Nucleus]
    d_IkBaNFkB  = d_IkBaNFkB - f_d_c_an
    d_IkBa      = d_IkBa + f_d_c_an
    d_NFkB      = d_NFkB + f_d_c_an

    d_IkBbNFkB  = d_IkBbNFkB - f_d_c_bn
    d_IkBb      = d_IkBb + f_d_c_bn
    d_NFkB      = d_NFkB + f_d_c_bn

    d_IkBeNFkB  = d_IkBeNFkB - f_d_c_en
    d_IkBe      = d_IkBe + f_d_c_en
    d_NFkB      = d_NFkB + f_d_c_en

    d_IkBaNFkBn = d_IkBaNFkBn - f_d_n_an
    d_IkBan     = d_IkBan + f_d_n_an
    d_NFkBn     = d_NFkBn + f_d_n_an

    d_IkBbNFkBn = d_IkBbNFkBn - f_d_n_bn
    d_IkBbn     = d_IkBbn + f_d_n_bn
    d_NFkBn     = d_NFkBn + f_d_n_bn

    d_IkBeNFkBn = d_IkBeNFkBn - f_d_n_en
    d_IkBen     = d_IkBen + f_d_n_en
    d_NFkBn     = d_NFkBn + f_d_n_en

    # Free IkB Degradation [Cytoplasm & Nucleus]
    d_IkBa      = d_IkBa  - f_pd_c_a
    d_IkBb      = d_IkBb  - f_pd_c_b
    d_IkBe      = d_IkBe  - f_pd_c_e
    d_IkBan     = d_IkBan - f_pd_n_a
    d_IkBbn     = d_IkBbn - f_pd_n_b
    d_IkBen     = d_IkBen - f_pd_n_e

    # IkB:NFkB Degradation [Cytoplasm & Nucleus]
    d_IkBaNFkB  = d_IkBaNFkB - f_pd_c_2an
    d_NFkB      = d_NFkB + f_pd_c_2an

    d_IkBbNFkB  = d_IkBbNFkB - f_pd_c_2bn
    d_NFkB      = d_NFkB + f_pd_c_2bn

    d_IkBeNFkB  = d_IkBeNFkB - f_pd_c_2en
    d_NFkB      = d_NFkB + f_pd_c_2en

    d_IkBaNFkBn = d_IkBaNFkBn - f_pd_n_2an
    d_NFkBn     = d_NFkBn + f_pd_n_2an

    d_IkBbNFkBn = d_IkBbNFkBn - f_pd_n_2bn
    d_NFkBn     = d_NFkBn + f_pd_n_2bn

    d_IkBeNFkBn = d_IkBeNFkBn - f_pd_n_2en
    d_NFkBn     = d_NFkBn + f_pd_n_2en

    # IKK Mediated IkB Degradation
    d_IkBa      = d_IkBa - f_ph_c_a
    d_IkBb      = d_IkBb - f_ph_c_b
    d_IkBe      = d_IkBe - f_ph_c_e

    # IKK Mediated IkB:NFkB Degradation
    d_IkBaNFkB  = d_IkBaNFkB - f_ph_c_an
    d_NFkB      = d_NFkB + f_ph_c_an

    d_IkBbNFkB  = d_IkBbNFkB - f_ph_c_bn
    d_NFkB      = d_NFkB + f_ph_c_bn

    d_IkBeNFkB  = d_IkBeNFkB - f_ph_c_en
    d_NFkB      = d_NFkB + f_ph_c_en

    # Upstream Pathway  [LPS] --> IKK activation

    # Ligand binding
    d_LPSo      = d_LPSo - f_b_LPS
    d_LPS       = d_LPS  + f_b_LPS

    d_LPS       = d_LPS       -   f_b_TLR4LPS
    d_TLR4       = d_TLR4       -   f_b_TLR4LPS
    d_TLR4LPS   = d_TLR4LPS   +   f_b_TLR4LPS
    d_LPS       = d_LPS       +   f_ub_TLR4LPS
    d_TLR4       = d_TLR4     +   f_ub_TLR4LPS
    d_TLR4LPS   = d_TLR4LPS   -   f_ub_TLR4LPS

    d_LPSen     = d_LPSen     -   f_b_TLR4LPSen
    d_TLR4en     = d_TLR4en     -   f_b_TLR4LPSen
    d_TLR4LPSen = d_TLR4LPSen +   f_b_TLR4LPSen
    d_LPSen     = d_LPSen     +   f_ub_TLR4LPSen
    d_TLR4en     = d_TLR4en   +   f_ub_TLR4LPSen
    d_TLR4LPSen = d_TLR4LPSen -   f_ub_TLR4LPSen

    # Activation
    d_MyD88     = d_MyD88     -   f_a_MyD88
    d_MyD88s    = d_MyD88s    +   f_a_MyD88
    d_MyD88     = d_MyD88     +   f_i_MyD88
    d_MyD88s    = d_MyD88s    -   f_i_MyD88

    d_TRIF      = d_TRIF      -   f_a_TRIF
    d_TRIFs     = d_TRIFs     +   f_a_TRIF
    d_TRIF      = d_TRIF      +   f_i_TRIF
    d_TRIFs     = d_TRIFs     -   f_i_TRIF

    # shuttling
    d_LPS       = d_LPS       -   f_in_LPS
    d_LPSen     = d_LPSen     +   f_in_LPS
    # d_LPSo       = d_LPSo     +   f_out_LPS
    d_LPS       = d_LPS       +   f_out_LPS
    d_LPSen     = d_LPSen     -   f_out_LPS

    d_TLR4      = d_TLR4      -   f_in_TLR4
    d_TLR4en    = d_TLR4en    +   f_in_TLR4
    d_TLR4      = d_TLR4      +   f_out_TLR4
    d_TLR4en    = d_TLR4en    -   f_out_TLR4

    d_TLR4LPS      = d_TLR4LPS      -   f_in_TLR4LPS
    d_TLR4LPSen    = d_TLR4LPSen    +   f_in_TLR4LPS
    d_TLR4LPS      = d_TLR4LPS      +   f_out_TLR4LPS
    d_TLR4LPSen     = d_TLR4LPSen   -   f_out_TLR4LPS

    # generation & degradation
    d_TLR4      = d_TLR4 + f_g_TLR4
    d_TLR4      = d_TLR4 - f_d_TLR4
    d_TLR4en     = d_TLR4en - f_d_TLR4en
    d_LPSen      = d_LPSen - f_d_LPSen
    d_TLR4LPSen      = d_TLR4LPSen - f_d_TLR4LPSen
    # -- modified 0515
    d_TLR4LPS      = d_TLR4LPS - f_d_TLR4LPS

    # TRAF6
    d_TRAF6       = d_TRAF6  - f_a_TRAF6_MyD88s
    d_TRAF6s      = d_TRAF6s + f_a_TRAF6_MyD88s
    d_TRAF6       = d_TRAF6  - f_a_TRAF6_TRIFs
    d_TRAF6s      = d_TRAF6s + f_a_TRAF6_TRIFs
    d_TRAF6       = d_TRAF6  + f_i_TRAF6
    d_TRAF6s      = d_TRAF6s - f_i_TRAF6

    # IKKK Regulation
    d_IKKK      = d_IKKK     + f_IKKK_on
    d_IKKK_off  = d_IKKK_off - f_IKKK_on

    d_IKKK      = d_IKKK     + f_IKKK_on_TRAF6s
    d_IKKK_off  = d_IKKK_off - f_IKKK_on_TRAF6s

    d_IKKK_off  = d_IKKK_off + f_IKKK_off
    d_IKKK      = d_IKKK     - f_IKKK_off

    # IKK Regulation
    d_IKK       = d_IKK     + f_IKK_on
    d_IKK_off   = d_IKK_off - f_IKK_on

    d_IKK       = d_IKK     + f_IKK_on_IKKK
    d_IKK_off   = d_IKK_off - f_IKK_on_IKKK

    d_IKK       = d_IKK     - f_IKK_off
    d_IKK_off   = d_IKK_off + f_IKK_off

    d_IKK       = d_IKK     - f_IKK_off_i
    d_IKK_i     = d_IKK_i   + f_IKK_off_i

    d_IKK_off   = d_IKK_off + f_IKKi
    d_IKK_i     = d_IKK_i   - f_IKKi

    ## TNR module [NEW] + A20 [NEW]
    # A20 Expression Transcrviption
    d_a20t      = d_a20t + f_rsu_a20
    d_a20t      = d_a20t + f_rsr_a20
    d_a20t      = d_a20t - f_rd_a20

    # Translation
    d_a20       = d_a20  + f_ps_c_a20

    # Degradation
    d_a20       = d_a20 - f_pd_c_a20

    # Upstream Pathway  [TNF] --> IKK activation

    # TNF-independent Activation

    # tnfrm trimerization 3tnfrm<--> TNFR
    d_tnfrm     = d_tnfrm - 3 * f_a_tnfrm
    d_TNFR      = d_TNFR +   f_a_tnfrm

    d_tnfrm     = d_tnfrm + 3 * f_d_TNFR
    d_TNFR      = d_TNFR -   f_d_TNFR

    # TNF Receptor Metabolism
    d_TNFR     = d_TNFR     - f_i_TNFR
    d_TNFRtnf  = d_TNFRtnf  - f_i_TNFRtnf

    d_tnfrm    = d_tnfrm    - f_pd_tnfrm
    d_tnfrm    = d_tnfrm    + f_syn_tnfrm

    # Complex I Generation  TTR + TNFR <--> C1_off; C1
    d_C1_off    = d_C1_off  + f_a_C1_off
    d_TTR       = d_TTR     - f_a_C1_off
    d_TNFR      = d_TNFR    - f_a_C1_off

    d_C1_off    = d_C1_off  - f_d_C1_off
    d_TTR       = d_TTR     + f_d_C1_off
    d_TNFR      = d_TNFR    + f_d_C1_off

    d_C1        = d_C1      - f_d_C1
    d_TTR       = d_TTR     + f_d_C1
    d_TNFR      = d_TNFR    + f_d_C1

    # C1_off  <--> C1
    d_C1_off    = d_C1_off  - f_a_C1
    d_C1        = d_C1      + f_a_C1

    d_C1_off    = d_C1_off  + f_C1_A20
    d_C1        = d_C1      - f_C1_A20

    d_C1_off    = d_C1_off  + f_C1_off
    d_C1        = d_C1      - f_C1_off

    # C1 & C1_off internalization [treated as deg]
    d_C1        = d_C1      - f_i_C1
    d_C1_off    = d_C1_off  - f_i_C1_off

    # TNF-Dependent Activation

    # TNF degradation
    d_TNF       = d_TNF     - f_pd_m_tnf

    # TNF:TNFR binding
    d_TNFRtnf      = d_TNFRtnf    +   f_a_tnfrmtnf
    d_tnfrm     = d_tnfrm   - 3 * f_a_tnfrmtnf
    d_TNF       = d_TNF     -   f_a_tnfrmtnf

    d_TNFR      = d_TNFR    - f_a_TNFRtnf
    d_TNF       = d_TNF     - f_a_TNFRtnf
    d_TNFRtnf   = d_TNFRtnf + f_a_TNFRtnf

    d_TNFR      = d_TNFR    + f_d_TNFRtnf
    d_TNF       = d_TNF     + f_d_TNFRtnf
    d_TNFRtnf   = d_TNFRtnf - f_d_TNFRtnf

    # Complex I Generation  TTR + TNFRtnf <--> C1tnf_off; C1tnf
    d_C1tnf_off = d_C1tnf_off   + f_a_C1tnf_off
    d_TTR       = d_TTR         - f_a_C1tnf_off
    d_TNFRtnf   = d_TNFRtnf     - f_a_C1tnf_off

    d_C1tnf_off = d_C1tnf_off   - f_d_C1tnf_off
    d_TTR       = d_TTR         + f_d_C1tnf_off
    d_TNFRtnf   = d_TNFRtnf     + f_d_C1tnf_off

    d_C1tnf     = d_C1tnf       - f_d_C1tnf
    d_TTR       = d_TTR         + f_d_C1tnf
    d_TNFRtnf   = d_TNFRtnf     + f_d_C1tnf

    # C1tnf_off <--> C1tnf
    d_C1tnf_off = d_C1tnf_off   - f_a_C1tnf
    d_C1tnf     = d_C1tnf       + f_a_C1tnf

    d_C1tnf_off = d_C1tnf_off   + f_C1tnf_off
    d_C1tnf     = d_C1tnf       - f_C1tnf_off

    d_C1tnf_off = d_C1tnf_off   + f_C1tnf_A20
    d_C1tnf     = d_C1tnf       - f_C1tnf_A20

    # C1tnf & C1tnf_off internalization [treated as deg]
    d_C1tnf_off = d_C1tnf_off   - f_i_C1tnf_off
    d_C1tnf     = d_C1tnf       - f_i_C1tnf

    # C1tnf_off;C1tnf <--> C1_off;C1 + TNF
    d_C1tnf_off = d_C1tnf_off   - f_d_tnf_C1_off
    d_C1_off    = d_C1_off      + f_d_tnf_C1_off
    d_TNF       = d_TNF         + f_d_tnf_C1_off

    d_C1tnf_off = d_C1tnf_off   + f_a_tnf_C1_off
    d_C1_off    = d_C1_off      - f_a_tnf_C1_off
    d_TNF       = d_TNF         - f_a_tnf_C1_off

    d_C1tnf     = d_C1tnf       - f_d_tnf_C1
    d_C1        = d_C1          + f_d_tnf_C1
    d_TNF       = d_TNF         + f_d_tnf_C1

    d_C1tnf     = d_C1tnf       + f_a_tnf_C1
    d_C1        = d_C1          - f_a_tnf_C1
    d_TNF       = d_TNF         - f_a_tnf_C1

    # IKKK Regulation

    d_IKKK      = d_IKKK     + f_IKKK_on_C1
    d_IKKK_off  = d_IKKK_off - f_IKKK_on_C1

    d_IKKK      = d_IKKK     + f_IKKK_on_C1tnf
    d_IKKK_off  = d_IKKK_off - f_IKKK_on_C1tnf

    ## ADDED
    # A20 mediated TRAF6 inactivation
    d_TRAF6    = d_TRAF6    + v.TP[15] * a20 * TRAF6s
    d_TRAF6s    = d_TRAF6s  - v.TP[15] * a20 * TRAF6s

    ## Save concentrations for next time step
    delta[1]      = d_IkBa
    delta[2]      = d_IkBan
    delta[3]      = d_IkBaNFkB
    delta[4]      = d_IkBaNFkBn
    delta[5]      = d_IkBat
    delta[6]      = d_IkBb
    delta[7]      = d_IkBbn
    delta[8]      = d_IkBbNFkB
    delta[9]      = d_IkBbNFkBn
    delta[10]     = d_IkBbt
    delta[11]     = d_IkBe
    delta[12]     = d_IkBen
    delta[13]     = d_IkBeNFkB
    delta[14]     = d_IkBeNFkBn
    delta[15]     = d_IkBet

    delta[16]     = d_NFkB
    delta[17]     = d_NFkBn
    #
    delta[18]     = d_LPSo
    delta[19]     = d_LPS
    delta[20]     = d_LPSen
    delta[21]     = d_TLR4
    delta[22]     = d_TLR4en
    delta[23]     = d_TLR4LPS
    delta[24]     = d_TLR4LPSen
    delta[25]     = d_MyD88
    delta[26]     = d_MyD88s
    delta[27]     = d_TRIF
    delta[28]     = d_TRIFs
    delta[29]     = d_TRAF6
    delta[30]     = d_TRAF6s
    delta[31]     = d_IKKK_off
    delta[32]     = d_IKKK
    delta[33]     = d_IKK_off
    delta[34]     = d_IKK
    delta[35]     = d_IKK_i

    # TNF component
    delta[36]     = d_TNFnas
    delta[37]     = d_TNFmRNA
    delta[38]     = d_TNFpro
    delta[39]     = d_TNF

    # TNFR part
    delta[40]     = d_tnfrm
    delta[41]     = d_TNFR
    delta[42]     = d_TNFRtnf
    delta[43]     = d_C1
    delta[44]     = d_C1_off
    delta[45]     = d_C1tnf
    delta[46]     = d_C1tnf_off
    delta[47]     = d_TTR

    # A20
    delta[48]     = d_a20
    delta[49]     = d_a20t

    delta
end

max_stars_repo_name: jiweiqi/NFkB

max_issues_repo_name: jiweiqi/NFkB

max_forks_repo_name: jiweiqi/NFkB


text:
# horizontal or vertical array of widgets
# nestable. But the top-most List handles the canvas navigation

# Every time a nesting layer is added, all the children's windows have to be redone,
# to reflect their locations in the main pad.
# This is bad!

using Lazy

function newTwList( scr::TwObj;
        height::Real = 25, width::Real = 80,
        posy::Any = :center, posx::Any = :center,
        canvasheight = 80,
        canvaswidth = 128,
        box=true,
        horizontal=false,
        title=utf8(""),
        showLineInfo=true)
    obj = TwObj( TwListData(), Val{:List } )
    obj.box = box
    obj.title = title
    obj.borderSizeV= box ? 1 : 0
    obj.borderSizeH= box ? 1 : 0
    obj.data.horizontal = horizontal
    obj.data.showLineInfo = showLineInfo
    obj.data.canvasheight = canvasheight
    obj.data.canvaswidth = canvaswidth

    link_parent_child( scr, obj, height,width,posy,posx )
    if objtype( scr ) == :Screen
        obj.data.pad = newpad( obj.data.canvasheight, obj.data.canvaswidth )
    end
    obj
end

# move a fully formed widget into this list. need more bookkeeping
function push_widget!( o::TwObj{TwListData}, w::TwObj )
    global rootwin
    # change the widget's window to reflect its location on the canvas
    begx = 0
    begy = 0
    if o.data.horizontal
        for sw in o.data.widgets
            begx += sw.width
        end
    else
        for sw in o.data.widgets
            begy += sw.height
        end
    end

    # This widget must have been previously registered to screen
    unregisterTwObj( o.screen.value, w )
    if typeof( w.window ) <: Ptr && w.window != rootwin
        delwin( w.window ) # so we don't leak memory
    end

    # by the time a list is being added, its contents must be fully populated
    if objtype( w ) == :List
        if typeof( w.data.pad ) <: Ptr && w.data.pad != o.data.pad # the list has its own pad, get rid of it.
            delwin( w.data.pad )
        end
        update_list_canvas( w )
        w.height = w.data.canvasheight
        w.width = w.data.canvaswidth
        w.data.pad = nothing
    end

    w.window = TwWindow( WeakRef( o ), begy, begx, w.height, w.width )
    push!( o.data.widgets, w )
    w.hasFocus = false
end

function update_list_canvas( o::TwObj{TwListData} )
    ws = o.data.widgets
    if isempty( ws )
        # TODO: reconsider this
        o.data.canvasheight = 80
        o.data.canvaswidth = 128
    else
        # TODO: is this necessary?
        for w in o.data.widgets
            if objtype( w ) == :List
                update_list_canvas(w)
            end
        end
        if o.data.horizontal
            o.data.canvasheight = maximum( map( _->objtype(_)==:List? _.data.canvasheight : _.height, ws ) )
            o.data.canvaswidth = sum( map( _->objtype(_)==:List? _.data.canvaswidth : _.width, ws ) )
        else
            o.data.canvasheight = sum( map( _->objtype(_)==:List? _.data.canvasheight: _.height, ws ) )
            o.data.canvaswidth = maximum( map( _->objtype(_)==:List? _.data.canvaswidth : _.width, ws ) )
        end
        if !(typeof( o.window ) <: Ptr)
            o.height = o.data.canvasheight + (o.box?2:0)
            o.width = o.data.canvaswidth + (o.box?2:0)
        end
    end
    if o.data.pad != nothing
        delwin( o.data.pad )
        o.data.pad = newpad( o.data.canvasheight, o.data.canvaswidth )
    end
end

function draw( o::TwObj{TwListData} )
    werase( o.window ) # this is important, or attributes on the pad may be lost

    if typeof( o.window ) <: Ptr
        set_default_focus( o )
        viewContentHeight = o.height - 2*o.borderSizeV
        viewContentWidth  = o.width - 2*o.borderSizeH
        if o.box
            if o.data.navigationmode
                wattron( o.window, COLOR_PAIR( 12 ) )
            end

            box( o.window, 0, 0 )

            if o.data.showLineInfo
                if o.data.canvasheight <= viewContentHeight
                    vscale = "v:all"
                else
                    vscale = @sprintf( "v:%d/%d", o.data.canvaslocy, o.data.canvasheight )
                end

                if o.data.canvaswidth <= viewContentWidth
                    hscale = "h:all"
                else
                    hscale = @sprintf( "h:%d/%d", o.data.canvaslocx, o.data.canvaswidth )
                end

                msg = vscale * " " * hscale
                mvwprintw( o.window, 0, o.width - length( msg ) - 3, "%s", msg )
            end

            if o.data.bottomText != ""
                mvwprintw( o.window, o.height-1, 2, "%s", o.data.bottomText )
            end

            if o.data.navigationmode
                mvwprintw( o.window, 0, 2, "%s", "Navigation mode" )
                wattroff( o.window, COLOR_PAIR( 12 ) )
            end
        end
    end

    for w in o.data.widgets
        # TODO: no need to draw widget outside visible range? or just draw everything?
        if w.isVisible
            draw( w )
        end
    end
    # handle pushing the canvas to the screen if it is the root list
    if typeof( o.window ) <: Ptr
        #TODO: how much of the canvas are we showing?
        borderSizeH = o.box ? 1 : 0
        borderSizeV = o.box ? 1 : 0
        contentheight = o.height - borderSizeV*2
        contentwidth  = o.width  - borderSizeH*2
        canvasheight = o.data.canvasheight
        canvaswidth  = o.data.canvaswidth
        copywin( o.data.pad, o.window, o.data.canvaslocy, o.data.canvaslocx,
            borderSizeV, borderSizeH,
            min( contentheight, canvasheight - o.data.canvaslocy),
            min( contentwidth , canvaswidth  - o.data.canvaslocx) )
    end
end

function lowest_widget( o::TwObj{ TwListData } )
    if o.data.focus == 0
        error( "cannot locate focused widget")
    end
    w = o.data.widgets[ o.data.focus ]
    if objtype( w ) == :List
        return lowest_widget( w )
    else
        return w
    end
end

function lowest_widget_location_area( o::TwObj{ TwListData }, y::Int=0, x::Int=0 )
    if o.data.focus == 0
        error( "cannot locate focused widget")
    end
    w = o.data.widgets[ o.data.focus ]
    if objtype( w ) == :List
        return lowest_widget_location_area( w, y+w.window.yloc, x+w.window.xloc )
    else
        return ( w, y+w.window.yloc, x+w.window.xloc, w.window.height, w.window.width )
    end
end

function ensure_visible_on_canvas( o::TwObj )
    h = o.height
    w = o.width
    y = o.window.yloc
    x = o.window.xloc
    log( @sprintf( "ensure %s is visible", string( o ) ) )
    log( @sprintf( "  init local y,x: %d %d", y, x ) )
    win = o.window
    par = win.parent.value
    while( !(typeof( win.parent.value.window) <: Ptr) )
        y += win.parent.value.window.yloc
        x += win.parent.value.window.xloc
        par = win.parent.value
        win = win.parent.value.window
    end
    par = win.parent.value
    log( @sprintf( "  actual coord y,x: %d %d", y, x ) )
    @assert objtype( par ) == :List
    contentwidth = par.width - (par.box?2:0)
    contentheight = par.height - (par.box?2:0)
    log( @sprintf( "  canvas size     : %d %d", par.data.canvasheight, par.data.canvaswidth ) )
    log( @sprintf( "  window geom     : %d %d", contentheight, contentwidth ) )
    log( @sprintf( "  canvas wind.orig: %d %d", par.data.canvaslocy, par.data.canvaslocx) )

    if par.data.canvaslocx > x
        par.data.canvaslocx = x
    end
    if par.data.canvaslocy > y
        par.data.canvaslocy = y
    end
    if x + w - par.data.canvaslocx > contentwidth
        par.data.canvaslocx = max( 0, x + w - contentwidth )
    end
    if y + h - par.data.canvaslocy > contentheight
        par.data.canvaslocy = max( 0, y + h - contentheight )
    end
    if par.data.canvaslocx > par.data.canvaswidth - contentwidth
        par.data.canvaslocx = max(0,par.data.canvaswidth - contentwidth)
    end
    if par.data.canvaslocy > par.data.canvasheight - contentheight
        par.data.canvaslocy = max(0,par.data.canvasheight- contentheight)
    end
    log( @sprintf( "  canvas  new orig: %d %d", par.data.canvaslocy, par.data.canvaslocx) )
end

function inject( o::TwObj{TwListData}, token::Any )
    retcode = :pass
    dorefresh = false
    isrootlist = typeof( o.window ) <: Ptr
    focus = o.data.focus
    if focus == 0
        return :pass
    end

    if token == :esc
        return :exit_nothing
    end

    if !o.data.navigationmode
        result = inject( o.data.widgets[ focus], token )
        if result != :pass
            refresh(o)
            return result
        end
    end

    function check_accept_focus(w::TwObj,stepsign::Int)
        if w.isVisible && w.acceptsFocus
            if objtype(w) == :List
                r = 1:length(w.data.widgets)
                if stepsign == -1
                    r = reverse(r)
                end
                for i in r
                    if check_accept_focus(w.data.widgets[i], stepsign )
                        return true
                    end
                end
                return false
            else
                return true
            end
        end
        return false
    end

    # TODO: what's the behavior of :esc
    # TODO: what's the behavior of :exit_ok
    if !o.data.navigationmode && token in [ :tab, :shift_tab ]
        prevw = o.data.widgets[focus]
        # note that if the widget is a list and can take a tab/shift tab as, we
        # wouldn't be here in the first place
        dorefresh = true
        if token == :tab
            stp = 1
        else
            stp = -1
        end
        if isrootlist # wrap around
            i = mod1(focus + stp,length(o.data.widgets) )
            while (i != focus)
                w = o.data.widgets[i]
                if check_accept_focus(w,stp)
                    deep_unfocus( prevw )
                    deep_focus( w, stp == -1 ) # 2nd arg is reverse
                    retcode = :got_it
                    break
                end
                i = mod1(i+ stp,length(o.data.widgets))
            end
        else
            i = focus + stp
            if stp == 1
                r = (focus+stp):length(o.data.widgets)
            else
                r = (focus+stp):-1:1
            end
            for i in r
                w = o.data.widgets[i]
                if check_accept_focus(w,stp)
                    deep_unfocus( prevw )
                    deep_focus( w, stp == -1 ) # 2nd arg is reverse
                    retcode = :got_it
                    break
                end
            end
        end
    elseif !o.data.navigationmode && token in [ :left, :right, :up, :down,
        :ctrl_left, :ctrl_right, :ctrl_up, :ctrl_down ] && isrootlist
        # find the location of the current focus
        (w, yloc,xloc, height, width) = lowest_widget_location_area( o )
        if token in [ :up, :ctrl_up ]
            distfunc = function( to::Tuple{Int,Int,Int,Int} )
                updown_arrow_distance( to, (yloc,xloc,height,width), 1 )
            end
        elseif token in [ :down, :ctrl_down ]
            distfunc = function( to::Tuple{Int,Int,Int,Int} )
                updown_arrow_distance( to, (yloc,xloc,height,width), -1 )
            end
        elseif token in [ :left, :ctrl_left ]
            distfunc = function( to::Tuple{Int,Int,Int,Int} )
                leftright_arrow_distance( to, (yloc,xloc,height,width), 1 )
            end
        else
            distfunc = function( to::Tuple{Int,Int,Int,Int} )
                leftright_arrow_distance( to, (yloc,xloc,height,width), -1 )
            end
        end
        wdists = Any[]
        geometric_filter( o, distfunc, 0, 0, wdists, true, 40 )
        candidate = nothing
        mindist = 999999
        for (cw,dist) in wdists
            if dist < mindist
                candidate = cw
                mindist = dist
            end
        end
        if candidate != nothing
            deep_unfocus( w )
            deep_focus( candidate )
            dorefresh = true
            retcode = :got_it
        end
    elseif o.data.navigationmode && token in [ :left, :right]
        if token == :left
            xstep = -1
        elseif token == :right
            xstep = 1
        end
        newlocx = o.data.canvaslocx + xstep
        if newlocx < 0
            o.data.canvaslocx = 0
        elseif newlocx + o.width - o.borderSizeH*2 >= o.data.canvaswidth
            o.data.canvaslocx = max(0, o.data.canvaswidth - o.width + o.borderSizeH*2 )
        else
            o.data.canvaslocx = newlocx
        end
        dorefresh = true
    elseif o.data.navigationmode && token in [ :up, :down, :pageup, :pagedown ]
        # the canvas location can move by one screen-size in either direction
        # but it's kept in by the maximum canvas sizes
        # also, focus would change to the new one closest to the current focused widget
        if token == :up
            ystep = -1
        elseif token == :down
            ystep = 1
        elseif token == :pageup
            ystep = -( o.height - o.borderSizeV*2 )
        elseif token == :pagedown
            ystep = ( o.height - o.borderSizeV*2 )
        end
        newlocy = o.data.canvaslocy + ystep
        if newlocy < 0
            o.data.canvaslocy = 0
        elseif newlocy + o.height - o.borderSizeV*2 >= o.data.canvasheight
            o.data.canvaslocy = max(0, o.data.canvasheight - o.height+ o.borderSizeV*2 )
        else
            o.data.canvaslocy = newlocy
        end
        dorefresh = true
    elseif token == :KEY_MOUSE && isrootlist
        (mstate, x, y, bs ) = getmouse()
        if mstate == :button1_pressed
            (rely, relx) = screen_to_relative( o.window, y, x )
            if 0<=relx<o.width && 0<=rely<o.height
                rely -= o.borderSizeV
                relx -= o.borderSizeH
                # find the closest widget
                distfunc = function( to::Tuple{Int,Int,Int,Int} )
                    point_from_area( rely, relx, to )
                end
                wdists = Any[]
                geometric_filter( o, distfunc, 0, 0, wdists, false, 40 )
                candidate = nothing
                mindist = 999999
                for (cw,dist) in wdists
                    if dist < mindist
                        candidate = cw
                        mindist = dist
                    end
                end
                if candidate != nothing
                    w = lowest_widget( o )
                    deep_unfocus( w )
                    deep_focus( candidate )
                    dorefresh = true
                    retcode = :got_it
                    o.data.navigationmode = false
                end
            else
                retcode = :pass
            end
        end
    elseif token == :ctrl_F4 && isrootlist
        o.data.navigationmode = !o.data.navigationmode
        # if no longer in navigationmode
        # switch focus to the closest visible widget (% overlap with visible window makes a big difference)
        # partially visible widgets suffer a substantial penalty
        if !o.data.navigationmode
            (w, yloc,xloc, height, width) = lowest_widget_location_area( o )
            canvaslocx = o.data.canvaslocx
            canvaslocy = o.data.canvaslocy
            canvaslocx2 = o.data.canvaslocx + o.width - o.borderSizeH*2
            canvaslocy2 = o.data.canvaslocy + o.height - o.borderSizeV*2
            distfunc = function( to::Tuple{Int,Int,Int,Int} )
                # a = area of the candidate box
                # aOverlap = area overlap between candidate box and the canvas window
                # d = distance between the candidate box and the current box's center
                # final distance = D * A / (AOverlap + epsilon)
                a = to[3]*to[4]
                aOverlap = max( 0, min( to[2] + to[4], canvaslocx2 )-max( to[2], canvaslocx ) ) *
                    max( 0, min( to[1] + to[3], canvaslocy2 )-max( to[1], canvaslocy ) )
                d = point_from_area( yloc + height >> 1, xloc + width >> 1, to )
                if aOverlap == 0 # this this the best way to do effective distance?
                    return d + 1000
                else
                    return d + a / aOverlap
                end
            end
            wdists = Any[]
            geometric_filter( o, distfunc, 0, 0, wdists, false, 2000 )
            candidate = nothing
            mindist = 999999
            for (cw,dist) in wdists
                if dist < mindist
                    candidate = cw
                    mindist = dist
                end
            end
            if candidate != nothing
                w = lowest_widget( o )
                deep_unfocus( w )
                deep_focus( candidate )
            end
        end
        dorefresh = true
        retcode = :got_it
    elseif token == :F1 && isrootlist
        helper = newTwViewer( o.screen.value, helptext( o ), posy=:center, posx=:center, showHelp=false, showLineInfo=false, bottomText = "Esc to continue" )
        activateTwObj( helper )
        unregisterTwObj( o.screen.value, helper )
        dorefresh = true
        retcode = :got_it
    end

    if dorefresh
        refresh( o )
    end
    return retcode
end

function deep_unfocus( w::TwObj )
    tmpw = w
    while typeof( tmpw ) <: TwObj{ TwListData }
        focus = tmpw.data.focus
        tmpw2 = tmpw.data.widgets[focus]
        tmpw.data.focus = 0
        tmpw2.hasFocus = false
        tmpw = tmpw2
    end

    local par::TwObj{TwListData}
    w.hasFocus = false
    inject(w, :focus_off)
    if tmpw != w
        inject( tmpw, :focus_off )
    end
    tmpw = w
    @oncethen while( !( typeof( tmpw.window ) <: Ptr ) )
        par = tmpw.window.parent.value
        par.data.focus = 0
        par.hasFocus = false
        tmpw = par
    end
end

function set_default_focus( w::TwObj{TwListData}, rev=false )
    if w.data.focus == 0 && !isempty( w.data.widgets )
        if rev
            w.data.focus = length(w.data.widgets)
        else
            w.data.focus = 1
        end

        subw = w.data.widgets[w.data.focus]
        subw.hasFocus = true
        if typeof( subw ) <: TwObj{TwListData}
            set_default_focus( subw, rev )
        else
            ensure_visible_on_canvas( subw )
        end
    end
end

function deep_focus( w::TwObj, rev=false )
    local par::TwObj{TwListData}
    w.hasFocus = true
    if objtype(w) == :List
        set_default_focus( w, rev )
    else
        ensure_visible_on_canvas( w )
    end
    tmpw = w
    @oncethen while( !( typeof( tmpw.window ) <: Ptr ) )
        par = tmpw.window.parent.value
        for (i,c) in enumerate( par.data.widgets )
            if c == tmpw
                par.data.focus = i
                par.hasFocus = true
                tmpw = par
                break
            end
        end
        @assert par.data.focus != 0
        @assert par.hasFocus
    end
end

function geometric_filter( o::TwObj{TwListData}, distfunc::Function,
        y::Int, x::Int, list::Array, excludezero::Bool=false, cutoff::Int=100 )
    for w in o.data.widgets
        if objtype( w ) == :List
            geometric_filter( w, distfunc, y+w.window.yloc, x+w.window.xloc,
                list, excludezero, cutoff )
        elseif w.isVisible && w.acceptsFocus
            to = ( y+w.window.yloc, x+w.window.xloc, w.window.height, w.window.width )
            dist = distfunc( to )
            if dist > cutoff || dist == 0 && excludezero
                continue
            end
            push!( list, (w, dist) )
        end
    end
end

function updown_arrow_distance( to::(@compat Tuple{Int,Int,Int,Int}), from::(@compat Tuple{Int,Int,Int,Int}), sgn::Int )
    tocentx = to[2] + to[4] >> 1
    tocenty = to[1] + to[3] >> 1

    fromcentx = from[2] + from[4] >> 1
    fromcenty = from[1] + from[3] >> 1

    if sgn*tocenty >= sgn*fromcenty
        return 99999
    end

    ret = (fromcenty - tocenty)*sgn

    if from[2] > tocentx
        ret += ( from[2] - tocentx ) * 5
    elseif tocentx >= from[2]+from[4]
        ret += (tocentx - (from[2]+from[4])+1) * 5
    end
    return ret
end

function leftright_arrow_distance( to::(@compat Tuple{Int,Int,Int,Int}), from::(@compat Tuple{Int,Int,Int,Int}),sgn::Int )
    tocentx = to[2] + to[4] >> 1
    tocenty = to[1] + to[3] >> 1

    fromcentx = from[2] + from[4] >> 1
    fromcenty = from[1] + from[3] >> 1

    if sgn*tocentx >= sgn*fromcentx
        return 99999
    end

    ret = (fromcentx - tocentx)*sgn

    if from[1] > tocenty
        ret += ( from[1] - tocenty ) * 5
    elseif tocenty >= from[1]+from[3]
        ret += (tocenty - (from[1]+from[3])+1) * 5
    end
    return ret
end

# used for mouse click
# inside the area it returns zero, otherwise it's manhattan distance to the boundary
function point_from_area( y::Int, x::Int, from::(@compat Tuple{Int,Int,Int,Int}) )
    xdist = 0
    if x < from[2]
        xdist += from[2]-x
    elseif x >= from[2]+from[4]
        xdist += x - (from[2]+from[4]-1)
    end

    ydist = 0
    if y < from[1]
        ydist += from[1]-y
    elseif y >= from[1]+from[3]
        ydist += y - (from[1]+from[3]-1)
    end

    return xdist + ydist
end

function helptext( o::TwObj{TwListData} )
    focus = o.data.focus
    isrootlist = typeof( o.window ) <: Ptr
    if focus == 0
        return utf8("")
    end
    s = helptext( o.data.widgets[ focus ] )
    if isrootlist
h = utf8("""
ctrl-F4 : toggle navigation mode
mouse-click: activate nearest widget
ctrl-arrows: directional focus movements
  (normal arrows work too if not consumed by the current widget)
tab/shift-tab: cycle through all widgets
""")
        if s == "" # just the navigation text
            s = h
        else # merge the help text into a single window
            s *= "\n"* ("—"^7) * " canvas navigation " * ("—"^7) * "\n"* h
        end
    end
    s
end