text:
# Draw actomyosin networks, force components, and stress
# Alex Tam, 12/10/2020

"Draw the network configuration"
function draw_network(s, af, mm, parN, parA, Force, Lxx, Lxy, Lyx, Lyy)
    gr(); plot(); # Load GR plotting back-end and clear previous plots
    default(titlefont = (18, "times"), guidefont = (26, "times"), tickfont = (18, "times"))
    # Actin Filaments
    fx = zeros(Float64, parA.nSeg+1, length(af));
    fy = zeros(Float64, parA.nSeg+1, length(af)); # Pre-allocate un-translated filament node positions
    n1x = Vector{Float64}(); n1y = Vector{Float64}();
    n2x = Vector{Float64}(); n2y = Vector{Float64}(); # Pre-allocate translated node positions (non-plus end segments)
    p1x = Vector{Float64}(); p1y = Vector{Float64}();
    p2x = Vector{Float64}(); p2y = Vector{Float64}(); # Pre-allocate translated node positions (plus end segments)
    for f in af
        nx = Vector{}(); ny = Vector{}(); # Pre-allocate vectors of current filament nodes
        for seg in f.segments
            for j = 1:length(seg.t)
                if j == 1
                    mx, my, px, py = get_segment_nodes(f, seg, s, seg.t[j], Lxx, Lxy, Lyx, Lyy); # Get segment nodes
                    push!(nx, mx); push!(nx, px); # Store unprojected nodes of a segment (x)
                    push!(ny, my); push!(ny, py); # Store unprojected nodes of a segment (y)
                else
                    mx, my, px, py = get_segment_nodes(f, seg, s, seg.t[j], Lxx, Lxy, Lyx, Lyy); # Get segment nodes
                    if seg.index == length(f.segments)
                        push!(p1x, mx); push!(p2x, px); push!(p1y, my); push!(p2y, py); # Store nodes (plus-end)
                    else
                        push!(n1x, mx); push!(n2x, px); push!(n1y, my); push!(n2y, py); # Store nodes (non-plus-end)
                    end
                end
            end
        end
        nx = unique(nx); ny = unique(ny); # Remove duplicates
        fx[:, f.index] = nx; fy[:, f.index] = ny; # Store filament nodes in matrix
    end
    plot!(fx, fy, legend = false, color = "red", aspect_ratio=:equal, xlabel = L"$x$", ylabel = L"$y$", xlims = [0, parN.lxx], ylims = [0, parN.lyy]) # Plot un-translated segments
    scatter!(fx[1:end-1,:], fy[1:end-1,:], color = "red"); # Plot nodes as points
    scatter!(fx[end,:], fy[end,:], color = "black"); # Plot plus end nodes
    plot!(transpose(hcat(n1x, n2x)), transpose(hcat(n1y, n2y)), color = "green") # Plot translated segments (non-plus-end)
    plot!(transpose(hcat(p1x, p2x)), transpose(hcat(p1y, p2y)), color = "green") # Plot translated segments (plus end)
    scatter!([n1x, n2x], [n1y, n2y], color = "green"); # Plot translated nodes (non-plus-ends)
    scatter!([p1x], [p1y], color = "green"); # Plot translated nodes
    scatter!([p2x], [p2y], color = "black"); # Plot translated nodes (plus ends)
    # Myosin Motors
    m1x = Vector{Float64}(); tm1x = Vector{Float64}(); # Pre-allocate outside loop to plot all motors at once
    m1y = Vector{Float64}(); tm1y = Vector{Float64}();
    m2x = Vector{Float64}(); tm2x = Vector{Float64}();
    m2y = Vector{Float64}(); tm2y = Vector{Float64}();
    if all(s.mp[1] .<= 1)
        for m in mm
            x1, y1, x2, y2 = get_motor_pos(m, s, Lxx, Lxy, Lyx, Lyy); # Obtain motor positions (dimensionless, un-translated)
            # Convert to dimensional, translated positions
            mx1 = (x1 - m.t1[1])*Lxx + (y1 - m.t1[2])*Lyx;
            my1 = (x1 - m.t1[1])*Lxy + (y1 - m.t1[2])*Lyy;
            mx2 = (x2 - m.t2[1])*Lxx + (y2 - m.t2[2])*Lyx;
            my2 = (x2 - m.t2[1])*Lxy + (y2 - m.t2[2])*Lyy;
            # Compute translations for visualisation only
            myosin_translations = periodic([mx1/parN.lxx, my1/parN.lyy], [mx2/parN.lxx, my2/parN.lyy]);
            for j = 1:length(myosin_translations)
                if j == 1
                    push!(m1x, mx1); push!(m1y, my1); push!(m2x, mx2); push!(m2y, my2);
                else
                    t = myosin_translations[j];
                    # Store dimensional, translated positions
                    push!(tm1x, mx1 - (t[1]*Lxx + t[2]*Lyx));
                    push!(tm1y, my1 - (t[1]*Lxy + t[2]*Lyy));
                    push!(tm2x, mx2 - (t[1]*Lxx + t[2]*Lyx));
                    push!(tm2y, my2 - (t[1]*Lxy + t[2]*Lyy));
                end
            end
        end
        scatter!([m1x, m2x], [m1y, m2y], color = "blue"); # Plot binding sites as points
        scatter!([tm1x, tm2x], [tm1y, tm2y], color = "orange"); # Plot binding sites as points
        plot!(transpose(hcat(m1x, m2x)), transpose(hcat(m1y, m2y)), color = "blue") # Draw motors
        plot!(transpose(hcat(tm1x, tm2x)), transpose(hcat(tm1y, tm2y)), color = "orange") # Draw projected motors
    end
    # Principal stress directions
    Stress = [ Force[1]/Lyy Force[2]/Lyy; Force[3]/Lxx Force[4]/Lxx ]; # Compute stress tensor
    evals = eigvals(Stress); evecs = eigvecs(Stress);# Compute eigenvalues and eigenvectors
    if evals[1] >= 0
        plot!([parN.lxx/2 , parN.lxx/2 + evals[1]*evecs[1,1]], [parN.lyy/2 , parN.lyy/2 + evals[1]*evecs[1,2]], arrow = :arrow, color = "orange")
    else
        plot!([parN.lxx/2 , parN.lxx/2 + abs(evals[1])*evecs[1,1]], [parN.lyy/2 , parN.lyy/2 + abs(evals[1])*evecs[1,2]], arrow = :arrow, color = "blue")
    end
    if evals[2] >= 0
        plot!([parN.lxx/2 , parN.lxx/2 + evals[2]*evecs[2,1]], [parN.lyy/2 , parN.lyy/2 + evals[2]*evecs[2,2]], arrow = :arrow, color = "orange")
    else
        plot!([parN.lxx/2 , parN.lxx/2 + abs(evals[2])*evecs[2,1]], [parN.lyy/2 , parN.lyy/2 + abs(evals[2])*evecs[2,2]], arrow = :arrow, color = "blue")
    end
end

"Plot stress components"
function draw_stress(parN, Force, time_step, Lxx, Lyy)
    gr(); plot(); # Load GR plotting back-end and clear previous plots
    default(titlefont = (18, "times"), guidefont = (26, "times"), tickfont = (18, "times"))
    # Pre-allocate
    sxx = Vector{Float64}(); # Pre-allocate vector for stress component (xx)
    sxy = Vector{Float64}(); # Pre-allocate vector for stress component (xy)
    syx = Vector{Float64}(); # Pre-allocate vector for stress component (yx)
    syy = Vector{Float64}(); # Pre-allocate vector for stress component (yy)
    sxx_tot = Vector{Float64}(); # Pre-allocate vector for stress component (xx)
    sxy_tot = Vector{Float64}(); # Pre-allocate vector for stress component (xy)
    syx_tot = Vector{Float64}(); # Pre-allocate vector for stress component (yx)
    syy_tot = Vector{Float64}(); # Pre-allocate vector for stress component (yy)
    times = (0:time_step-1).*parN.dt; # Vector of times at which we obtain measurements
    for i = 1:time_step
        push!(sxx, Force[i][1]/Lyy); # Construct vector for stress component (xx)
        push!(sxy, Force[i][2]/Lyy); # Construct vector for stress component (xy)
        push!(syx, Force[i][3]/Lxx); # Construct vector for stress component (yx)
        push!(syy, Force[i][4]/Lxx); # Construct vector for stress component (yy)
        push!(sxx_tot, trapz(sxx, times)); # Integrate stress over time
        push!(sxy_tot, trapz(sxy, times)); # Integrate stress over time
        push!(syx_tot, trapz(syx, times)); # Integrate stress over time
        push!(syy_tot, trapz(syy, times)); # Integrate stress over time
    end
    # Plot stress
    p1 = plot(times, sxx, color = "orange", linewidth = 2, label=:false, xlabel = L"$t$", ylabel = "Stress (pN/μm)", legend=:bottomleft, size = (1000, 500));
    plot!(times, sxy, color = "purple", linewidth = 2, label=:false);
    plot!(times, syx, color = "green", linewidth = 2, label=:false);
    plot!(times, syy, color = "blue", linewidth = 2, label=:false);
    plot!(times, movingaverage(sxx, 100), color = "orange", linewidth = 2, linestyle = :dash, label = L"\sigma_{xx}");
    plot!(times, movingaverage(sxy, 100), color = "purple", linewidth = 2, linestyle = :dash, label = L"\sigma_{xy}");
    plot!(times, movingaverage(syx, 100), color = "green", linewidth = 2, linestyle = :dash, label = L"\sigma_{yx}");
    plot!(times, movingaverage(syy, 100), color = "blue", linewidth = 2, linestyle = :dash, label = L"\sigma_{yy}");
    # Plot integrated stress
    p2 = plot(times, sxx_tot, color = "orange", linewidth=2, label = L"\int \sigma_{xx} \; dt", xlabel = L"$t$", ylabel = "Integrated Stree (pNs/μm)", legend=:bottomleft, size = (1000, 500));
    plot!(times, sxy_tot, color = "purple", linewidth=2, label = L"\int \sigma_{xy} \;dt");
    plot!(times, syx_tot, color = "green", linewidth=2, label = L"\int \sigma_{yx} \;dt");
    plot!(times, syy_tot, color = "blue", linewidth=2, label = L"\int \sigma_{yy} \;dt");
    plot(p1, p2, layout = (1,2))
end

"Plot stress"
function draw_bulk_stress(parN, Force, time_step, Lxx, Lyy)
    gr(); plot(); # Load GR plotting back-end and clear previous plots
    default(titlefont = (18, "times"), guidefont = (26, "times"), tickfont = (18, "times"))
    Bulk_Stress = Vector{Float64}(); Bulk_Stress_Int = Vector{Float64}(); # Pre-allocate
    times = (0:parN.nT-1).*parN.dt; # Vector of times at which we obtain measurements
    for i = 1:parN.nT
        Stress = [ Force[i][1]/Lyy Force[i][2]/Lyy ; Force[i][3]/Lxx Force[i][4]/Lxx ]; # Compute stress tensor
        push!(Bulk_Stress, sum(eigvals(Stress))/2);
        push!(Bulk_Stress_Int, trapz(Bulk_Stress, times)); # Integrate force over time
    end
    Bulk_Stress_Moving = movingaverage(Bulk_Stress[1:time_step], 100);
    plot(times[1:time_step], Bulk_Stress[1:time_step], linewidth = 2, label = "Bulk Stress", xlabel = L"$t$", ylabel = "Bulk Stress (pN/μm)", legend=:topleft);
    plot!(times[1:time_step], Bulk_Stress_Moving[1:time_step], linewidth = 2, linestyle = :dot, label = "Moving Average");
    return times, Bulk_Stress, Bulk_Stress_Int
end

"Plot angle"
function draw_angle(parN, Theta, time_step)
    gr(); plot(); # Load GR plotting back-end and clear previous plots
    default(titlefont = (18, "times"), guidefont = (26, "times"), tickfont = (18, "times"))
    times = (0:parN.nT-1).*parN.dt; # Vector of times at which we obtain measurements
    plot(times[1:time_step], Theta[1:time_step], linewidth = 2, label =:false, xlabel = L"$t$", ylabel = L"$\theta$", ylims = [0, pi]);
end

"Integrate vector data using the trapezoidal rule"
function trapz(y, x)
    int = 0.0
    for i = 1:length(y)-1
        int += 0.5*(y[i] + y[i+1])*(x[i+1]-x[i]); # Composite trapezoidal rule
    end
    return int
end

