text:
### Extra Functions for Route Planning ###
### Ted Steiner, MIT AeroAstro ###
### August 2014 - February 2015  ###

# Convenience function to get a list of routable waypoints in an environment.
function getWaypoints( network::OpenStreetMap.Network,  # Routing network
                       candidates::Vector{Int},         # List of candidate waypoints (vertices in network)
                       min_dist::Real,                  # Minimum distance between waypoints
                       num_waypoints::Int,              # Desired waypoint count
                       avoid::Vector{Int}=Int[],        # Optional list of candidates to avoid selecting
                       min_distance::Real=0 )         # Optional minimum route distance

    waypoints = Int[]
    valid, waypoints, route_dist = randomRouteRecursive(network,candidates,waypoints,min_distance=min_distance,min_waypoints=num_waypoints,min_waypoint_dist=min_dist,avoid=avoid)
    return waypoints
end

### Generate a list of routable waypoints
# NOTE: If network.w is times (for fastest routing), distance below refers
# instead to time.
function randomRouteRecursive( network::OpenStreetMap.Network,  # Routing network
                               candidates::Vector{Int},         # List of candidate waypoints
                               route::Vector{Int},              # List of waypoints
                               route_distance::Float64=0.0;     # Accumulated route distance
                               avoid::Array{Int,1}=Int[],       # Optional list of candidates to skip
                               min_waypoints::Int=2,            # Minimum number of waypoints in route
                               min_waypoint_dist::Real=0,       # Minimum distance between waypoints
                               min_distance::Real=0 )           # Minimum route distance

    if length(route) >= min_waypoints && route_distance >= min_distance
        return true, route, route_distance
    end

    order = randperm(length(candidates))

    for index in order
        if !(candidates[index] in avoid)
            if length(route) < 1
                push!(route,candidates[index])
                valid_future, route_new, route_distance = randomRouteRecursive(network,candidates,route,route_distance,avoid=avoid,min_waypoints=min_waypoints,min_waypoint_dist=min_waypoint_dist,min_distance=min_distance)
                if !valid_future
                    pop!(route)
                end
            else
                valid_this, next_dist = existsValidRoute(network,route[end],candidates[index],min_waypoint_dist)
                if valid_this
                    if route[end] != candidates[index]
                        # To prevent getting stuck or cycling, make sure next point is not one of the previous 5 points
                        if length(route) < 3 || (route[end-1] != candidates[index])
                            if length(route) < 4 || (route[end-2] != candidates[index])
                                if length(route) < 5 || (route[end-3] != candidates[index])
                                    if length(route) < 6 || (route[end-4] != candidates[index])
                                        push!(route,candidates[index])
                                        valid_future, route_new, route_distance_new = randomRouteRecursive(network,candidates,route,route_distance+next_dist,avoid=avoid,min_waypoints=min_waypoints,min_waypoint_dist=min_waypoint_dist,min_distance=min_distance)
                                        if valid_future
                                            return true, route_new, route_distance_new
                                        else
                                            pop!(route)
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    return false, route, route_distance
end

# Compute a route along a list of waypoints
function computeRoute( network::OpenStreetMap.Network,  # Routing network
                       waypoints::Vector{Int},          # List of waypoint node IDs
                       fastest::Bool )                  # True: fastest time, False: shortest distance

    route = Int[]
    metric = 0.0 # Time if fastest=true, else distance

    for k = 2:length(waypoints)
        s = waypoints[k-1]
        t = waypoints[k]
        if fastest
            route_nodes, route_metric = OpenStreetMap.fastestRoute(network,s,t)
        else
            route_nodes, route_metric = OpenStreetMap.shortestRoute(network,s,t)
        end

        if k == 2
            push!(route, route_nodes...)
        else
            push!(route, route_nodes[2:end]...)
        end

        metric += route_metric
    end

    return route, metric
end

# Get distances from route start to all subsequent route distances
function distanceFromStart( nodes::Dict{Int,Geodesy.ENU},
                            route::Vector{Int} )
    dist = 0.0
    dists = zeros(length(route))
    prev_point = nodes[route[1]]
    for i = 2:length(route)
        point = nodes[route[i]]
        dist += OpenStreetMap.distance(prev_point,point)
        dists[i] = deepcopy(dist)
        prev_point = point
    end
    return dists
end

# Add nodes along a route according to some fixed-size spacing
function interpolateRoute!( nodes::Dict{Int,Geodesy.ENU},   # Node locations
                            route::Vector{Int},             # List of node IDs
                            d_spacing::Real )               # Spacing between interpolated nodes (meters)
    interp_route = Int[]
    push!(interp_route,route[1])
    for k = 2:length(route)
        push!(interp_route,interpolateNodes!(nodes,route[k-1],route[k],d_spacing)...)
        push!(interp_route,route[k])
    end

    return interp_route
end

# Helper function for interpolateRoute! to place the nodes
function interpolateNodes!( nodes, s, t, d_spacing )
    if s == t
        return Int[]
    end

    x1 = nodes[s].east
    x2 = nodes[t].east
    y1 = nodes[s].north
    y2 = nodes[t].north
    d_tot = sqrt((x2-x1)^2 + (y2-y1)^2)
    d_frac = d_spacing / d_tot

    num_interp = int(floor(d_tot/d_spacing))-1
    if num_interp < 1
        return Int[]
    end

    interp = zeros(Int,num_interp)
    for k = 1:length(interp)
        x = x1 + k*d_frac*(x2-x1)
        y = y1 + k*d_frac*(y2-y1)
        loc = OpenStreetMap.ENU(x,y)
        interp[k] = OpenStreetMap.addNewNode!(nodes,loc)
    end

    return interp
end

# Helper function for randomRouteRecursive()
function existsValidRoute(network, node0, node1, min_dist=0.0)
    start_vertex = network.v[node0]
    finish_vertex = network.v[node1]
    finish_index = network.v[node1].index

    dijkstra_result = Graphs.dijkstra_shortest_paths(network.g,network.w,start_vertex,visitor=EndWhenAnyNode([finish_vertex]))

    dist = dijkstra_result.dists[finish_index]
    if dist < Inf && dist > min_dist
        return true, dist
    end

    return false, -1
end


text:
"""
	Constraint(A, b, convex_set_type, dim = 0, indices = 0:0)

Creates a COSMO constraint: `Ax + b ∈ convex_set`.

By default the following convex set types are supported: `ZeroSet`, `Nonnegatives`, `SecondOrderCone`, `PsdCone`, `PsdConeTriangle`.

# Examples
```jldoctest; setup = :(using COSMO)
julia> COSMO.Constraint([1 0;0 1], zeros(2), COSMO.Nonnegatives)
Constraint
Size of A: (2, 2)
ConvexSet: COSMO.Nonnegatives{Float64}
```

For convex sets that require their own data, it is possible to pass the pass the instantiated object directly rather than the type name.
# Examples
```jldoctest; setup = :(using COSMO)
julia> COSMO.Constraint([1 0;0 1], zeros(2), COSMO.Box([-1.;-1.],[1.;1.]))
Constraint
Size of A: (2, 2)
ConvexSet: COSMO.Box{Float64}
```

---
The optional arguments `dim` and `indices` can be used to specify A and b for subparts of variable `x`. If `x` has dimension `dim = 4`,
then x[2] and x[3] can be constrained to the zero cone in the following way:

# Examples
```jldoctest; setup = :(using COSMO)
julia> c = COSMO.Constraint([1 0;0 1], zeros(2), COSMO.ZeroSet, 4, 2:3)
Constraint
Size of A: (2, 4)
ConvexSet: COSMO.ZeroSet{Float64}
```
Notice that extra columns of A have been added automatically.
```
julia>Matrix(c.A)
2×4 Array{Float64,2}:
0.0  1.0  0.0  0.0
0.0  0.0  1.0  0.0
```
"""
struct Constraint{T <: AbstractFloat}
	A::Union{AbstractMatrix{T}, AbstractVector{T}}
	b::AbstractVector{T}
	convex_set::AbstractConvexSet{T}

	# constructor
	function Constraint{T}(
		A::AbstractMatrix{T},
		b::AbstractVector{T},
		convex_set::AbstractConvexSet{T},
		dim::Integer = 0,
		indices::UnitRange = 0:0) where{T}

		size(A, 1) != length(b) && throw(DimensionMismatch("The dimensions of matrix A and vector b don't match."))
		size(A,1)  != convex_set.dim && throw(DimensionMismatch("The row dimension of A doesn't match the dimension of the constraint set."))
		size(b, 2) != 1 && throw(DimensionMismatch("Input b must be a vector or a scalar."))

		if indices != 0:0
			(indices.start < 1 || indices.stop < indices.start) && throw(DomainError("The index range for x has to be increasing and nonnegative."))
			dim < indices.stop && throw(DomainError("The dimension of x: $(dim) must be equal or higher than the the stop value of indices: $(indices.stop)."))
			Ac = spzeros(size(A, 1), dim)
			Ac[:, indices] = A
			A = Ac
		end
		new(A, b, convex_set)
	end
end

#We always choose to internally convert whatever datatype
#we are given to a floating point type either specified by
#the user or to DefaultFloat
Constraint(args...) = Constraint{DefaultFloat}(args...)

#support the case where only a Type is specified for the set, and we need to
#create an instance of the appropriate size
function Constraint{T}(
	A::AbstractMatrix{T},
	b::AbstractVector{T},
	set_type::Type{ <: AbstractConvexSet},args...) where{T}

	# this constructor doesnt work with cones that need special arguments like the power cone
	set_type <: ArgumentCones && error("You can't create a constraint by passing the convex set as a type, if your convex set is a $(set_type). Please pass an object.")
	# call the appropriate set constructor
	convex_set = set_type{T}(size(A, 1))
	Constraint{T}(A, b, convex_set, args...)
end

#allow various ways of passing A and b and convert  to AbstractMatrix and AbstractVector types
function Constraint{T}(A::AbstractMatrix, b::AbstractVector, args...) where{T}
	Constraint{T}(AbstractMatrix{T}(A), AbstractVector{T}(b), args...)
end

#all others convert first o matrix / vector args
function Constraint{T}(A::Real,b::Real,args...) where{T}
	Constraint{T}(reshape([A], 1, 1), [b], args...)
end

function Constraint{T}(A::AbstractMatrix,b::AbstractMatrix, args...) where {T}
	Constraint{T}(A, vec(b), args...)
end

