text:
@testset "scan char" begin
    @test FWF.scan("testfiles/test1.txt") == [1:1, 3:19, 21:48]
    @test FWF.scan("testfiles/test1.txt", skip=1) == [3:8, 11:16, 21:24, 27:32, 36:40, 44:48]
    @test FWF.scan("testfiles/test1.txt", [' ';'a':'f']) == [3:8, 11:16, 21:24, 27:32, 36:40, 44:48]
    @test FWF.scan("testfiles/test2.txt") == [1:1, 3:5]
    @test FWF.scan("testfiles/test2.txt", nrow=2) == [1:1, 3:3, 5:5]
    @test FWF.scan("testfiles/test2.txt", nrow=3) == [1:1, 3:5]
end

@testset "scan byte" begin
    @test FWF.scan("testfiles/test1.txt", UInt8(' ')) == [1:1, 3:19, 21:48]
    @test FWF.scan("testfiles/test1.txt", UInt8(' '), skip=1) == [3:8, 11:16, 21:24, 27:32, 36:40, 44:48]
    @test FWF.scan("testfiles/test2.txt", UInt8(' ')) == [1:1, 3:5]
    @test FWF.scan("testfiles/test2.txt", UInt8(' '), nrow=2) == [1:1, 3:3, 5:5]
    @test FWF.scan("testfiles/test2.txt", UInt8(' '), nrow=3) == [1:1, 3:5]
end


text:
const DEPS_FILE = joinpath(@__DIR__, "deps.jl")

if isfile(DEPS_FILE)
    rm(DEPS_FILE)
end

using Libdl

function get_paths(path)
    extension = if Sys.isapple()
        ".dylib"
    elseif Sys.iswindows()
        Libdl.dlext
    elseif Sys.isunix()
        ".so"
    else
        error("System not supported.")
    end
    libcbcsolver = escape_string(joinpath(path, "libCbcSolver" * extension))
    libcbc       = escape_string(joinpath(path, "libCbc" * extension))
    return libcbcsolver, libcbc
end

function write_depsfile(libcbcsolver, libcbc)
    open(DEPS_FILE, "w") do io
        println(io, "const libcbcsolver = \"$(libcbcsolver)\"")
        println(io, "const libCbc = \"$(libcbc)\"")
    end
end

path = get(ENV, "JULIA_CBC_LIBRARY_PATH", nothing)
@show path
if path === nothing
    if VERSION >= v"1.3"
        # Use Cbc_jll instead.
        exit()
    else
        # Use version from BinaryProvider.jl
        include(joinpath(@__DIR__, "build_binary_provider.jl"))
    end
else
    # use the version from JULIA_CBC_LIBRARY_PATH!
    libcbcsolver, libcbc = get_paths(path)
    d = Libdl.dlopen_e(libcbcsolver)
    if d != C_NULL
        write_depsfile(libcbcsolver, libcbc)
    else
        error("""
        Unable to build Cbc.jl because your environment variable `JULIA_CBC_LIBRARY_PATH`
        points to an invalid location. It is currently:

        $(ENV["JULIA_CBC_LIBRARY_PATH"])

        Set it to point to the directory containing `libCbcSolver`, or to install the
        default binaries, run:

        delete!(ENV, "JULIA_CBC_LIBRARY_PATH")
        import Pkg
        Pkg.build("Cbc")

        then restart Julia for the changes to take effect.
        """)
    end
end

ext: jl


text:
# Input argument x specifies the initial iterate of the optimization
# algorithm. When x = -1, or when any of the entries of x are
# negative, the default setting for x is used.
#
# Since the MOSEK solver for the quadratic subproblem sometimes does
# not tolerate iterates in which most of the entries are nonzero
# ("dense" vectors), the default initial estimate for qpsubprob =
# "mosek" is a sparse vector in which only the first two entries are
# nonzero.
function mixsqp(L; x = -1,
                convergence_tol = 1e-8, eps_tol = 1e-8, sparse_tol = 1e-3,
                lowrank_method = "svd", lowrank_criterion = "by_rank",
                lowrank_rtol = 1e-8, lowrank_rank = minimum(size(L)),
                maxiter = 100, maxqpiter = 100,
                nullprior = 0, gradient_calc = "approx",
                linesearch = true, verbose = true)