"Moving average"
function movingaverage(x::Vector, numofele::Int)
    BackDelta = div(numofele,2)
    ForwardDelta = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1
    len = length(x)
    y = similar(x)
    for n = 1:len
        lo = max(1, n-BackDelta)
        hi = min(len, n+ForwardDelta)
        y[n] = mean(x[lo:hi])
    end
    return y
end


text:
# Issue is that type of err value depends on
# gradient method used
# Solutions
# Specialize: indomain/applytrackerr based on type of Omega

conjoinerror!(err::Ref{<:Real}, b) = err[] &= b

"Is `ω` in the domain of `x`?"
function applytrackerr(x, ω, errinit = softtrue())
  ω_ = tagerror(ω, errinit) # zt: surprisingly expensive
  fx = apl(x, ω_)
  (fx = fx, err = ω_.tags.err.x)
end

tagnotrackerr(ω) = tag(ω, (donttrack = true,))
applynotrackerr(x, ω) = apl(x, tagnotrackerr(ω))

"Soft `indomain`: distance from `ω` to the domain of `x`"
indomainₛ(x, ω, errinit = softtrue()) = applytrackerr(x, ω, errinit).err
indomainₛ(x::RandVar) = ciid(ω -> indomainₛ(x, ω))

"Is `ω` in the domain of `x`?"
indomain(x, ω, errinit = true) = applytrackerr(x, ω, errinit).err
indomain(x::RandVar) = ciid(ω -> indomain(x, ω))

function Omega.condf(tω::TaggedΩ, x, y)
  if !haskey(tω.tags, :donttrack) && haskey(tω.tags, :err)
    conjoinerror!(tω.tags.err, apl(y, tω))
  end
  apl(x, tω)
end

function Omega.cond(tω::TaggedΩ, bool)
  if !haskey(tω.tags, :donttrack) && haskey(tω.tags, :err)
    conjoinerror!(tω.tags.err, bool)
  end
end

"Tag `ω` with `err`"
# tagerror(ω, errinit::AbstractBool) = tag(ω, (err = Ref(errinit),)) #FIXME UNCOMMENT THIS AND SPECIALISE BELOW TO OMEGA
tagerror(ω, errinit::AbstractBool) = tag(ω, (err = Ref{Real}(errinit),))

"Tag tω with error"
function tagerror(tω::TaggedΩ, errinit::AbstractBool)
  # If error already there, use that!
  haskey(tω.tags, :err) ? tω : tag(tω, (err = Ref{Real}(errinit),))
end

ext: jl


text:
# Analyze the blow-up time of numerical simulations of the 3D inviscid
# Taylor-Green vortex for different dissipation-free SBP discretizations
# implemented in Trixi.jl:
# - DGSEM
# - Gauss collocation
# - Periodic FD SBP methods
# - TODO: One-block FD SBP methods?
#
# For optimal performance, run this script with
#   julia --threads=1 --check-bounds=no taylor_green_vortex.jl
# to let Julia use half of your `Sys.CPU_THREADS` cores. You can also specify
# the number of processes to use via
#   julia --threads=1 --procs=N --check-bounds=no taylor_green_vortex.jl
# where `N` is the number of processes you want to use.

# Package setup on one processor to avoid duplicating precompilation work etc.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.precompile()

using Distributed

# Add worker processes if none are present
if nprocs() == 1
  addprocs(Sys.CPU_THREADS ÷ 2)
end

# Package setup required on each processor
@everywhere begin
  using Pkg
  Pkg.activate(@__DIR__)
end
@everywhere using SummationByPartsOperators
@everywhere using Trixi
@everywhere using OrdinaryDiffEq

# Setup on the driver process
using CSV, DataFrames

# Further setup required on each processor
@everywhere begin
  # classical inviscid Taylor-Green vortex
  function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
    A  = 1.0 # magnitude of speed
    Ms = 0.1 # maximum Mach number

    rho = 1.0
    v1  =  A * sin(x[1]) * cos(x[2]) * cos(x[3])
    v2  = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
    v3  = 0.0
    p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
    p   = p + A^2 / 16 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3]))

    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
  end

  function blow_up_time(args...; abstol=1.0e-7, reltol=1.0e-7)

    # 3D compressible Euler equations with ratio of specific heats 7/5
    equations = CompressibleEulerEquations3D(1.4)

    # dissipation-free setup using EC fluxes
    surface_flux = flux_ranocha
    volume_flux  = flux_ranocha

    # use a try catch block since not all accuracy order of FD operators
    # work with all numbers of nodes
    local mesh, solver
    try
      solver, mesh = setup_solver_mesh(args..., surface_flux, volume_flux)
    catch e
      @info "solver setup faild" args e
      return NaN
    end

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

    tspan = (0.0, 20.0) # We used (0.0, 50.0) before but crashes are earlier
    ode = semidiscretize(semi, tspan)

    alive_callback = AliveCallback(alive_interval=1000)

    integrator = init(ode, RDPK3SpFSAL49(), abstol=abstol, reltol=reltol,
                      save_everystep=false, callback=alive_callback)

    try
      solve!(integrator)
    catch e
      @info "Blow-up" integrator.t args e
    end

    return integrator.t
  end

  function setup_solver_mesh(::Type{DGSEM}, polydeg, refinement_level,
                             surface_flux, volume_flux)

    solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux,
                   volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0, -1.0) .* π
    coordinates_max = ( 1.0,  1.0,  1.0) .* π
    mesh = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=refinement_level,
                    periodicity=(true, true, true),
                    n_cells_max=100_000)

    return solver, mesh
  end

  function setup_solver_mesh(::Type{GaussSBP}, polydeg, refinement_level,
                             surface_flux, volume_flux)

    solver = DGMulti(polydeg=polydeg, element_type=Hex(),
                     approximation_type=GaussSBP(),
                     surface_integral=SurfaceIntegralWeakForm(surface_flux),
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0, -1.0) .* π
    coordinates_max = ( 1.0,  1.0,  1.0) .* π
    cells_per_dimension = 2^refinement_level .* (1, 1, 1)
    mesh = DGMultiMesh(solver;
                       coordinates_min, coordinates_max, cells_per_dimension,
                       periodicity=(true, true, true))

    return solver, mesh
  end

  # CGSEM
  function setup_solver_mesh(::typeof(couple_continuously),
                             polydeg, refinement_level,
                             surface_flux, volume_flux)

    approximation_type = couple_continuously(
      legendre_derivative_operator(xmin=0.0, xmax=1.0, N=polydeg+1),
      UniformPeriodicMesh1D(xmin=0.0, xmax=1.0, Nx=2^refinement_level))
    solver = DGMulti(element_type=Hex(), approximation_type=approximation_type,
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0, -1.0) .* π
    coordinates_max = ( 1.0,  1.0,  1.0) .* π
    mesh = DGMultiMesh(solver; coordinates_min, coordinates_max)

    return solver, mesh
  end

  # Periodic FD
  function setup_solver_mesh(::typeof(periodic_derivative_operator),
                             accuracy_order, nnodes_per_dimension,
                             surface_flux, volume_flux)

    approximation_type = periodic_derivative_operator(
      derivative_order=1, accuracy_order=accuracy_order,
      xmin=0.0, xmax=1.0, N=nnodes_per_dimension)
    solver = DGMulti(element_type=Hex(), approximation_type=approximation_type,
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0, -1.0) .* π
    coordinates_max = ( 1.0,  1.0,  1.0) .* π
    mesh = DGMultiMesh(solver; coordinates_min, coordinates_max)

    return solver, mesh
  end

  # SBP FD
  function setup_solver_mesh(::typeof(derivative_operator),
                             accuracy_order, nnodes_per_dimension,
                             surface_flux, volume_flux)

    approximation_type = derivative_operator(
      MattssonNordström2004(),
      # MattssonAlmquistCarpenter2014Extended(),
      # MattssonAlmquistCarpenter2014Optimal(),
      # MattssonAlmquistVanDerWeide2018Minimal(),
      # MattssonAlmquistVanDerWeide2018Accurate(),
      derivative_order=1, accuracy_order=accuracy_order,
      xmin=0.0, xmax=1.0, N=nnodes_per_dimension)
    solver = DGMulti(element_type=Hex(), approximation_type=approximation_type,
                     surface_integral=SurfaceIntegralWeakForm(surface_flux),
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0, -1.0) .* π
    coordinates_max = ( 1.0,  1.0,  1.0) .* π
    cells_per_dimension = (1, 1, 1)
    mesh = DGMultiMesh(solver; coordinates_min, coordinates_max, cells_per_dimension)

    return solver, mesh
  end
end # @everywhere

# Prepare parameters
## DG
polydegs = 1:7
refinement_levels = 1:3
## FD
accuracy_orders = 2:2:10
nnodes_per_dimension = 4:2:12

# Run simulations and save data
dgsem_parameters = [(DGSEM, polydeg, refinement_level) for refinement_level in refinement_levels for polydeg in polydegs]
gauss_parameters = [(GaussSBP, polydeg, refinement_level) for refinement_level in refinement_levels for polydeg in polydegs]
cgsem_parameters = [(couple_continuously, polydeg, refinement_level) for refinement_level in refinement_levels for polydeg in polydegs]
fd_periodic_parameters = [(periodic_derivative_operator, accuracy_order, nnodes) for accuracy_order in accuracy_orders for nnodes in nnodes_per_dimension]

parameters = vcat(dgsem_parameters, gauss_parameters, cgsem_parameters,
                  fd_periodic_parameters)
blow_up_times = pmap(args -> blow_up_time(args...), parameters)
data = DataFrame(variant = map(x -> x[1], parameters),
                 polydeg_or_accuracy_order = map(x -> x[2], parameters),
                 refinement_level_or_nnodes = map(x -> x[3], parameters),
                 blow_up_time = blow_up_times)
CSV.write(joinpath(@__DIR__, "taylor_green_vortex.txt"), data)


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using ADCMEKit
Random.seed!(233)

reset_default_graph()
include("eikonal_op.jl")

m = 60
n = 30
h = 0.1

f = ones(n+1, m+1)
f[12:18, :] .= 10.
srcx = 5
srcy = 15

u = PyObject[]
for (k,(x,y)) in enumerate(zip(5*ones(Int64, length(1:n)), 1:n))
    push!(u,eikonal(f,x,y,h))
end

for (k,(x,y)) in enumerate( zip(1:m,5*ones(Int64, length(1:m))))
    push!(u,eikonal(f,x,y,h))
end

for (k,(x,y)) in enumerate( zip(1:m,25*ones(Int64, length(1:m))))
    push!(u,eikonal(f,x,y,h))
end

sess = Session()
uobs = run(sess, u)

F = Variable(ones(n+1, m+1))
# pl = placeholder(F0'[:])
# F = reshape(pl, n+1, m+1)

# F = Variable(ones(n+1, m+1))
u = PyObject[]
for (k,(x,y)) in enumerate(zip(5*ones(Int64, length(1:n)), 1:n))
    push!(u,eikonal(F,x,y,h))
end

for (k,(x,y)) in enumerate(zip(1:m,5*ones(Int64, length(1:m))))
    push!(u,eikonal(F,x,y,h))
end

for (k,(x,y)) in enumerate(zip(1:m,25*ones(Int64, length(1:m))))
    push!(u,eikonal(F,x,y,h))
end

# loss = sum([sum((uobs[i][5:5:end,55] - u[i][5:5:end,55])^2) for i = 1:length(u)])
loss = sum([sum((uobs[i][1:end,55] - u[i][1:end,55])^2) for i = 1:length(u)])

init(sess)
@show run(sess, loss)

# lineview(sess, F, loss, f, ones(n+1, m+1))
# gradview(sess, pl, loss, ones((n+1)* (m+1)))

# meshview(sess, pl, loss, F0'[:])

BFGS!(sess, loss, 1000,var_to_bounds=Dict(F=>(0.5,100.0)))

pcolormesh(run(sess, F))
colorbar()


text:
"""
    $(SIGNATURES)
Plots the structure of the factor graph. GraphPlot must be imported before DistributedFactoGraphs for these functions to be available.
Returns the plot context.

E.g.
```
using GraphPlot
using DistributedFactorGraphs, DistributedFactorGraphs.DFGPlots
# ... Make graph...
# Using GraphViz plotting
dfgplot(fg)
# Save to PDFG
using Compose
draw(PDF("/tmp/graph.pdf", 16cm, 16cm), dfgplot(fg))
```

More information at [GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl)
"""
function dfgplot(dfg::DistributedFactorGraphs.MetaGraphsDFG, p::DFGPlotProps = DFGPlotProps())
    @info "Deprecated, please use GraphsDFG or LightDFG."
    nodesize = [has_prop(dfg.g,i,:factor) ?  p.nodesize.fac : p.nodesize.var for i=vertices(dfg.g)]
    if p.drawlabels
        nodelabel = [has_prop(dfg.g,i,:factor) ? "" : string(get_prop(dfg.g,i,:label)) for i=vertices(dfg.g)]
    else
        nodelabel = nothing
    end
    nodefillc = [has_prop(dfg.g,i,:factor) ? p.nodefillc.fac : p.nodefillc.var for i=vertices(dfg.g)]

    gplot(dfg.g, nodelabel=nodelabel, nodesize=nodesize, nodefillc=nodefillc, layout=p.layout)