function Constraint{T}(A::Union{AbstractVector,AbstractMatrix}, b::Real, args...) where{T}
	Constraint{T}(reshape(A, 1, length(A)), [b], args...)
end

function Constraint{T}(A::AbstractVector,b::AbstractVector,args...) where{T}
	Constraint{T}(reshape(A, length(A), 1), b, args...)
end

function Base.show(io::IO, obj::COSMO.Constraint)
	print(io,"Constraint\nSize of A: $(size(obj.A))\nConvexSet: $(typeof(obj.convex_set))")
end

ext: jl


text:
module InteractiveViz

#using AbstractPlotting
using GLMakie
using Printf
using Colors
using Statistics

include("types.jl")
include("core.jl")
include("samplers.jl")
include("simple.jl")
include("demo.jl")

end # module


text:
@testset "92.reverse-linked-list-ii.jl" begin
    @test reverse_between(ListNode{Int}([1, 2, 3, 4, 5]), 2, 4) ==
          ListNode{Int}([1, 4, 3, 2, 5])
    @test reverse_between(ListNode{Int}([1, 2, 3, 4, 5]), 1, 3) ==
          ListNode{Int}([3, 2, 1, 4, 5])
end


text:
macro inbounds(ex)
    esc(ex)
end

include(joinpath(@__DIR__, "../../../../examples/reduce_partition_by.jl"))

using Test
@testset begin
    @test sum(c_xs) == length(xs)
    @test sum(m_xs .* c_xs) ≈ sum(xs)
end


text:
import Base: IteratorSize, HasLength, IsInfinite, length

"""
`AbstractTrigger` is an abstract type for Triggers

A Trigger define when a job should be run.

`AbstractTrigger` should implement `get_next_dt_fire` function which returns instant at which a job should be run (given current instant `dt_now` and instant when job was previously run (`dt_previous_fire`)
"""
abstract type AbstractTrigger end
abstract type AbstractFiniteTrigger <: AbstractTrigger end
abstract type AbstractInfiniteTrigger <: AbstractTrigger end

"""
NoTrigger define a trigger that never trigger.

It's a useful struct for triggers operations such as applying offset or jitter to a trigger.
"""
struct NoTrigger <: AbstractTrigger
end

function IteratorSize(trigger::AbstractFiniteTrigger)
    HasLength()
end

function IteratorSize(trigger::AbstractInfiniteTrigger)
    IsInfinite()
end

function length(trigger::AbstractFiniteTrigger)
    trigger.n
end

include("instant.jl")
include("time.jl")
include("period.jl")
include("timeframe.jl")
include("offset.jl")
include("jitter.jl")
include("iterator.jl")
include("custom.jl")


text:
"""
    ScaledKernel(k::Kernel, σ²::Real=1.0)

Scaled kernel derived from `k` by multiplication with variance `σ²`.

# Definition

For inputs ``x, x'``, the scaled kernel ``\\widetilde{k}`` derived from kernel ``k`` by
multiplication with variance ``\\sigma^2 > 0`` is defined as
```math
\\widetilde{k}(x, x'; k, \\sigma^2) = \\sigma^2 k(x, x').
```
"""
struct ScaledKernel{Tk<:Kernel,Tσ²<:Real} <: Kernel
    kernel::Tk
    σ²::Vector{Tσ²}
end

function ScaledKernel(kernel::Tk, σ²::Tσ²=1.0) where {Tk<:Kernel,Tσ²<:Real}
    @check_args(ScaledKernel, σ², σ² > zero(Tσ²), "σ² > 0")
    return ScaledKernel{Tk,Tσ²}(kernel, [σ²])
end

@functor ScaledKernel

(k::ScaledKernel)(x, y) = first(k.σ²) * k.kernel(x, y)

function kernelmatrix(κ::ScaledKernel, x::AbstractVector, y::AbstractVector)
    return κ.σ² .* kernelmatrix(κ.kernel, x, y)
end

function kernelmatrix(κ::ScaledKernel, x::AbstractVector)
    return κ.σ² .* kernelmatrix(κ.kernel, x)
end

function kernelmatrix_diag(κ::ScaledKernel, x::AbstractVector)
    return κ.σ² .* kernelmatrix_diag(κ.kernel, x)
end

function kernelmatrix_diag(κ::ScaledKernel, x::AbstractVector, y::AbstractVector)
    return κ.σ² .* kernelmatrix_diag(κ.kernel, x, y)
end

function kernelmatrix!(
    K::AbstractMatrix, κ::ScaledKernel, x::AbstractVector, y::AbstractVector
)
    kernelmatrix!(K, κ.kernel, x, y)
    K .*= κ.σ²
    return K
end

function kernelmatrix!(K::AbstractMatrix, κ::ScaledKernel, x::AbstractVector)
    kernelmatrix!(K, κ.kernel, x)
    K .*= κ.σ²
    return K
end

function kernelmatrix_diag!(K::AbstractVector, κ::ScaledKernel, x::AbstractVector)
    kernelmatrix_diag!(K, κ.kernel, x)
    K .*= κ.σ²
    return K
end

Base.:*(w::Real, k::Kernel) = ScaledKernel(k, w)

Base.show(io::IO, κ::ScaledKernel) = printshifted(io, κ, 0)

function printshifted(io::IO, κ::ScaledKernel, shift::Int)
    printshifted(io, κ.kernel, shift)
    return print(io, "\n" * ("\t"^(shift + 1)) * "- σ² = $(first(κ.σ²))")
end


text:
# Define shorthand functions depending on the returns to scale
FDH_CRS(X,Y,input) = FDH{CRS}(X,Y,input)
FDH_VRS(X,Y,input) = FDH{VRS}(X,Y,input)
FDH_NIRS(X,Y,input) = FDH{NIRS}(X,Y,input)
FDH_NDRS(X,Y,input) = FDH{NDRS}(X,Y,input)

immutable FDH{T<:RS} <: AbstractDEA{FreeDisposal,T}
  Data::DEAData
  input::Bool

  function FDH{T}(X::Array,Y::Array,input::Bool) where T<:RS
    new(DEAData(X,Y),input)
  end
end

function getdata(F::FDH)
  return F.Data
end

# Solve FDH program under VRS with (Xk,Yk) as evaluation point
function (F::FDH{VRS})(Xk::Array,Yk::Array)
  X,Y = F[1:end]
  dompeer = NaN
  dom = find(all(Y .>= Yk,2) & all(X .<= Xk,2))
  if(F.input) #Input-oriented
    theta = Inf
    for k in dom
      Ixk = find(X[k,:] .> 0.0)
      curmin = maximum(X[k,Ixk]./Xk[Ixk])
      if(curmin < theta)
        theta = curmin
        dompeer = k
      end
    end
  else #Output-oriented
    theta = -Inf
    for k in dom
      Iyk = find(Y[k,:] .> 0.0)
      curmax = minimum(Y[k,Iyk]./Yk[Iyk])
      if(curmax > theta)
        theta = curmax
        dompeer = k
      end
    end
  end
  return DEAResult(theta,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:size(X,1)])
end

# Solve FDH program under NDRS with (Xk,Yk) as evaluation point
function (F::FDH{NDRS})(Xk::Array,Yk::Array)
  X,Y = F[1:end]
  K = size(X,1)
  dompeer = NaN
  domx = find(all(X .<= Xk,2))
  if(F.input) #Input-oriented
    theta = Inf
    for k in domx
      Ixk = find(X[k,:] .> 0.0)
      Jyk = find(Y[k,:] .> 0.0)
      curmin = max(maximum(Yk[Jyk]./Y[k,Jyk]),1.0)*maximum(X[k,Ixk]./Xk[Ixk])
      if(curmin < theta)
        theta = curmin
        dompeer = k
      end
    end
  else #Output-oriented
    theta = -Inf
    for k in domx
      Ixk = find(X[k,:] .> 0.0)
      Jyk = find(Y[k,:] .> 0.0)
      curmax = minimum(Y[k,Jyk]./Yk[Jyk])*minimum(Xk[Ixk]./X[k,Ixk])
      if(curmax > theta)
        theta = curmax
        dompeer = k
      end
    end
  end
  return DEAResult(theta,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:K])
end

# Solve FDH program under CRS with (Xk,Yk) as evaluation point
function (F::FDH{CRS})(Xk::Array,Yk::Array)
  X,Y = F[1:end]
  K = size(X,1)
  dompeer = NaN
  dom = find(maximum(reshape(Yk,1,:)./Y,2) .<= minimum(reshape(Xk,1,:)./X,2))
  if(F.input) #Input-oriented
    theta = Inf
    for k in dom
      Ixk = find(X[k,:] .> 0.0)
      Jyk = find(Y[k,:] .> 0.0)
      curmin = maximum(Yk[Jyk]./Y[k,Jyk])*maximum(X[k,Ixk]./Xk[Ixk])
      if(curmin < theta)
        theta = curmin
        dompeer = k
      end
    end
  else #Output-oriented
    theta = -Inf
    for k in dom
      Ixk = find(X[k,:] .> 0.0)
      Jyk = find(Y[k,:] .> 0.0)
      curmax = minimum(Y[k,Jyk]./Yk[Jyk])*minimum(Xk[Ixk]./X[k,Ixk])
      if(curmax > theta)
        theta = curmax
        dompeer = k
      end
    end
  end
  return DEAResult(theta,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:K])
end

# Solve FDH program under NIRS with (Xk,Yk) as evaluation point
function (F::FDH{NIRS})(Xk::Array,Yk::Array)
  X,Y = F[1:end]
  K = size(X,1)
  dompeer = NaN
  domy = find(all(Y .>= Yk,2))
  if(F.input) #Input-oriented
    theta = Inf
    for k in domy
      Ixk = X[k,:] .> 0.0
      Jyk = Y[k,:] .> 0.0
      curmin = maximum(Yk[:,Jyk]./Y[k,Jyk])*maximum(X[k,Ixk]./Xk[:,Ixk])
      if(curmin < theta)
        theta = curmin
        dompeer = k
      end
    end
  else #Output-oriented
    theta = -Inf
    for k in domy
      Ixk = X[k,:] .> 0.0
      Jyk = Y[k,:] .> 0.0
      curmax = minimum(Y[k,Jyk]./Yk[:,Jyk])*min(minimum(Xk[:,Ixk]./X[k,Ixk]),1.0)
      if(curmax > theta)
        theta = curmax
        dompeer = k
      end
    end
  end
  return DEAResult(theta,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:K])
end