# Get the number of rows (n) and columns (k) of L.
  n = size(L,1);
  k = size(L,2);

  # If any entries of input x are negative, set the initial estimate
  # to the default setting.
  #
  # When the MOSEK solver is used, the initial estimate is set to a
  # sparse vector in which all but two of the entries are zero.
  if any(x .< 0)
    x = ones(k)/k;
  end

  # If requested (i.e., if lowrank_rtol > 0), compute the partial QR
  # decomposition with relative precision "tol", then retrieve the
  # permutation matrix, P. For details on the partial QR, see
  # https://github.com/klho/LowRankApprox.jl.

  # Start timing for low-rank approximation of L.

  tic();

    if lowrank_method == "qr"
      if lowrank_criterion == "by_rank"
        F = pqrfact(L, rank = lowrank_rank);
        P = sparse(F[:P]);
      elseif lowrank_criterion == "by_rtol"
        F = pqrfact(L, rtol = lowrank_rtol);
        P = sparse(F[:P]);
      else
        error("lowrank_criterion must be either by_rank or by_rtol");
      end

    elseif lowrank_method == "svd"
      if lowrank_criterion == "by_rank"
        F = psvdfact(L, rank = lowrank_rank);
        S = Diagonal(F[:S]);
      elseif lowrank_criterion == "by_rtol"
        F = psvdfact(L, rtol = lowrank_rtol);
        S = Diagonal(F[:S]);
      else
        error("lowrank_criterion must be either by_rank or by_rtol");
      end
    end

  lowranktime = toq();

  # Summarize the analysis here.
  if verbose
    @printf("Running SQP algorithm with the following settings:\n")
    @printf("- %d x %d data matrix\n",n,k)
    @printf("- convergence tolerance = %0.2e\n",convergence_tol)
    @printf("- zero threshold        = %0.2e\n",sparse_tol)
    if lowrank_method == "qr"
      err = maximum(abs.(F[:Q]*F[:R]*P' - L));
      @printf("- partial QR tolerance  = %0.2e\n",lowrank_rtol)
      @printf("- partial QR max. error = %0.2e\n",err)
    elseif lowrank_method == "svd"
      err = maximum(abs.(F[:U]*S*F[:Vt] - L));
      @printf("- partial SVD tolerance  = %0.2e\n",lowrank_rtol)
      @printf("- partial SVD max. error = %0.2e\n",err)
    else
      @printf("- Exact derivative computation (partial QR not used).\n")
    end
  end

  # Initialize storage for the outputs obj, gmin, nnz and nqp.
  obj      = zeros(maxiter);
  gmin     = zeros(maxiter);
  nnz      = zeros(Int,maxiter);
  nqp      = zeros(Int,maxiter);
  timing   = zeros(maxiter);
  qptiming = zeros(maxiter);

  # Initialize loop variables used in the loops below so that they
  # are available outside the scope of the loop.
  i     = 0;
  j     = 0;
  D     = 0;
  t     = 0;
  numls = -1;

  # Print the column labels for reporting the algorithm's progress.
  if verbose
    @printf("iter      objective -min(g+1)  #nz #qp #ls\n")
  end

  # Repeat until we reach the maximum number of iterations, or until
  # convergence is reached.
  for i = 1:maxiter

    # Start timing the iteration.
    tic();

    # Compute the gradient and Hessian, optionally using the partial
    # QR decomposition to increase the speed of these computations.
    # gradient and Hessian computation -- Rank reduction method
        if lowrank_method == "qr"
          D = 1 ./(F[:Q]*(F[:R]*(P'*x)) .+ eps_tol);
          g = -P * F[:R]' * (F[:Q]'*D)/n;
          H = P * F[:R]' * (F[:Q]'*Diagonal(D.^2)*F[:Q])*F[:R]*P'/n + eps_tol * Diagonal(ones(k));
        elseif lowrank_method == "svd"
          D = 1 ./(F[:U]*(S*(F[:Vt]*x)) .+ eps_tol);
          g = -F[:Vt]'*(S * (F[:U]'*D))/n;
          H = (F[:V]*S*(F[:U]'*Diagonal(D.^2)*F[:U])* S*F[:Vt])/n + eps_tol * Diagonal(ones(k));
        else
          D = 1 ./(L*x .+ eps_tol);
          g = -L'*D/n;
          H = L'*Diagonal(D.^2)*L/n + eps_tol * Diagonal(ones(k));
        end

    if nullprior > 0
      g[1] -= nullprior/x[1]/n;
      H[1,1] += nullprior/x[1]^2/n;
    end

    # Report on the algorithm's progress.
    if lowrank_method == "qr"
      obj[i] = -sum(log.(F[:Q]*(F[:R]*(P'*x)) .+ eps_tol));
    elseif lowrank_method == "svd"
      obj[i] = -sum(log.(F[:U]*(S*(F[:Vt]*x)) .+ eps_tol));
    else
      obj[i] = mixobjective(L,x,eps_tol);
    end
    gmin[i] = minimum(g) .+ 1;
    nnz[i]  = length(find(x .> sparse_tol));
    nqp[i]  = j;
    if verbose
      if i == 1
          @printf("%4d %0.8e %+0.2e %4d\n",
              i,obj[i],-gmin[i],nnz[i]);
      else
          @printf("%4d %0.8e %+0.2e %4d %3d %3d\n",
              i,obj[i],-gmin[i],nnz[i],nqp[i-1],numls);
      end
    end

    # Check convergence of outer loop
    if gmin[i] >= -convergence_tol
      break
    end

    # Solve the QP subproblem using either the active-set or
    # interior-point (MOSEK) method.
    out, qptiming[i], bytes, gctime,
    memallocs = @timed if true
      y,nqp[i] = qpactiveset(x,g,H,convergence_tol = convergence_tol,sparse_tol = sparse_tol,
                      maxiter = maxqpiter);
    end

    # Perform backtracking line search
    if linesearch == true
        for t = 1:10
          if lowrank_method == "qr"
            D_new = 1 ./(F[:Q]*(F[:R]*(P'*y)) .+ eps_tol);
          elseif lowrank_method == "svd"
            D_new = 1 ./(F[:U]*(S*(F[:Vt]*y)) .+ eps_tol);
          else
            D_new = 1 ./(L*y .+ eps_tol);
          end
          if all(D_new .> 0)
            if sum(log.(D)) - sum(log.(D_new)) > sum((x-y) .* g) / 2
              break
            end
          end
          y = (y-x)/2 + x;
        end
        numls = t;
    else
        numls = -1;
    end

    # Update the solution to the original optimization problem.
    x = copy(y);

    # Get the elapsed time for the ith iteration.
    timing[i] = toq();
  end

  # Return: (1) the solution (after zeroing out any values below the
  # tolerance); (2) the value of the objective at each iteration; (3)
  # the minimum gradient value of the modified objective at each
  # iteration; (4) the number of nonzero entries in the vector at each
  # iteration; and (5) the number of inner iterations taken to solve
  # the QP subproblem at each outer iteration.
  x[x .< sparse_tol] .= 0;
  x             = x/sum(x);
  totaltime     = lowranktime + sum(timing[1:i]);
  if verbose
    @printf("Optimization took %d iterations and %0.4f seconds.\n",i,totaltime)
  end
  return Dict([("x",Vector(x)), ("totaltime",totaltime),
               ("lowranktime",lowranktime), ("obj",obj[1:i]),
               ("gmin",gmin[1:i]), ("nnz",nnz[1:i]),
               ("nqp",nqp[1:i]), ("timing",timing[1:i]),
               ("qptiming",qptiming[1:i])])
end

# Solve the QP subproblem for the mix-SQP algorithm using an active
# set method.
function qpactiveset(x, g, H; convergence_tol = 1e-8, sparse_tol = 1e-3, maxiter = 100)

  # Get the number of degrees of freedom in the optimization problem.
  k = length(x);

  # Initialize the solution to the QP subproblem.
  y      = sparse(zeros(k));
  ind    = find(x .> sparse_tol);
  y[ind] .= 1/length(ind);
  i = 0;

  # Repeat until we reach the maximum number of iterations, or until
  # convergence is reached.
  for i = 1:maxiter

    # Define the smaller QP subproblem.
    s   = length(ind);
    H_s = H[ind,ind];
    d   = H*y + 2*g .+ 1;
    d_s = d[ind];

    # Solve the smaller problem.
    p      = sparse(zeros(k));
    p_s    = -H_s\d_s;
    p[ind] = p_s;

    # Check convergence using KKT
    if norm(p_s) < convergence_tol

      # Compute the Lagrange multiplier.
      z = d;
      if all(z .>= -convergence_tol)
        break
      elseif length(ind) < k
        notind  = setdiff(1:k,ind);
        ind_min = notind[findmin(z[notind])[2]];
        ind     = sort([ind; ind_min]);
      end
    else

      # Find a feasible step length.
      alpha     = 1;
      alpha0    = -y[ind] ./ p_s;
      ind_block = find(p_s .< 0);
      alpha0    = alpha0[ind_block];
      if ~isempty(ind_block)
        v, t = findmin(alpha0);
        if v < 1

          # Blocking constraint.
          ind_block = ind[ind_block[t]];
          alpha     = v;

          # Update working set if there is a blocking constraint.
          deleteat!(ind,find(ind .- ind_block .== 0));
        end
      end

      # Move to the new "inner loop" iterate (y) along the search
      # direction.
      y = y + alpha * p;
    end
  end

  # Return the solution to the quadratic program.
  return y, Int(i)
end

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper


text:
mutable struct VecModel{Tv <: AbstractVector} <: AbstractModel
    objective::Union{Nothing, Objective}
    eq_constraints::VectorOfFunctions
    ineq_constraints::VectorOfFunctions
    sd_constraints::VectorOfFunctions
    box_min::Tv
    box_max::Tv
    init::Tv
    integer::BitVector
end
function VecModel(
        objective::Union{Nothing, Objective},
        eq_constraints::VectorOfFunctions,
        ineq_constraints::VectorOfFunctions,
        box_min,
        box_max,
        init,
        integer::BitVector)
    return VecModel(objective, eq_constraints, ineq_constraints, VectorOfFunctions(SDConstraint[]), box_min, box_max, init, integer)
end

function isfeasible(model::VecModel, x::AbstractVector; ctol = 1e-4)
    return all(getmin(model) .<= x .<= getmax(model)) &&
        all(getineqconstraints(model)(x) .<= ctol) &&
        all(-ctol .<= geteqconstraints(model)(x) .<= ctol)
end

function addvar!(m::VecModel, lb::Real, ub::Real; init::Real = lb, integer = false)
    push!(getmin(m), lb)
    push!(getmax(m), ub)
    push!(m.init, init)
    push!(m.integer, integer)
    return m
end
function addvar!(m::VecModel, lb::Vector{<:Real}, ub::Vector{<:Real}; init::Vector{<:Real} = copy(lb), integer = falses(length(lb)))
    append!(getmin(m), lb)
    append!(getmax(m), ub)
    append!(m.init, init)
    append!(m.integer, integer)
    return m
end

function getinit(m::VecModel)
    ma = getmax(m)
    mi = getmin(m)
    init = m.init
    return map(1:length(mi)) do i
        if isfinite(init[i])
            return init[i]
        else
            _ma = ma[i]
            _mi = mi[i]
            _ma == Inf && _mi == -Inf && return 0.0
            _ma == Inf && return _mi + 1.0
            _mi == -Inf && return _ma - 1.0
            return (_ma + _mi) / 2
        end
    end
end

get_objective_multiple(model::VecModel) = getobjective(model).multiple[]

function set_objective_multiple!(model::VecModel, m)
    getobjective(model).multiple[] = m
    return model
end

"""
Generic `optimize` for VecModel
"""
function optimize(model::VecModel, optimizer::AbstractOptimizer, x0, args...; kwargs...)
    workspace = Workspace(model, optimizer, copy(x0), args...; kwargs...)
    return optimize!(workspace)
end

"""
 Workspace constructor without x0
"""
function optimize(model::VecModel, optimizer::AbstractOptimizer, args...; kwargs...)
    workspace = Workspace(model, optimizer, args...; kwargs...)
    return optimize!(workspace)
end

function tovecmodel(m::AbstractModel, x0 = getmin(m))
    v, _unflatten = flatten(x0)
    unflatten = Unflatten(x0, _unflatten)
    return VecModel(
        # objective
        Objective(x -> m.objective(unflatten(x)), flags = m.objective.flags),
        # eq_constraints
        length(m.eq_constraints.fs) != 0 ? VectorOfFunctions(map(m.eq_constraints.fs) do c
            EqConstraint(x -> maybeflatten(c.f(unflatten(x)))[1], maybeflatten(c.rhs)[1], c.dim, c.flags)
        end) : VectorOfFunctions(EqConstraint[]),
        # ineq_constraints
        length(m.ineq_constraints.fs) != 0 ? VectorOfFunctions(map(m.ineq_constraints.fs) do c
            IneqConstraint(x -> maybeflatten(c.f(unflatten(x)))[1], maybeflatten(c.rhs)[1], c.dim, c.flags)
        end) : VectorOfFunctions(IneqConstraint[]),
        # sd_constraints
        length(m.sd_constraints.fs) != 0 ? VectorOfFunctions(map(m.sd_constraints.fs) do c
            SDConstraint(x -> c.f(unflatten(x)), c.dim)
        end) : VectorOfFunctions(SDConstraint[]),
        # box_min
        float.(flatten(m.box_min)[1]),
        # box_max
        float.(flatten(m.box_max)[1]),
        # init
        float.(flatten(m.init)[1]),
        # integer
        convert(BitVector, flatten(m.integer)[1]),
    ), float.(v), unflatten
end


text:
# set working directory to directory of this file
cd(@__DIR__)
# load package manager
using Pkg
# activate environment in parent folder.
# call Pkg.instantiate() if error, then continue
Pkg.activate("..")
# load dependency in Project.toml, and common functions in mytest
using Distributions, mytest
# load model specific functions
include("model_specific_functions.jl")
# from model_specific_functions.jl
model_function()

# add a package to your Project.toml
# ] add Distributions

# Package versions are added in Project.toml  under [compat]

# A function from mytest. Put common functions there.
greet()

# tests can be ran with ] test mytest


text:
# Resampling schemes for particle filters
# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

# default resampling scheme
function resample(w::Vector{Float64},
                  num_particles::Int = length(w))
    resampleSystematic(w, num_particles)
end

# More stable, faster version of rand(Categorical)
function randcat(p::Vector{Float64})
    # if(any(p .< 0)) error("Negative probabilities not allowed"); end
    r, s = rand(), one(Int)
    for j = 1:length(p)
        r -= p[j]
        if(r <= 0.0) s = j; break; end
    end

    s
end

function resampleMultinomial(w::Vector{Float64},
                             num_particles::Int)

    s = Distributions.sampler(Categorical(w))
    indx = rand(s, num_particles)

end

function resampleResidual(w::Vector{Float64},
                          num_particles::Int)

    M = length( w )

    # "Repetition counts" (plus the random part, later on):
    Ns = floor.(length(w) .* w)

    # The "remainder" or "residual" count:
    R = Int(sum( Ns ))

    # The number of particles which will be drawn stocastically:
    M_rdn = num_particles-R;

    # The modified weights:
    Ws = (M .* w - floor.(M .* w))/M_rdn;

    # Draw the deterministic part:
    indx1 = Array{Int}(undef, R)
    i=1
    for j=1:M
        for k=1:Ns[j]
            indx1[i]=j
            i = i + 1
        end
    end

    # And now draw the stocastic (Multinomial) part:
    s = Distributions.sampler(Categorical(w))
    indx2 = rand(s, M_rdn)
    indx = append!(indx1, indx2)

    return indx

end

function resampleStratified(w::Vector{Float64},
                            num_particles::Int)

    N = num_particles
    Q = cumsum(w)

    T = Array{Float64}(undef, N+1)
    for i=1:N,
        T[i] = rand()/N + (i-1)/N
    end
    T[N+1] = 1

    i=1
    j=1

    indx = Array{Int}(undef, N)
    while i<=N
        if T[i]<Q[j]
            indx[i]=j
            i=i+1
        else
            j=j+1
        end
    end

    return indx

end

function resampleSystematic(w::Vector{Float64},
                            num_particles::Int)

    N = num_particles
    Q = cumsum(w)

    T = collect(range(0, stop = maximum(Q)-1/N, length = N)) .+ rand()/N
    push!(T, 1)

    i=1
    j=1

    indx = Array{Int}(undef, N)
    while i<=N
        if (T[i]<Q[j])
            indx[i]=j
            i=i+1
        else
            j=j+1
        end
        # if j == length(w)
        #   indx[i] = j
        #   break;
        # end
    end

    return indx

end


text:
module Bayes
export StochMapDom, OpenBayesNet, FreeBayesFunctions,
  Ob, Hom, dom, codom, compose, id, otimes, mzero, braid,
  mcopy, delete

import Base: +
using AutoHashEquals
using LinearMaps
import LinearMaps: adjoint
const LMs = LinearMaps
using LinearOperators
import LinearOperators:
  adjoint, opEye, opExtension, opRestriction, opZeros

using BayesNets
import BayesNets:
  BayesNet

using LightGraphs
import LightGraphs:
  indegree, outdegree, nv

using ...Catlab, ...Doctrines
import ...Doctrines:
  Ob, Hom, dom, codom, compose, ⋅, ∘, id, oplus, ⊕, mzero, braid,
  dagger, dunit, dcounit, mcopy, Δ, delete, ◊, mmerge, ∇, create, □,
  plus, zero, coplus, cozero, meet, top, join, bottom
using ...Programs
import ...Programs: evaluate_hom

@syntax FreeBayesFunctions(ObExpr,HomExpr) MonoidalCategoryWithDiagonals begin
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  otimes(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = new(f,g; strict=true) # No normalization!
end

# BayesMaps instance
#--------------------

@auto_hash_equals struct StochMapDom
  N::Int #Number of nodes in StochMapDom
end

struct OpenBayesNet
  Net::BayesNet
  dom::Array{Int}
  codom::Array{Int}
end

@instance MonoidalCategoryWithDiagonals(StochMapDom,OpenBayesNet) begin
#   @import adjoint, +

  dom(f::OpenBayesNet) = begin
#     f_dag = f.Net.dag
# #     vIDtoSymbol_dict = Dict()
# #     for entry in f.name_to_index;
# #       vIDtoSymbol_dict[f.name_to_index[entry[1]]] = entry[1]
# #     end
#     indeg_bayes = indegree(f_dag)
#     names = []
#     for i in 1:size(indeg_bayes,1)
#         if indeg_bayes[i] == 0
#             push!(names,i)
#         end
#     end
    StochMapDom(length(f.dom))
  end
  codom(f::OpenBayesNet) = begin
#     f_dag = f.dag
#     indeg_bayes = indegree(f_dag)
#     names = []
#     for i in 1:size(indeg_bayes,1)
#         if indeg_bayes[i] != 0
#             push!(names,i)
#         end
#     end
    StochMapDom(length(f.codom))
  end

  compose(f::OpenBayesNet, g::OpenBayesNet) = begin
    dom = f.dom
    codom = [length(f.dom)+ elem for elem in g.codom]
    combined_net = deepcopy(f.Net)
    for elem in g.codom
      push!(combined_net, g.cpds[elem])
    end
    combined_dict = merge(f.Net.name_to_index, g.Net.name_to_index)
    for (key,value) in g.Net.name_to_index
      if value <= length(g.dom)
        !delete(combined_dict, key)
      else
        combined_dict[key] = value + length(f.dom)
      end
    end
    OpenBayesNet(BayesNet(combined_net.dag, combined_net.cpds, combined_dict), dom, codom)
  end

  #id(V::LinearMapDom) = LMs.UniformScalingMap(1, V.N)

  otimes(V::StochMapDom, W::StochMapDom) = begin
    StochMapDom(V.N + W.N)
  end

  otimes(f::OpenBayesNet, g::OpenBayesNet) = begin
    dom = vcat(f.dom, [nv(f.Net.dag) + elem for elem in g.dom])
    codom = vcat(f.codom, [nv(f.Net.dag) + elem for elem in g.codom])
    combined_dag = blockdiag(f.Net.dag, g.Net.dag)
    combined_cpd = vcat(f.Net.cpds, g.Net.cpds)
    combined_dict = merge(f.Net.name_to_index, g.Net.name_to_index)
    for (key,value) in g.Net.name_to_index
      combined_dict[key] = value + nv(f.Net.dag)
    end
    OpenBayesNet(BayesNet(combined_dag, combined_cpd, combined_dict), dom, codom)
  end

  id(V::StochMapDom) = begin
    OpenBayesNet(BayesNet(),[],[])
  end
  munit(::Type{StochMapDom}) = OpenBayesNet(BayesNet(),[],[])
#   mzero(::Type{LinearMapDom}) = LinearMapDom(0)
  braid(V::StochMapDom, W::StochMapDom) = begin
    OpenBayesNet(BayesNet(),[],[])
  end
  mcopy(V::StochMapDom) = OpenBayesNet(BayesNet(),[],[])
  delete(V::StochMapDom) = OpenBayesNet(BayesNet(),[],[])
#   plus(V::LinearMapDom) = LinearMap(plus_lm, mcopy_lm, V.N, 2*V.N)
#   zero(V::LinearMapDom) = LinearMap(zero_lm(V.N), delete_lm, V.N, 0)

#   plus(f::LinearMap, g::LinearMap) = f+g
#   scalar(V::LinearMapDom, c::Number) = LMs.UniformScalingMap(c, V.N)
#   antipode(V::LinearMapDom) = LMs.UniformScalingMap(-1, V.N)
end

# braid_lm(n::Int) = x::AbstractVector -> vcat(x[n+1:end], x[1:n])
# mcopy_lm(x::AbstractVector) = vcat(x, x)
# delete_lm(x::AbstractVector) = eltype(x)[]
# plus_lm(x::AbstractVector) = begin
#   n = length(x) ÷ 2
#   x[1:n] + x[n+1:end]
# end
# zero_lm(n::Int) = x::AbstractVector -> zeros(eltype(x), n)

end


text:
using KinematicRepresentations
using Test

@testset "KinematicRepresentations.jl" begin
    # Write your tests here.
end


text:
using LinearAlgebra, Interpolations, Plots
using HA_solver

Model = HA_solver.NCS_growth_params()
HA_solver.Find_eq_NCS_Smolyak(Model)


text:
module EnvelopeAnalysis
using DSP, RecipesBase, Distributions, StatsBase

include("probmodels.jl")
include("recipes.jl")
include("envstats.jl")

end


text:
struct DimsSelector{N}
    dims::NTuple{N, Int}
end
dims(args...) = DimsSelector(args)

function (d::DimsSelector)(c::CartesianIndex{N}) where N
    t = ntuple(N) do n
        return n in d.dims ? c[n] : 1
    end
    return CartesianIndex(t)
end

select(data, d::DimsSelector) = (d,) => identity => ""

function select(data, name::StringLike)
    v = getcolumn(data, Symbol(name))
    return (v,) => identity => string(name)
end

function select(data, idx::Integer)
    name = columnnames(data)[idx]
    return select(data, name)
end

select(::Nothing, v::AbstractArray) = (v,) => identity => ""

function select(data, x::Pair{<:Any, <:StringLike})
    name, label = x
    vs, _ = select(data, name)
    return vs => identity => string(label)
end

function select(data, x::Pair{<:Any, <:Any})
    name, transformation = x
    if name isa Tuple
        vs = map(n -> only(first(select(data, n))), name)
        label = ""
    else
        vs, (_, label) = select(data, name)
    end
    return vs => transformation => label
end

function select(data, x::Pair{<:Any, <:Pair})
    name, transformation_label = x
    transformation, label = transformation_label
    names = name isa Tuple ? name : (name,)
    vs = map(n -> only(first(select(data, n))), names)
    return vs => transformation => label
end


text:
function get_params(d::DataEntry)
    sort(cat(d.folder_entries..., d.group_entries..., dims=1))
end

function build_new_registry(d::DataEntry)
    params = get_params(d)
    DataFrame(vcat([x=>String[] for x in params], ["file"=>String[], "group"=>String[]]))
end

function register_entry(registry, d::DataEntry, vals, file_name, data_name)
    file_path = replace(get_folder_path(d, vals) * "/" * file_name, "\\"=>"/")
    group_path = get_group_path(d, vals) * "/" * data_name
    params_template = get_param_set()
    params = get_params(d)
    res_dict = Dict(vcat([Symbol(p)=>cfmt(params_template[p], vals[p]) for p in params], [:file=>file_path, :group=>group_path])...)
    push!(registry, res_dict)
end

function write_registry(registry, file)
    CSV.write(file, registry)
end

function load_registry(file)
    CSV.File(file, type=String) |> DataFrame
end

function load_registry(file, entry::DataEntry)
    try
        load_registry(file)
    catch
        @warn "No valid registery file found. Build a new registery."
        build_new_registry(entry)
    end
end

function query_registry(reg, vals)
    params_template = get_param_set()
    vals_str = Dict(Symbol(k)=>cfmt(params_template[k], v) for (k, v) in vals)
    vals_keys = collect(keys(vals_str))
    @from i in reg begin
        @where all((x)->getproperty(i, x)==vals_str[x], vals_keys)
        @select i
        @collect DataFrame
    end
end

function query_registry(reg, vals, condition)
    params_template = get_param_set()
    vals_str = Dict(Symbol(k)=>cfmt(params_template[k], v) for (k, v) in vals)
    vals_keys = collect(keys(vals_str))
    @from i in reg begin
        @where all((x)->getproperty(i, x)==vals_str[x], vals_keys) && condition(i)
        @select i
        @collect DataFrame
    end
end


text:
# This data structure is a singly linked list (i.e. each node has a reference to the item after).
# The value lies in efficient insertion/deletion from the middle of the list, at the cost of more storage/complexity at each insert/delete.
# We represent this structure using:
#     HeadNode --> first --> second --> ... --> last --> TailNode

# Note: this implementation is mostly copied from DoublyLinkedList

type SLLNode{T}
  ishead::Bool
  istail::Bool
  nextnode::SLLNode{T}
  data::T

  SLLNode{T}(::Type{T}, ishead::Bool, istail::Bool) = new(ishead, istail)
  SLLNode{T}(ishead::Bool, istail::Bool, nextnode::SLLNode{T}, data::T) = new(ishead, istail, nextnode, data)
end

Base.string{T}(x::SLLNode{T}) = "SLLNode{ishead=$(x.ishead), istail=$(x.istail), data=$(isdefined(x,:data) ? x.data : "")}"
Base.print{T}(io::IO, x::SLLNode{T}) = print(io, string(x))
Base.show{T}(io::IO, x::SLLNode{T}) = print(io, string(x))

type SinglyLinkedList{T}
  n::Int
  headnode::SLLNode{T}
  tailnode::SLLNode{T}

  # constructor creates an empty list with head/tail pointing to each other
  function SinglyLinkedList{T}(::Type{T})
    headnode = SLLNode{T}(T, true, false)
    tailnode = SLLNode{T}(T, false, true)
    headnode.nextnode = tailnode
    new(0, headnode, tailnode)
  end
end

SinglyLinkedList{T}(::Type{T}) = SinglyLinkedList{T}(T)

head{T}(l::SinglyLinkedList{T}) = l.headnode
tail{T}(l::SinglyLinkedList{T}) = l.tailnode
peek{T}(l::SinglyLinkedList{T}) = l.headnode.nextnode.data

# Note: iterating like "for x in X; (do something); end" translates to:
# state = start(X)
# while !done(X, state)
#   (x, state) = next(X, state)
#   (do something)
# end

# iterate with the syntax "for x in list"
Base.start{T}(l::SinglyLinkedList{T}) = l.headnode.nextnode # start with the first item (will be the tail if list is empty)
Base.done{T}(l::SinglyLinkedList{T}, state::SLLNode{T}) = state.istail
Base.next{T}(l::SinglyLinkedList{T}, state::SLLNode{T}) = (state, state.nextnode)

function insertAfter!{T}(l::SinglyLinkedList{T}, node::SLLNode{T}, newnode::SLLNode{T})
  if node.istail
    error("SinglyLinkedList: Can't insert after the tail node")
  end

  after = node.nextnode
  node.nextnode = newnode
  newnode.nextnode = after
  l.n += 1
  newnode
end

# insert data after node and return new node
function insertAfter!{T}(l::SinglyLinkedList{T}, node::SLLNode{T}, data::T)

  if node.istail
    error("SinglyLinkedList: Can't insert after the tail node")
  end

  insertAfter!(l, node, SLLNode{T}(false, false, node.nextnode, data))
end

function insertSorted!{T}(l::SinglyLinkedList{T}, data::T, compareFunction::Function, ascending::Bool = true)
  insertSorted!(l, SLLNode{T}(false, false, l.headnode, data), compareFunction, ascending)
end

function insertSorted!{T}(l::SinglyLinkedList{T}, newnode::SLLNode{T}, compareFunction::Function, ascending::Bool = true)

  # NOTE: it is expected that compareFunction is defined for 2 arguments of type T, and that
  # if the first arg is <, ==, > than the second, then it returns an int: neg, 0, pos

  # in this function, we search for the node which comes directly before the correct position, then insertAfter

  node = head(l)
  while !node.nextnode.istail
    nextnode = node.nextnode
    cmpResult = compareFunction(newnode.data, nextnode.data)
    if cmpResult != 0 && (cmpResult < 0) == ascending
      return insertAfter!(l, node, newnode)
    end
    node = nextnode
  end
  insertAfter!(l, node, newnode)
end

function Base.deleteat!{T}(l::SinglyLinkedList{T}, node::SLLNode{T})

  if node.ishead || node.istail
    return
  end

  tmp = head(l)
  while !tmp.nextnode.istail && tmp.nextnode !== node
    tmp = tmp.nextnode
  end

  if tmp.nextnode === node
    tmp.nextnode = node.nextnode
    l.n -= 1
  else
    error("Couldn't find $node in $l")
  end

  # clear this node's pointers
  node.nextnode = node
  return
end

# # add to front or back and return the node
Base.unshift!{T}(l::SinglyLinkedList{T}, data::T) = insertAfter!(l, l.headnode, data)

function Base.shift!{T}(l::SinglyLinkedList{T})
  if isempty(l)
    error("trying to shift from an empty list")
  end
  firstnode = l.headnode.nextnode
  l.headnode.nextnode = firstnode.nextnode
  l.n -= 1
  firstnode.data
end
# getfirst!{T}(l::SinglyLinkedList{T}) = shift!(l)

# # add to front or back and return the node
Base.unshift!{T}(l::SinglyLinkedList{T}, newnode::SLLNode{T}) = insertAfter!(l, l.headnode, newnode)

Base.length(l::SinglyLinkedList) = l.n
Base.isempty(l::SinglyLinkedList) = length(l) == 0

function Base.find{T}(l::SinglyLinkedList{T}, data::T)
  for x in l
    if x.data == data
      return x
    end
  end
  return tail(l)
end

function Base.print{T}(io::IO, l::SinglyLinkedList{T})
  print(io, "{")
  print(io, join([x.data for x in l], ", "))
  print(io, "}")
end

function Base.show{T}(io::IO, l::SinglyLinkedList{T})
  print(io, "SinglyLinkedList(n=$(l.n))")
  print(io, l)
end

toarray{T}(l::SinglyLinkedList{T}) = T[node.data for node in l]


text:
function setupDivSigGrad(redoSetup,redoData,pFor,Ainv,Mesh2Mesh,sigtrue,sigref)
println("------- setup or load experiment from disk ---------")
srcSpacing = [2,2]
srcPad     = [1,1,1,1]
setupFile = "testDivSigGradSEGsetupMUMPS.mat"
if isfile(setupFile)
	mf = matread(setupFile)
	if any(srcSpacing .!= mf["srcSpacing"]) || any(srcPad .!= mf["srcPad"]) || any(Mfwd.n.!=mf["nfwd"])|| any(Minv.n.!=mf["n"])
		println("source spacing or padding has changed. need to setup experiment")
		redoSetup = true
	end
else
	redoSetup = true
end
if redoSetup
	Sources,Receivers,Wd = getDCResistivitySourcesAndRecAll(nfwd,Mfwd,srcSpacing=srcSpacing,srcPad=srcPad)
	results = Dict("sigtrue"=>sigtrue,"Sources"=>Sources,"Receivers"=>Receivers,"Wd"=>Wd,"srcSpacing"=>srcSpacing,
					"srcPad"=>srcPad,"n"=>Minv.n,"nfwd"=>Mfwd.n)
	matwrite(setupFile,results)
else
	mf        = matread(setupFile)
	Sources   = mf["Sources"]
	Receivers = mf["Receivers"]
	Wd        = mf["Wd"]
end

nsrc        = size(Sources,2)
println("number of sources $nsrc")

println("------- generate or load data ----------")
# if worker == 0
	pFor = DivSigGradParam(Mfwd,Sources,Receivers,[],spzeros(0,0),Ainv);
# else
	# pFor = remotecall(worker,DivSigGradParam,Mfwd,Sources,Receivers,[],Ainv);
	# pFor = [pFor]
# end
dataFile = "testDivSigGradSEGdataMUMPS.mat"
if isfile(dataFile)
	mf = matread(dataFile)
	if size(Sources)!=size(mf["Sources"]) || norm(Sources-mf["Sources"],Inf)>1e-10
		println("sources have changed. need to recompute data.")
		redoData = true;
	end
	if size(sigtrue)!=size(mf["sigtrue"]) || norm(sigtrue-mf["sigtrue"],Inf)>1e-8
		println("sigma changed. need to recompute data.")
		redoData = true;
	end
        if length(mf["dobs"])!=size(Receivers,2)*nsrc
		println("receivers have changed. need to recompute data.")
		redoData = true;
	end
end

if redoData || !isfile(dataFile)
	println("get DC data for true model")

	@time dobsDC,     = getData(Mesh2Mesh'*vec(sigtrue),pFor)
	println("get DC data for homogeneous model")

	@time dobs0,     = getData(Mesh2Mesh'*vec(sigref),pFor)
	# if worker!=0
		# dobs0 = fetch(dobs0[1]);
		# dobsDC = fetch(dobsDC[1]);
	# end
	results = Dict("sigtrue"=>sigtrue,"sigref"=>sigref,"dobs"=>dobsDC,"dobs0"=>dobs0,"Sources"=>Sources)
	matwrite(dataFile,results)
else
	println("Reading data from file")
	mf     = matread(dataFile)
	dobsDC = mf["dobs"]
	dobs0  = mf["dobs0"]
end

return pFor,Sources,Receivers,Wd,dobsDC,dobs0
end


text:
###### TODO This file helped to refactor Bf.jl in the dir src but it doesn't executed !! ##
## <summary>
## This is the brainfuck interpreter.
##
## > 	Increment the pointer.
## < 	Decrement the pointer.
## + 	Increment the byte at the pointer.
## - 	Decrement the byte at the pointer.
## . 	Output the byte at the pointer.
## , 	Input a byte and store it in the byte at the pointer.
## [ 	Jump forward past the matching ] if the byte at the pointer is zero.
## ] 	Jump backward to the matching [ unless the byte at the pointer is zero.
##
## Extended commands, included in BrainPlus.
## @   Exits the program or if inside a function, return to prior position in main program and restore state.
## $   Overwrites the byte in storage with the byte at the pointer.
## !   Overwrites the byte at the pointer with the byte in storage.
## a,b Call function a - z.
## 0-F Sets the value of the current memory pointer to a multiple of 16.
## *   Sets the return value of a function to the value at the current memory pointer; parent storage will get return value.
## </summary>
module Interpreter
    #region Private Members

    ## <summary>
    ## Object used to swap state for a function call. This data is restored when the function terminates.
    ## </summary>
    mutable struct FunctionCallObj
        InstructionPointer :: Int
        DataPointer :: Int
        FunctionInputPointer :: Int
        #public Stack<int> CallStack #TODO
        ExitLoop :: Bool
        ExitLoopInstructionPointer :: Int
        Ticks :: Int
        char Instruction :: Char
        byte Storage :: UInt8
        byte? ReturnValue :: Char
        MaxIterationCount :: Int
    end

    struct This

        ## <summary>
        ## The "call stack"
        ## </summary>
        #private readonly Stack<int> m_CallStack = new Stack<int>();
        #TODO

        ## <summary>
        ## The input function
        ## </summary>
        m_Input :: String;

        ## <summary>
        ## The instruction set
        ## </summary>
        #private readonly IDictionary<char, Action> m_InstructionSet = new Dictionary<char, Action>();
        m_InstructionSet :: Dict() # = Dict()

        ## <summary>
        ## The memory of the program
        ## </summary>
        m_Memory :: Array{Int,1} # = [0 for i=1:32768]

        ## <summary>
        ## The output function
        ## </summary>
        m_Output :: Array{Char,1} # = []

        ## <summary>
        ## The program code
        ## </summary>
        m_Source :: Array{Char,1} # = []

        ## <summary>
        ## The data pointer
        ## </summary>
        m_DataPointer:: Int # = 0

        ## <summary>
        ## The instruction pointer
        ## </summary>
        m_InstructionPointer :: Int # = 1

        ## <summary>
        ## Boolean flag to indicate if we should skip the loop and continue execution at the next valid instruction. Used if the pointer is zero and a begin loop [ instruction is read, in which case we jump forward past the matching ].
        ## </summary>
        m_ExitLoop :: Bool # =  false

        ## <summary>
        ## Holds the instruction pointer for the start of the loop. Used to bypass all inner-loops when searching for the end of the current loop.
        ## </summary>
        m_ExitLoopInstructionPointer :: Int # = 0

        ## <summary>
        ## The list of functions and their starting instruction index.
        ## </summary>
        #private readonly Dictionary<char, FunctionInst> m_Functions = new Dictionary<char, FunctionInst>();
        m_Functions :: Dict # = Dict()

        ## <summary>
        ## Identifier for next function. Will serve as the instruction to call this function.
        ## </summary>
        m_NextFunctionCharacter :: Char

        ## <summary>
        ## The function "call stack".
        ## </summary>
        #private readonly Stack<FunctionCallObj> m_FunctionCallStack = new Stack<FunctionCallObj>();
        #TODO
        m_FunctionCallStack ::Array

        ## <summary>
        ## Pointer to the current call stack (m_FunctionCallStack or m_CallStack).
        ## </summary>
        #private Stack<int> m_CurrentCallStack;
        # TODO
        m_CurrentCallStack :: Int

        ## <summary>
        ## Pointer to a function's parent memory. When an input (,) command is executed from within a function, the function's current memory cell gets a copy of the value of the parent memory at this pointer. This allows passing multiple values as input to a function.
        ## For example: ++>++++>+<<a!.@,>,-[-<+>]<+$@
        ## Parent memory contains: 2, 4, 1. Function will contain: 2, 4 and store a value of 6 in storage. Resulting parent memory remains: 2, 4, 1. Upon next command !, parent memory will contain: 6, 4, 1. The value 6 is then displayed as output.
        ## </summary>
        m_FunctionInputPointer :: Int

        ## <summary>
        ## Number of cells available to functions. When a function is executed, an array of cells are allocated in upper-addresses (eg., 1000-1999, 2000-2999, etc.) for usage.
        ## </summary>
        _functionSize :: Int

        ## <summary>
        ## Max number of iterations for a program or function to run. Can be custom specified within a function using the syntax: @maxit=1234|function_code_here
        ## </summary>
        m_MaxIterationCount :: Int

        ## <summary>
        ## Storage memory value. Usually used to hold return values from function calls.
        ## </summary>
        m_Storage :: Int

        ## <summary>
        ## Function return value. Set by using * command (instead of print . command).
        ## </summary>
        m_ReturnValue :: UInt8

        ## <summary>
        ## Options for function behavior in the interpreter.
        ## </summary>
        #private Function[] m_Options = null;
        #TODO
        m_Options :: Array

        #endregion

        #region Public Members

        ## <summary>
        ## Number of instructions executed within the main program or the current function.
        ## </summary>
        m_Ticks :: Int

        ## <summary>
        ## Number of total instructions executed, including within functions.
        ## </summary>
        m_TotalTicks :: Int

        ## <summary>
        ## Flag to stop execution of the program.
        ## </summary>
        m_Stop :: Bool

        ## <summary>
        ## List of executed functions in the main program. Used for reference purposes by the GA to determine which functions were executed in the program (not functions calling other functions).
        ## </summary>
        #public Dictionary<char, int> m_ExecutedFunctions = new Dictionary<char, int>();
        m_ExecutedFunctions :: Dict

        This() = new("", Dict(), [0 for i=1:32768], [], [], 0, 1, false, 0,
        Dict(), 'a', [], 0, 0, 300, 0, 0, 0, [], 0, 0, false, Dict() )
    end

    ## <summary>
    ## Read-only access to the current data pointer index in memory.
    ## </summary>
    m_CurrentDataPointer() = m_DataPointer # { get { return m_DataPointer; } }

    ## <summary>
    ## Read-only access to the current instruction pointer index.
    ## </summary>
    #m_CurrentInstructionPointer{ get { return m_InstructionPointer; } }
    m_CurrentInstructionPointer() = m_InstructionPointer

    ## <summary>
    ## True if a function is currently running. False if the main program is running.
    ## </summary>
    #public bool IsInsideFunction { get { return m_FunctionCallStack.Count > 0; } }
    # TODO IsInsideFunction() = m_FunctionCallStack.Count > 0

    ## <summary>
    ## True if currently inside a loop []. False otherwise. Note, check IsInsideFunction to tell if this is a loop within a function or the main program.
    ## </summary>
    #public bool IsInsideLoop { get { return m_CurrentCallStack.Count > 0; } }
    # TODO

    ## <summary>
    ## The name of the currently executing function or null.
    ## </summary>
    #public char? m_CurrentFunction { get { if (IsInsideFunction) return m_FunctionCallStack.Peek().Instruction; else return null; } }
    #TODO

    #endregion

    ## <summary>
    ## Constructor
    ## </summary>
    ## <param name="programCode"></param>
    ## <param name="input">Function to call when input command (,) is executed.</param>
    ## <param name="output">Function to call when output command (.) is executed.</param>
    ## <param name="function">Callback handler to notify that a function is being executed: callback(instruction).</param>
    ## <param name="options">Additional interpreter options.</param>
    function Interpreter(programCode :: String, input :: Array{UInt8,1}, Action<byte> output; Action<char> function = null, Function[] options = null)
        // Save the program code
        this = This()
        this.m_Source = collect(programCode) #programCode.ToCharArray();

        // Store the i/o delegates
        this.m_Input = input;
        this.m_Output = output;

        // Set any additional options.
        if (options != null)
        {
            this.m_Options = options;
        }

        m_CurrentCallStack = m_CallStack;

        // Create the instruction set for Basic Brainfuck.
        this.m_InstructionSet['+'] = () -> ( (!this.m_ExitLoop) ? this.m_Memory[this.m_DataPointer] += 1 )
        this.m_InstructionSet['+'] = () -> ( (!this.m_ExitLoop) ? this.m_Memory[this.m_DataPointer] -= 1 )

        this.m_InstructionSet.Add('+', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer]++; });
        this.m_InstructionSet.Add('-', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer]--; });

        this.m_InstructionSet.Add('>', () => { if (!m_ExitLoop) this.m_DataPointer++; });
        this.m_InstructionSet.Add('<', () => { if (!m_ExitLoop) this.m_DataPointer--; });

        this.m_InstructionSet.Add('.', () => { if (!m_ExitLoop) this.m_Output(this.m_Memory[this.m_DataPointer]); });

        // Prompt for input. If inside a function, pull input from parent memory, using the current FunctionInputPointer. Each call for input advances the parent memory cell that gets read from, allowing the passing of multiple values as input to a function.
        this.m_InstructionSet.Add(',', () => { if (!m_ExitLoop) m_Memory[this.m_DataPointer] = IsInsideFunction ? this.m_Memory[this.m_FunctionInputPointer++] : this.m_Input(); });

        this.m_InstructionSet.Add('[', () =>
        {
            if (!m_ExitLoop && this.m_Memory[this.m_DataPointer] == 0)
            {
                // Jump forward to the matching ] and exit this loop (skip over all inner loops).
                m_ExitLoop = true;

                // Remember this instruction pointer, so when we get past all inner loops and finally pop this one off the stack, we know we're done.
                m_ExitLoopInstructionPointer = this.m_InstructionPointer;
            }

            this.m_CurrentCallStack.Push(this.m_InstructionPointer);
        });
        this.m_InstructionSet.Add(']', () =>
        {
            var temp = this.m_CurrentCallStack.Pop();

            if (!m_ExitLoop)
            {
                this.m_InstructionPointer = this.m_Memory[this.m_DataPointer] != 0
                    ? temp - 1
                    : this.m_InstructionPointer;
            }
            else
            {
                // Continue executing after loop.
                if (temp == m_ExitLoopInstructionPointer)
                {
                    // We've finally exited the loop.
                    m_ExitLoop = false;
                    m_ExitLoopInstructionPointer = 0;
                }
            }
        });

        // Create the instruction set for Brainfuck Extended Type 3.
        this.m_InstructionSet.Add('0', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 0; });
        this.m_InstructionSet.Add('1', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 16; });
        this.m_InstructionSet.Add('2', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 32; });
        this.m_InstructionSet.Add('3', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 48; });
        this.m_InstructionSet.Add('4', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 64; });
        this.m_InstructionSet.Add('5', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 80; });
        this.m_InstructionSet.Add('6', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 96; });
        this.m_InstructionSet.Add('7', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 112; });
        this.m_InstructionSet.Add('8', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 128; });
        this.m_InstructionSet.Add('9', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 144; });
        this.m_InstructionSet.Add('A', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 160; });
        this.m_InstructionSet.Add('B', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 176; });
        this.m_InstructionSet.Add('C', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 192; });
        this.m_InstructionSet.Add('D', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 208; });
        this.m_InstructionSet.Add('E', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 224; });
        this.m_InstructionSet.Add('F', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = 240; });
        this.m_InstructionSet.Add('*', () => { if (!m_ExitLoop) this.m_ReturnValue = this.m_Memory[this.m_DataPointer]; });
        this.m_InstructionSet.Add('@', () =>
        {
            if (IsInsideFunction)
            {
                // Exit function.
                var temp = m_FunctionCallStack.Pop();

                // Restore the data pointer.
                this.m_DataPointer = temp.DataPointer;

                /*if (this.m_ReturnValue.HasValue)
                {
                    this.m_Memory[this.m_DataPointer] = this.m_ReturnValue.Value;
                }*/

                // Restore the call stack.
                this.m_CurrentCallStack = temp.CallStack;
                // Restore exit loop status.
                this.m_ExitLoop = temp.ExitLoop;
                // Restore exit loop instruction pointer.
                this.m_ExitLoopInstructionPointer = temp.ExitLoopInstructionPointer;
                // Restore ticks.
                this.m_Ticks = temp.Ticks;
                // Restore global storage.
                this.m_Storage = this.m_ReturnValue.HasValue ? this.m_ReturnValue.Value : temp.Storage;
                // Restore parent return value.
                this.m_ReturnValue = temp.ReturnValue;
                // Restore max iteraction count.
                this.m_MaxIterationCount = temp.MaxIterationCount;
                // Restore the instruction pointer.
                this.m_InstructionPointer = temp.InstructionPointer;
                // Restore function input pointer.
                this.m_FunctionInputPointer = temp.FunctionInputPointer;
            }
            else
            {
                // Exit program.
                this.m_Stop = true;
            }
        });
        this.m_InstructionSet.Add('$', () =>
        {
            if (!m_ExitLoop)
            {
                // If we're inside a function, use the function's own global storage (separate from the main program).
                // However, if this is the last storage command in the function code, then use the main/calling-function storage, to allow returning a value.
                if (IsInsideFunction && this.m_Source[m_InstructionPointer + 1] == '@')
                {
                    // Set function return value.
                    this.m_ReturnValue = this.m_Memory[this.m_DataPointer];
                }
                else
                {
                    // Set global storage for this main program or function.
                    this.m_Storage = this.m_Memory[this.m_DataPointer];
                }
            }
        });

        this.m_InstructionSet.Add('!', () => { if (!m_ExitLoop) this.m_Memory[this.m_DataPointer] = this.m_Storage; });

        // Scan code for function definitions and store their starting memory addresses.
        ScanFunctions(programCode);

        // If we found any functions, create the instruction set for them.
        for (char inst = 'a'; inst < m_NextFunctionCharacter; inst++)
        {
            char instruction = inst; // closure
            this.m_InstructionSet.Add(instruction, () =>
            {
                if (!m_ExitLoop)
                {
                    // Record a list of executed function names from the main program (not a function calling another function).
                    if (!IsInsideFunction)
                    {
                        if (m_ExecutedFunctions.ContainsKey(instruction))
                        {
                            m_ExecutedFunctions[instruction]++;
                        }
                        else
                        {
                            m_ExecutedFunctions.Add(instruction, 1);
                        }
                    }

                    if (function != null)
                    {
                        // Notify caller of a function being executed.
                        function(instruction);
                    }

                    // Store the current instruction pointer and data pointer before we move to the function.
                    var functionCallObj = new FunctionCallObj { InstructionPointer = this.m_InstructionPointer, DataPointer = this.m_DataPointer, FunctionInputPointer = this.m_FunctionInputPointer, CallStack = this.m_CurrentCallStack, ExitLoop = this.m_ExitLoop, ExitLoopInstructionPointer = this.m_ExitLoopInstructionPointer, Ticks = this.m_Ticks, Instruction = instruction, Storage = this.m_Storage, ReturnValue = this.m_ReturnValue, MaxIterationCount = this.m_MaxIterationCount };
                    this.m_FunctionCallStack.Push(functionCallObj);

                    // Give the function a fresh call stack.
                    this.m_CurrentCallStack = new Stack<int>();
                    this.m_ExitLoop = false;
                    this.m_ExitLoopInstructionPointer = 0;

                    // Initialize the function global storage.
                    this.m_Storage = 0;
                    this.m_ReturnValue = null;

                    // Load options for this function.
                    FunctionInst functionOptions = m_Functions[instruction];

                    // Set the function input pointer to the parent's starting memory. Calls for input (,) from within the function will read from parent's memory, each call advances the parent memory cell that gets read from. This allows passing multiple values to a function.
                    // Note, if we set the starting m_FunctionInputPointer to 0, functions will read from the first input position (0).
                    // If we set it to m_DataPointer, functions will read input from the current position in the parent memory (n). This is trickier for the GA to figure out, because it may have to downshift the memory back to 0 before calling the function so that the function gets all input. Setting this to 0 makes it easier for the function to get the input.
                    this.m_FunctionInputPointer = functionOptions.ReadInputAtMemoryStart ? 0 : this.m_DataPointer;

                    // Set the data pointer to the functions starting memory address.
                    this.m_DataPointer = _functionSize * (instruction - 96); // each function gets a space of 1000 memory slots.

                    // Clear function memory.
                    Array.Clear(this.m_Memory, this.m_DataPointer, _functionSize);

                    // Set ticks to 0.
                    this.m_Ticks = 0;

                    // Set the max iteration count for this function, if one was specified.
                    this.m_MaxIterationCount = functionOptions.MaxIterationCount > 0 ? functionOptions.MaxIterationCount : this.m_MaxIterationCount;

                    // Set the instruction pointer to the beginning of the function.
                    this.m_InstructionPointer = functionOptions.InstructionPointer;
                }
            });
        }
    end

    ## <summary>
    ## Run the program
    ## </summary>
    public void Run(int maxInstructions = 0)
    {
        m_Ticks = 0;
        m_TotalTicks = 0;
        m_Stop = false;

        if (maxInstructions > 0)
        {
            RunLimited(maxInstructions);
        }
        else
        {
            RunUnlimited();
        }
    }

    ## <summary>
    ## Run the program with a maximum number of instructions before throwing an exception. Avoids infinite loops.
    ## </summary>
    ## <param name="maxInstructions">Max number of instructions to execute</param>
    function runLimited(int maxInstructions)
        m_MaxIterationCount = maxInstructions;

        # Iterate through the whole program source
        while (this.m_InstructionPointer < this.m_Source.Length && !m_Stop)
            # Fetch the next instruction
            char instruction = this.m_Source[this.m_InstructionPointer];

            # See if that IS an instruction and execute it if so
            Action action;
            if (this.m_InstructionSet.TryGetValue(instruction, out action))
                # Yes, it was - execute
                action();
            end

            # Next instruction
            this.m_InstructionPointer++;

            # Have we exceeded the max instruction count?
            if (m_MaxIterationCount > 0 && m_Ticks >= m_MaxIterationCount)
                if (IsInsideFunction)
                    # We're inside a function, but ran out of instructions. Exit the function, but continue.
                    if (this.m_InstructionSet.TryGetValue('@', out action))
                        action();
                        this.m_InstructionPointer++;
                    end
                else
                    break;
                end
            end

            m_Ticks++;
            m_TotalTicks++;
        end
    end

    ## <summary>
    ## Run the program
    ## </summary>
    private void RunUnlimited()
    {
        // Iterate through the whole program source
        while (this.m_InstructionPointer < this.m_Source.Length && !m_Stop)
        {
            // Fetch the next instruction
            char instruction = this.m_Source[this.m_InstructionPointer];

            // See if that IS an instruction and execute it if so
            Action action;
            if (this.m_InstructionSet.TryGetValue(instruction, out action))
            {
                // Yes, it was - execute
                action();
            }

            // Next instruction
            this.m_InstructionPointer++;

            m_Ticks++;
            m_TotalTicks++;
        }
    }

    ## <summary>
    ## Pre-scan the program code to record function instruction pointers.
    ## </summary>
    private void ScanFunctions(string source)
    {
        this.m_InstructionPointer = source.IndexOf('@');
        while (this.m_InstructionPointer > -1 && this.m_InstructionPointer < source.Length - 1 && !m_Stop)
        {
            // Retrieve any settings for this function.
            Function functionDetail = m_Options != null && m_Options.Length > m_Functions.Count ? m_Options[m_Functions.Count] : null;

            // Store the function.
            m_Functions.Add(m_NextFunctionCharacter++, new FunctionInst(this.m_InstructionPointer, functionDetail));

            this.m_InstructionPointer = source.IndexOf('@', this.m_InstructionPointer + 1);
        }

        this.m_InstructionPointer = 0;
    }