end

function gplot(dfg::DistributedFactorGraphs.MetaGraphsDFG; keyargs...)
    @info "Deprecated, please use GraphsDFG or LightDFG."
    gplot(dfg.g; keyargs...)
end


text:
function declare_old_executable_product(product_name)
    path_name = Symbol(string(product_name, "_path"))
    return quote
        # This is the old-style `withenv()`-based function
        function $(product_name)(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
            # We sub off to a shared function to avoid compiling the same thing over and over again
            return Base.invokelatest(
                JLLWrappers.withenv_executable_wrapper,
                f,
                $(Symbol("$(product_name)_path")),
                PATH[],
                LIBPATH[],
                adjust_PATH,
                adjust_LIBPATH,
            )
        end

        # This will eventually be replaced with a `Ref{String}`
        $(path_name) = ""
        function $(Symbol(string("get_", product_name, "_path")))()
            return $(path_name)::String
        end
    end
end

function declare_new_executable_product(product_name)
    @static if VERSION < v"1.6.0-DEV"
        return nothing
    else
        path_name = Symbol(string(product_name, "_path"))
        return quote
            # This is the new-style `addenv()`-based function
            function $(product_name)(; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
                env = Base.invokelatest(
                    JLLWrappers.adjust_ENV!,
                    copy(ENV),
                    PATH[],
                    LIBPATH[],
                    adjust_PATH,
                    adjust_LIBPATH,
                )
                return Cmd(Cmd([$(path_name)]); env)
            end

            # Signal to concerned parties that they should use the new version, eventually.
            #@deprecate $(product_name)(func) $(product_name)()
        end
    end
end

macro declare_executable_product(product_name)
    return excat(
        # We will continue to support `withenv`-style for as long as we must
        declare_old_executable_product(product_name),
        # We will, however, urge users to move to the thread-safe `addenv`-style on Julia 1.6+
        declare_new_executable_product(product_name),
    )
end

macro init_executable_product(product_name, product_path)
    path_name = Symbol(string(product_name, "_path"))
    return esc(quote
        # Locate the executable on-disk, store into $(path_name)
        global $(path_name) = joinpath(artifact_dir, $(product_path))

        # Add this executable's directory onto the list of PATH's that we'll need to expose to dependents
        push!(PATH_list, joinpath(artifact_dir, $(dirname(product_path))))
    end)
end


text:
# MLBase transformers.
module MLBaseWrapper
using DataFrames
using LinearAlgebra

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils
import AutoMLPipeline.AbsTypes: fit!, transform!

export fit!,transform!

export Standardize
export standardize, standardize!, estimate, transform,StandardScaler

"""
    StandardScaler(
       Dict(
          :center => true,
          :scale => true
       )
    )

Standardizes each feature using (X - mean) / stddev.
Will produce NaN if standard deviation is zero.
"""
mutable struct StandardScaler <: Transformer
  model
  args

  function StandardScaler(args=Dict())
    default_args = Dict(
      :center => true,
      :scale => true
    )
    new(nothing, mergedict(default_args, args))
  end
end

"""
    fit!(st::StandardScaler, features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}

Compute the parameters to center and scale.
"""
function fit!(st::StandardScaler, features::DataFrame, labels::Vector=[])
  mfeatures = convert(Matrix{Float64},features)
  pfeatures = mfeatures' |> collect |> Matrix{Float64}
  st_transform = estimate(Standardize, Array(mfeatures'); st.args...)
  st.model = Dict(
    :standardize_transform => st_transform
  )
end

"""
    transform!(st::StandardScaler, features::T)  where {T<:Union{Vector,Matrix,DataFrame}}

Apply the computed parameters for centering and scaling to new data.
"""
function transform!(st::StandardScaler, features::DataFrame)
  mfeatures = convert(Matrix{Float64},features)
  st_transform = st.model[:standardize_transform]
  pfeatures = mfeatures' |> collect |> Matrix{Float64}
  transposed_instances = Array(pfeatures)
  pres = transform(st_transform, transposed_instances)
  return (pres' |> collect |> Array{Float64}) |> DataFrame
end

### Standardization

"""
    Standardize(d::Int, m::Vector{Float64}, s::Vector{Float64})

Standardization type.
"""
struct Standardize
    dim::Int
    mean::Vector{Float64}
    scale::Vector{Float64}

    function Standardize(d::Int, m::Vector{Float64}, s::Vector{Float64})
        lenm = length(m)
        lens = length(s)
        lenm == d || lenm == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lens == d || lens == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new(d, m, s)
    end
end

indim(t::Standardize) = t.dim
outdim(t::Standardize) = t.dim

function transform!(y::DenseArray{YT,1}, t::Standardize, x::DenseArray{XT,1}) where {YT<:Real,XT<:Real}
    d = t.dim
    length(x) == length(y) == d || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(x, y)
                copy!(y, x)
            end
        else
            for i = 1:d
                @inbounds y[i] = x[i] * s[i]
            end
        end
    else
        if isempty(s)
            for i = 1:d
                @inbounds y[i] = x[i] - m[i]
            end
        else
            for i = 1:d
                @inbounds y[i] = s[i] * (x[i] - m[i])
            end
        end
    end
    return y
end

function transform!(y::DenseArray{YT,2}, t::Standardize, x::DenseArray{XT,2}) where {YT<:Real,XT<:Real}
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(x,2)
    size(y,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(x, y)
                copy!(y, x)
            end
        else
            for j = 1:n
                xj = view(x, :, j)
                yj = view(y, :, j)
                for i = 1:d
                    @inbounds yj[i] = xj[i] * s[i]
                end
            end
        end
    else
        if isempty(s)
            for j = 1:n
                xj = view(x, :, j)
                yj = view(y, :, j)
                for i = 1:d
                    @inbounds yj[i] = xj[i] - m[i]
                end
            end
        else
            for j = 1:n
                xj = view(x, :, j)
                yj = view(y, :, j)
                for i = 1:d
                    @inbounds yj[i] = s[i] * (xj[i] - m[i])
                end
            end
        end
    end
    return y
end

transform!(t::Standardize, x::DenseArray{T,1}) where {T<:AbstractFloat} = transform!(x, t, x)
transform!(t::Standardize, x::DenseArray{T,2}) where {T<:AbstractFloat} = transform!(x, t, x)

transform(t::Standardize, x::DenseArray{T,1}) where {T<:Real} = transform!(Array{Float64}(undef,size(x)), t, x)
transform(t::Standardize, x::DenseArray{T,2}) where {T<:Real} = transform!(Array{Float64}(undef,size(x)), t, x)

# estimate a standardize transform

function estimate(::Type{Standardize}, X::DenseArray{T,2}; center::Bool=true, scale::Bool=true) where T<:Real
    d, n = size(X)
    n >= 2 || error("X must contain at least two columns.")

    m = Array{Float64}(undef,ifelse(center, d, 0))
    s = Array{Float64}(undef,ifelse(scale, d, 0))

    if center
        fill!(m, 0.0)
        for j = 1:n
            xj = view(X, :, j)
            for i = 1:d
                @inbounds m[i] += xj[i]
            end
        end
        rmul!(m, 1.0 / n)
    end

    if scale
        fill!(s, 0.0)
        if center
            for j = 1:n
                xj = view(X, :, j)
                for i = 1:d
                    @inbounds s[i] += abs2(xj[i] - m[i])
                end
            end
        else
            for j = 1:n
                xj = view(X, :, j)
                for i = 1:d
                    @inbounds s[i] += abs2(xj[i])
                end
            end
        end
        for i = 1:d
            @inbounds s[i] = sqrt((n - 1) / s[i])
        end
    end

    return Standardize(d, m, s)
end

# standardize

function standardize(X::DenseArray{T,2}; center::Bool=true, scale::Bool=true) where T<:Real
    t = estimate(Standardize, X; center=center, scale=scale)
    Y = transform(t, X)
    return (Y, t)
end

function standardize!(X::DenseArray{T,2}; center::Bool=true, scale::Bool=true) where T<:AbstractFloat
    t = estimate(Standardize, X; center=center, scale=scale)
    Y = transform!(t, X)
    return (Y, t)
end

end


text:
import Gtk.GConstants.GdkModifierType, Gtk.GConstants.GdkModifierType.SHIFT

if Sys.isapple()
    global const PrimaryModifier = GdkModifierType.MOD2
    global const SecondaryModifer = GdkModifierType.CONTROL
end
if Sys.iswindows()
    global const PrimaryModifier = GdkModifierType.CONTROL
    global const SecondaryModifer = GdkModifierType.MOD1 #alt key
end
if Sys.islinux()
    global const PrimaryModifier = GdkModifierType.CONTROL
    global const SecondaryModifer = GdkModifierType.MOD1
end
global const NoModifier  = zero(typeof(PrimaryModifier))

mutable struct Action
    keyval::Integer
    state::Integer
    description::AbstractString

    Action(k::Integer, s::Integer) = new(k, s, "")
    Action(k::Integer, s::Integer, d::AbstractString) = new(k, s, d)
    Action(k::Integer, d::AbstractString) = new(k, NoModifier, d)
    Action(k::AbstractString, s::Integer, d::AbstractString) = new(keyval(k), s, d)
    Action(k::AbstractString, s::Integer) = new(keyval(k), s, "")
end

#https://developer.gnome.org/gtk3/unstable/checklist-modifiers.html
function doing(a::Action, event::Gtk.GdkEvent)

    mod = get_default_mod_mask()
    #on os x, the command key is also the meta key
    if Sys.isapple()
        if a.state == NoModifier && event.state == NoModifier
             return event.keyval == a.keyval
        end
        if (event.keyval == a.keyval) && (event.state & mod == a.state)
            return true
        end
        return (event.keyval == a.keyval) &&
               (event.state & mod == a.state + GdkModifierType.META)
    end

    return (event.keyval == a.keyval) && (event.state & mod == a.state)
end

function rightclick(event)
    mod = get_default_mod_mask()
    return event.button == 3 || (event.button == 1 && event.state & mod == SecondaryModifer)
end

global const Actions = Dict{AbstractString}{Action}()

Actions["save"]     = Action("s", PrimaryModifier, "Save file")
Actions["open"]     = Action("o", PrimaryModifier, "Open file")
Actions["closetab"] = Action("w", PrimaryModifier, "Close current tab")
Actions["newtab"]   = Action("n", PrimaryModifier, "New tab")
Actions["datahint"] = Action("D", PrimaryModifier+GdkModifierType.SHIFT, "Show data hint")
Actions["search"]   = Action(keyval("f"), PrimaryModifier, "Search")

Actions["runline"]  = Action(GdkKeySyms.Return, PrimaryModifier + GdkModifierType.SHIFT, "Execute current line")
Actions["runcode"]  = Action(GdkKeySyms.Return, PrimaryModifier, "Execute code")
Actions["runfile"]  = Action(GdkKeySyms.F5, "Run current file")

#TODO I should allow for several shortcuts
Actions["runline_kp"]  = Action(GdkKeySyms.KP_Enter, PrimaryModifier + GdkModifierType.SHIFT, "Execute current line")
Actions["runcode_kp"]  = Action(GdkKeySyms.KP_Enter, PrimaryModifier, "Execute code")

Actions["copy"]     = Action("c", PrimaryModifier, "Copy")
Actions["paste"]    = Action("v", PrimaryModifier, "Paste")
Actions["cut"]      = Action("x", PrimaryModifier, "Cut")
Actions["move_to_line_start"] = Action("a", SecondaryModifer, "Move cursor to line start")
Actions["move_to_line_end"]   = Action("e", SecondaryModifer, "Move cursor to line end")
Actions["delete_line"] = Action("k", PrimaryModifier, "Delete line")
Actions["duplicate_line"] = Action("d", PrimaryModifier, "Duplicate line")
Actions["toggle_comment"] = Action(0x02f, PrimaryModifier, "Toggle comment")
Actions["undo"] = Action("z", PrimaryModifier, "Undo")
Actions["redo"] = Action("Z", PrimaryModifier + GdkModifierType.SHIFT, "Redo")
Actions["goto_line"] = Action("g", PrimaryModifier, "Go to line")

Actions["select_all"] = Action("a", PrimaryModifier, "Select all")

##
Actions["extract_method"] = Action("e", PrimaryModifier, "Extract method")

#console
Actions["interrupt_run"] = Action("x", SecondaryModifer, "Interrupt current task")
Actions["clear_console"] = Action("k", PrimaryModifier, "Clear Console")

#main window

# § key when pressing cmd
Actions["console_editor_switch"] = Action(167, PrimaryModifier, "Switch focus between editor and console")


text:
using Random
using LinearAlgebra

export spqr_ssi

#==
SPQR_SSI block power method or subspace iteration applied to inv(R)
 to estimate rank, smallest singular values, and left/right singular vectors.

 [U, S, V, stats] = spqr_ssi(R, opts...)

 Uses the block power method or subspace iteration applied to the inverse of
 the triangular matrix R (implicitly) to estimate the numerical rank, smallest
 singular values, and left and right singular vectors R.  The algorithm will
 be efficient in the case the nullity of R is relatively low and inefficient
 otherwise.

 Input:

   R -- an n by n nonsingular triangular matrix
   opts (optional) -- see 'help spqr_rank_opts' for details.

 Output:

   U  -- n by k matrix containing estimates of the left singular vectors of R
       corresponding to the singular values in S. When stats.flag is 0 (and
       when the optional parameter opts.nsvals_large has its default value 1),
       r = n - k + 1 is the estimated rank of R. Also NT=U(:,2:k) is an
       orthonormal basis for the numerical null space of R'.
   S -- A k by k diagonal matrix whose diagonal entries are the estimated
       smallest k singular values of R, with k as described above.  For i=1:k,
       S(i,i) is an an estimate of singular value (r + i - 1) of R.
       Note that, unless stats.flag = 3 (see below), S(1,1)) > tol and for
       i =2:k, S(i,i) <= tol.
   V  -- n by k matrix containing estimates of the right singular vectors of R
       corresponding to the singular values in S.   When stats.flag is 0,
       N=V(:,2:k) is an orthonormal basis for the numerical null space of R.
   stats -- statistics, type 'help spqr_rank_stats' for details.

 Note that U' * R = S * V' (in exact arithmetic) and that R * V is
 approximately equal to U * S.

 output (for one or two output arguments): [s,stats] = spqr_ssi (...)
    s -- diag(S), with S as above.

 Example:
    R = sparse(matrixdepot("kahan", 100))
    U,S,V = spqr_ssi(R)
    norm_of_residual = opnorm( U' * R - diagm(S) * V' ) # should be near zero
     or
    U,S,V,stats = spqr_ssi(R, nothing, stats)
    N = V(:,2:end)   # orthonormal basis for numerical null space of R
    NT = U(:,2:end)  # orthonormal basis for numerical null space of R'
    norm_R_times_N = opnorm(R*N)
    norm_R_transpose_times_NT = opnorm(R'*NT)
     or
    opts = struct('tol',1.e-5,'nsvals_large',3)
    s,stats = spqr_ssi(R, opts, stats)
    stats  # information about several singular values

 See also spqr_basic, spqr_null, spqr_pinv, spqr_cod.

 Copyright 2012, Leslie Foster and Timothy A Davis.

 Outline of algorithm:
    let b = initial block size
    U = n by b random matrix with orthonormal columns
    repeat
       V1 = R \ U
       V,D1,X1 = svd(V1) = the compact svd of V1
       U1 = R' \ V
       U,D2,X2 = svd(U1) = the compact svd of U1
       sdot = diag(D2)
       s = 1 ./ sdot = estimates of singular values of R
       if s(end) <= tol
          increase the number of columns of U and repeat loop
       else
          repeat loop until stopping criteria (see code) is met
       end
    end
    k = smallest i with  sdot(i) > tol
    V = V*X2 (see code)
    V = first k columns of V in reverse order
    U = first k columns of U in reverse order
    s = 1 ./ (first k entries in sdot in reverse order)
    S = diag(s)

    Note for i = 1, 2, ..., k, S(i,i) is an upper bound on singular
    value (r + i - 1) of R in exact arithmetic but not necessarily
    in finite precision arithmetic.
==#

"""
start block power method to estimate the smallest singular
     values and the corresponding right singular vectors (in the
     r by nblock matrix V) and left singular vectors (in the r by
     nblock matrix U) of the n by n matrix R
"""
function spqr_ssi(R::AbstractMatrix{T}; nargout=4, opts...) where T<:Number
        ##      get_details::Int = 1,   # 0: basic statistics (default)
        # 1: detailed statistics:  basic stats, plus input options, time taken by
        #    various phases, statistics from spqr and spqr_rank subfunctions called,
        #    and other details.  Normally of interest only to the developers.
        # 2: basic statistics and a few additional statistics.  Used internally
        #    by some routines to pass needed information.
        ##      repeatable::Bool = true, # by default, results are repeatable
        ##      tol::Real = -1.0,        # a negative number means the default
        # tolerance should be computed
        ##      tol_norm_type::Int = 2, # 1: use opnorm(A, 1) to compute the default tol
                                        # 2: use normest(A, 0.01)
        ##      nsvals_large::Int = 1,  # default number of large singular values to estimate
        ##      min_block::Int = 3,     #
        ##      max_block::Int = 10,    #
        ##      min_iters::Int = 3,     #
        ##      max_iters::Int = 100,   #
        ##      nblock_increment = 5,   #
        ##      convergence_factor::Real = 0.1  #

    start_tic = time_ns()
    opts, get_details, repeatable, min_iters, max_iters, convergence_factor,
        nsvals_large, min_block, max_block, nblock_increment =
    get_opts(opts, :get_details, :repeatable, :ssi_min_iters, :ssi_max_iters,
                 :ssi_convergence_factor, :nsvals_large,
                 :ssi_min_block, :ssi_max_block, :ssi_nblock_increment)

    opts, tol, normest_R = get_tol_norm(opts, R)

    println("ssi(max_block=$max_block)")

    m, n = size(R)
    if m != n
        error("R must be square")
    end

    private_stream = repeatable ? MersenneTwister(1) : Random.GLOBAL_RNG

    #-------------------------------------------------------------------------------
    # adjust the options to reasonable ranges
    #-------------------------------------------------------------------------------

    # cannot compute more than n large singular values, where R is n-by-n
    nsvals_large = min(nsvals_large, n)

    # make the block size large enough so that there is a possiblity of
    # calculating all nsvals_large singular values
    max_block = max(max_block, nsvals_large)

    # max_block cannot be larger than n
    max_block = min(max_block, n)

    # min_block cannot be larger than n
    min_block = min(min_block, n)

    # start with nblock = min_block
    nblock = min_block
    opts = merge_opts(opts, nsvals_large=nsvals_large, ssi_max_block=max_block, ssi_min_block=min_block)
    stats = Statistics(real(T))
    stats.normest_A = normest_R
    stats.opts_used = opts.data
    stats.time_initialize = time_ns() - start_tic

    #-------------------------------------------------------------------------------
    # initializations
    #-------------------------------------------------------------------------------

    # stats.flag and stats.rank initialized in spqr_rank_get_inputs
    # set the order of the remaining stats fields
    stats.tol = tol
    stats.tol_alt = -1.0   # removed later if remains -1

    if get_details == 2
        stats.ssi_max_block_used = -1
        stats.ssi_min_block_used = -1
    end
    if (get_details == 1)
        stats.norm_R_times_N = -1
        stats.norm_R_transpose_times_NT = -1
        stats.iters = -1
        stats.nsvals_large_found = 0
        stats.final_blocksize = -1
        stats.ssi_max_block_used = -1
        stats.ssi_min_block_used = -1
        stats.time = 0
        stats.time_iters = 0
        stats.time_est_error_bounds = 0
        stats.time_svd = 0
    end

    stats.rank = 0
    if get_details == 1 || get_details == 2
        stats.ssi_max_block_used = max_block
        stats.ssi_min_block_used = nblock
    end

    get_details == 1 && ( start_iters_tic = time_ns() )

    U = randn(private_stream, n, nblock)
    U = qr(U).Q * Matrix{T}(I, n, nblock) # random orthogonal matrix with nblock cols
    # est_error_bound_calculated = 0      # set to 1 later if bound calculated
    flag_overflow = 0                     # set to 1 later if overflow occurs

    #-------------------------------------------------------------------------------
    # iterations
    #-------------------------------------------------------------------------------
    D2 = T[]
    X2 = T[]
    est_error_bound = T(0)

    if R isa SubArray{<:Number,2,<:AbstractSparseMatrix}
        R = sparse(R)
    end

    iters = 0
    while iters < max_iters
        iters += 1

        U0 = U
        V1 = R \ U
        if !all(isfinite.(V1))
            flag_overflow = 1
            break   # *************>>>>> early exit from for loop
        end

        get_details == 1 && ( time_svd = time_ns() )
        V, D1, X1 = svd(V1)
        get_details == 1 && ( stats.time_svd += time_ns() - time_svd )

        U1 = R' \ V
        # Note: with the inverse power method overflow is a potential concern
        if !all(isfinite.(U1))
            # *************>>>>> early exit from for loop
            # We know of no matrix that triggers this condition, so the next
            # two lines are untested.
            flag_overflow = 1      # untested
            break                   # untested
        end

        if get_details == 1
            time_svd = time_ns()
        end

        U, D2, X2 = svd(U1)

        if (get_details == 1)
            stats.time_svd += time_ns() - time_svd
        end

        k = findfirst(D2 .< tol^-1)
        if k == nothing && nblock == n
            # success since the numerical rank of R is zero
            break   # *************>>>>> early exit from for loop
        end

        errbound(k::Int) = norm(U0 * (X1 * (X2[:,k] ./ D1)) - U[:,k] ./ D2[k]) / sqrt(2)

        if k != nothing
            # k equals the calculated nullity + 1, k corresponds to
            # the first singular value of R that is greater than tol
            # reduce nsvals_large if necessary due to max_block limit:
            nsvals_large_old = nsvals_large
            nsvals_large = min(nsvals_large, max_block - k + 1)

            if nblock >= k + nsvals_large - 1
                # estimate error bound for singular value n - k + 1 of R:
                k2 = k + nsvals_large - 1
                est_error_bound = errbound(k)
                # When nsvals_large is 1, k and k2 are the same.  When
                # nsvals_large > 1, k2 corresponds to the largest singular
                # value of R that is returned.
                est_error_bound2 = k == k2 ? est_error_bound : errbound(k2)
            end

            # Note that
            #     [ 0     R ] [  U   ]  =    [ U0 * X1 * ( D1 \ X2 ) ]
            #     [ R'    0 ] [ V*X2 ]       [     V * ( X2 / D2 )   ]
            # It follows that, by Demmel's Applied Numerical Linear Algebra,
            # Theorem 5.5, some eigenvalue of
            #     B  =  [ 0    R ]
            #           [ R'   0 ]
            # will be within a distance of est_error_bound of 1 / d2(k), where s(1)
            # = 1 / d2(k) is our estimate of singular value n - k + 1 of R.
            # Typically the eigenvalue of B is singular value number n - k + 1 of
            # R, although this is not guaranteed.  est_error_bound is our estimate
            # of a bound on the error in using s(1) to approximate singular value
            # number n - k + 1 of R.

            if nblock >= k + nsvals_large - 1 &&
               est_error_bound <= convergence_factor * abs(1 / D2[k] - tol) &&
               est_error_bound2 <= convergence_factor * abs( 1/D2[k2] ) &&
               iters >= min_iters

                # Motivation for the tests:
                # The first test in the if statement is an attempt to insure
                # that nsvals_large singular values of R larger than tol are
                # calculated.
                #
                # Goal of the power method is to increase nblock until we find
                # sigma = (singular value n - k + 1 of R) is > tol.  If
                # this is true it is guaranteed that n - k + 1 is the
                # correct numerical rank of R.  If we let s(1) = 1 / d2(k)
                # then s(1) > tol.  However s(1) is only an estimate of sigma.
                # However, in most cases
                #     | s(1) - sigma | <= est_error_bound                   (1)
                # and to be conservative assume only that
                #  |s(1) - sigma|<= (1/convergence_factor)*est_error_bound  (2)
                # where convergence_factor<=1. By the second test in the if
                # statement
                #  est_error_bound <= convergence_factor * | s(1) - tol |   (3)
                # Equations (2) and (3) imply that
                #      | s(1) - sigma | <= | s(1) - tol |
                # This result and s(1) > tol imply that sigma > tol, as
                # desired.  Thus the second test in the if statement attempts
                # to carry out enough iterations to insure that the calculated
                # numerical rank is correct.
                #
                # The third test in the if statement checks on the accuracy of
                # the estimate for singular values n - k2 + 1.  Let sigma2 be
                # singular value n - k2 + 1 of R.  Usually it is true
                # that
                #     | s( k2 ) - sigma2 | <= est_error_bound2.             (4)
                # Assuming (4) and the third test in the if statement it
                # follows that the estimated relative
                # error in s(k2),  as measured by est_error_bound2 / s( k2) ,
                # is less that or equal to convergence_factor.  Therefore
                # the third test in the if statement attempts to insure that
                # the largest singular value returned by ssi has a relative
                # error <= convergence_factor.
                #
                # SUCCESS!!!, found singular value or R larger than tol
                break  # *************>>>>> early exit from for loop
            end
            nsvals_large = nsvals_large_old  # restore original value
        end

        if nblock == max_block && iters >= min_iters && k == nothing
            # reached max_block block size without encountering any
            # singular values of R larger than tolerance
            break    # *************>>>>> early exit from for loop
        end

        if 1 <= iters && iters < max_iters &&
            ( k == nothing  || ( k != nothing && nblock < k + nsvals_large - 1) )

            # increase block size
            nblock_prev = nblock
            if k == nothing
                nblock = min(nblock + nblock_increment, max_block)
            else
                nblock = min(k + nsvals_large - 1, max_block )
            end
            if nblock > nblock_prev
                Y = randn(private_stream, n, nblock-nblock_prev)
                Y = Y - U * ( U' * Y )
                Y = qr(Y).Q * Matrix{T}(I, n, nblock-nblock_prev)
                U = [U Y]
            end
        end
    end

    if get_details == 1
        stats.final_blocksize = nblock    # final block size
        stats.time_iters = time_ns() - start_iters_tic
        stats.iters = iters               # number of iterations taken in ssi
    end

    #-------------------------------------------------------------------------------
    # check for early return
    #-------------------------------------------------------------------------------

    if flag_overflow == 1
        warning("spqr_rank:overflow", "overflow in inverse power method")
        stats, U, S, V = spqr_failure(4, stats, get_details, start_tic)
        return nargout == 1 ? S : nargout == 2 ? (S, stats) : (U, S, V, stats)
    end

    #-------------------------------------------------------------------------------
    # determine estimated singular values of R
    #-------------------------------------------------------------------------------

    nkeep = k != nothing ? min(nblock, k + nsvals_large - 1) : nblock
    est_error_bounds = Vector{T}(undef, nkeep)

    if k != nothing
        # Note: in this case nullity = k - 1 and rank = n - k + 1
        est_error_bounds[nsvals_large] = est_error_bound
        numerical_rank = n - k + 1
    else
        k = nblock
        if nblock == n
            numerical_rank = 0
            # success since numerical rank is 0
        else
            # in this case rank not successfully determined
            # Note: In th/is case k is a lower bound on the nullity and
            #       n - k is an upper bound on the rank
            numerical_rank = n - k  #upper bound on numerical rank
            nsvals_large = 0 # calculated no singular values > tol
        end
    end

    stats.rank = numerical_rank

    if get_details == 1
        stats.nsvals_large_found = nsvals_large
    end

    #-------------------------------------------------------------------------------
    # adjust V so that R'*U = V*S (in exact arithmetic)
    #-------------------------------------------------------------------------------

    V = V * X2
    # reverse order of S, U and V and keep only nkeep singular vectors
    S = T(1) ./ D2[nkeep:-1:1]
    V = V[:,nkeep:-1:1]
    U = U[:,nkeep:-1:1]

    if get_details == 1
        t = time_ns()
    end

    # this recalculates est_error_bounds(nsvals_large)
    U0 = [ R * V - U * Diagonal(S); R' * U - V * Diagonal(S) ]
    est_error_bounds .= vec(sqrt.(sum(abs2, U0, dims=1) / 2))
    # Note that
    #    [ 0      R ]  [ U ]   -    [ U ] * S = [ R * V - U * S ]
    #    [ R'     0 ]  [ V ]   -    [ V ]       [ R'* U - V * S ].
    # It follows, by Demmel's Applied Numerical Linear Algebra,
    # Theorem 5.5, for i = 2, . . .,  k, some eigenvalue of
    #     B  =  [ 0    R ]
    #           [ R'   0 ]
    # will be within a distance of the norm of the ith column of
    # [ R * V - U * S R'*U - V*S] / sqrt(2) from S(i,i). Typically this
    # eigenvalue will be singular value number n - k + i of R,
    # although it is not guaranteed that the singular value number
    # is correct.  est_error_bounds(i) is the norm of the ith column
    # of [ R * V - U * S; R' * U - V * S ] / sqrt(2).

    if get_details == 1
        # Note that stats.time_est_error_bounds includes the time for the error
        #   bound calculations done outside of the subspace iterations loop
        stats.time_est_error_bounds = time_ns() - t
    end

    nr1 = numerical_rank - nsvals_large
    stats.sval_numbers_for_bounds = nr1 + 1 : nr1 + nkeep

    #-------------------------------------------------------------------------------
    # compute norm R*N and R'*NT
    #-------------------------------------------------------------------------------

    # compute norm R*N where N = V(:,nsvals_large+1:end) is the approximate
    #         null space of R and
    #         norm R'*NT where NT = U(:,nsvals_large+1:end) is the approximate
    #         null space of R'

    if get_details == 1
        t = time_ns()
    end

    norm_R_times_N = opnorm(R * V)
    # svals_R_transpose_times_NT = svd(R'*U)'
    norm_R_transpose_times_NT = opnorm(R' * U)

    if get_details == 1
        stats.norm_R_times_N = norm_R_times_N
        stats.norm_R_transpose_times_NT = norm_R_transpose_times_NT
        stats.time_svd += time_ns() - t
    end

    # Note: norm_R_times_N is an upper bound on sing. val. rank1+1 of R
    # and norm_R_transpose_times_NT is also an upper bound on sing. val.
    #      rank1+1 of R
    max_norm_RN_RTNT = max(norm_R_times_N, norm_R_transpose_times_NT)
    # choose max here to insure that both null spaces are good

    #-------------------s------------------------------------------------------------
    # determine flag indicating the accuracy of the rank calculation
    #-------------------------------------------------------------------------------

    if numerical_rank == 0
        # numerical rank is 0 in this case
        stats.flag = 0
    elseif (numerical_rank == n  || max_norm_RN_RTNT <= tol) &&
           (nsvals_large > 0 && S[nsvals_large] - est_error_bounds[nsvals_large] > tol )
        # in this case, assuming est_error_bounds(nsvals_large) is a true
        # error bound, then the numerical rank is correct.  Also
        # N = V(:,nsvals_large+1:end) and NT = U(:,nsvals_large+1:end)
        # are bases for the numerical null space or R and R', respectively
        stats.flag = 0
    elseif ( nsvals_large > 0 && numerical_rank == n ) ||
            ( nsvals_large > 0 &&
             S[nsvals_large] - est_error_bounds[nsvals_large] > max_norm_RN_RTNT )
        # in this case, assuming est_error_bounds(nsvals_large) is a true
        # error bound, then the numerical rank is correct with a modified
        # tolerance.  This is a rare case.
        stats.flag = 1
        tol_alt = ( S[nsvals_large] - est_error_bounds[nsvals_large] )
        tol_alt -= eps(tol_alt)         # so that tol_alt satisfies the >
                                        # portion of the inequality below
        # tol_alt = max_norm_RN_RTNT
        # Note: satisfactory values of tol_alt are in the range
        #    S[nsvals_large] - est_error_bounds[nsvals_large] > tol_alt
        #    >= max_norm_RN_RTNT
        stats.tol_alt = tol_alt
    elseif  nsvals_large > 0 && S[nsvals_large] > tol &&
            max_norm_RN_RTNT <= tol &&
            S[nsvals_large] - est_error_bounds[nsvals_large] <= max_norm_RN_RTNT
        # in this case, assuming est_error_bounds(nsvals_large) is a true
        # error bound, the error bound is too large to determine the
        # numerical rank.  This case is very rare.
        stats.flag = 2
    else
        # in this case all the calculated singular values are
        # smaller than tol or either N is not a basis for the numerical
        # null space of R with tolerance tol or NT is not such a basis for R'.
        # stats.rank is an upper bound on the numerical rank.
        stats.flag = 3
    end

    #-------------------------------------------------------------------------------
    # return results
    #-------------------------------------------------------------------------------
    stats.est_svals_of_R = S
    stats.est_error_bounds = est_error_bounds

    if get_details == 1
        stats.time = time_ns() - start_tic
    end

    nargout == 1 ? S : nargout == 2 ? (S, stats) : (U, S, V, stats)
end


text:
# # potentially create subtype of SimpleDiGraph that has fadjlist sorted

function sort_by_deg_or!(DG::SimpleDiGraph)
    for i in vertices(DG) # for each node
        # # order its out-neighbor list from lowest degree to highest degree
        sort!(DG.fadjlist[i], lt=(x, y) -> (degree(DG, x) < degree(DG, y)) ||(degree(DG, x) == degree(DG, y) && x < y))
    end
end

"""
    toDAG

Returns degree oriented DAG
"""
function toDAG(G::SimpleGraph)
    n = nv(G)

    # # SimpleDiGraphs contain both out- and in-degree information
    DG = SimpleDiGraph(n)

    for i in vertices(G) # for each node
        for j in neighbors(G, i) # for each neighbor

            # # if edge has already been seen, skip processing
            if i > j
                continue
            end

            degi = degree(G, i)
            degj = degree(G, j)

            # # conditions for adding a directed edge
            if (degi < degj) || (degi == degj && i < j)
                add_edge!(DG, i, j)
            else
                add_edge!(DG, j, i)
            end
        end
    end

    # for i in vertices(DG) # for each node
    #     # # order its out-neighbor list from lowest degree to highest degree
    #     sort!(DG.fadjlist[i], lt=(x, y) -> (degree(DG, x) < degree(DG, y)) ||(degree(DG, x) == degree(DG, y) && x < y))
    # end

    sort_by_deg_or!(DG)
    # println(DG.fadjlist)
    # println(DG.badjlist)

    return DG

end


text:
using ChemometricsData
using Test

@testset "ChemometricsData.jl" begin
    #unfortunately there are no unit tests at this time...
end


text:
"""
SamplerState for standard SGD

Fields
    - x: the current state
    - teration counter
    - stepsize: s.stepsize(s.t) gives the stepsize to be used in iteration s.t

Constructors
    construct using `SGDState(x::Array{Float64},stepsize::Function)` or `SGDState(x::Array{Float64},stepsize::Float64)`
"""
type SGDState <: SamplerState
    x::Array{Float64}
    t::Int64
    stepsize::Function
    function SGDState(x::Array{Float64},stepsize::Function)
        nparams = size(x)
        new(x,0,stepsize)
    end
    function SGDState(x::Array{Float64},stepsize::Float64)
        eps = stepsize
        stepsize = t::Int64 -> eps
        new(x,0,stepsize)
    end
end

function sample!(s::SGDState, grad::Function)
  gt = grad(s.x)
  s.t += 1
  Base.LinAlg.axpy!(s.stepsize(s.t),grad(s.x),s.x)
  s
end


text:
type DCOut
	df::DataFrame
	dims::Vector{Dimension}
	groups::Vector{Group}
	charts::Vector{DCChart}
	widgets::Vector{DCWidget}
	output_id::Integer
	elastic_height::Bool
	html_debug::Bool

	function DCOut(df::DataFrame)
		retval = new()
		retval.df = df
		retval.dims = Dimension[]
		retval.groups = Group[]
		retval.charts = DCChart[]
		retval.widgets = DCWidget[]
		retval.output_id = rand(0:999999)
		retval.elastic_height = false
		retval.html_debug = false
		retval
	end
end

"""
	get_group_by_name

Returns the group inside the given DCCout instance with the given name.
"""
function get_group_by_name(dcout::DCOut, name::ASCIIString)
	results = find(x -> x.name == name, dcout.groups)
	if length(results) == 0
		error("group \"", name, "\" not found")
	elseif length(results) == 1
		return dcout.groups[results[1]]
	else
		error("dcout in invalid state: two groups have the same name")
	end
end

"""
	get_groups_by_col

Returns the groups inside the given DCCout instance associated with
the given column.
"""
function get_groups_by_col(dcout::DCOut, col::Symbol)
	idxs = find(x -> x.dim.name == col, dcout.groups)
	result = Group[]
	for idx in idxs
		push!(result, dcout.groups[idx])
	end
	result
end

"""
	get_group_by_name

Returns the dimension inside the given DCCout instance created from the given column.
"""
function get_dim_by_col(dcout::DCOut, col::Symbol)
	results = find(x -> x.name == col, dcout.dims)
	if length(results) == 0
		error("dimension from column \"", col, "\" not found")
	elseif length(results) == 1
		return dcout.dims[results[1]]
	else
		error("dcout in invalid state: two dimensions have the same name")
	end
end

"""
	get_charts

Returns all charts constructed from the given column of the given type.
`chart_type` can be: piechart, barchart, linechart, rowchart, bubblechart
"""
function get_charts(dcout::DCOut, col::Symbol, chart_type::ASCIIString)
	idxs = find(x -> (x.group.dim.name == col && uppercase(chart_type) == uppercase(x.typ.concreteName)), dcout.charts)
	result = DCChart[]
	for idx in idxs
		push!(result, dcout.charts[idx])
	end
	result
end

"""
	add_dimension!

Append the Dimension to the list of dimensions in the DCOut object.
"""
function add_dimension!(dcout::DCOut, dim::Dimension)
	if findfirst(x -> x.name == dim.name, dcout.dims) > 0
		error(string("attempt to add dimension \"", dim.name, "\" failed: a dimension with that name already exists"))
	end
	push!(dcout.dims, dim)
	dcout
end

"""
	add_chart!

Append the DCChart to the list of charts in the DCOut object.
"""
function add_chart!(dcout::DCOut, chart::DCChart)
	push!(dcout.charts, chart)
	dcout
end

"""
	add_widget!

Append the DCWidget to the list of widgets in the DCOut object.
"""
function add_widget!(dcout::DCOut, widget::DCWidget)
	push!(dcout.widgets, widget)
	dcout
end

"""
	add_group!

Append the Group to the list of groups in the DCOut object.
"""
function add_group!(dcout::DCOut, group::Group)
	if length(find(x -> x.name == group.name, dcout.groups)) > 0
		error(string("attempt to add group \"", group.name, "\" failed: a group with that name already exists"))
	end
	push!(dcout.groups, group)
	dcout
end

"""
	remove_group!

Removes the given Group from the list of groups in the DCOut object.
"""
function remove_group!(dcout::DCOut, group::Group)
	targets = find(x -> x.name == group.name, dcout.groups)
	if length(targets) == 0
		error("group \"", group.name, "\" does not exist")
	elseif length(targets) > 1
		# This case should technically not happen, but if it is the case,
		# we remove all groups that match `group`
		deleteat!(dcout.groups, targets[1])
		remove_group!(dcout, group)
	else
		deleteat!(dcout.groups, targets[1])
	end
	dcout
end

"""
	remove_chart!

Removes the given DCChart from the list of charts in the DCOut object.
"""
function remove_chart!(dcout::DCOut, chart::DCChart)
	targets = find(x -> x.parent == chart.parent, dcout.charts)
	if length(targets) == 0
		error("chart \"", chart.parent, "\" does not exist")
	elseif length(targets) > 1
		deleteat!(dcout.charts, targets[1])
		remove_group!(dcout, group)
	else
		deleteat!(dcout.charts, targets[1])
	end
	dcout
end

"""
	quick_add!

A utility function for quickly building a chart and adding it.
Works with linechart, barchart, rowchart, piechart.
Requires previously a constructed group.
"""
function quick_add!(dcout::DCOut, group::Group, chart_constructor::Function)
	new_chart = chart_constructor(dcout.df[group.dim.name], group)
	add_chart!(dcout, new_chart)
end
"""
	quick_add!

Quickly build a chart using the group constructed from `column`.
If multiple groups are found, an error is thrown unless `use_first`
is set, in which case the first group found is used.
"""
function quick_add!(dcout::DCOut, column::Symbol, chart_constructor::Function, use_first::Bool = false)
	groups = get_groups_by_col(dcout, column)
	if length(groups) == 0
		error("group from column \"", column, "\" not found")
	elseif length(groups) > 1
		if use_first
			quick_add!(dcout, groups[1], chart_constructor)
		else
			error("unable to infer group: multiple groups use column \"", column, "\"")
		end
	else
		quick_add!(dcout, groups[1], chart_constructor)
	end
end
"""
	quick_add!

Quickly build a chart using the group with the given name.
"""
function quick_add!(dcout::DCOut, group_name::ASCIIString, chart_constructor::Function)
	group = get_group_by_name(dcout, group_name)
	quick_add!(dcout, group, chart_constructor)
end

"""
	add_datacountwidget!

Construct and append a Data Count Widget to the DCOut object.
"""
function add_datacountwidget!(dcout::DCOut)
	add_widget!(dcout, datacountwidget())
	dcout
end

"""
	add_datatablewidget!

Construct and append a Data Table Widget to the DCOut object. Requires
a previously constructed dimension for the first column.
"""
function add_datatablewidget!(dcout::DCOut)
	columns = Symbol[]
	for dim in dcout.dims
		push!(columns, dim.name)
	end
	add_widget!(dcout, datatablewidget(columns))
	dcout
end

"""
	add_bubblechart!

Construct and append a Bubble Chart to the DCOut object. Requires a
previously created dimension. x_col, y_col, and r_col denote the
DataFrame fields whose sums will determine the x position, y position, and
radius of the bubbles in the final chart.

Use `:DCCount` to access the count field.
"""
function add_bubblechart!(dcout::DCOut, dim::Dimension, x_col::Symbol, y_col::Symbol, r_col::Symbol)
	cols = [x_col, y_col, r_col]
	for col in cols
		if col != :DCCount
			get_dim_by_col(dcout, col)
		end
	end
	group = reduce_master(dim, [x_col, y_col, r_col])
	add_group!(dcout, group)
	add_chart!(dcout, bubblechart(group, x_col, y_col, r_col))
	dcout
end
function add_bubblechart!(dcout::DCOut, dim_col::Symbol, x_col::Symbol, y_col::Symbol, r_col::Symbol)
	add_bubblechart!(dcout, get_dim_by_col(dcout, dim_col), x_col, y_col, r_col)
end

"""
	clear_charts!

Remove all charts from the DCOut object.
"""
function clear_charts!(dcout::DCOut)
	dcout.charts = DCChart[]
	dcout
end

"""
	clear_widgets

Remove all widgets from the DCOut object.
"""
function clear_widgets!(dcout::DCOut)
	dcout.widgets = DCWidget[]
	dcout
end

"""
	set_elastic_height!

Changes the elasticity of the output div height. A value of true
indicates that the div will resize to fit all elements without
needing a scrollbar.
"""
function set_elastic_height!(dcout::DCOut, elastic_height::Bool)
	dcout.elastic_height = elastic_height
	dcout
end

"""
	randomize_ids!

Randomly re-initialize all dcout random ids.
These are used when exporting to HTML to prevent charts from
referencing one another across IJulia cells.
"""
function randomize_ids!(dcout::DCOut)
	dcout.output_id = rand(0:999999)

  for chart in dcout.charts
    randomize_parent(chart)
  end
  for widget in dcout.widgets
    randomize_parent(widget)
  end
end

"""
	infer_dimensions!

Infer a dimension for each column in the DCOut DataFrame
"""
function infer_dimensions!(dcout::DCOut)
	for name in names(dcout.df)
		arr = dcout.df[name]
		if can_infer_chart(arr)
			dim = infer_dimension(arr, name)
			add_dimension!(dcout, dim)
		end
	end
	dcout
end

"""
	infer_groups!

Infer a group for each dimension in the DCOut DataFrame
"""
function infer_groups!(dcout::DCOut)
	for dim in dcout.dims
		arr = dcout.df[dim.name]
		if can_infer_chart(arr)
			group = infer_group(arr, dim)
			add_group!(dcout, group)
		end
	end
	dcout
end

"""
	dc(df::DataFrame)

Construct a DC.js visualization based on the columns in the given DataFrame.
This is the easiest and most straightforward way to use DC.jl.
"""
function dc(df::DataFrame)
	dcout = DCOut(df)

	for name in names(dcout.df)
		arr = dcout.df[name]
		if can_infer_chart(arr)
			dim = infer_dimension(arr, name)
			add_dimension!(dcout, dim)

			group = infer_group(arr, dim)
			add_group!(dcout, group)

			chart = infer_chart(arr, group)
			add_chart!(dcout, chart)
		end
	end

	dcout
end


text:
using Base64

"""
    run(doc::WeaveDoc; doctype = :auto,
        mod::Union{Module, Symbol} = :sandbox, out_path=:doc,
        args=Dict(), fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache = :off, throw_errors=false)

Run code chunks and capture output from parsed document.

* `doctype`: :auto = set based on file extension or specify one of the supported formats.
  See `list_out_formats()`
* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document, `:pwd`: Julia working directory,
  `"somepath"`: Path as a AbstractString e.g `"/home/mpastell/weaveout"`
* `args`: dictionary of arguments to pass to document. Available as WEAVE_ARGS.
* `mod`: Module where Weave `eval`s code. Defaults to `:sandbox`
   to create new sandbox module, you can also pass a module e.g. `Main`.
* `fig_path`: where figures will be generated, relative to out_path
* `fig_ext`: Extension for saved figures e.g. `".pdf"`, `".png"`. Default setting depends on `doctype`.
* `cache_path`: where of cached output will be saved.
* `cache`: controls caching of code: `:off` = no caching, `:all` = cache everything,
  `:user` = cache based on chunk options, `:refresh`, run all code chunks and save new cache.

**Note:** Run command from terminal and not using IJulia, Juno or ESS, they tend to mess with capturing output.
"""
function Base.run(doc::WeaveDoc; doctype = :auto,
        mod::Union{Module, Symbol} = :sandbox, out_path=:doc,
        args=Dict(), fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache = :off, throw_errors=false)
    #cache :all, :user, :off, :refresh

    doc.cwd = get_cwd(doc, out_path)
    doctype == :auto && (doctype = detect_doctype(doc.source))
    doc.doctype = doctype
    doc.format = formats[doctype]
    isdir(doc.cwd) || mkpath(doc.cwd)

    if occursin("2pdf", doctype) && cache == :off
        fig_path = mktempdir(abspath(doc.cwd))
    elseif occursin("2html", doctype)
        fig_path = mktempdir(abspath(doc.cwd))
    end

    cache == :off || @eval import Serialization

    #This is needed for latex and should work on all output formats
    Sys.iswindows() && (fig_path = replace(fig_path, "\\" => "/"))

    doc.fig_path = fig_path
    set_rc_params(doc.format.formatdict, fig_path, fig_ext)

    #New sandbox for each document with args exposed
    if mod == :sandbox
        sandbox = "WeaveSandBox$(rcParams[:doc_number])"
        mod = Core.eval(Main, Meta.parse("module $sandbox\nend"))
    end
    @eval mod WEAVE_ARGS = Dict()
    merge!(mod.WEAVE_ARGS, args)

    rcParams[:doc_number] += 1

    if haskey(doc.format.formatdict, :mimetypes)
      mimetypes = doc.format.formatdict[:mimetypes]
    else
      mimetypes = default_mime_types
    end

    report = Report(doc.cwd, doc.basename, doc.format.formatdict, mimetypes, throw_errors)
    pushdisplay(report)

    if cache != :off && cache != :refresh
        cached = read_cache(doc, cache_path)
        cached == nothing && @info("No cached results found, running code")
    else
        cached = nothing
    end

    executed = Any[]
    n = length(doc.chunks)

    for i = 1:n
        chunk = doc.chunks[i]

        if isa(chunk, CodeChunk)
            options = merge(rcParams[:chunk_defaults], chunk.options)
            merge!(chunk.options, options)
        end

        restore = (cache ==:user && typeof(chunk) == CodeChunk && chunk.options[:cache])

        if cached != nothing && (cache == :all || restore)
            result_chunks = restore_chunk(chunk, cached)
        else
            result_chunks = run_chunk(chunk, report, mod)
        end

        executed = [executed; result_chunks]
    end

    doc.header_script = report.header_script

    popdisplay(report)

    #Clear variables from used sandbox
    mod == :sandbox && clear_sandbox(SandBox)
    doc.chunks = executed

    if cache != :off
        write_cache(doc, cache_path)
    end

    return doc
end

"""Detect the output format based on file extension"""
function detect_doctype(source::AbstractString)
  ext = lowercase(splitext(source)[2])
  ext == ".jl" && return "md2html"
  occursin("md", ext) && return "md2html"
  occursin("rst", ext) && return "rst"
  occursin("tex", ext) && return "texminted"
  occursin("txt", ext) && return "asciidoc"

  return "pandoc"
end

function run_chunk(chunk::CodeChunk, report::Report, SandBox::Module)
    @info("Weaving chunk $(chunk.number) from line $(chunk.start_line)")
    result_chunks = eval_chunk(chunk, report, SandBox)
    occursin("2html", report.formatdict[:doctype]) && (result_chunks = embed_figures(result_chunks, report.cwd))
    return result_chunks
end

function embed_figures(chunk::CodeChunk, cwd)
    chunk.figures = [img2base64(fig, cwd) for fig in chunk.figures]
    return chunk
end

function embed_figures(result_chunks, cwd)
    for i in 1:length(result_chunks)
        figs = result_chunks[i].figures
        if !isempty(figs)
            result_chunks[i].figures = [img2base64(fig, cwd) for fig in figs]
        end
    end
    return result_chunks
end

function img2base64(fig, cwd)
  ext = splitext(fig)[2]
  f = open(joinpath(cwd, fig), "r")
    raw = read(f)
  close(f)
  if ext == ".png"
    return "data:image/png;base64," * stringmime(MIME("image/png"), raw)
  elseif ext == ".svg"
    return "data:image/svg+xml;base64," * stringmime(MIME("image/svg"), raw)
  elseif ext == ".gif"
    return "data:image/gif;base64," * stringmime(MIME("image/gif"), raw)
  else
    return(fig)
  end
end

function run_chunk(chunk::DocChunk, report::Report, SandBox::Module)
    chunk.content =  [run_inline(c, report, SandBox) for c in chunk.content]
    return chunk
end

function run_inline(inline::InlineText, report::Report, SandBox::Module)
    return inline
end

function run_inline(inline::InlineCode, report::Report, SandBox::Module)
    #Make a temporary CodeChunk for running code. Collect results and don't wrap
    chunk = CodeChunk(inline.content, 0, 0, "", Dict(:hold => true, :wrap => false))
    options = merge(rcParams[:chunk_defaults], chunk.options)
    merge!(chunk.options, options)

    chunks = eval_chunk(chunk, report, SandBox)
    occursin("2html", report.formatdict[:doctype]) && (chunks = embed_figures(chunks, report.cwd))

    output = chunks[1].output
    endswith(output, "\n") && (output = output[1:end-1])
    inline.output = output
    inline.rich_output = chunks[1].rich_output
    inline.figures = chunks[1].figures
    return inline
end

function reset_report(report::Report)
    report.cur_result = ""
    report.figures = AbstractString[]
    report.term_state = :text
end

function run_code(chunk::CodeChunk, report::Report, SandBox::Module)
    expressions = parse_input(chunk.content)
    N = length(expressions)
    #@show expressions
    result_no = 1
    results = ChunkOutput[ ]

    for (str_expr, expr) = expressions
        reset_report(report)
        lastline = (result_no == N)
        (obj, out) = capture_output(expr, SandBox, chunk.options[:term],
                      chunk.options[:display], lastline, report.throw_errors)
        figures = report.figures #Captured figures
        result = ChunkOutput(str_expr, out, report.cur_result, report.rich_output, figures)
        report.rich_output = ""
        push!(results, result)
        result_no += 1
    end
    return results
end

getstdout() = stdout

function capture_output(expr, SandBox::Module, term, disp,
                        lastline, throw_errors=false)
    #oldSTDOUT = STDOUT
    oldSTDOUT = getstdout()
    out = nothing
    obj = nothing
    rw, wr = redirect_stdout()
    reader = @async read(rw, String)
    try
        obj = Core.eval(SandBox, expr)
        if (term || disp) && typeof(expr) == Expr && expr.head != :toplevel
            obj != nothing && display(obj)
        elseif typeof(expr) == Symbol
            display(obj)
        #This shows images and lone variables, result can
        #Handle last line sepately
        elseif lastline && obj != nothing
            (expr.head != :toplevel && expr.head != :(=)) && display(obj)
        end
    catch E
        throw_errors && throw(E)
        display(E)
        @warn("ERROR: $(typeof(E)) occurred, including output in Weaved document")
    finally
        redirect_stdout(oldSTDOUT)
        close(wr)
        out = fetch(reader)
        close(rw)
    end
    return (obj, out)
end

#Parse chunk input to array of expressions
function parse_input(input::AbstractString)
    parsed = Tuple{AbstractString, Any}[]
    input = lstrip(input)
    n = length(input)
    pos = 1 #The first character is extra line end
    while pos ≤ n
        oldpos = pos
        code,  pos = Meta.parse(input, pos)
        push!(parsed, (input[oldpos:pos-1] , code ))
    end
    parsed
end

function eval_chunk(chunk::CodeChunk, report::Report, SandBox::Module)

    if !chunk.options[:eval]
        chunk.output = ""
        chunk.options[:fig] = false
        return chunk
    end

    #Run preexecute_hooks
    for hook in preexecute_hooks
      chunk = Compat.invokelatest(hook, chunk)
    end

    report.fignum = 1
    report.cur_chunk = chunk

    if haskey(report.formatdict, :out_width) && chunk.options[:out_width] == nothing
        chunk.options[:out_width] = report.formatdict[:out_width]
    end

    chunk.result = run_code(chunk, report, SandBox)

    #Run post_execute chunks
    for hook in postexecute_hooks
      chunk = Compat.invokelatest(hook, chunk)
    end

    if chunk.options[:term]
        chunks = collect_results(chunk, TermResult())
    elseif chunk.options[:hold]
        chunks = collect_results(chunk, CollectResult())
    else
        chunks = collect_results(chunk, ScriptResult())
    end

      #else
     #   chunk.options[:fig] && (chunk.figures = copy(report.figures))
    #end

    chunks
end

#function eval_chunk(chunk::DocChunk, report::Report, SandBox)
#    chunk
#end

#Set all variables to nothing
function clear_sandbox(SandBox::Module)
    for name = names(SandBox, all=true)
        if name != :eval && name != names(SandBox)[1]
            try eval(SandBox, Meta.parse(AbstractString(AbstractString(name), "=nothing"))) catch; end
        end
    end
end

function get_figname(report::Report, chunk; fignum = nothing, ext = nothing)
    figpath = joinpath(report.cwd, chunk.options[:fig_path])
    isdir(figpath) || mkpath(figpath)
    ext == nothing && (ext = chunk.options[:fig_ext])
    fignum == nothing && (fignum = report.fignum)

    chunkid = (chunk.options[:label] == nothing) ? chunk.number : chunk.options[:label]
    full_name = joinpath(report.cwd, chunk.options[:fig_path],
        "$(report.basename)_$(chunkid)_$(fignum)$ext")
    rel_name = "$(chunk.options[:fig_path])/$(report.basename)_$(chunkid)_$(fignum)$ext" #Relative path is used in output
    return full_name, rel_name
end

function get_cwd(doc::WeaveDoc, out_path)
    #Set the output directory
    if out_path == :doc
        cwd = doc.path
    elseif out_path == :pwd
        cwd = pwd()
    else
        #If there is no extension, use as path
        splitted = splitext(out_path)
        if splitted[2] == ""
            cwd = expanduser(out_path)
        else
            cwd = splitdir(expanduser(out_path))[1]
        end
    end
    return cwd
end

"""Get output file name based on out_path"""
function get_outname(out_path::Symbol, doc::WeaveDoc; ext = nothing)
    ext == nothing && (ext = doc.format.formatdict[:extension])
    outname = "$(doc.cwd)/$(doc.basename).$ext"
end

"""Get output file name based on out_path"""
function get_outname(out_path::AbstractString, doc::WeaveDoc; ext = nothing)
    ext == nothing && (ext = doc.format.formatdict[:extension])
    splitted = splitext(out_path)
    if (splitted[2]) == ""
        outname = "$(doc.cwd)/$(doc.basename).$ext"
    else
        outname = expanduser(out_path)
    end
end

function set_rc_params(formatdict, fig_path, fig_ext)
    if fig_ext == nothing
        rcParams[:chunk_defaults][:fig_ext] = formatdict[:fig_ext]
        docParams[:fig_ext] = formatdict[:fig_ext]
    else
        rcParams[:chunk_defaults][:fig_ext] = fig_ext
        docParams[:fig_ext] = fig_ext
    end
        rcParams[:chunk_defaults][:fig_path] = fig_path
        docParams[:fig_path] = fig_path
    return nothing
end

function collect_results(chunk::CodeChunk, fmt::ScriptResult)
    content = ""
    result_no = 1
    result_chunks = CodeChunk[ ]
    for r = chunk.result
        #Check if there is any output from chunk
        if strip(r.stdout) == "" && isempty(r.figures)  && strip(r.rich_output) == ""
            content *= r.code
        else
            content = "\n" * content * r.code
            rchunk = CodeChunk(content, chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
            content = ""
            rchunk.result_no = result_no
            result_no *=1
            rchunk.figures = r.figures
            rchunk.output = r.stdout * r.displayed
            rchunk.rich_output = r.rich_output
            push!(result_chunks, rchunk)
        end
    end
    if content != ""
        startswith(content, "\n") || (content = "\n" * content)
        rchunk = CodeChunk(content, chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_results(chunk::CodeChunk, fmt::TermResult)
    output = ""
    prompt = chunk.options[:prompt]
    result_no = 1
    result_chunks = CodeChunk[ ]
    for r = chunk.result
        output *= prompt * r.code
        output *= r.displayed * r.stdout
        if !isempty(r.figures)
            rchunk = CodeChunk("", chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
            rchunk.output = output
            output = ""
            rchunk.figures = r.figures
            push!(result_chunks, rchunk)
        end
    end
    if output != ""
         rchunk = CodeChunk("", chunk.number, chunk.start_line, chunk.optionstring, copy(chunk.options))
         rchunk.output = output
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_results(chunk::CodeChunk, fmt::CollectResult)
    result_no = 1
    for r =chunk.result
        chunk.output *=  r.stdout
        chunk.rich_output *= r.rich_output
        chunk.figures = [chunk.figures; r.figures]
    end
    return [chunk]
end


text:
"""
    unsafe_createcolortable(palette::GDALPaletteInterp)

Construct a new color table.
"""
unsafe_createcolortable(palette::GDALPaletteInterp) =
    ColorTable(GDAL.gdalcreatecolortable(palette))

"Destroys a color table."
function destroy(ct::ColorTable)
    GDAL.gdaldestroycolortable(ct.ptr)
    ct.ptr = C_NULL
end

"""
    unsafe_clone(ct::ColorTable)

Make a copy of a color table.
"""
unsafe_clone(ct::ColorTable) = ColorTable(GDAL.gdalclonecolortable(ct.ptr))

"""
    paletteinterp(ct::ColorTable)

Fetch palette interpretation.

### Returns
palette interpretation enumeration value, usually `GPI_RGB`.
"""
paletteinterp(ct::ColorTable) = GDAL.gdalgetpaletteinterpretation(ct.ptr)

"""
    ncolorentry(ct::ColorTable)

Get number of color entries in table.
"""
ncolorentry(ct::ColorTable) = GDAL.gdalgetcolorentrycount(ct.ptr)

"""
    getcolorentry(ct::ColorTable, i::Integer)

Fetch a color entry from table.
"""
getcolorentry(ct::ColorTable, i::Integer) =
    unsafe_load(GDAL.gdalgetcolorentry(ct.ptr, i))

"""
    getcolorentryasrgb(ct::ColorTable, i::Integer)

Fetch a table entry in RGB format.

In theory this method should support translation of color palettes in non-RGB
color spaces into RGB on the fly, but currently it only works on RGB color
tables.

### Parameters
* `i`   entry offset from zero to GetColorEntryCount()-1.

### Returns
`true` on success, or `false` if the conversion isn't supported.
"""
function getcolorentryasrgb(ct::ColorTable, i::Integer)
    colorentry = Ref{GDAL.GDALColorEntry}(GDAL.GDALColorEntry(0, 0, 0, 0))
    result = Bool(GDAL.gdalgetcolorentryasrgb(ct.ptr, i, colorentry))
    result || @warn("The conversion to RGB isn't supported.")
    return colorentry[]
end

"""
    setcolorentry!(ct::ColorTable, i::Integer, entry::GDAL.GDALColorEntry)

Set entry in color table.

Note that the passed in color entry is copied, and no internal reference to it
is maintained. Also, the passed in entry must match the color interpretation of
the table to which it is being assigned.

The table is grown as needed to hold the supplied offset.

### Parameters
* `i`     entry offset from `0` to `ncolorentry()-1`.
* `entry` value to assign to table.
"""
function setcolorentry!(ct::ColorTable, i::Integer, entry::GDAL.GDALColorEntry)
    GDAL.gdalsetcolorentry(ct.ptr, i, Ref{GDAL.GDALColorEntry}(entry))
    return ct
end

"""
    createcolorramp!(ct::ColorTable, startindex, startcolor::GDAL.GDALColorEntry,
        endindex, endcolor::GDAL.GDALColorEntry)

Create color ramp.

Automatically creates a color ramp from one color entry to another. It can be
called several times to create multiples ramps in the same color table.

### Parameters
* `startindex` index to start the ramp on the color table [0..255]
* `startcolor` a color entry value to start the ramp
* `endindex`   index to end the ramp on the color table [0..255]
* `endcolor`   a color entry value to end the ramp

### Returns
total number of entries, -1 to report error
"""
function createcolorramp!(
        ct::ColorTable,
        startindex::Integer,
        startcolor::GDAL.GDALColorEntry,
        endindex::Integer,
        endcolor::GDAL.GDALColorEntry
    )
    return GDAL.gdalcreatecolorramp(
        ct.ptr,
        startindex,
        Ref{GDAL.GDALColorEntry}(startcolor),
        endindex,
        Ref{GDAL.GDALColorEntry}(endcolor)
    )
end


text:
__precompile__(true)

module PerronFrobenius

using Reexport
@reexport using CausalityToolsBase
using StateSpaceReconstruction

import StateSpaceReconstruction:
    GroupSlices,
    Embeddings.AbstractEmbedding,
    DelaunayTriangulation,
    cembed,
    invariantize

import Simplices:
    subsample_coeffs,
    simplexintersection

using Distributed
using SharedArrays
using StaticArrays
using SparseArrays
using InplaceOps
using RecipesBase
using Printf
using LinearAlgebra
using DelayEmbeddings

# Abstract estimator type
abstract type TransferOperatorEstimator end

include("BinningSchemes/BinningSchemes.jl")
include("TransferOperators/TransferOperators.jl")
include("InvariantMeasures/InvariantMeasures.jl")

include("is_markov.jl")

end


text:
## --- Load required pacages, install Chron if required

    try
        using Chron
    catch
        using Pkg
        Pkg.add("Chron")
        using Chron
    end

    using Statistics, StatsBase, SpecialFunctions
    using Plots; gr();

## --- Define sample properties

    # # # # # # # # # # # Enter sample information here! # # # # # # # # # # # #
    # Input the number of samples we wish to model (must match below)
    nSamples = 4
    # Make an instance of a ChronSection object for nSamples
    smpl = NewChronAgeData(nSamples)
    smpl.Name           = ("Sample 1", "Sample 2", "Sample 3", "Sample 4") # Et cetera
    smpl.Age_14C       .= [ 6991,  7088,  7230,  7540,] # Measured ages
    smpl.Age_14C_sigma .= [   30,    70,    50,    50,] # Measured 1-σ uncertainties
    smpl.Height        .= [ -355,  -380,-397.0,-411.5,] # Depths below surface should be negative
    smpl.Height_sigma  .= fill(0.01, nSamples) # Usually assume little or no sample height uncertainty
    smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
    smpl.Age_Unit = "Years BP" # Unit of measurement for ages
    smpl.Height_Unit = "m" # Unit of measurement for Height and Height_sigma

    # IMPORTANT: smpl.Height must increase with increasing stratigraphic height
    # -- i.e., stratigraphically younger samples must be more positive. For this
    # reason, it is convenient to represent depths below surface as negative
    # numbers.
    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

## --- Calculate and plot calendar age PDFs for each sample

    smpl.Params = fill(NaN, length(intcal13["Age_Calendar"]), nSamples)
    for i = 1:nSamples
        # The likelihood that a measured 14C age could result from a sample of
        # a given calendar age is proportional to the intergral of the product
        # of the two respective distributions
        likelihood = normproduct.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], intcal13["Age_14C"], intcal13["Age_sigma"])
        likelihood ./= sum(likelihood) * intcal13["dt"] # Normalize

        # Estimate mean and standard deviation by drawing samples from distribution
        samples = draw_from_distribution(likelihood, 10^6) .* maximum(intcal13["Age_Calendar"])
        smpl.Age[i] = mean(samples)
        smpl.Age_sigma[i] = std(samples)
        smpl.Age_025CI[i] = percentile(samples,2.5)
        smpl.Age_975CI[i] = percentile(samples,97.5)

        # Populate smpl.Params with log likelihood for each sample
        smpl.Params[:,i] = normproduct_ll.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], intcal13["Age_14C"], intcal13["Age_sigma"])

        # Plot likelihood vector for each sample
        t = (intcal13["Age_Calendar"] .> smpl.Age[i] - 5*smpl.Age_sigma[i]) .& (intcal13["Age_Calendar"] .< smpl.Age[i] + 5*smpl.Age_sigma[i])
        plot(intcal13["Age_Calendar"][t], likelihood[t], label=smpl.Name[i], xlabel="Calendar Age", ylabel="Likelihood")
        savefig("$(smpl.Name[i])_CalendarAge.pdf")
    end

## --- Run stratigraphic model

    # # # # # # # # # # Configure stratigraphic model here! # # # # # # # # # #
    # Configure the stratigraphic Monte Carlo model
    config = NewStratAgeModelConfiguration()
    # If you in doubt, you can probably leave these parameters as-is
    config.resolution = 0.2 # Same units as sample height. Smaller is slower!
    config.bounding = 0.5 # how far away do we place runaway bounds, as a fraction of total section height
    (bottom, top) = extrema(smpl.Height)
    npoints_approx = round(Int,length(bottom:config.resolution:top) * (1 + 2*config.bounding))
    config.nsteps = 15000 # Number of steps to run in distribution MCMC
    config.burnin = 10000*npoints_approx # Number to discard
    config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

    # Run the stratigraphic MCMC model
    @time (mdl, agedist, lldist) = StratMetropolis14C(smpl, config)

    # Plot results (mean and 95% confidence interval for both model and data)
    hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
    plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="")
    plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black)
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))")
    savefig(hdl,"AgeDepthModel.pdf");
    display(hdl)

## --- Interpolate results at a specific height

    # Stratigraphic height at which to interpolate
    height = -404

    age_at_height = linterp1s(mdl.Height,mdl.Age,height)
    age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,height)
    age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,height)
    print("Interpolated age at height=$height: $age_at_height +$(age_at_height_max-age_at_height)/-$(age_at_height-age_at_height_min) $(smpl.Age_Unit)")

    # Optional: interpolate full age distribution
    interpolated_distribution = Array{Float64}(undef,size(agedist,2))
    for i=1:size(agedist,2)
        interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],height)
    end
    hdl = histogram(interpolated_distribution, nbins=50, label="")
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit)) at height=$height", ylabel="Likelihood (unnormalized)")
    savefig(hdl, "Interpolated age distribution.pdf")
    display(hdl)