function (F::FDH{T})() where T<:RS
  	Data = getdata(F)
  	res = Array{DEAResult}(getnrdmu(Data))

  	@sync @parallel for k in eachindex(Data)
  		res[k] = F(Data[k]...)
  	end

  	return res
end

max_stars_repo_name: kepiej/NEPA


text:
using Test, FinancialMonteCarlo;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

underlying_name = "ENJ"
underlying_name_2 = "ABPL"

underlying_ = underlying_name * "_" * underlying_name_2

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model_enj = BlackScholesProcess(sigma, Underlying(S0, d));
Model_abpl = BlackScholesProcess(sigma, Underlying(S0, d));
Model_n = BlackScholesProcess(sigma, Underlying(S0, d));

Model1 = GaussianCopulaNVariateLogProcess(Model_enj, Model_abpl)
Model = GaussianCopulaNVariateLogProcess(Model_enj, Model_abpl, 0.0)
@test Model1.rho == Model.rho

display(Model)

mktdataset = underlying_ → Model
mktdataset_2 = underlying_name → Model_enj
mktdataset += "notneeded" → Model_n

portfolio = underlying_name → 1.0 * EUData
portfolio += underlying_name → 1.0 * AMData
portfolio += underlying_name → -1.0 * (-1.0) * BarrierData

price_mkt = pricer(mktdataset, rfCurve, mc, portfolio)
price_old = pricer(mktdataset_2, rfCurve, mc, portfolio)

@test abs(price_mkt - price_old) < 1e-8


text:
using HomotopyLasso
using Test

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


text:
# rna.jl
#
# Fit G state models (generalized telegraph models) to RNA abundance data
# For single cell RNA (scRNA) technical noise is included as a yieldfactor
# single molecule FISH (smFISH) is treated as loss less
#

"""
    fit_rna(nchains::Int,gene::String,cell::String,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,fish::Bool,runcycle::Bool,inlabel,label,nsets::Int,cv=0.,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/",yieldprior = 0.05,ejectprior = 1.0)
    fit_rna(nchains::Int,data::AbstractRNAData,gene::String,cell::String,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,fish::Bool,runcycle,inlabel,label,nsets,cv=0.,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/",yieldprior = 0.05,ejectprior = 1.0)

Fit steady state or transient GM model to RNA data for a single gene, write the result (through function finalize), and return nothing.

# Arguments
- `nchains`: number of MCMC chains
- `gene`: gene name
- `cell`: cell type
- `datacond`: condition, if more than one condition use vector of strings e.g. ["DMSO","AUXIN"]
- `maxtime`: float maximum time for entire run
- `infolder`: folder pointing to results used as initial conditions
- `resultfolder`: folder where results go
- `datafolder`: folder for data, string or array of strings
- `inlabel`: name of input files (not including gene name but including condition)
- `label`: = name of output files
- `nsets`: int number of rate sets
- `runcycle`: if true, cycle through all parameters sequentially in MCMC
- `cyclesteps`: int number of steps per cycle
- `samplesteps`: int number of samples
- `warmupsteps`: int number of warmup steps
- `annealsteps`: in number of annealing steps
- `temp`: MCMC temperature
- `tempanneal`: starting temperature for annealing
- `root`: root folder of data and Results folders
- `fittedparam`: vector of rate indices,  indices of parameters to be fit (input as string of ints separated by "-")
- `fixedeffects`: (tuple of vectors of rate indices) string indicating which rate is fixed, e.g. "eject"
- `data`: data structure

"""

function fit_rna(nchains::Int,gene::String,cell::String,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,fish::Bool,runcycle::Bool,inlabel,label,nsets::Int,cv=0.,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = ".",yieldprior=0.05)
    gene = check_genename(gene,"[")
    datafolder = joinpath("data",datafolder)
    if occursin("-",datafolder)
        datafolder = string.(split(datafolder,"-"))
    end
    if occursin("-",datacond)
        datacond = string.(split(datacond,"-"))
    end
    if transient
        data = data_rna(gene,datacond,datafolder,fish,label,root,["T0","T30","T120"],[0.,30.,120.])
    else
        data = data_rna(gene,datacond,datafolder,fish,label,root)
    end
    fit_rna(nchains,data,gene,cell,fittedparam,fixedeffects,datacond,G,maxtime,infolder,resultfolder,datafolder,fish,runcycle,inlabel,label,nsets,cv,transient,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root,yieldprior)
end
function fit_rna(nchains::Int,data::AbstractRNAData,gene::String,cell::String,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,fish::Bool,runcycle,inlabel,label,nsets,cv=0.,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = ".",yieldprior=0.05)
    println(now())
    printinfo(gene,G,datacond,datafolder,infolder,resultfolder,maxtime)
    println("size of histogram: ",data.nRNA)

    resultfolder = joinpath("results",resultfolder)
    infolder = joinpath("results",infolder)

    if fish
        yieldprior = 1.
    end

    model = model_rna(gene,cell,G,cv,fittedparam,fixedeffects,inlabel,infolder,nsets,root,data,yieldprior)
    options = MHOptions(samplesteps,0,warmupsteps,annealsteps,maxtime*.9,temp,tempanneal)
    if runcycle > 0
        model = cycle(nchains,fish,fixedeffects,model,data,options)
    end
    print_ll(data,model)
    fit,stats,measures = run_mh(data,model,options,nchains);
    finalize(data,model,fit,stats,measures,temp,resultfolder,root)
    println(now())
    nothing
end

"""
cycle(nchains,fish,fixedeffects,model,data,options)

run_mh by cycling through individual parameters sequentially
"""
function cycle(nchains,fish,fixedeffects,model,data,options)
    maxtime = options.maxtime/10
    options = MHOptions(100,0,0,0,maxtime,options.temp,options.tempanneal)
    print_ll(data,model,"pre-cycle ll: ")
    t0 = time()
    nsets = length(data.nRNA)
    r = model.rates
    nalleles = model.nalleles
    G = model.G
    fittedparam = model.fittedparam
    cv = 0.02
    rateprior = model.rateprior
    while (time() - t0 < maxtime)
        for i in eachindex(fittedparam)
            model = model_rna(r,[rateprior[i]],G,nalleles,cv,[fittedparam[i]],fixedeffects,0)
            fit,_,_ = run_mh(data,model,options,nchains);
            r = get_rates(fit.parml,model)
        end
    end
    return model_rna(r,rateprior,G,nalleles,cv,fittedparam,fixedeffects,0)
end
"""
check_genename(gene,p1)

Check genename for p1
if p1 = "[" change to "("
(since swarm cannot parse "(")

"""
function check_genename(gene,p1)
    if occursin(p1,gene)
        if p1 == "["
            gene = replace(gene,"[" => "(")
            gene = replace(gene,"]" => ")")
        elseif p1 == "("
            gene = replace(gene,"(" => "]")
            gene = replace(gene,")" => "]")
        end
    end
    return gene
end

"""
print_ll(param,data,model,message="initial ll:")

compute and print initial loglikelihood
"""
function print_ll(param,data,model,message)
    ll,_ = loglikelihood(param,data,model)
    println(message,ll)
end
function print_ll(data,model,message="initial ll: ")
    ll,_ = loglikelihood(get_param(model),data,model)
    println(message,ll)
end

"""
printinfo(gene,G,cond,datafolder,infolder,resultfolder,maxtime)

print out run information
"""
function printinfo(gene,G,cond,datafolder,infolder,resultfolder,maxtime)
    println("Gene: ",gene," G: ",G," Treatment:  ",cond)
    println("data: ",datafolder)
    println("in: ", infolder," out: ",resultfolder)
    println("maxtime: ",maxtime)
end

"""
finalize(data,model,fit,stats,waic,temp,resultfolder,root)

write out run results and print out final loglikelihood and deviance
"""
function finalize(data,model,fit,stats,measures,temp,resultfolder,root)
    writefile = joinpath(root,resultfolder)
    writeall(writefile,fit,stats,measures,data,temp,model)
    println("final max ll: ",fit.llml)
    print_ll(vec(stats.meanparam),data,model,"mean ll: ")
    println("Mean fitted rates: ",stats.meanparam[:,1])
    println("Acceptance: ",fit.accept,"/",fit.total)
    println("Deviance: ",deviance(fit,data,model))
    println("rhat: ",maximum(measures.rhat))
end

#Prepare data structures
"""
    data_rna(gene::String,cond::String,datafolder,fish,label,root)
    data_rna(path,time,gene,time)
    data_rna(path,time,gene)

Load data structure
"""

function data_rna(gene::String,cond::String,datafolder::String,fish::Bool,label,root)
    if cond == "null"
        cond = ""
    end
    datafile = fish ? FISHpath(gene,cond,datafolder,root) : scRNApath(gene,cond,datafolder,root)
    data_rna(datafile,label,gene,fish)
end
function data_rna(gene::String,cond::Array,datafolder::String,fish::Bool,label,root)
    datafile = Array{String,1}(undef,length(cond))
    for i in eachindex(cond)
        datafile[i] = fish ? FISHpath(gene,cond[i],datafolder,root) : scRNApath(gene,cond[i],datafolder,root)
    end
    data_rna(datafile,label,gene,fish)
end
function data_rna(gene::String,cond::Array,datafolder::Array,fish::Bool,label,root)
    datafile = Array{String,1}(undef,length(cond))
    for i in eachindex(cond)
        datafile[i] = fish ? FISHpath(gene,cond[i],datafolder,root) : scRNApath(gene,cond[i],datafolder,root)
    end
    println(datafile)
    data_rna(datafile,label,gene,fish)
end
function data_rna(gene::String,cond::String,datafolder::String,fish::Bool,label,root,sets::Vector,time::Vector)
    if cond == "null"
        cond = ""
    end
    datafile =[]
    for set in sets
        folder = joinpath(datafolder,set)
        path = fish ? FISHpath(gene,cond,datafolder,root) : scRNApath(gene,cond,datafolder,root)
        datafile = vcat(datafile,path)
    end
    data_rna(datafile,label,times,gene,false)
end
function data_rna(path,name,time,gene::String,fish::Bool)
    len,h = histograms_rna(path,gene,fish)
    TransientRNAData(name,gene,len,time,h)
end
function data_rna(path,name,gene::String,fish::Bool)
    len,h = histograms_rna(path,gene,fish)
    RNAData(name,gene,len,h)
end
function data_rna(path,name,gene::String,fish::Array{Bool,1})
    len,h = histograms_rna(path,gene,fish)
    RNAMixedData(name,gene,len,fish,h)
end
"""
histograms_rna(path,gene)

prepare mRNA histograms for gene given data folder path
"""
function histograms_rna(path::Array,gene::String,fish::Bool)
    n = length(path)
    h = Array{Array,1}(undef,n)
    lengths = Array{Int,1}(undef,n)
    for i in eachindex(path)
        lengths[i], h[i] = histograms_rna(path[i],gene,fish)
    end
    return lengths,h
end
function histograms_rna(path::String,gene::String,fish::Bool)
    if fish
        h = read_fish(path,gene,.98)
    else
        h = read_scrna(path,.99)
        if length(h) == 0
            throw("data not found")
        end
    end
    return length(h),h
end
function histograms_rna(path::Array,gene::String,fish::Array{Bool,1})
    n = length(path)
    h = Array{Array,1}(undef,n)
    lengths = Array{Int,1}(undef,n)
    for i in eachindex(path)
        lengths[i], h[i] = histograms_rna(path[i],gene,fish[i])
    end
    return lengths,h
end

# Prepare model structures
"""
model_rna(gene,cell,G,fish,fittedparam,fixedeffects,inlabel,infolder,nsets,root,data,verbose=true)
model_rna(r::Vector,d::Vector,G::Int,nalleles::Int,propcv,fittedparam,fixedeffects,fish::Bool,method=0)
model_rna(r,G,nalleles,nsets,propcv,fittedparam,decayprior,yieldprior,method)
model_rna(r,G,nalleles,nsets,propcv,fittedparam,decayprior,method)
model_rna(r,G,nalleles,nsets,propcv,fittedparam,decayprior,noisepriors,method)

make model structure

"""
function model_rna(gene::String,cell::String,G::Int,cv,fittedparam,fixedeffects,inlabel,infolder,nsets,root,data,yieldprior,verbose=true)
    decayrate = get_decay(gene,cell,root)
    nalleles = alleles(gene,cell,root)
    if verbose
        println("alleles: ",nalleles)
        if decayrate < 0
            throw("decayrate < 0")
        else
            println("decay rate: ",decayrate)
        end
    end
    if cv <= 0
        cv = getcv(gene,G,nalleles,fittedparam,inlabel,infolder,root,verbose)
    end
    if G == 1
        ejectrate = mean_histogram(data.histRNA) * decayrate/nalleles * yieldprior
    else
        ejectrate = yieldprior
    end
    r = getr(gene,G,nalleles,decayrate,ejectrate,inlabel,infolder,nsets,root,verbose)
    model = model_rna(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,ejectrate)
    return model
end
# function model_rna2(gene::String,cell::String,G::Int,fish::Bool,cv,fittedparam,fixedeffects,inlabel,infolder,nsets,root,data,yield,ejectrate,verbose=true)
#     decayrate = get_decay(gene,cell,root)
#     nalleles = alleles(gene,cell,root)
#     if verbose
#         println("alleles: ",nalleles)
#         if decayrate < 0
#             throw("decayrate < 0")
#         else
#             println("decay rate: ",decayrate)
#         end
#     end
#     if cv <= 0
#         cv = getcv(gene,G,nalleles,fittedparam,inlabel,infolder,root,verbose)
#     end
#     if G == 1
#         ejectrate = mean_histogram(data.histRNA) * decayrate/nalleles
#     end
#     r = getr2(gene,G,nalleles,decayrate,ejectrate,yield,inlabel,infolder,nsets,root,fish,verbose)
#     model = model_rna2(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,ejectrate,yield,fish)
#     return model
# end

function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam,fixedeffects,decayprior,ejectprior)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,ejectprior)
    model_rna(r::Vector,d,G::Int,nalleles,propcv,fittedparam,fixedeffects)
end
function model_rna(r::Vector,d,G::Int,nalleles,propcv,fittedparam,fixedeffects,method=0)
    if length(fixedeffects) > 0
        model = GMfixedeffectsmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,fixedeffects,method)
    else
        model = GMmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
    end
    return model
end
# function model_rna2(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam,fixedeffects,decayprior,ejectprior,yieldprior,fish::Bool,method=0)
#     if fish
#         d = prior_rna(r,G,nsets,fittedparam,decayprior,ejectprior)
#     else
#         d = prior_rna2(r,G,nsets,fittedparam,decayprior,ejectprior,yieldprior)
#     end
#     model_rna2(r,d,G,nalleles,propcv,fittedparam,fixedeffects,fish,method)
# end
# function model_rna2(r::Vector,d,G::Int,nalleles,propcv,fittedparam,fixedeffects,fish,method)
#         if length(fixedeffects) > 0
#             model = GMfixedeffectsmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,fixedeffects,method)
#         else
#             model = GMmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
#         end
#     return model
# end

function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,decayprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,LogNormal)
    GMmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,decayprior::Float64,yieldprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,yieldprior)
    GMlossmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,fixedeffects::Tuple,decayprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior)
    GMfixedeffectsodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,fixedeffects,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,fixedeffects::Tuple,decayprior::Float64,yieldprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,yieldprior)
    GMfixedeffectslossmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,fixedeffects,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,propcv,fittedparam::Array,decayprior::Float64,noisepriors::Array,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G::Int,1,fittedparam,decayprior,noisepriors)
    if method == 1
        GMrescaledmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
    else
        GMmultimodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
    end
end
function model_delay_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,decayprior,delayprior)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,delayprior)
    GMdelaymodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),Int64}(G,nalleles,r,d,propcv,fittedparam,1)
end

"""
transient_rna(nchains,gene::String,fittedparam,cond::Vector,G::Int,maxtime::Float64,infolder::String,resultfolder,datafolder,inlabel,label,nsets,runcycle::Bool=false,samplesteps::Int=40000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

make structures for transient model fit
"""
function transient_rna(nchains,gene::String,cell,fittedparam,cond::Vector,G::Int,maxtime::Float64,infolder::String,resultfolder,datafolder,inlabel,label,nsets,runcycle::Bool=false,samplesteps::Int=40000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    data = make_data(gene,cond,G,datafolder,label,nsets,root)
    model = make_model(gene,cell,G,fittedparam,inlabel,infolder,nsets,root,data)
    param,_ = initial_proposal(model)
    return param, data, model
end

"""
getr(gene,G,nalleles,decayrate,ejectrate,inlabel,infolder,nsets::Int,root,verbose)

"""
function getr(gene,G,nalleles,decayrate,ejectrate,inlabel,infolder,nsets::Int,root,verbose)
    r = getr(gene,G,nalleles,inlabel,infolder,root,verbose)
    if ~isnothing(r)
        if length(r) == 2*G*nsets + 1
            for n in nsets
                r[2*G*n-1] *= clamp(r[2*G*nsets + 1],eps(Float64),1-eps(Float64))
            end
            r = r[1:2*G*nsets]
        end
        if length(r) == 2*G*nsets
            if verbose
                println("init rates: ",r)
            end
            return r
        end
    end
    println("No r")
    setr(G,nsets,decayrate,ejectrate)
end
# function getr2(gene,G,nalleles,decayrate,ejectrate,yield,inlabel,infolder,nsets::Int,root,fish::Bool,verbose)
#     r = getr(gene,G,nalleles,inlabel,infolder,root,verbose)
#     if ~isnothing(r)
#         if length(r) == 2*G*nsets + 1
#             for i in 1:nsets
#                 r[2*G*i-1] *= clamp(r[2*G*nsets + 1],eps(Float64),1-eps(Float64))
#             end
#             r = r[1:2*G*nsets]
#         end
#         if length(r) == 2*G*nsets
#             if verbose
#                 println(r)
#             end
#             return r
#         end
#     end
#     println("No r")
#     setr2(G,nsets,decayrate,ejectrate,yield,fish)
# end
function getr(gene,G,nalleles,inlabel,infolder,root,verbose)
    ratefile = path_Gmodel("rates",gene,G,nalleles,inlabel,infolder,root)
    if verbose
        println("rate file: ",ratefile)
    end
    if isfile(ratefile)
        r = readrates(ratefile,2)
    else
        return nothing
    end
end
function getcv(gene,G,nalleles,fittedparam,inlabel,infolder,root,verbose = true)
    paramfile = path_Gmodel("param-stats",gene,G,nalleles,inlabel,infolder,root)
    if isfile(paramfile)
        cv = read_covlogparam(paramfile)
        cv = float.(cv)
        if ~ isposdef(cv) || size(cv)[1] != length(fittedparam)
            cv = .02
        end
    else
        cv = .02
    end
    if verbose
        println("cv: ",cv)
    end
    return cv
end
"""
    get_decay(gene::String,cell::String,root::String,col::Int=2)
    get_decay(gene::String,path::String,col::Int)

    Get decay rate for gene and cell

"""
function get_decay(gene::String,cell::String,root::String,col::Int=2)
    path = get_file(root,"data/halflives",cell,"csv")
    if isnothing(path)
        println(gene," has no decay time")
        return -1.
    else
        get_decay(gene,path,col)
    end
end
function get_decay(gene::String,path::String,col::Int)
    a = nothing
    in = readdlm(path,',')
    ind = findfirst(in[:,1] .== gene)
    if ~isnothing(ind)
        a = in[ind,col]
    end
    get_decay(a,gene)
end
function get_decay(a,gene::String)
    if typeof(a) <: Number
        return get_decay(a)
    else
        println(gene," has no decay time")
        return -1.
    end
end
get_decay(a::Float64) = log(2)/a/60.

"""
    alleles(gene::String,cell::String,root::String,col::Int=3)
    alleles(gene::String,path::String,col::Int=3)

    Get allele number for gene and cell
"""
function alleles(gene::String,cell::String,root::String,col::Int=3)
    path = get_file(root,"data/alleles",cell,"csv")
    if isnothing(path)
        return 2
    else
        alleles(gene,path,col)
    end
end

function alleles(gene::String,path::String,col::Int=3)
    a = nothing
    in,h = readdlm(path,',',header=true)
    ind = findfirst(in[:,1] .== gene)
    if isnothing(ind)
        return 2
    else
        a = in[ind,col]
        return isnothing(a) ? 2 : Int(a)
    end
end

"""
prior_rna(r::Vector,G::Int,nsets::Int,propcv,fittedparam::Array,decayprior,yieldprior)
prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,noisepriors::Array)

prepare prior distribution
r[mod(1:2*G,nsets)] = rates for each model set  (i.e. rates for each set are stacked)
r[2*G*nsets + 1] == yield factor (i.e remaining after loss due to technical noise)

prior for multithresholded smFISH data
r[1:2G] = model rates
r[2G+1] = additive noise mean
r[2G + 1 + 1:length(noisepriors)] = remaining fraction after thresholding (i.e. yield)
"""
# function prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,ejectprior::Float64,f=LogNormal)
#         ind = 2*G * nsets + 1
#         if length(r) == ind
#             rm,rcv = setrate(G,nsets,decayprior,ejectprior,yieldprior)
#             if in(ind,fittedparam)
#                 return distributionBeta_array(rm[fittedparam],rcv[fittedparam],findfirst(ind.==fittedparam),f)
#             else
#                 return distribution_array(rm[fittedparam],rcv[fittedparam],f)
#             end
#         else
#             throw("rates have wrong length")
#         end
# end
# function prior_rna2(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,ejectprior::Float64,yieldprior::Float64,f=LogNormal)
#         ind = 2*G * nsets
#         if length(r) == ind
#             rm,rcv = setrate2(G,nsets,decayprior,ejectprior,yieldprior)
#             return distribution_array(rm[fittedparam],rcv[fittedparam],f)
#         else
#             throw("rates have wrong length")
#         end
# end
function prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,ejectprior::Float64,f=LogNormal)
        if length(r) == 2*G*nsets
            rm,rcv = setrate(G,nsets,decayprior,ejectprior)
            return distribution_array(rm[fittedparam],rcv[fittedparam],f)
        else
            throw("rates have wrong length")
        end
end
function prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,ejectprior::Float64,noisepriors::Array,f=LogNormal)
        if length(r) == 2*G * nsets + length(noisepriors)
            rm,rcv = setrate(G,nsets,decayprior,ejectprior,noisepriors)
            return distribution_array(rm[fittedparam],rcv[fittedparam],f)
        else
            throw("rates have wrong length")
        end
end

"""
    setrate(G::Int,nsets::Int,decayrate::Float64,ejectrate::Float64)
    setrate(G::Int,nsets::Int,decayrate,ejectrate)
    setrate(G::Int,nsets::Int,decayrate,ejectrate,yield)
    setrate(G::Int,nsets::Int,decayrate::Float64,ejectrate::Float64,noisepriors::Array)

    Set mean and cv of rates

"""

function setr(G,nsets,decayrate,ejectrate)
    r = setrate(G,nsets,decayrate,ejectrate)[1]
    println("init rates: ",r)
    return r
end
# function setr2(G,nsets,decayrate,ejectrate,yield,fish)
#     if fish
#         r = setrate(G,nsets,decayrate,ejectrate)[1]
#     else
#         r = setrate2(G,nsets,decayrate,ejectrate,yield)[1]
#     end
#     println(r)
#     return r
# end
function setrate(G::Int,nsets::Int,decayrate::Float64,ejectrate::Float64)
    rm = Vector{Float64}(undef,0)
    rc = similar(rm)
    for i in 1:nsets
        rm = vcat(rm,[repeat([0.01;0.1],(G-1));ejectrate;decayrate])
        rc = vcat(rc,[.4*ones(2*(G-1));1.;0.01])
    end
    return rm,rc
end
# function setrate(G::Int,nsets::Int,decayrate::Float64,ejectrate::Float64,yield::Float64)
#     rm,rcv = setrate(G,nsets,decayrate,ejectrate)
#     return [rm;yield],[rcv;.25]
# end
# function setrate2(G::Int,nsets::Int,decayrate::Float64,ejectprior::Float64,yieldprior::Float64)
#     rm = Vector{Float64}(undef,0)
#     rc = similar(rm)
#     for i in 1:nsets
#         rm = vcat(rm,[repeat([0.01;0.1],(G-1));ejectprior*yieldprior;decayrate])
#         rc = vcat(rc,[.4*ones(2*(G-1));1.14;0.01])
#     end
#     return rm,rc
# end
function setrate(G::Int,nsets::Int,decayrate::Float64,ejectrate::Float64,noisepriors::Array)
    rm,rc = setrate(G,nsets,decayrate,ejectrate)
    for nm in noisepriors
        rm = vcat(rm,nm)
        rc = vcat(rc,.25)
    end
    return rm,rc
end

"""
proposal_cv_rna(cv)
proposal_cv_rna(propcv,fittedparam)

set proposal cv to a vector or matrix
"""
proposal_cv_rna(cv) = typeof(cv) == Float64 ? propcv*ones(length(fittedparam)) : propcv

proposal_cv_rna(propcv,fittedparam) = typeof(propcv) == Float64 ? propcv*ones(length(fittedparam)) : propcv

"""
fittedparam_rna(G,nsets,loss)
fittedparam_rna(G,nsets)

select all parameters to be fitted except decay rates
yieldfactor is last parameter in loss models
"""
function fittedparam_rna(G,nsets,loss)
    fp = fittedparam_rna(G,nsets)
    if loss == 1
        return vcat(fp,nsets*2*G+1)
    else
        return fp
    end
end
function fittedparam_rna(G,nsets)
    nrates = 2*G  # total number of rates in GM model
    k = nrates - 1  # adjust all rate parameters except decay time
    fp = Array{Int,1}(undef,nsets*k)
    for j in 0:nsets-1
        for i in 1:k
            fp[k*j + i] = nrates*j + i
        end
    end
    return fp
end

"""
rescale_rate_rna(r,G,decayrate::Float64)

set new decay rate and rescale
transition rates such that steady state distribution is the same
"""
function rescale_rate_rna(r,G,decayrate::Float64)
    rnew = copy(r)
    if mod(length(r),2*G) == 0
        rnew *= decayrate/r[2*G]
    else
        stride = fld(length(r),fld(length(r),2*G))
        for i in 0:stride:length(r)-1
            rnew[i+1:i+2*G] *= decayrate/r[2*G]
        end
    end
    return rnew
end

# Read in data and construct histograms
"""
read_scrna(filename::String,yield::Float64=.99,nhistmax::Int=1000)
Construct mRNA count per cell histogram array of a gene
"""
function read_scrna(filename::String,threshold::Float64=.99,nhistmax::Int=500)
    if isfile(filename) && filesize(filename) > 0
        x = readdlm(filename)[:,1]
        x = truncate_histogram(x,threshold,nhistmax)
        if x == 0
            dataFISH = Array{Int,1}(undef,0)
        else
            dataFISH = x
        end
        return dataFISH
    else
        # println("data file not found")
        return Array{Int,1}(undef,0)
    end
end

"""
read_fish(path,gene,threshold)
Read in FISH data from 7timepoint type folders

"""
function read_fish(path::String,cond::String,threshold::Float64=.98,maxlength = 1800)
    xr = zeros(maxlength)
    lx = 0
    for (root,dirs,files) in walkdir(path)
        for file in files
            target = joinpath(root, file)
            if occursin(cond,target) && occursin("cellular",target)
                # println(target)
                x1 = readdlm(target)[:,1]
                x1 = truncate_histogram(x1,threshold,maxlength)
                lx = length(x1)
                # println(lx)
                xr[1:min(lx,maxlength)] += x1[1:min(lx,maxlength)]
            end
        end
    end
    return truncate_histogram(xr,1.0,maxlength)
end

function read_fish(path1::String,cond1::String,path2::String,cond2::String,threshold::Float64=.98)
    x1 = read_fish(path1,cond1,threshold)
    x2 = read_fish(path2,cond2,threshold)
    combine_histogram(x1,x2)
end

# functions to build paths to RNA/FISH data and results
"""
scRNApath(gene::String,cond::String,datapath::String,root::String)

generate path to scRNA data for given gene and condition cond
    data files must have format genename_cond.txt
"""
function scRNApath(gene::String,cond::String,datapath::String,root::String)
    datapath = joinpath(root,datapath)
    if cond == ""
        joinpath(datapath,gene * ".txt")
    else
        joinpath(datapath,gene * "_" * cond * ".txt")
    end
end
scRNApath(gene,cond,datapath) = joinpath(datapath,gene * "_" * cond * ".txt")

"""
FISHpath(gene,cond,datapath,root)

generate path to FISH data
"""
# FISHpath(gene,cond,datapath,root) = joinpath(joinpath(joinpath(root,datapath),gene),cond)
FISHpath(gene,cond,datapath,root) = joinpath(root,datapath,gene,cond)
FISHpath(gene,cond,datapath) = joinpath(datapath,gene,cond)

"""
ratepath_Gmodel(gene::String,cond::String,G::Int,nalleles::Int,label,folder,root)

"""
function ratepath_Gmodel(gene::String,cond::String,G::Int,nalleles::Int,label,folder,root)
    path_Gmodel("rates",gene,G,nalleles,label * "-" * cond,folder,root)
end

"""
path_Gmodel(type,gene::String,G::Int,nalleles::Int,label::String,folder,root)
"""
function path_Gmodel(type,gene::String,G::Int,nalleles::Int,label::String,folder,root)
    filelabel = label  * "_" * gene *  "_" * "$G" * "_" * "$nalleles" * ".txt"
    ratefile = type * "_" * filelabel
    joinpath(root,folder,ratefile)
end

"""
stats_rna(genes::Vector,conds,datapath,threshold=.99)
"""
function stats_rna(genes::Vector,conds,datapath,threshold=.99)
    g = Vector{String}(undef,0)
    z = Vector{Float64}(undef,0)
    h1 = Vector{Float64}(undef,0)
    h2 = Vector{Float64}(undef,0)
    for gene in genes
        t,m1,m2 = expression_rna(gene,conds,datapath,threshold)
        push!(g,gene)
        push!(z,t)
        push!(h1,m2)
        push!(h2,m1)
    end
    return g,z,h1,h2
end

"""
expression_rna(gene,conds::Vector,folder::String,threshold=.99)
"""
function expression_rna(gene,conds::Vector,folder::String,threshold=.99)
    h = Vector{Vector}(undef,2)
    for i in eachindex(conds)
        datapath = scRNApath(gene,conds[i],folder)
        h[i] = read_scrna(datapath,threshold)
    end
    if length(h[1]) > 0 && length(h[2]) > 0
        return log_2sample(h[1],h[2]), mean_histogram(h[1]), mean_histogram(h[2])
    else
        return 0,0,0
    end
end

"""
expression_rna(gene,conds::Vector,folder::Vector,threshold=.99)
"""
function expression_rna(gene,conds::Vector,folder::Vector,threshold=.99)
    h = Vector{Vector}(undef,2)
    for i in eachindex(conds)
        datapath = scRNApath(gene,conds[i],folder[i])
        h[i] = read_scrna(datapath,threshold)
    end
    if length(h[1]) > 0 && length(h[2]) > 0
        return log_2sample(h[1],h[2]), mean_histogram(h[1]), mean_histogram(h[2])
    else
        return 0,0,0
    end
end

"""
expression_rna(genes::Vector,cond::String,folder,threshold=.99)

"""
function expression_rna(genes::Vector,cond::String,folder,threshold=.99)
    h1 = Array{Any,2}(undef,0,2)
    for gene in genes
        datapath = scRNApath(gene,cond,folder)
        h = read_scrna(datapath,threshold)
        if length(h) > 0
            h1 = vcat(h1,[gene mean_histogram(h)])
        end
    end
    return h1
end

"""
rna_fish(gene,cond,fishfolder,rnafolder,yield,root)

output RNA histogram and downsampled FISH histogram with loss
"""
function rna_fish(gene,cond,fishfolder,rnafolder,yield,root)
    datarna = histograms_rna(scRNApath(gene,cond,rnafolder,root),gene,false)
    f=reduce_fish(gene,cond,datarna[1],fishfolder,yield,root)
    return datarna[2],f
end

"""
reduce_fish(gene,cond,nhist,fishfolder,yield,root)

sample fish histogram with loss (1-yield)
"""
function reduce_fish(gene,cond,nhist,fishfolder,yield,root)
    fish = histograms_rna(FISHpath(gene,cond,fishfolder,root),gene,true)
    technical_loss(fish[2],yield,nhist)
    fish[2]
end


text:
##==============================================================================
## Overloading show
##==============================================================================
# Base.show_default(io, v)
Base.show(io::IO, ::MIME"text/plain", v::DFGVariable) = show(IOContext(io, :limit=>true, :compact=>true), v)

Base.show(io::IO, ::MIME"text/plain", f::DFGFactor) = show(IOContext(io, :limit=>true, :compact=>true), f)

function Base.show(io::IO, ::MIME"text/plain", dfg::AbstractDFG)
    summary(io, dfg)
    println("\n  UserId: ", dfg.userId)
    println("  RobotId: ", dfg.robotId)
    println("  SessionId: ", dfg.sessionId)
    println("  Description: ", dfg.description)
    println("  Nr variables: ", length(ls(dfg)))
    println("  Nr factors: ",length(lsf(dfg)))
    println("  User Data: ", keys(dfg.userData))
    println("  Robot Data: ", keys(dfg.robotData))
    println("  Session Data: ", keys(dfg.sessionData))
end

#default for Atom/Juno
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::Union{AbstractDFG, DFGVariable, DFGFactor}) = x

##==============================================================================
## Printing Variables and Factors
##==============================================================================

printVariable(vert::DFGVariable; kwargs...) = printVariable(stdout::IO, vert; kwargs...)

function printVariable(io::IO, vert::DFGVariable;
                       short::Bool=false,
                       compact::Bool=true,
                       limit::Bool=true,
                       skipfields::Vector{Symbol}=Symbol[],
                       solveKeys::Vector{Symbol}=Symbol[])

    ioc = IOContext(io, :limit=>limit, :compact=>compact)

    if short
        printstyled(ioc, summary(vert),"\n", bold=true)
        vnd = getSolverData(vert)
        println(ioc, "label: $(vert.label)")
        println(ioc, "tags: $(getTags(vert))")
        println(ioc, "size marginal samples: $(size(vnd.val))")
        println(ioc, "kde bandwidths: $((vnd.bw)[:,1])")
        if 0 < length(getPPEDict(vert))
            println(ioc, "PPE.suggested: $(round.(getPPE(vert).suggested,digits=4))")
        else
            println(ioc, "No PPEs")
        end
        # println(ioc, "kde max: $(round.(getKDEMax(getKDE(vnd)),digits=4))")
        # println(ioc, "kde max: $(round.(getKDEMax(getKDE(vnd)),digits=4))")
    else

        printstyled(ioc, summary(vert),"\n", bold=true, color=:blue)

        :solver in skipfields && push!(skipfields, :solverDataDict)
        :ppe in skipfields && push!(skipfields, :ppeDict)

        t = typeof(vert)
        fields = setdiff(fieldnames(t), skipfields)
        nf = nfields(vert)

        for f in fields
            printstyled(ioc, f,":\n", color=:blue)
            show(ioc, getproperty(vert, f))
            println(ioc)
        end
    end
end

printFactor(vert::DFGFactor; kwargs...) = printFactor(stdout::IO, vert; kwargs...)
function printFactor(io::IO, vert::DFGFactor;
                     short::Bool=false,
                     compact::Bool=true,
                     limit::Bool=true,
                     skipfields::Vector{Symbol}=Symbol[])

    ioc = IOContext(io, :limit=>limit, :compact=>compact)

    if short
        printstyled(ioc, summary(vert),"\n", bold=true)
        println(ioc, "label: ", vert.label)
        println(ioc, "timestamp: ", vert.timestamp)
        println(ioc, "tags: ", vert.tags)
        println(ioc, "solvable: ", vert.solvable)
        println(ioc, "VariableOrder: ", vert._variableOrderSymbols)
    else

        printstyled(ioc, summary(vert),"\n", bold=true, color=:blue)

        :solver in skipfields && push!(skipfields, :solverData)

        t = typeof(vert)
        fields = setdiff(fieldnames(t), skipfields)
        nf = nfields(vert)

        for f in fields
            printstyled(ioc, f,":\n", color=:blue)
            show(ioc, getproperty(vert, f))
            println(ioc)
        end
    end
end

"""
    $SIGNATURES

Display and return to console the user factor identified by tag name.
"""
printFactor(dfg::AbstractDFG, sym::Symbol) = print(getFactor(dfg, sym))

"""
   $SIGNATURES

Display the content of `VariableNodeData` to console for a given factor graph and variable tag`::Symbol`.

Dev Notes
- TODO split as two show macros between AMP and DFG
"""
printVariable(dfg::AbstractDFG, sym::Symbol) = print(getVariable(dfg, sym))

printNode(dfg::AbstractDFG, sym::Symbol) = isVariable(dfg,sym) ? printVariable(dfg, sym) : printFactor(dfg, sym)


text:
@testset "Linked read datastores" begin

    function ucdavis_tag(rec::FASTQ.Record)
        id = FASTQ.identifier(rec)
        tag = zero(UInt32)
        @inbounds for i in 1:16
            tag <<= 2
            letter = id[i]
            if letter == 'C'
                tag = tag + 1
            elseif letter == 'G'
                tag = tag + 2
            elseif letter == 'T'
                tag = tag + 3
            elseif letter != 'A'
                tag = zero(UInt32)
                break
            end
        end
        return tag
    end

    function get_fastq_seqs(::Type{A}, r1, r2, maxlen) where {A<:DNAAlphabet}
        r1s = open(FASTQ.Reader, r1) do rdr
            collect(rdr)
        end
        r2s = open(FASTQ.Reader, r2) do rdr
            collect(rdr)
        end
        tags = map(r1s) do rec
            ucdavis_tag(rec)
        end
        keep = tags .!= zero(UInt32)
        tags = tags[keep]
        r1s = FASTQ.sequence.(LongSequence{A}, r1s[keep])
        r2s = FASTQ.sequence.(LongSequence{A}, r2s[keep])
        for seq in r1s
            if length(seq) > maxlen
                resize!(seq, maxlen)
            end
        end
        for seq in r1s
            if length(seq) > maxlen
                resize!(seq, maxlen)
            end
        end
        order = sortperm(tags)
        return r1s[order], r2s[order]
    end

    function check_round_trip(::Type{A}, R1, R2, maxlen, chunksize = 1000000) where {A<:DNAAlphabet}
        r1_seqs, r2_seqs = get_fastq_seqs(A, R1, R2, maxlen)

        fqa = open(FASTQ.Reader, R1)
        fqb = open(FASTQ.Reader, R2)

        ds = LinkedReads{A}(fqa, fqb, "10xtest", Symbol("ucdavis-test"), UCDavis10x, maxlen, chunksize)
        ds2 = open(LinkedReads{A}, "10xtest.lrseq")
        @test ReadDatastores.deduce_datastore_type("10xtest.lrseq") == LinkedReads{A}
        ds_seqs = collect(ds)
        ds2_seqs = collect(ds2)

        return ds_seqs[1:2:end] == ds2_seqs[1:2:end] == r1_seqs && ds_seqs[2:2:end] == ds2_seqs[2:2:end] == r2_seqs
    end

    function check_show(ds, msg)
        buf = IOBuffer()
        show(buf, ds)
        return String(take!(buf)) == msg
    end

    @test check_round_trip(DNAAlphabet{2}, "10x_tester_R1.fastq", "10x_tester_R2.fastq", 250, 10)
    @test check_round_trip(DNAAlphabet{4}, "10x_tester_R1.fastq", "10x_tester_R2.fastq", 250, 10)

    ds = open(LinkedReads{DNAAlphabet{4}}, "10xtest.lrseq")
    @test ReadDatastores.name(ds) == Symbol("ucdavis-test")
    @test ReadDatastores.max_read_length(ds) == 250
    @test check_show(ds, "Linked Read Datastore 'ucdavis-test': 166 reads (83 pairs)")
    @test firstindex(ds) == 1
    @test lastindex(ds) == 166
    @test Base.IteratorSize(ds) == Base.HasLength()
    @test Base.IteratorEltype(ds) == Base.HasEltype()
    @test Base.eltype(ds) == LongSequence{DNAAlphabet{4}}
    @test_throws BoundsError ds[200]
    @test_throws BoundsError buffer(ds)[200]
    @test_throws ErrorException buffer(ds, 1)
    @test_throws ErrorException buffer(ds, ReadDatastores._bytes_per_read(ds) - 1)
    @test_throws BoundsError load_sequence!(ds, 200, dna"")
    @test collect(ds) == collect(buffer(ds)) == open(LinkedReads{DNAAlphabet{4}}, "10xtest.lrseq") do ds
        collect(ds)
    end

    @test_throws ReadDatastores.DatastoreVersionError{LinkedReads{DNAAlphabet{2}}} open(LinkedReads{DNAAlphabet{2}}, "10xtest-old.lrseq")
    @test_throws ReadDatastores.DatastoreEncodingError{LinkedReads{DNAAlphabet{2}}} open(LinkedReads{DNAAlphabet{2}}, "10xtest.lrseq")
end


text:
"""
    mrca(tree::AbstractTree, target)

Returns the node within `tree` that is the Most Recent Common Ancestor of the
leaves (or internal nodes) defined by `target`. `target` can be an iterator
over nodes or an AbstractArray of nodes. The return value has the same node
type as the elements of `target`.

### Examples
≡≡≡≡≡≡≡≡≡≡≡

julia> tree = open(parsenewick, Phylo.path("H1N1.newick"))
RootedTree with 507 tips, 1013 nodes and 1012 branches.
Leaf names are 227, 294, 295, 110, 390, ... [501 omitted] ... and 418

julia> tips = rand(collect(nodefilter(isleaf, tree)), 3)
3-element Vector{LinkNode{OneRoot, String, Dict{String, Any}, LinkBranch{OneRoot, String, Dict{String, Any}, Float64}}}:
 LinkNode 153, a tip of the tree with an incoming connection (branch 888).

 LinkNode 120, a tip of the tree with an incoming connection (branch 195).

 LinkNode 504, a tip of the tree with an incoming connection (branch 44).

julia> mrca(tree, tips)
LinkNode Node 1003, an internal node with 1 inbound and 2 outbound connections (branches 1001 and 999, 1000)
"""
function mrca(tree::AbstractTree, target)
    ancestors = getancestors(tree, first(target))
    ranks = Dict(j => i for (i,j) in enumerate(ancestors))
    checked = Set(ancestors)
    oldest = 1
    for species in target
        while !(species ∈ checked)
            push!(checked, species)
            species = getparent(tree, species)
        end
        oldest = max(oldest, get(ranks, species, 0))
    end
    ancestors[oldest]
end

"""
    nodeheights(tree::Phylo.AbstractTree; onlyleaves = false, noleaves = false)

Returns an AxisArray of the height of all nodes in `tree` over the root node.
`onlyleaves` and `noleaves` filter the nodes for leaves and internal nodes,
respectively
"""
function nodeheights(tree::Phylo.AbstractTree; onlyleaves = false, noleaves = false)
    function findheights!(clade::String, parentheight::Float64 = 0.0)
        myheight = parentheight
        if hasinbound(tree, clade)
             myheight += getlength(tree, getinbound(tree, clade))
        end
        height[clade] = myheight
        for ch in getchildren(tree, clade)
            findheights!(ch, myheight)
        end
    end
    names = getnodename.((tree, ), traversal(tree, preorder))
    height = AxisArray(Vector{Float64}(undef, nnodes(tree)), x = names)
    root = getnodename(tree, getroot(tree))
    findheights!(root)
    onlyleaves && return height[filter(t->isleaf(tree, t), names)]
    noleaves && return height[filter(t->!isleaf(tree, t), names)]
    height
end

"""
    distance(tree::AbstractTree, node1, node2)

Distance between two nodes on a tree
"""
function distance(tree::AbstractTree, node1, node2)
    branches = branchroute(tree, getnode(tree, node1), getnode(tree, node2))
    return mapreduce(branch -> getlength(tree, branch), +, branches; init = 0.0)
end

"""
    distances(tree::AbstractTree)

Pairwise distances between all leaf nodes on a tree
"""
function distances(tree::AbstractTree)
    leaves = getleaves(tree)
    names = getnodename.(tree, leaves)
    return AxisArray([distance(tree, li, lj) for li in leaves, lj in leaves],
                     Axis{:x}(names), Axis{:y}(names))
end

"""
    heighttoroot(tree::AbstractTree, node)

Height of a node of the tree above the root
"""
function heighttoroot(tree::AbstractTree{OneTree, <:Rooted}, node)
    return mapreduce(branch -> getlength(tree, branch), +,
                     branchhistory(tree, node); init = 0.0)
end

"""
    heightstoroot(tree::AbstractTree)

Height of all of the leaves of the tree above the root
"""
function heightstoroot(tree::AbstractTree)
    return nodeheights(tree, onlyleaves = true)
end


text:
push!(LOAD_PATH,"../src")
using Documenter: Documenter, makedocs, deploydocs, doctest, DocMeta
using ReversedStrings
using Test

docdir = joinpath(dirname(pathof(ReversedStrings)),"../docs/src/")
mandir = joinpath(docdir,"man")
DocMeta.setdocmeta!(ReversedStrings, :DocTestSetup, quote
    using TextParse
    using ReversedStrings
end; recursive=true)

makedocs(;
         source=docdir,
         modules=[ReversedStrings],
         authors="Gregor Kappler",
         repo="https://github.com/gkappler/ReversedStrings.jl/blob/{commit}{path}#L{line}",
         sitename="ReversedStrings.jl",
         format=Documenter.HTML(;
                                prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://gkappler.github.io/ReversedStrings.jl",
                                assets=String[],
                                ),
         pages=[
             "Home" => "index.md"
             # "Developer Guide" => "developer.md"
         ],
         )

deploydocs(;
    repo="github.com/gkappler/ReversedStrings.jl",
)


text:
using FileIO
using Parameters

include("random_matrices.jl")
include("network_sims_gpu.jl")
include("LS_calculations.jl")

#------------------------------------------------------------------------------
# This is the part where we find which eta and network realization to

J_seeds = collect(5000:5019)
etas = collect(-1.0f0:0.1f0:1.0f0)

ind = parse(Int64, ARGS[1])
ps = [(eta, J_seed) for eta in etas, J_seed in J_seeds]
eta, J_seed = ps[ind]
#------------------------------------------------------------------------------

# Where the results will be saved
dir = "results1"

N = 3000
g = 1.5f0
s = 1f0

h0_seed = 6000
dh_seed = 7000
Q0_seed = 8000

J = generate_J(N, eta, seed=J_seed)
h = generate_nv(N, Float32, seed=h0_seed)
dh = generate_nv(N, Float32, seed=dh_seed)
Q = generate_ONB(N, Float32, seed=Q0_seed)

dt = 1f-2
t_bound = 2f3
t_ONS = 1f0
wait = 2f2

@time Lspec, max_LE = sim_lyapunov_gpu(J, g, h, dh, Q; s=s, dt=dt,
		t_bound=t_bound, t_ONS=t_ONS, wait=wait, return_max_LE=true,
		printstuff=false, p_factor=1e-3, scale_limit=100)

dimensionality = KY_attractor_dim(Lspec)
entropy_rate = KS_entropy_rate(Lspec)

fn = "Lspec_E$(ind)_N=$N,g=$g,eta=$eta,J_seed=$J_seed.jld2"
out = Dict{String,Any}()
@pack! out = N, g, eta, s, J_seed, h0_seed, dh_seed, Q0_seed, dt, t_bound,
		t_ONS, wait, Lspec, max_LE, dimensionality, entropy_rate

save(joinpath(dir, fn), out)


text:
###########################################
### Building the meta-rule Jazz Grammar ###
###########################################

function make_meta_rule{T<:Tuple}(name::String, catslist::Vector{T}, f=x->x)
    rules = [ContextFreeRule(cats[1], map(f, cats[2:end])) for cats in catslist]
    eval( :( $(Symbol(name)) = MetaRule($name, $rules) ) )
end

const all_category_rule_symbols = [
    :prolongation;
    :diatonic_prep;
    :true_modulation;
    :parallel_modulation;
    :full_mode_change;
    :partial_mode_change;
    :plagal_prep;
    :octa_dom_subs3;
    :octa_dom_subs6;
    :octa_dom_subs9;
    :diatonic_subs;
    :min_M6_subs
]

# back-off order is from left to right
dependent_components(cat::Category) = (scaledeg(cat),)

function make_jazz_grammar(Score :: DataType = LogProb, category_rules_symbols=all_category_rule_symbols)
    K_M = [Key(r,"major") for r in 1:12]
    K_m = [Key(r,"minor") for r in 1:12]
    K = [K_M; K_m]
    S = [ScaleDeg(deg) for deg in 1:7]
    B = [true, false]
    C = [Category(s,k,b) for s in S for k in K for b in B]

    make_meta_rule("prolongation",
        [(c,c,c) for c in C], notdomsubst)
    make_meta_rule("diatonic_prep",
        [(c,transpose(c, 4), c) for c in C], notdomsubst)
    make_meta_rule("true_modulation",
        [(c, Category("V", modulate(c), false), c) for c in C if numeral(scaledeg(c)) != "I"], notdomsubst)
    make_meta_rule("parallel_modulation",
        [(c, Category("V", change_mode(modulate(c)), false), c) for c in C], notdomsubst)
    make_meta_rule("full_mode_change",
        [(Category("I", k, false), Category("V", change_mode(k), false), Category("I", change_mode(k), false)) for k in K])
    make_meta_rule("partial_mode_change",
        [(Category("I", k, false), Category("V", k, false), Category("I", change_mode(k), false)) for k in K])
    make_meta_rule("plagal_prep",
        [(Category("I", k, false), Category("IV", k, false), Category("I", k, false)) for k in K])
    make_meta_rule("octa_dom_subs3",
        [(Category("V", k, false), Category("V", k+3, true)) for k in K])
    make_meta_rule("octa_dom_subs6",
        [(Category("V", k, false), Category("V", k+6, true)) for k in K])
    make_meta_rule("octa_dom_subs9",
        [(Category("V", k, false), Category("V", k+9, true)) for k in K])
    make_meta_rule("diatonic_subs",
        [[(Category("I", k, false), Category("VI", k, false)) for k in K_M];
         [(Category("I", k, false), Category("III", k, false)) for k in K_m];
         [(Category("II", k, false), Category("IV", k, false)) for k in K];
         [(Category("V", k, b), Category("VII", k, false)) for k in K_M for b in B]])
    make_meta_rule("min_M6_subs",
        [(Category("V", k, b), Category("II", k, true)) for k in K_M for b in B])

    category_rules = map(eval, category_rules_symbols)

    make_meta_rule("triadic_termination", [(c, Terminal(c)) for c in C])
    make_meta_rule("alt_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "alt")) for k in K for b in B])
    make_meta_rule("aug_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "aug")) for k in K for b in B])
    make_meta_rule("sus_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "sus4")) for k in K for b in B])
    make_meta_rule("min_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "maj")) for k in K_m for b in B])
    # make_meta_rule("min_M6_dom", ...)

    terminal_rules = [triadic_termination,
                      alt_dominant,
                      aug_dominant,
                      sus_dominant,
                      min_dominant]

    startsymbols = [Category("I",k,false) for k in K]

    Grammar(category_rules, terminal_rules, startsymbols, Score, dependent_components)
end

####################################################
### Building the plain context-free Jazz Grammar ###
####################################################

function make_PCFG_rule_list{T<:Tuple}(name::String, catslist::Vector{T}, f=x->x)
    rules = [ContextFreeRule(cats[1], map(f, cats[2:end])) for cats in catslist]
    eval( :( $(Symbol(name)) = $rules) )
end

function make_plain_PCFG_jazz_grammar(Score :: DataType = LogProb, category_rules_symbols=all_category_rule_symbols)
    K_M = [Key(r,"major") for r in 1:12]
    K_m = [Key(r,"minor") for r in 1:12]
    K = [K_M; K_m]
    S = [ScaleDeg(deg) for deg in 1:7]
    B = [true, false]
    C = [Category(s,k,b) for s in S for k in K for b in B]

    make_PCFG_rule_list("prolongation",
        [(c,c,c) for c in C], notdomsubst)
    make_PCFG_rule_list("diatonic_prep",
        [(c,transpose(c, 4), c) for c in C], notdomsubst)
    make_PCFG_rule_list("true_modulation",
        [(c, Category("V", modulate(c), false), c) for c in C if numeral(scaledeg(c)) != "I"], notdomsubst)
    make_PCFG_rule_list("parallel_modulation",
        [(c, Category("V", change_mode(modulate(c)), false), c) for c in C], notdomsubst)
    make_PCFG_rule_list("full_mode_change",
        [(Category("I", k, false), Category("V", change_mode(k), false), Category("I", change_mode(k), false)) for k in K])
    make_PCFG_rule_list("partial_mode_change",
        [(Category("I", k, false), Category("V", k, false), Category("I", change_mode(k), false)) for k in K])
    make_PCFG_rule_list("plagal_prep",
        [(Category("I", k, false), Category("IV", k, false), Category("I", k, false)) for k in K])
    make_PCFG_rule_list("octa_dom_subs3",
        [(Category("V", k, false), Category("V", k+3, true)) for k in K])
    make_PCFG_rule_list("octa_dom_subs6",
        [(Category("V", k, false), Category("V", k+6, true)) for k in K])
    make_PCFG_rule_list("octa_dom_subs9",
        [(Category("V", k, false), Category("V", k+9, true)) for k in K])
    make_PCFG_rule_list("diatonic_subs",
        [[(Category("I", k, false), Category("VI", k, false)) for k in K_M];
         [(Category("I", k, false), Category("III", k, false)) for k in K_m];
         [(Category("II", k, false), Category("IV", k, false)) for k in K];
         [(Category("V", k, b), Category("VII", k, false)) for k in K_M for b in B]])
    make_PCFG_rule_list("min_M6_subs",
        [(Category("V", k, b), Category("II", k, true)) for k in K_M for b in B])

    category_rules = vcat(map(eval, category_rules_symbols)...)

    make_PCFG_rule_list("triadic_termination", [(c, Terminal(c)) for c in C])
    make_PCFG_rule_list("alt_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "alt")) for k in K for b in B])
    make_PCFG_rule_list("aug_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "aug")) for k in K for b in B])
    make_PCFG_rule_list("sus_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "sus4")) for k in K for b in B])
    make_PCFG_rule_list("min_dominant", [(Category("V",k,b), Terminal(Category("V",k,b), "maj")) for k in K_m for b in B])
    # make_meta_rule("min_M6_dom", ...)

    terminal_rules = [triadic_termination;
                      alt_dominant;
                      aug_dominant;
                      sus_dominant;
                      min_dominant]

    startsymbols = [Category("I",k,false) for k in K]

    Grammar(category_rules, terminal_rules, startsymbols, Score, dependent_components)
end


text:
"""
    isboundary!(img::AbstractArray; background = 0, dims = coords_spatial(A), kwargs...)

Finds the boundaries that are just within each object, replacing the original image.
`background` is the scalar value of the background pixels which will not be marked as boundaries.
Keyword arguments are passed to `extreme_filter` which include `dims` indicating the dimension(s)
over which to discover boundaries.

See out-of-place version, [`isboundary`](@ref), for examples.
"""
function isboundary!(img::AbstractArray{Bool,N}; background=false, kwargs...) where {N}
    # Find regions where there is at least one background pixel
    # centered on a foreground pixel
    if background == true
        # background_img = img
        # foreground_img = .!img
        img .= .!img .& extreme_filter(|, img; kwargs...)
        # Alternatively, with ImageFiltering.jl
        # img .= .!img .& mapwindow(any, img, (3,3))
    elseif background == false
        # background_img = .!img
        # foreground_img = img
        img .&= extreme_filter(|, .!img; kwargs...)
        # Alternatively, with ImageFiltering.jl
        # img .&= mapwindow(any, .!img, (3,3))
    else
        # This should be the same as isboundary_thick
        # background is neither true or false, use more generic algorithm
        background = Int8(-1) # Any value other than true or false will do
        #allequal(x,y) = ifelse(x == y, x, background)
        allequal(x, y) = x == y ? x : background
        # the entire image is the foreground
        img .= (extreme_filter!(allequal, similar(img, Int8), img; kwargs...) .== background)
    end
    return img
end
function isboundary!(img::AbstractArray{T}; background=zero(T), kwargs...) where {T}
    # Find regions where that are not homogeneous (all equal)
    # centered on a foreground pixel
    # background_img = img .== background
    # foreground_img = img .!= background
    #allequal(x,y) = ifelse(x == y, x, background)
    allequal(x, y) = x == y ? x : background
    m = extreme_filter(allequal, img; kwargs...)
    img .= (m .== background) .& (img .!= background)
    # Alternatively, with ImageFiltering.jl
    # allequal(x) = all(==(first(x)), @view(x[2:end]))
    # .!mapwindow(allequal, A, (3,3)) .& (A .!= background)
    return img
end

"""
    isboundary(img::AbstractArray; background = 0, dims = coords_spatial(A), kwargs...)

Finds the boundaries that are just within each object.
`background` is the scalar value of the background pixels which will not be marked as boundaries.
Keyword arguments are passed to `extremefilt!` which include `dims` indicating the dimension(s)
over which to discover boundaries.

See also its in-place version [`isboundary!`](@ref) and the alternative version that finds thick boundaries, [`isboundary_thick`](@ref ImageMorphology.isboundary_thick).

# Examples

```@meta
DocTestSetup = quote
    import ImageMorphology: isboundary
end
```

```jldoctest
julia> A = zeros($Int, 16, 16); A[4:8, 4:8] .= 5; A[4:8, 9:12] .= 6; A[10:12,13:15] .= 3; A[10:12,3:6] .= 9; A
16×16 Matrix{$Int}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> isboundary(A)
16×16 Matrix{$Int}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> isboundary(A .!= 0)
16×16 BitMatrix:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> isboundary(A .!= 0; dims = 1)
16×16 BitMatrix:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> isboundary(A .!= 0; dims = 2)
16×16 BitMatrix:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
```
"""
function isboundary(img::AbstractArray{T}; background=zero(T), kwargs...) where {T}
    return isboundary!(copy(img); background=background, kwargs...)
end

## Alternate implementations using dilate and erode
# The implementations below are based on scikit-image

#=
Copyright (C) 2021, Mark Kittisopikul and JuliaImages Github Organization
Copyright (C) 2019, the scikit-image team
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

 1. Redistributions of source code must retain the above copyright
    notice, this list of conditions and the following disclaimer.
 2. Redistributions in binary form must reproduce the above copyright
    notice, this list of conditions and the following disclaimer in
    the documentation and/or other materials provided with the
    distribution.
 3. Neither the name of skimage nor the names of its contributors may be
    used to endorse or promote products derived from this software without
    specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
=#

# Specific implementation of `isboundary_thick` using dilate and erode similar to scikit-image, uses extreme_filter twice
function isboundary_thick_dilate_erode(img::AbstractArray; kwargs...)
    # https://github.com/scikit-image/scikit-image/blob/d44ceda6241cb23a22dc8abf09a05090ed14da7f/skimage/segmentation/boundaries.py#L165-L167
    return dilate(img; kwargs...) .!= erode(img; kwargs...)
end
function isboundary_dilate_erode(img::AbstractArray{T}; background=zero(T), kwargs...) where {T}
    # https://github.com/scikit-image/scikit-image/blob/d44ceda6241cb23a22dc8abf09a05090ed14da7f/skimage/segmentation/boundaries.py#L165-L170
    thick_boundaries = isboundary_thick_dilate_erode(img; kwargs...)
    foreground_img = img .!= background
    return thick_boundaries .& foreground_img
end
"""
    isboundary_thick(img::AbstractArray; dims = coords_spatial(img), kwargs...)

Find thick boundaries that are just outside and just inside the objects.
This is a union of the inner and outer boundaries.
Keyword `dims` indicates over which dimensions to look for boundaries. This `dims` and additional keywords `kwargs` are passed to [`extreme_filter`](@ref).

See also [`isboundary`](@ref) which just yields the inner boundaries.

# Examples

```@meta
DocTestSetup = quote
    import ImageMorphology: isboundary_thick
end
```

```jldoctest
julia> A = zeros($Int, 16, 16); A[4:8, 4:8] .= 5; A[4:8, 9:12] .= 6; A[10:12,13:15] .= 3; A[10:12,3:6] .= 9; A
16×16 Matrix{$Int}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> isboundary_thick(A)
16×16 BitMatrix:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0
 0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0
 0  0  1  1  0  0  0  1  1  0  0  1  1  0  0  0
 0  0  1  1  0  0  0  1  1  0  0  1  1  0  0  0
 0  0  1  1  0  0  0  1  1  0  0  1  1  0  0  0
 0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0
 0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
 0  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1
 0  1  1  0  0  1  1  0  0  0  0  1  1  0  1  1
 0  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1
 0  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

 julia> isboundary_thick(A) .& (A .!= 0)
16×16 BitMatrix:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

 julia> isboundary_thick(A) == isboundary(A; background = -1)
true

julia> isboundary_thick(A) .& (A .!= 0) == isboundary(A) # inner boundaries
true

julia> isboundary_thick(A .!= 0) .& (A .== 0)  == isboundary(A .== 0) # outer boundaries
true
 ```
"""
const isboundary_thick = isboundary_thick_dilate_erode