end


text:
using Test, SciMLBase

showtest(x, with_color=true) = sprint() do io
    Base.show(IOContext(io, :color => with_color), MIME"text/plain"(), x)
end

@testset "Respect :color setting in IOContext" begin
    local prob = ODEProblem(identity, 1.0, (0.0, 1.0), ones(1,1))

    # This test might be too specific
    # @test showtest(prob, true) == "\e[36mODEProblem\e[0m with uType \e[36mFloat64\e[0m and tType \e[36mFloat64\e[0m. In-place: \e[36mfalse\e[0m\ntimespan: (0.0, 1.0)\nu0: 1.0"

    # Only test whether the colorizers are printed:
    @test occursin(SciMLBase.TYPE_COLOR, showtest(prob, true))
    @test occursin(SciMLBase.NO_COLOR, showtest(prob, true))

    @test !occursin(SciMLBase.TYPE_COLOR, showtest(prob, false))
    @test !occursin(SciMLBase.NO_COLOR, showtest(prob, false))
end


text:
using InfinitesimalGenerators, Test, Statistics, LinearAlgebra,  Expokit

xbar = 0.0
κ = 0.1
σ = 0.02
X = OrnsteinUhlenbeck(; xbar = xbar, κ = κ, σ = σ, length = 1000)

## Feynman-Kac
ψ = X.x.^2
t = range(0, stop = 100, step = 1/10)
u = feynman_kac(generator(X); t = t, ψ = ψ, direction = :forward)
@test maximum(abs, u[:, 50] .- expmv(t[50], generator(X), ψ)) <= 1e-3
@test maximum(abs, u[:, 200] .- expmv(t[200], generator(X), ψ)) <= 1e-3
@test maximum(abs, u[:, end] .- expmv(t[end], generator(X), ψ)) <= 1e-5
@test maximum(abs, feynman_kac(generator(X); t = t, ψ = ψ, direction = :forward) .- feynman_kac(generator(X); t = collect(t), ψ = ψ, direction = :forward)) <= 1e-5

## Multiplicative Functional dM/M = x dt
m = AdditiveFunctionalDiffusion(X, X.x, zeros(length(X.x)))
η, r = cgf(m)(1)
@test η ≈ xbar + 0.5 * σ^2 / κ^2 atol = 1e-2
r_analytic = exp.(X.x ./ κ)
@test norm(r ./ sum(r) .- r_analytic ./ sum(r_analytic)) <= 2 * 1e-3
t = range(0, stop = 200, step = 1/10)
u = feynman_kac(generator(m)(1); t = t, direction = :forward)
@test log.(stationary_distribution(X)' * u[:, end]) ./ t[end] ≈ η atol = 1e-2

# test speed
μm = - 0.06
m = AdditiveFunctionalDiffusion(X, X.x .+ μm, zeros(length(X.x)))
ζ = tail_index(m)
@test μm * ζ + 0.5 * ζ^2 * (σ^2 / κ^2) ≈ 0.0 atol = 1e-2
η, r = cgf(m ; eigenvector = :right)(ζ)
η, l = cgf(m ; eigenvector = :left)(ζ)
f =  exp.(ζ .* X.x ./ κ)
norm(f ./ sum(f) .- r ./ sum(r)) <= 1e-2
ψ_reaching = r.* l ./ sum(r .* l)
speed = sum(ψ_reaching .* m.μm)
@test speed ≈ μm  +  ζ * (σ^2 / κ^2) atol = 1e-2

# test transformation with a funciton # m2 = p * M.
# This does not work very well. Note that it works only if the distribution of p has a thinner tail than the distirbuiton of M
m = AdditiveFunctionalDiffusion(X, X.x .- 0.06, zeros(length(X.x)))
ζ = tail_index(m)
η, r = cgf(m ; eigenvector = :right)(ζ)
η, l = cgf(m ; eigenvector = :left)(ζ)
p =  exp.(5 .* X.x)
m2 = AdditiveFunctionalDiffusion(X, m.μm .+ (generator(m.X) * log.(p)),  (InfinitesimalGenerators.∂(X) * log.(p)) .* X.σx, ρ = 1)
η2, r2 = cgf(m2; eigenvector = :right)(ζ)
η2, l2 = cgf(m2; eigenvector = :left)(ζ)
r3 = (r ./ p.^ζ) ./ sum(r ./ p.^ζ)
r2 = r2 ./ sum(r2)
#@test r2 ≈ r3 rtol = 1e-2
l3 = (l .* p.^ζ) ./ sum(l .* p.^ζ)
#@test l2 ≈ l3 rtol = 1e-1

l' * (generator(m.X) * log.(p))
ψ_reaching = l.*r ./ sum(l.* r)
ψ_reaching' * (generator(m.X) * log.(p))

## test left and right eigenvector with correlation
m = AdditiveFunctionalDiffusion(X, X.x, 0.01 * ones(length(X.x)); ρ = 1)
η, r = cgf(m; eigenvector = :right)(1)
η, l = cgf(m; eigenvector = :left)(1)

ψ_tilde = stationary_distribution(DiffusionProcess(X.x, X.μx .+ m.ρ .* m.σm .* X.σx, X.σx))
@test (r .* ψ_tilde) ./ sum(r .* ψ_tilde) ≈ l rtol = 1e-3

# Test Multiplicative with ρ = 0.0
μm = -0.01
σm = 0.1
m = AdditiveFunctionalDiffusion(X, μm .+ X.x, σm .* ones(length(X.x)))
ζ = tail_index(m)
ζ_analytic = 2 * (-μm) / (σm^2 + (σ / κ)^2)
@test ζ ≈ ζ_analytic atol = 1e-2
η, r = cgf(m; eigenvector = :right)(ζ)
η, l = cgf(m; eigenvector = :left)(ζ)
@test η ≈ 0.0 atol = 1e-4
ψ = stationary_distribution(X)
@test (r .* ψ) ./ sum(r .* ψ) ≈ l rtol = 1e-3

# Test that the modified process μ + σ^2 ∂ ln(r) has a stationary distribution given by $r^2ψ$
X = OrnsteinUhlenbeck(;κ =κ, σ = σ, length = 1000)
m = AdditiveFunctionalDiffusion(X, μm .+ X.x .- 0.02, σm .* ones(length(X.x)))
ψ = stationary_distribution(X)
ζ = tail_index(m)
η, r = cgf(m; eigenvector = :right)(ζ)
η, l = cgf(m; eigenvector = :left)(ζ)
ψ_cond = stationary_distribution(DiffusionProcess(X.x, X.μx .+ X.σx.^2 .* (InfinitesimalGenerators.∂(X) * log.(r)), X.σx))
@test (r.^2 .* ψ) ./ sum(r.^2 .* ψ) ≈ ψ_cond rtol = 1e-1

# Test Multiplicative with ρ = 1.0
m = AdditiveFunctionalDiffusion(X, m.μm, m.σm; ρ = 1.0)
ζ = tail_index(m)
η, r = cgf(m; eigenvector = :right)(ζ)
η, l = cgf(m; eigenvector = :left)(ζ)
@test η ≈ 0.0 atol = 1e-3
ψ_tilde = stationary_distribution(DiffusionProcess(X.x, X.μx .+ ζ .* m.σm .* m.ρ .* X.σx , X.σx))
@test (r .* ψ_tilde) ./ sum(r .* ψ_tilde) ≈ l rtol = 1e-3

# Test CIR
gbar = 0.03
σ = 0.01
X = CoxIngersollRoss(xbar = gbar, κ = κ, σ = σ)
m = AdditiveFunctionalDiffusion(X, X.x, zeros(length(X.x)))
η_analytic = gbar * κ^2 / σ^2 * (1 - sqrt(1 - 2 * σ^2 / κ^2))
@test cgf(m)(1.0)[1] ≈ η_analytic rtol = 1e-2

# for CIR the speed is given by
# speed_analytic = - (g .- 0.009) + xbar * κ / sqrt(κ^2 - 2 * σ^2 * ζ)
# η, r = cgf(m, eigenvector = :right)(ζ)


text:
# TODO(kykim): Try a lock-free approach using Atomics.
mutable struct PISActionNode{S,A}
    id::Int
    a_label::A
    n::Int
    q::Float64
    transitions::Vector{Tuple{S,Float64}}
    unique_transitions::Set{S}
    n_a_children::Int
    c_cdf_est::RunningCDFEstimator  # Conditional CDF.
    a_lock::ReentrantLock
end
PISActionNode(id::Int, a::A, n::Int, q::Float64, transitions::Vector{Tuple{S,Float64}}) where {S,A} =
    PISActionNode{S,A}(id, a, n, q, transitions, Set{S}(), 0, RunningCDFEstimator([0.0], [1e-5]), ReentrantLock())

@inline n(n::PISActionNode) = n.n
@inline q(n::PISActionNode) = n.q
@inline n_a_children(n::PISActionNode) = n.n_a_children

mutable struct PISStateNode{S,A}
    id::Int
    s_label::S
    total_n::Int
    children::Vector{A}
    s_lock::ReentrantLock
    # Action nodes currently being evaluated. Used for applying virtual loss.
    a_selected::Set{A}
end
PISStateNode(id::Int, s::S, total_n::Int, children::Vector{A}) where {S,A} =
    PISStateNode{S,A}(id, s, total_n, children, ReentrantLock(), Set{A}())

@inline total_n(n::PISStateNode) = n.total_n
@inline children(n::PISStateNode) = n.children
@inline n_children(n::PISStateNode) = length(children(n))
@inline isroot(n::PISStateNode) = (n.id == 1)

mutable struct PISTree{S,A}
    root::Union{Nothing,S}

    state_nodes::Dict{S, PISStateNode}
    state_action_nodes::Dict{Tuple{S,A}, PISActionNode}

    cdf_est::RunningCDFEstimator
    costs::Vector{Float64}
    weights::Vector{Float64}

    state_nodes_lock::ReentrantLock
    state_action_nodes_lock::ReentrantLock
    costs_weights_lock::ReentrantLock

    # Used to assign unique ids to nodes.
    _s_id_counter::Threads.Atomic{Int}
    _a_id_counter::Threads.Atomic{Int}

    function PISTree{S,A}(root::Union{Nothing,S}=nothing) where {S,A}
        return new(root,
                   Dict{S, PISStateNode}(),
                   Dict{Tuple{S,A}, PISActionNode}(),

                   RunningCDFEstimator([0.0], [1e-7]),
                   sizehint!(Float64[], 10_000),
                   sizehint!(Float64[], 10_000),

                   ReentrantLock(),
                   ReentrantLock(),
                   ReentrantLock(),

                   Threads.Atomic{Int}(1),
                   Threads.Atomic{Int}(1))
    end
end

Base.isempty(tree::PISTree) = isempty(tree.state_nodes)

function insert_state_node!(tree::PISTree{S,A}, s::S) where {S,A}
    snode = nothing
    Base.@lock tree.state_nodes_lock begin
        haskey(tree.state_nodes, s) && return tree.state_nodes[s]
        id = Threads.atomic_add!(tree._s_id_counter, 1)
        snode = PISStateNode(id, s, 0, Vector{A}())
        tree.state_nodes[s] = snode
    end
    return snode
end

function sanode_key(s_label::S, a_label::A) where {S,A}
    return tuple(s_label, a_label)
end

function insert_action_node!(tree::PISTree{S,A}, snode::PISStateNode{S,A},
                             a::Union{S,A}, n0::Int, q0::Float64) where {S,A}
    sanode = nothing
    Base.@lock tree.state_action_nodes_lock begin
        state_action_key = sanode_key(snode.s_label, a)
        haskey(tree.state_action_nodes, state_action_key) && return tree.state_action_nodes[state_action_key]
        id = Threads.atomic_add!(tree._a_id_counter, 1)
        sanode = PISActionNode(id, a, n0, q0, Vector{Tuple{S,Float64}}())
        tree.state_action_nodes[state_action_key] = sanode
    end
    Base.@lock snode.s_lock begin
        push!(snode.children, a)
        snode.total_n += n0
    end
    return sanode
end


text:
struct ParameterContainer
    update_ref::UpdateRef
    parameter_array::JuMP.Containers.DenseAxisArray
    multiplier_array::JuMP.Containers.DenseAxisArray
end

get_parameter_array(c::ParameterContainer) = c.parameter_array
get_multiplier_array(c::ParameterContainer) = c.multiplier_array
Base.length(c::ParameterContainer) = length(c.parameter_array)
Base.size(c::ParameterContainer) = size(c.parameter_array)

const ParametersContainer = Dict{Symbol, ParameterContainer}


text:
module V

const NAMES = [
    "EGF",
    "ErbB1",
    "EGF_ErbB1",
    "HRG",
    "ErbB3",
    "HRG_ErbB3",
    "ErbB4",
    "HRG_ErbB4",
    "E11",
    "ErbB2",
    "E12",
    "E23",
    "E34",
    "E24",
    "E44",
    "E11P",
    "E12P",
    "E23P",
    "E24P",
    "E34P",
    "E44P",
    "Grb2",
    "E11P_Grb2",
    "Shc",
    "E11P_Shc",
    "RasGAP",
    "E11P_RasGAP",
    "E12P_Grb2",
    "E12P_Shc",
    "E12P_RasGAP",
    "E23P_Grb2",
    "E23P_Shc",
    "PI3K",
    "E23P_PI3K",
    "E23P_RasGAP",
    "E24P_Grb2",
    "E24P_Shc",
    "E24P_PI3K",
    "E24P_RasGAP",
    "E34P_Grb2",
    "E34P_Shc",
    "E34P_PI3K",
    "E34P_RasGAP",
    "E44P_Grb2",
    "E44P_Shc",
    "E44P_PI3K",
    "E44P_RasGAP",
    "E11P_Shc_Grb2",
    "E12P_Shc_Grb2",
    "E23P_Shc_Grb2",
    "E24P_Shc_Grb2",
    "E34P_Shc_Grb2",
    "E44P_Shc_Grb2",
    "E11P_ShcP_Grb2",
    "E12P_ShcP_Grb2",
    "E23P_ShcP_Grb2",
    "E24P_ShcP_Grb2",
    "E34P_ShcP_Grb2",
    "E44P_ShcP_Grb2",
    "PTP1B",
    "E11P_ShcP_Grb2_PTP1B",
    "E12P_ShcP_Grb2_PTP1B",
    "E23P_ShcP_Grb2_PTP1B",
    "E24P_ShcP_Grb2_PTP1B",
    "E34P_ShcP_Grb2_PTP1B",
    "E44P_ShcP_Grb2_PTP1B",
    "SOS",
    "E11P_Grb2_SOS",
    "E12P_Grb2_SOS",
    "E23P_Grb2_SOS",
    "E24P_Grb2_SOS",
    "E34P_Grb2_SOS",
    "E44P_Grb2_SOS",
    "E11P_ShcP_Grb2_SOS",
    "E12P_ShcP_Grb2_SOS",
    "E23P_ShcP_Grb2_SOS",
    "E24P_ShcP_Grb2_SOS",
    "E34P_ShcP_Grb2_SOS",
    "E44P_ShcP_Grb2_SOS",
    "Gab1",
    "E11P_Grb2_Gab1",
    "E12P_Grb2_Gab1",
    "E23P_Grb2_Gab1",
    "E24P_Grb2_Gab1",
    "E34P_Grb2_Gab1",
    "E44P_Grb2_Gab1",
    "E11P_ShcP_Grb2_Gab1",
    "E12P_ShcP_Grb2_Gab1",
    "E23P_ShcP_Grb2_Gab1",
    "E24P_ShcP_Grb2_Gab1",
    "E34P_ShcP_Grb2_Gab1",
    "E44P_ShcP_Grb2_Gab1",
    "RasGDP",
    "E11P_Grb2_SOS_RasGDP",
    "E12P_Grb2_SOS_RasGDP",
    "E23P_Grb2_SOS_RasGDP",
    "E24P_Grb2_SOS_RasGDP",
    "E34P_Grb2_SOS_RasGDP",
    "E44P_Grb2_SOS_RasGDP",
    "E11P_ShcP_Grb2_SOS_RasGDP",
    "E12P_ShcP_Grb2_SOS_RasGDP",
    "E23P_ShcP_Grb2_SOS_RasGDP",
    "E24P_ShcP_Grb2_SOS_RasGDP",
    "E34P_ShcP_Grb2_SOS_RasGDP",
    "E44P_ShcP_Grb2_SOS_RasGDP",
    "RasGTP",
    "Raf",
    "RafP",
    "MEK",
    "MEKP",
    "ERK",
    "ERKP",
    "ERKPP",
    "SOSP",
    "Gab1P",
    "E11P_Grb2_Gab1P",
    "E12P_Grb2_Gab1P",
    "E23P_Grb2_Gab1P",
    "E24P_Grb2_Gab1P",
    "E34P_Grb2_Gab1P",
    "E44P_Grb2_Gab1P",
    "E11P_ShcP_Grb2_Gab1P",
    "E12P_ShcP_Grb2_Gab1P",
    "E23P_ShcP_Grb2_Gab1P",
    "E24P_ShcP_Grb2_Gab1P",
    "E34P_ShcP_Grb2_Gab1P",
    "E44P_ShcP_Grb2_Gab1P",
    "PIP3_Gab1",
    "PIP3_Gab1P",
    "E11P_Grb2_Gab1P_PTP1B",
    "E12P_Grb2_Gab1P_PTP1B",
    "E23P_Grb2_Gab1P_PTP1B",
    "E24P_Grb2_Gab1P_PTP1B",
    "E34P_Grb2_Gab1P_PTP1B",
    "E44P_Grb2_Gab1P_PTP1B",
    "E11P_ShcP_Grb2_Gab1P_PTP1B",
    "E12P_ShcP_Grb2_Gab1P_PTP1B",
    "E23P_ShcP_Grb2_Gab1P_PTP1B",
    "E24P_ShcP_Grb2_Gab1P_PTP1B",
    "E34P_ShcP_Grb2_Gab1P_PTP1B",
    "E44P_ShcP_Grb2_Gab1P_PTP1B",
    "PIP3_Gab1P_PTP1B",
    "E11P_Grb2_Gab1P_RasGAP",
    "E12P_Grb2_Gab1P_RasGAP",
    "E23P_Grb2_Gab1P_RasGAP",
    "E24P_Grb2_Gab1P_RasGAP",
    "E34P_Grb2_Gab1P_RasGAP",
    "E44P_Grb2_Gab1P_RasGAP",
    "E11P_ShcP_Grb2_Gab1P_RasGAP",
    "E12P_ShcP_Grb2_Gab1P_RasGAP",
    "E23P_ShcP_Grb2_Gab1P_RasGAP",
    "E24P_ShcP_Grb2_Gab1P_RasGAP",
    "E34P_ShcP_Grb2_Gab1P_RasGAP",
    "E44P_ShcP_Grb2_Gab1P_RasGAP",
    "PIP3_Gab1P_RasGAP",
    "E11P_RasGAP_RasGTP",
    "E12P_RasGAP_RasGTP",
    "E23P_RasGAP_RasGTP",
    "E24P_RasGAP_RasGTP",
    "E34P_RasGAP_RasGTP",
    "E44P_RasGAP_RasGTP",
    "E11P_Grb2_Gab1P_RasGAP_RasGTP",
    "E12P_Grb2_Gab1P_RasGAP_RasGTP",
    "E23P_Grb2_Gab1P_RasGAP_RasGTP",
    "E24P_Grb2_Gab1P_RasGAP_RasGTP",
    "E34P_Grb2_Gab1P_RasGAP_RasGTP",
    "E44P_Grb2_Gab1P_RasGAP_RasGTP",
    "E11P_ShcP_Grb2_Gab1P_RasGAP_RasGTP",
    "E12P_ShcP_Grb2_Gab1P_RasGAP_RasGTP",
    "E23P_ShcP_Grb2_Gab1P_RasGAP_RasGTP",
    "E24P_ShcP_Grb2_Gab1P_RasGAP_RasGTP",
    "E34P_ShcP_Grb2_Gab1P_RasGAP_RasGTP",
    "E44P_ShcP_Grb2_Gab1P_RasGAP_RasGTP",
    "PIP3_Gab1P_RasGAP_RasGTP",
    "E11P_Grb2_Gab1P_PI3K",
    "E12P_Grb2_Gab1P_PI3K",
    "E23P_Grb2_Gab1P_PI3K",
    "E24P_Grb2_Gab1P_PI3K",
    "E34P_Grb2_Gab1P_PI3K",
    "E44P_Grb2_Gab1P_PI3K",
    "E11P_ShcP_Grb2_Gab1P_PI3K",
    "E12P_ShcP_Grb2_Gab1P_PI3K",
    "E23P_ShcP_Grb2_Gab1P_PI3K",
    "E24P_ShcP_Grb2_Gab1P_PI3K",
    "E34P_ShcP_Grb2_Gab1P_PI3K",
    "E44P_ShcP_Grb2_Gab1P_PI3K",
    "PIP3_Gab1P_PI3K",
    "PIP2",
    "E23P_PI3K_PIP2",
    "E24P_PI3K_PIP2",
    "E34P_PI3K_PIP2",
    "E44P_PI3K_PIP2",
    "E11P_Grb2_Gab1P_PI3K_PIP2",
    "E12P_Grb2_Gab1P_PI3K_PIP2",
    "E23P_Grb2_Gab1P_PI3K_PIP2",
    "E24P_Grb2_Gab1P_PI3K_PIP2",
    "E34P_Grb2_Gab1P_PI3K_PIP2",
    "E44P_Grb2_Gab1P_PI3K_PIP2",
    "E11P_ShcP_Grb2_Gab1P_PI3K_PIP2",
    "E12P_ShcP_Grb2_Gab1P_PI3K_PIP2",
    "E23P_ShcP_Grb2_Gab1P_PI3K_PIP2",
    "E24P_ShcP_Grb2_Gab1P_PI3K_PIP2",
    "E34P_ShcP_Grb2_Gab1P_PI3K_PIP2",
    "E44P_ShcP_Grb2_Gab1P_PI3K_PIP2",
    "PIP3_Gab1P_PI3K_PIP2",
    "PIP3",
    "PTEN",
    "Akt",
    "AktP",
    "PP2A",
    "RafP_inactive",
    "RafPP_inactive",
    "PreduspmRNA",
    "duspmRNA",
    "DUSP",
    "DUSPP",
    "DUSP_ERKPP",
    "DUSP_ERKP",
    "DUSP_ERK",
    "DUSPP_ERKPP",
    "DUSPP_ERKP",
    "DUSPP_ERK",
    "GSK3b",
    "GSK3bP",
    "PrecmycmRNA",
    "cMycP",
    "cmycmRNA",
    "cMyc",
    "cMycPP",
]

for (idx, name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module


text:
module APipe

export @>, @>>, @.>>

include("manipulate.jl")
include("macros.jl")

end  # module


text:
# Contains structures and functions for simulating decays
struct DecaySimulationSettings
    # Configuration for simulating decay
    irf::TcspcTransient
    bck::TcspcTransient
    model::TcspcFitModel
    function DecaySimulationSettings(irf::TcspcTransient,
        bck::TcspcTransient,model::TcspcFitModel)
        # Ensure that the irf and background have the same axis
        if irf.t != bck.t
            error("irf and bck must have same time axis")
        end
        # Ensure that the linear model components are nonnegative
        if !isnonnegative(model)
            error("Linear coefficients must be nonnegative")
        end
        new(irf,bck,model)
    end
end

function sim_component(τ::Real,nphots::Integer,t::TcspcTimeAxis)
    # Simulates single component

    # Construct Exponential distribution
    d = Exponential(τ)
    phots = rand(d,nphots)
    # Fit data to histogram
    bins = time_bins(t)
    h = fit(Histogram,phots,bins)
    return h.weights
end

function sim_noise(nphots::Integer,t::TcspcTimeAxis)
    # Simulates random noise photons
    bins = time_bins(t)
    d = Uniform(bins[1],bins[end])
    phots = rand(d,nphots)
    h = fit(Histogram,phots,bins)
    return h.weights
end

function simulate_dk(settings::DecaySimulationSettings)
    # Generates a multi-exponential decay w/ IRF convolution

    # First, shift our IRF
    irf = shift(settings.irf,settings.model.Q.val)
    irf = TcspcTransient(settings.irf.t,irf)
    # Now, determine the number of photons we need per component
    M = model_vec(settings.model,irf,settings.bck)
    N = sum(M)
    # Multiply N by the fraction of the contribution to determine number of photons
    # Only the offset and the exponential components need this
    # To account for scattering and background effects, we have the IRF and background
    # already provided. Hence, we can just multiply the vector by the fraction, and add to the
    # binned histogram prior to Poisson weighting.
    τₘ = τₛ(settings.model)
    αₑ = τₘ[1:2:end,1]
    τₑ = τₘ[2:2:end,1]
    if nacomp(settings.model) > 0
        θₘ = θₛ(settings.model)
        αₑₘ = θₘ[1:2:end,[1, 4]]
        αₑ = vcat(αₑ,αₑₘ[:,1] .* αₑₘ[:,2])
        τₑₘ = θₘ[2:2:end,[1, 4]]
        τₑ = vcat(τₑ,((1 ./ τₑₘ[:,1]) + (1 ./ τₑₘ[:,2])) .^ (-1))
    end
    ncomp = length(τₑ)
    A = zeros(irf.t.nbins,ncomp)
    for i = 1:ncomp
        τtemp = τₑ[i]
        nphots = round(Int, N * αₑ[i])
        dkᵢ = sim_component(τtemp,nphots,irf.t)
        dkᵢ /= maximum(dkᵢ)
        A[:,i] = dkᵢ
    end
    # Now, convolve matrix in place
    convmat!(A,irf)
    # Now, multiply by effective weights
    dk = A[:,1:nlcomp(settings.model)] * αₑ[1:nlcomp(settings.model)]
    if nacomp(settings.model) > 0
        dk += A[:,(nlcomp(settings.model) + 1):end] *
            αₑ[(nlcomp(settings.model) + 1):end]
    end
    # Now simulate noise
    nphot = round(Int,settings.model.Z.val * N)
    dk += sim_noise(nphot,irf.t)
    # Now, account for scatter and background
    dk += irf.data * settings.model.S.val
    dk += settings.bck.data * settings.model.V.val
    # Finally, weight everything using Poisson
    sim_poisson(round.(Int,dk))
end

function simulate_dk(irf::TcspcTransient,bck::TcspcTransient,
    model::TcspcFitModel)
    # Wrapper for simulation
    sett = DecaySimulationSettings(irf,bck,model)
    simulate_dk(sett)
end

function simulate_dk(start::AbstractFloat,
    stop::AbstractFloat,nbins::Integer,
    ψ::AbstractFloat, irf::AbstractVector{<:Integer},
    bck::AbstractVector{<:Integer},
    param::AbstractVector{<:AbstractFloat},
    ntau::Integer,nani::Integer)

    # Simulates decays in public API

    # First, build time axis
    tax = TcspcTimeAxis(start,stop,nbins,ψ)
    # Next, build irf and bck
    irf = TcspcTransient(tax,irf)
    bck = TcspcTransient(tax,bck)
    # Finally, build generic model from params
    # Keep all parameters free - no fitting here
    fxparam = falses(size(param))
    gparam = falses(size(param))
    mod = TcspcFitModel(param,fxparam,gparam,ntau,nani,true)
    simulate_dk(irf,bck,mod)
end

function simulate_dk(tcal::AbstractFloat,nbins::Integer,
    ψ::AbstractFloat, irf::AbstractVector{<:Integer},
    bck::AbstractVector{<:Integer},
    param::AbstractVector{<:AbstractFloat},
    ntau::Integer,nani::Integer)

    # Simulates decays in public API

    # First, build time axis
    tax = TcspcTimeAxis(tcal,nbins,ψ)
    # Next, build irf and bck
    irf = TcspcTransient(tax,irf)
    bck = TcspcTransient(tax,bck)
    # Finally, build generic model from params
    # Keep all parameters free - no fitting here
    fxparam = falses(size(param))
    gparam = falses(size(param))
    mod = TcspcFitModel(param,fxparam,gparam,ntau,nani,true)
    simulate_dk(irf,bck,mod)
end


text:
# The content of this file was directly copied from Documenter.jl

# REPL doctest splitter.

const PROMPT_REGEX = r"^julia> (.*)$"
const SOURCE_REGEX = r"^       (.*)$"

function repl_splitter(code)
    lines = split(string(code, "\n"), '\n')
    input = String[]
    output = String[]
    buffer = IOBuffer() # temporary buffer for doctest inputs and outputs
    found_first_prompt = false
    while !isempty(lines)
        line = popfirst!(lines)
        prompt = match(PROMPT_REGEX, line)
        # We allow comments before the first julia> prompt
        !found_first_prompt && startswith(line, '#') && continue
        if prompt === nothing
            source = match(SOURCE_REGEX, line)
            if source === nothing
                savebuffer!(input, buffer)
                println(buffer, line)
                takeuntil!(PROMPT_REGEX, buffer, lines)
            else
                println(buffer, source[1])
            end
        else
            found_first_prompt = true
            savebuffer!(output, buffer)
            println(buffer, prompt[1])
        end
    end
    savebuffer!(output, buffer)
    zip(input, output)
end

function savebuffer!(out, buf)
    n = bytesavailable(seekstart(buf))
    n > 0 ? push!(out, rstrip(String(take!(buf)))) : out
end

function takeuntil!(r, buf, lines)
    while !isempty(lines)
        line = lines[1]
        if !occursin(r, line)
            println(buf, popfirst!(lines))
        else
            break
        end
    end
end


text:
module KissThreading

using Random: MersenneTwister
using Future: randjump
using Core.Compiler: return_type

export trandjump, TRNG, tmap, tmap!, tmapreduce, getrange

"""
    trandjump(rng = MersenneTwister(0); jump=big(10)^20)

Return a vector of copies of `rng`, which are advanced by different multiples
of `jump`. Effectively this produces statistically independent copies of `rng`
suitable for multi threading. See also [`Random.randjump`](@ref).
"""
function trandjump end

_randjump(rng, n, jump=big(10)^20) = accumulate(randjump, [jump for i in 1:n], init = rng)

function trandjump(rng = MersenneTwister(0); jump=big(10)^20)
    n = Threads.nthreads()
    rngjmp = Vector{MersenneTwister}(undef, n)
    for i in 1:n
        rngjmp[i] = randjump(rng, jump*i)
    end
    rngjmp
end

const TRNG = trandjump()

"""
    TRNG

A vector of statistically independent random number generators. Useful of threaded code:
```julia
rng = TRNG[Threads.threadid()]
rand(rng)
```
"""
TRNG

default_batch_size(n) = min(n, round(Int, 10*sqrt(n)))

struct Mapper
    atomic::Threads.Atomic{Int}
    len::Int
end

@inline function (mapper::Mapper)(batch_size, f, dst, src...)
    ld = mapper.len
    atomic = mapper.atomic
    Threads.@threads for j in 1:Threads.nthreads()
        while true
            k = Threads.atomic_add!(atomic, 1)
            batch_start = 1 + (k-1) * batch_size
            batch_end = min(k * batch_size, ld)
            batch_start > ld && break
            batch_map!(batch_start:batch_end, f, dst, src...)
        end
    end
    dst
end

@inline function batch_map!(range, f, dst, src...)
    @inbounds for j in range
        dst[j] = f(getindex.(src, j)...)
    end
end

function _doc_threaded_version(f)
    """Threaded version of [`$f`](@ref). The workload is divided into chunks of length `batch_size`
    and processed by the threads. For very cheap `f` it can be advantageous to increase `batch_size`."""
end

"""
    tmap!(f, dst::AbstractArray, src::AbstractArray...; batch_size=1)

$(_doc_threaded_version(map!))
"""
function tmap!(f, dst::AbstractArray, src::AbstractArray...; batch_size=1)
    ld = length(dst)
    if (ld, ld) != extrema(length.(src))
        throw(ArgumentError("src and dst vectors must have the same length"))
    end
    atomic = Threads.Atomic{Int}(1)
    mapper = Mapper(atomic, ld)
    mapper(batch_size, f, dst, src...)
end

"""
    tmap(f, src::AbstractArray...; batch_size=1)

$(_doc_threaded_version(map))
"""
function tmap(f, src::AbstractArray...; batch_size=1)
    g = Base.Generator(f,src...)
    T = Base.@default_eltype(g)
    dst = similar(first(src), T)
    tmap!(f, dst, src..., batch_size=batch_size)
end

struct MapReducer{T}
    r::Base.RefValue{T}
    atomic::Threads.Atomic{Int}
    lock::Threads.SpinLock
    len::Int
end

@inline function (mapreducer::MapReducer{T})(batch_size, f, op, src...) where T
    atomic = mapreducer.atomic
    lock = mapreducer.lock
    len = mapreducer.len
    Threads.@threads for j in 1:Threads.nthreads()
        k = Threads.atomic_add!(atomic, batch_size)
        k > len && continue
        y = f(getindex.(src, k)...)
        r = convert(T, y)
        range = (k + 1) : min(k + batch_size - 1, len)
        r = batch_mapreduce(r, range, f, op, src...)
        k = Threads.atomic_add!(atomic, batch_size)
        while k ≤ len
            range = k : min(k + batch_size - 1, len)
            r = batch_mapreduce(r, range, f, op, src...)
            k = Threads.atomic_add!(atomic, batch_size)
        end
        Threads.lock(lock)
        mapreducer.r[] = op(mapreducer.r[], r)
        Threads.unlock(lock)
    end
    mapreducer.r[]
end

"""
    tmapreduce(f, op, src::AbstractArray...; init, batch_size=default_batch_size(length(src[1])))

$(_doc_threaded_version(mapreduce))

Warning: In contrast to `Base.mapreduce` it is assumed that `op` must be commutative. Otherwise
the result is undefined.
"""
function tmapreduce end

function tmapreduce(f, op, src::AbstractArray...; init, batch_size=default_batch_size(length(src[1])))
    T = get_reduction_type(init, f, op, src...)
    _tmapreduce(T, init, batch_size, f, op, src...)
end

function tmapreduce(::Type{T}, f, op, src::AbstractArray...; init, batch_size=default_batch_size(length(src[1]))) where T
    _tmapreduce(T, init, batch_size, f, op, src...)
end

@inline function _tmapreduce(::Type{T}, init, batch_size, f, op, src...) where T
    lss = extrema(length.(src))
    lss[1] == lss[2] || throw(ArgumentError("src vectors must have the same length"))

    atomic = Threads.Atomic{Int}(1)
    lock = Threads.SpinLock()
    len = lss[1]
    mapreducer = MapReducer{T}(Base.RefValue{T}(init), atomic, lock, len)
    return mapreducer(batch_size, f, op, src...)
end

@inline function get_reduction_type(init, f, op, src...)
    Tx = return_type(f, Tuple{eltype.(src)...})
    Trinit = return_type(op, Tuple{typeof(init), Tx})
    Tr = return_type(op, Tuple{Trinit, Tx})
    Tr === Union{} ? typeof(init) : Tr
end

@inline function batch_mapreduce(r, range, f, op, src...)
    @inbounds for i in range
        r = op(r, f(getindex.(src, i)...))
    end
    r
end

"""
    getrange(n)

Partition the range `1:n` into `Threads.nthreads()` subranges and return the one corresponding to `Threads.threadid()`.
Useful for splitting a workload among multiple threads. See also the `TiledIteration` package for more advanced variants.
"""
function getrange(n)
    tid = Threads.threadid()
    nt = Threads.nthreads()
    d , r = divrem(n, nt)
    from = (tid - 1) * d + min(r, tid - 1) + 1
    to = from + d - 1 + (tid ≤ r ? 1 : 0)
    from:to
end

end


text:
using GridViewer
using Test

@testset "GridViewer.jl" begin
    @test 1 == 1
    @test 1 == 0
end

@testset "Some other tests" begin
    @test 2 == 2
    @test 2 == 0
end


text:
### functions for plotting gene expression data

### Load GP libraries for psuedotime
using GaussianProcesses
using Distributions

### include
include("gtex.jl")

#### islet sc rnaseq
function loadisletscrna(metafile="scrnaseq_human_islets_meta.tsv", countfile="scrnaseq_human_islets_normcounts.tsv.gz"; projdir=getprojectdir())
    meta = CSV.read(joinpath(projdir, "data", metafile), DataFrame)
    normcounts = CSV.read(joinpath(projdir, "data", countfile), DataFrame)
    (; meta, normcounts)
end

### plot boxplot of islet sc data
function scboxplot(gene, data ; kwargs...)
    ind = findfirst(occursin.(gene, data.normcounts.Gene))
    plot(title=string("GSE101207\n", gene), leg=false, xlabel="Islet Cell Type", ylabel="Norm. counts", fontfamily="helvetica" , xrotation=45; kwargs...)
    y = [v for v in data.normcounts[ind, data.meta.Label]]

    celllabels = Dict("pp" => "Pancreatic polypeptide", "delta" => "Delta", "psc" => "Pancreatic stellate", "duct" => "Duct")
    celltype = get.(Ref(celllabels), data.meta.CellType, data.meta.CellType)
    boxplot!(celltype, y, outliers=false, c=:grey, linecolor=:black)
    dotplot!(celltype, y, c=:steelblue, marker=(stroke(:white)))
end

function scboxplotgenes(genes, data ; kwargs...)
    phs = [scboxplot(g, islet_scrnaseq, size=(1200, 400)) for g in genes]
    for p in phs[2:end]
        plot!(p, ylabel="", left_margin=-2mm)
    end
    plot!(phs[1], ylabel="Norm. Counts", left_margin=5mm)

    p = plot(phs..., layout=(1, length(phs)), size=(1100, 300),  bottom_margin=10mm, top_margin=5mm, fontfamily="helvetica", titlefont=font(12, "helvetica"); kwargs...)
end

### beta cell maturation scRNA-seq
function load_betadiff_scrnaseq(metafile="scrnaseq_balboa_xin_cell_stage_meta.tsv", countfile="scrnaseq_balboa_xin_cell_stage_exp.tsv.gz"; projdir=getprojectdir())
    meta = CSV.read(joinpath(projdir, "data", metafile), DataFrame)
    normcounts = CSV.read(joinpath(projdir, "data", countfile), DataFrame)

    (; meta, normcounts)
end

## bar plot
function barbetadiff(gene, data ; kwargs...)
    ind = findfirst(data.normcounts.GeneName .== gene)
    meta = @subset(data.meta, :cell_type .∈ Ref(["Adult_Alpha", "Adult_Beta"]))
    y = [data.normcounts[ind, l] for l in meta.Label]

    stages = ["S5",  "S6",  "S7",  "M1",  "M3",  "M6",  "A"]
    stageorder = Dict(s => i for (i, s) in enumerate(stages))
    cc = cgrad(:cividis, 9, categorical=true)[2:end]

    groupedbar(meta.cell_type, y, group=getindex.(Ref(stageorder), meta.developmental_stage), lab=permutedims(stages), c=permutedims(cc), title=string("scRNA-seq endocrine cell maturation\n$(gene)"),
        ylabel="Norm. Counts", xlabel="Cell Type", fontfamily="helvetica"; kwargs...)
end

######### beta cell differentiation Pseudotime

function load_betadiff_sc(; projdir=getprojectdir())

    labelsfile = joinpath(projdir, "data", "singlecell_betadiff_labels.csv")
    tpmfile    = joinpath(projdir, "data", "singlecell_betadiff_tpm.csv.gz")

    data = CSV.read(tpmfile, DataFrame)
    rename!(data, ["Gene" ; string.("t", names(data)[2:end])])

    labels = CSV.read(labelsfile, DataFrame)
    labels.mid = (labels.start .+ labels.stop)/2;
    labels

    (data=data, labels=labels)

end

### beta cell differnentiation by stages, not currently used in figures
function load_beta_diff_geusz_xie(; projdir=getprojectdir())

    metafile = joinpath(projdir, "data", "rnaseq_betadiff_geusz_xie_meta.tsv")
    tpmfile  = joinpath(projdir, "data", "rnaseq_betadiff_geusz_xie_tpm.tsv.gz")

    meta = CSV.read(metafile, DataFrame)
    tpm = CSV.read(tpmfile, DataFrame)

    meta.Cell = replace.(meta.Cell, Ref("endocrine" => "EN"))

    (meta=meta, tpm=tpm)
end

function plotgenestage(gene, data ; kwargs...)

    ind = findall(occursin.(gene, data.tpm.Gene))

    if isempty(ind)
        error("Gene $gene not found")
    else
        plotgenestage(ind, data; kwargs...)
    end
end
function plotgenestage(ind::Vector{T}, data; maxnorm=false, vert=false, cc=1:length(ind), kwargs...) where {T}

    meta = data.meta
    tpm = data.tpm
    tdf = unique(sort(meta[!, [:CellOrder, :Cell]], :CellOrder))
    p = plot()

    if  length(ind) > 1
        ofv = range(.1, -.1, length=length(ind))
        s = 2
    else
        ofv = [0]
        s = 1
    end

    for (k, i) in enumerate(ind)
        x = [tpm[i, l] for l in meta.Label]
        if maxnorm
            x ./= maximum(x)
        end
        mg = combine(groupby(DataFrame(CellOrder=meta.CellOrder, X=x), :CellOrder), :X => mean => :Mean, :X => std => :Std)

        if vert

            scatter!(x, meta.CellOrder .+ meta.Offset/s .+ ofv[k], yticks=(tdf.CellOrder, tdf.Cell), c=cc[k], linecolor=cc[k], lab=tpm.GeneName[i], leg=:outertopright, marker=(1, stroke(0));  kwargs...)
            plot!(mg.Mean, mg.CellOrder .+ ofv[k], marker=:hline, xerror=mg.Std, c=cc[k], lab="")
        else
            scatter!(meta.CellOrder .+ meta.Offset, x, xticks=(tdf.CellOrder, tdf.Cell), c=cc[k], lab="", leg=:outertopright, marker=(stroke(0));  kwargs...)
            plot!(mg.CellOrder, mg.Mean, marker=:hline, yerror=mg.Std, lab=tpm.GeneName[i], c=cc[k])
        end

    end
    p
end

########## GaussianProcess Regression for sc beta diff psuedotime

datatrans(x, α=100, β = 1) = log(α*x + β)
invtrans(y, α=100, β=1) = (exp(y)-β)/α
sqd(x) = ifelse(x > 0, x*x, 0);
function gp_reg(t, y, st; initparams = (1.0, 0.0, 4.0), datatrans=datatrans, invtrans=invtrans, f_prior = Normal(1.4, 4.0), ℓ_prior = Normal(2., 2.0), n_prior=Normal(1.0, .75), label="")
    σf, ℓ, σn = initparams

    kern = Mat52Iso(ℓ, σf)
    set_priors!(kern, [ℓ_prior, f_prior])
    gp = GP(t, datatrans.(y), MeanZero(), kern, σn)
    set_priors!(gp.logNoise, [n_prior] )
    optimize!(gp)
    μ, v = predict_y(gp, st)
    sf  = invtrans.(μ)
    cil = invtrans.(μ .- 1.96*sqrt.(v))
    ciu = invtrans.(μ .+ 1.96*sqrt.(v))

    (gp=gp, μ=μ, v=v, sf=sf, cil=cil, ciu=ciu, t=t, y=y, st=st, label=label)
end

function gpgene(gene, data, dx=0.5)
    x = collect(@subset(data, :Gene .== "PseudoTime")[1, 2:end])
    y = collect(@subset(data, :Gene .== gene)[1, 2:end]);
    y[end] = mean(y[end-3:end-1]) ## final datapoint appears outlying, replace with mean of previous 3
    st = minimum(x):dx:maximum(x)
    gp_reg(x, y, st, label=gene)
end

function plotgpr(gpr ; c=:steelblue)
    plot()
    scatter!(gpr.t, gpr.y, c=c, lab="", marker=(stroke(0), 2))
    plot!(gpr.st, gpr.sf, ribbon=[gpr.sf - gpr.cil gpr.ciu - gpr.sf], lab=gpr.label, c=c)
end

function plotgpr!(gpr ; c=:steelblue, maxnorm=false, vert=true)

    if maxnorm
        m = maximum(gpr.sf)
    else
        m = 1
    end

    if vert

        scatter!(gpr.y/m, gpr.t, c=c, lab="", marker=(stroke(0), 2))
        plot!(gpr.sf/m, gpr.st, lab="", c=c)

        plot!(Shape([gpr.cil/m ; reverse(gpr.ciu/m)], [gpr.st ; reverse(gpr.st)]), fillalpha=0.2 ,line=stroke(0), c=c, lab=gpr.label)
    else

        scatter!(gpr.t, gpr.y/m, c=c, lab="", marker=(stroke(0), 2))
        plot!(gpr.st, gpr.sf/m, ribbon=[gpr.sf/m - gpr.cil/m gpr.ciu/m - gpr.sf/m], lab=gpr.label, c=c)
    end
end

function plotgprgenes(genes, betadiffsc ; kwargs...)
    phs = Plots.Plot[]
    for (i, g) in enumerate(genes)
        p = plot()
        plotgpr!(gpgene(g, betadiffsc.data), c=i, vert=false)
        p = xticks!(betadiffsc.labels.mid, betadiffsc.labels.cell)
        plot!(title=g)
        push!(phs, p)
    end
    plot!(phs[1], ylabel="Norm counts", left_margin=5mm)
    for i in 2:length(phs)
        plot!(phs[i], yformatter=y->"", left_margin=-3mm)
    end
    p = plot(phs..., layout=(1, length(phs)), size=((200*length(phs)) + 70, 195), link=:y, leg=false, grid=false, right_margin=-2mm, framestyle=:box, fontfamily="helvetica"; kwargs...)

end

function gp_hk1_gck(betadiffsc)
    plot(size=(325, 175), xticks=0:2, ylabel="Expression (au)", leg=:top, grid=false, xlabel="Differentiation Pseudotime")
    plotgpr!(gpgene("HK1", betadiffsc.data), maxnorm=true, vert=false)
    plotgpr!(gpgene("GCK", betadiffsc.data), c=:orange, maxnorm=true, vert=false)
    xticks!(betadiffsc.labels.mid, betadiffsc.labels.cell)
    plot!(fontfamily="helvetica", yticks=0:2)
end