## --- Calculate deposition rate binned by age  - - - - - - - - - - - - - - - -

    # Set bin width and spacing
    binwidth = round(nanrange(mdl.Age)/10,sigdigits=1) # Can also set manually, commented out below
    # binwidth = 100 # Same units as smpl.Age
    binoverlap = 10
    ages = collect(minimum(mdl.Age):binwidth/binoverlap:maximum(mdl.Age))
    bincenters = ages[1+Int(binoverlap/2):end-Int(binoverlap/2)]
    spacing = binoverlap

    # Calculate rates for the stratigraphy of each markov chain step
    dhdt_dist = Array{Float64}(undef, length(ages)-binoverlap, config.nsteps)
    @time for i=1:config.nsteps
        heights = linterp1(reverse(agedist[:,i]), reverse(mdl.Height), ages)
        dhdt_dist[:,i] .= abs.(heights[1:end-spacing] - heights[spacing+1:end]) ./ binwidth
    end

    # Find mean and 1-sigma (68%) CI
    dhdt = nanmean(dhdt_dist,dim=2)
    dhdt_50p = nanmedian(dhdt_dist,dim=2)
    dhdt_16p = nanpctile(dhdt_dist,15.865,dim=2) # Minus 1-sigma (15.865th percentile)
    dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
    # Other confidence intervals (10:10:50)
    dhdt_20p = nanpctile(dhdt_dist,20,dim=2)
    dhdt_80p = nanpctile(dhdt_dist,80,dim=2)
    dhdt_25p = nanpctile(dhdt_dist,25,dim=2)
    dhdt_75p = nanpctile(dhdt_dist,75,dim=2)
    dhdt_30p = nanpctile(dhdt_dist,30,dim=2)
    dhdt_70p = nanpctile(dhdt_dist,70,dim=2)
    dhdt_35p = nanpctile(dhdt_dist,35,dim=2)
    dhdt_65p = nanpctile(dhdt_dist,65,dim=2)
    dhdt_40p = nanpctile(dhdt_dist,40,dim=2)
    dhdt_60p = nanpctile(dhdt_dist,60,dim=2)
    dhdt_45p = nanpctile(dhdt_dist,45,dim=2)
    dhdt_55p = nanpctile(dhdt_dist,55,dim=2)

    # Plot results
    hdl = plot(bincenters,dhdt, label="Mean", color=:black, linewidth=2)
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_20p; reverse(dhdt_80p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_25p; reverse(dhdt_75p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_30p; reverse(dhdt_70p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_35p; reverse(dhdt_65p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_40p; reverse(dhdt_60p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_45p; reverse(dhdt_55p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,bincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
    # savefig(hdl,"DepositionRateModelCI.pdf")
    display(hdl)

## --- Optional: Stratigraphic model including hiatuses

    # Data about hiatuses
    nHiatuses = 2 # The number of hiatuses you have data for
    hiatus = NewHiatusData(nHiatuses) # Struct to hold data
    hiatus.Height         = [-371.5, -405.0 ]
    hiatus.Height_sigma   = [   0.0,    0.0 ]
    hiatus.Duration       = [ 100.0,   123.0]
    hiatus.Duration_sigma = [  30.5,    20.0]

    # Run the model. Note the additional `hiatus` arguments
    @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis14C(smpl, hiatus, config)

    # Plot results (mean and 95% confidence interval for both model and data)
    hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
    plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white)
    plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black)
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))")

## --- End of File


text:
#
# Query output data.
#

immutable OutputFlow <: AbstractVector{Any}
    oty::Output
    col::Column
end

OutputFlow(oty, offs, vals) = OutputFlow(oty, Column(offs, vals))

# Array interface.

size(flow::OutputFlow) = (flow.col.len,)

length(flow::OutputFlow) = flow.col.len

function getindex(flow::OutputFlow, i)
    data = getdata(Any, flow.col, i)
    if isentity(flow.oty)
        N = classname(domain(flow.oty))
        E = Entity{N}
        data =
            isplural(flow.oty) ?
                EntityVector{N,typeof(data)}(data) :
            isoptional(flow.oty) ?
                !isnull(data) ? Nullable{E}(get(data)) : Nullable{E}() :
                E(data)
    end
    return data
end

getindex(flow::OutputFlow, idxs::AbstractVector{Int}) =
    OutputFlow(flow.oty, flow.col[idxs])

Base.linearindexing(flow::OutputFlow) = Base.LinearFast()

vcat(flow1::OutputFlow, flow2::OutputFlow) =
    OutputFlow(obound(flow1.oty, flow2.oty), vcat(flow1.col, flow2.col))

Base.array_eltype_show_how(::OutputFlow) = (true, "")

summary(flow::OutputFlow) = "OutputFlow[$(length(flow.col)) \ud7 $(flow.oty)]"

# Data components.

convert(::Type{Column}, flow::OutputFlow) =
    flow.col

column(flow::OutputFlow) = flow.col
offsets(flow::OutputFlow) = offsets(flow.col)
values(flow::OutputFlow) = values(flow.col)

# Output signature and its properties.

output(flow::OutputFlow) = flow.oty
domain(flow::OutputFlow) = domain(flow.oty)
mode(flow::OutputFlow) = mode(flow.oty)
decorations(flow::OutputFlow) = decorations(flow.oty)

isplain(flow::OutputFlow) = isregular(flow.oty)
isoptional(flow::OutputFlow) = isoptional(flow.oty)
isplural(flow::OutputFlow) = isplural(flow.oty)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto