text:
using StaticArrays
using Simplices.delaunayn

"""
A triangulation of a cloud of embedded points into disjoint simplices.

```julia
# Triangulate a set of random points in 3D space.
t = triangulate(rand(20, 3))

# Refine triangulation until all simplices are below the mean radius of the original
# triangulation.
target_radius = mean(t.radii)
refine_variable_k!(t, target_radius)
```

"""
abstract type AbstractTriangulation{D, N} <: Partition end

"""
    Triangulation{D <:Int, N <:Int}

A triangulation type. D is the dimension of the triangulated space. N = D + 1,
which equals the number of vertices each simplex in the triangulation has.
"""
struct Triangulation{D<:Int, N <:Int} <: AbstractTriangulation{D, N}
    E::AbstractEmbedding
    d::DelaunayTriangulation
    centroids::Array{Float64, 2}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
end

radii(t::AbstractTriangulation) = t.radii[1:end-1]
centroids(t::AbstractTriangulation) = t.centroids
orientations(t::AbstractTriangulation) = t.orientations
volumes(t::AbstractTriangulation) t.volumes

export dimension, npoints, radii, volumes, orientations

"""
    LinearlyInvariantTriangulation <: Triangulation`
A triangulation for which we have made sure the point corresponding to the last time
    index falls within the convex hull of the other points.
"""
struct LinearlyInvariantTriangulation{D <:Int, N <:Int} <: AbstractTriangulation{D, N}
    points::Array{Float64, 2}
    simplex_inds::Array{Int32, 2}
end

function triangulate(E::LinearlyInvariantEmbedding)
    simplex_inds = delaunay(E)
    D = length(simplex_inds[1]) - 1; N = D + 1
    N = D + 1
    points = E.points[1:end-1, :]
    LinearlyInvariantTriangulation{D, N}(points, simplex_inds)
end

function triangulate(E::AbstractEmbedding)
    simplex_inds = delaunayn(E)
    D = length(simplex_inds[1]) - 1; N = D + 1
    points = E.points[1:end-1, :]
    Triangulation{D, N}(points, simplex_inds)
end

triangulate(pts::AbstractArray{Float64, 2}) = triangulate(embed(pts))

function Base.summary(t::AbstractTriangulation)
    npts = size(t.embedding.points, 1)
    nsimplices = size(t.simplex_inds, 1)
    dim = size(t.embedding.points, 2)
    embeddingtype_tri = typeof(t)
    embeddingtype_emb = typeof(t.embedding)
    return """$dim-dimensional $(embeddingtype_tri) with $nsimplices simplices constructed
            from a $npts-pt $embeddingtype_emb"""
end

function matstring(t::AbstractTriangulation)
    fields = fieldnames(t)
    fields_str = String.(fields)
    maxlength = maximum([length(str) for str in fields_str]) + 2
    fields_str = [fields_str[i] *
                repeat(" ", maxlength - length(fields_str[i])) for i = 1:length(fields_str)]

    summaries = [join(":"*String(fields_str[i])*summary(getfield(t, fields[i]))*"\n") for i = 1:length(fields_str)] |> join
    infoline = "The following fields are available:\n"

    return summary(t)#*"\n\n"*infoline*summaries
end

Base.show(io::IO, t::AbstractTriangulation) = println(io, matstring(t))

"""
Find the indices of the simplices in the original triangulation that potentially
intersect with the image simplex with index `image_i`.
"""
function maybeintersecting_simplices(t::AbstractTriangulation, image_i::Int)
    inds_potential_simplices = Int[]

    n_simplices = length(t.radii)

    @inbounds for i = 1:n_simplices
        dist_difference = (transpose(t.centroids_im[image_i] - t.centroids[i]) *
                            (t.centroids_im[image_i] - t.centroids[i]) -
                                (t.radii_im[image_i] + t.radii[i])^2)[1]
        if dist_difference < 0
            push!(inds_potential_simplices, i)
        end
    end
    return inds_potential_simplices
end

"""
Find the indices of the image simplices in `t` that potentially intersect with
the original simplex with index `orig_i`.
"""
function maybeintersecting_imsimplices(t::AbstractTriangulation, orig_i::Int)
    inds_potential_simplices = Int[]

    n_simplices = length(t.radii)

    @inbounds for i = 1:n_simplices
        dist_difference = (transpose(t.centroids[orig_i] - t.centroids_im[i]) *
                            (t.centroids[orig_i] - t.centroids_im[i]) -
                                (t.radii[orig_i] + t.radii_im[i])^2)[1]
        if dist_difference < 0
            push!(inds_potential_simplices, i)
        end
    end
    return inds_potential_simplices
end

"""
Draw point representatives from the simplices of a triangulation `t`.

Precedure:
1) Generate one point per simplex.
2) Points are generated from the interior or on the boundary of each simplex.
3) Points are drawn according to a uniform distribution.
"""
function point_representatives(t::AbstractTriangulation)
    dim = size(t.points, 2)
    n_simplices = size(t.simplex_inds, 1)

    # Pre-allocate array to hold the points
    point_representatives = zeros(Float64, n_simplices, dim)

    # Loop over the rows of the simplex_inds array to access all the simplices.
    for i = 1:n_simplices
        simplex = t.points[t.simplex_inds[i], :]
        point_representatives[i, :] = Delaunay.childpoint(simplex)
    end

    return point_representatives
end


text:
export nca

using mlpack._Internal.io

import mlpack_jll
const ncaLibrary = mlpack_jll.libmlpack_julia_nca

# Call the C binding of the mlpack nca binding.
function nca_mlpackMain()
  success = ccall((:nca, ncaLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module nca_internal
  import ..ncaLibrary

end # module

"""
    nca(input; [armijo_constant, batch_size, labels, linear_scan, max_iterations, max_line_search_trials, max_step, min_step, normalize, num_basis, optimizer, seed, step_size, tolerance, verbose, wolfe])

This program implements Neighborhood Components Analysis, both a linear
dimensionality reduction technique and a distance learning technique.  The
method seeks to improve k-nearest-neighbor classification on a dataset by
scaling the dimensions.  The method is nonparametric, and does not require a
value of k.  It works by using stochastic ("soft") neighbor assignments and
using optimization techniques over the gradient of the accuracy of the neighbor
assignments.

To work, this algorithm needs labeled data.  It can be given as the last row of
the input dataset (specified with `input`), or alternatively as a separate
matrix (specified with `labels`).

This implementation of NCA uses stochastic gradient descent, mini-batch
stochastic gradient descent, or the L_BFGS optimizer.  These optimizers do not
guarantee global convergence for a nonconvex objective function (NCA's objective
function is nonconvex), so the final results could depend on the random seed or
other optimizer parameters.

Stochastic gradient descent, specified by the value 'sgd' for the parameter
`optimizer`, depends primarily on three parameters: the step size (specified
with `step_size`), the batch size (specified with `batch_size`), and the maximum
number of iterations (specified with `max_iterations`).  In addition, a
normalized starting point can be used by specifying the `normalize` parameter,
which is necessary if many warnings of the form 'Denominator of p_i is 0!' are
given.  Tuning the step size can be a tedious affair.  In general, the step size
is too large if the objective is not mostly uniformly decreasing, or if
zero-valued denominator warnings are being issued.  The step size is too small
if the objective is changing very slowly.  Setting the termination condition can
be done easily once a good step size parameter is found; either increase the
maximum iterations to a large number and allow SGD to find a minimum, or set the
maximum iterations to 0 (allowing infinite iterations) and set the tolerance
(specified by `tolerance`) to define the maximum allowed difference between
objectives for SGD to terminate.  Be careful---setting the tolerance instead of
the maximum iterations can take a very long time and may actually never converge
due to the properties of the SGD optimizer. Note that a single iteration of SGD
refers to a single point, so to take a single pass over the dataset, set the
value of the `max_iterations` parameter equal to the number of points in the
dataset.

The L-BFGS optimizer, specified by the value 'lbfgs' for the parameter
`optimizer`, uses a back-tracking line search algorithm to minimize a function.
The following parameters are used by L-BFGS: `num_basis` (specifies the number
of memory points used by L-BFGS), `max_iterations`, `armijo_constant`, `wolfe`,
`tolerance` (the optimization is terminated when the gradient norm is below this
value), `max_line_search_trials`, `min_step`, and `max_step` (which both refer
to the line search routine).  For more details on the L-BFGS optimizer, consult
either the mlpack L-BFGS documentation (in lbfgs.hpp) or the vast set of
published literature on L-BFGS.

By default, the SGD optimizer is used.

# Arguments

 - `input::Array{Float64, 2}`: Input dataset to run NCA on.
 - `armijo_constant::Float64`: Armijo constant for L-BFGS.  Default value
      `0.0001`.

 - `batch_size::Int`: Batch size for mini-batch SGD.  Default value `50`.

 - `labels::Array{Int, 1}`: Labels for input dataset.
 - `linear_scan::Bool`: Don't shuffle the order in which data points are
      visited for SGD or mini-batch SGD.  Default value `false`.

 - `max_iterations::Int`: Maximum number of iterations for SGD or L-BFGS
      (0 indicates no limit).  Default value `500000`.

 - `max_line_search_trials::Int`: Maximum number of line search trials for
      L-BFGS.  Default value `50`.

 - `max_step::Float64`: Maximum step of line search for L-BFGS.  Default
      value `1e+20`.

 - `min_step::Float64`: Minimum step of line search for L-BFGS.  Default
      value `1e-20`.

 - `normalize::Bool`: Use a normalized starting point for optimization.
      This is useful for when points are far apart, or when SGD is returning
      NaN.  Default value `false`.

 - `num_basis::Int`: Number of memory points to be stored for L-BFGS.
      Default value `5`.

 - `optimizer::String`: Optimizer to use; 'sgd' or 'lbfgs'.  Default value
      `sgd`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `step_size::Float64`: Step size for stochastic gradient descent
      (alpha).  Default value `0.01`.

 - `tolerance::Float64`: Maximum tolerance for termination of SGD or
      L-BFGS.  Default value `1e-07`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

 - `wolfe::Float64`: Wolfe condition parameter for L-BFGS.  Default value
      `0.9`.

# Return values

 - `output::Array{Float64, 2}`: Output matrix for learned distance
      matrix.

"""
function nca(input;
             armijo_constant::Union{Float64, Missing} = missing,
             batch_size::Union{Int, Missing} = missing,
             labels = missing,
             linear_scan::Union{Bool, Missing} = missing,
             max_iterations::Union{Int, Missing} = missing,
             max_line_search_trials::Union{Int, Missing} = missing,
             max_step::Union{Float64, Missing} = missing,
             min_step::Union{Float64, Missing} = missing,
             normalize::Union{Bool, Missing} = missing,
             num_basis::Union{Int, Missing} = missing,
             optimizer::Union{String, Missing} = missing,
             seed::Union{Int, Missing} = missing,
             step_size::Union{Float64, Missing} = missing,
             tolerance::Union{Float64, Missing} = missing,
             verbose::Union{Bool, Missing} = missing,
             wolfe::Union{Float64, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, ncaLibrary), Nothing, ());

  IORestoreSettings("Neighborhood Components Analysis (NCA)")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(armijo_constant)
    IOSetParam("armijo_constant", convert(Float64, armijo_constant))
  end
  if !ismissing(batch_size)
    IOSetParam("batch_size", convert(Int, batch_size))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(linear_scan)
    IOSetParam("linear_scan", convert(Bool, linear_scan))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(max_line_search_trials)
    IOSetParam("max_line_search_trials", convert(Int, max_line_search_trials))
  end
  if !ismissing(max_step)
    IOSetParam("max_step", convert(Float64, max_step))
  end
  if !ismissing(min_step)
    IOSetParam("min_step", convert(Float64, min_step))
  end
  if !ismissing(normalize)
    IOSetParam("normalize", convert(Bool, normalize))
  end
  if !ismissing(num_basis)
    IOSetParam("num_basis", convert(Int, num_basis))
  end
  if !ismissing(optimizer)
    IOSetParam("optimizer", convert(String, optimizer))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(step_size)
    IOSetParam("step_size", convert(Float64, step_size))
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if !ismissing(wolfe)
    IOSetParam("wolfe", convert(Float64, wolfe))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  # Call the program.
  nca_mlpackMain()

  return IOGetParamMat("output", points_are_rows)
end

ext: jl


text:
function generate_model_parameters_dictionary()
	# load stoichiometry
	stoichiometric_matrix = readdlm("./stoichiometry.dat")

	# initial species concentration
	initial_condition = [
		0.0,  # 1  BIOMASS
		0.0,  # 2  m_A_e
		0.0,  # 3  m_B_e
		0.0,  # 4  m_C_e
		0.0,  # 5  m_A_c
		0.0,  # 6  m_B_c
		0.0,  # 7  m_C_c
		0.0,  # 8  p_E1_c
		0.0,  # 9  p_E2_c
		0.0,  # 10  p_E3_c
		0.0,  # 11  p_E4_c
		0.0,  # 12  p_TA_c
		0.0,  # 13  p_TB_c
		0.0,  # 14  p_TC_c
		0.0,  # 15  mRNA_E1_c
		0.0,  # 16  mRNA_E2_c
		0.0,  # 17  mRNA_E3_c
		0.0,  # 18  mRNA_E4_c
		0.0,  # 19  mRNA_TA_c
		0.0,  # 20  mRNA_TB_c
		0.0   # 21 mRNA_TC_c
	]

	# all species dictionary
	all_species_dict = Dict{String, Int64}()  # String --> Int
	all_species_dict["m_A_e"] = 2
	all_species_dict["m_C_c"] = 7
	all_species_dict["mRNA_E2_c"] = 16
	all_species_dict["mRNA_E4_c"] = 18
	all_species_dict["BIOMASS"] = 1
	all_species_dict["p_TC_c"] = 14
	all_species_dict["p_E2_c"] = 9
	all_species_dict["mRNA_TB_c"] = 20
	all_species_dict["p_E3_c"] = 10
	all_species_dict["mRNA_TA_c"] = 19
	all_species_dict["p_E1_c"] = 8
	all_species_dict["m_C_e"] = 4
	all_species_dict["m_B_c"] = 6
	all_species_dict["p_TA_c"] = 12
	all_species_dict["p_TB_c"] = 13
	all_species_dict["mRNA_E1_c"] = 15
	all_species_dict["m_A_c"] = 5
	all_species_dict["mRNA_E3_c"] = 17
	all_species_dict["m_B_e"] = 3
	all_species_dict["p_E4_c"] = 11
	all_species_dict["mRNA_TC_c"] = 21

	# all species dictionary_reversed
	all_species_reversed_dict = Dict{Int64, String}()  # Int --> String
	for (key,val) in all_species_dict
		all_species_reversed_dict[val] = key
	end

	# all rnx species array
	rnx_species_array::Array{String,1} = [
		"m_A_e",
		"m_B_e",
		"m_C_e",
		"m_A_c",
		"m_B_c",
		"m_C_c",
		"p_E1_c",
		"p_E2_c",
		"p_E3_c",
		"p_E4_c",
		"p_TA_c",
		"p_TB_c",
		"p_TC_c"
	]

	# all mRNA species array
	mRNA_species_array::Array{String,1} = [
		"mRNA_E4_c",
		"mRNA_E2_c",
		"mRNA_E3_c",
		"mRNA_TA_c",
		"mRNA_TB_c",
		"mRNA_E1_c",
		"mRNA_TC_c"
	]

	# all protein species array
	protein_species_array::Array{String,1} = [
		"p_TA_c",
		"p_TC_c",
		"p_E2_c",
		"p_TB_c",
		"p_E1_c",
		"p_E4_c",
		"p_E3_c"
	]

	# signaling reaction kinetic constants, reaction name: reactant(s)<rnx type:catalyst(s)>product(s)
	kcat_signaling = ones(7)  # kcat[#reaction]: reaction name
	kcat_signaling[1] = 1.1e-3  # kcat: 1.0*m_A_e<uptake:1.0*p_TA_c>1.0*m_A_c
	kcat_signaling[2] = 8e-4  # kcat: 1.0*m_B_c<secrete:1.0*p_TB_c>1.0*m_B_e
	kcat_signaling[3] = 9e-4  # kcat: 1.0*m_C_c<secrete:1.0*p_TC_c>1.0*m_C_e
	kcat_signaling[4] = 1.8e-4  # kcat: 1.0*m_A_c<catalyze:1.0*p_E1_c>1.0*m_B_c
	kcat_signaling[5] = 1e-4  # kcat: 1.0*m_B_c<catalyze:1.0*p_E4_c>1.0*m_A_c
	kcat_signaling[6] = 1e-4  # kcat: 1.0*m_A_c<catalyze:1.0*p_E2_c>1.0*m_C_c
	kcat_signaling[7] = 1e-4  # kcat: 1.0*m_C_c<catalyze:1.0*p_E3_c>1.0*m_B_c

	# Monod affinity constant in signaling reaction, nomenclature: MonodK_#reaction_Substrate
	MonodAffinityConstantDict = Dict{String, Float64}()
	MonodAffinityConstantDict["MonodK~rnx1~m_A_e"] = 0.1
	MonodAffinityConstantDict["MonodK~rnx2~m_B_c"] = 0.1
	MonodAffinityConstantDict["MonodK~rnx3~m_C_c"] = 0.1
	MonodAffinityConstantDict["MonodK~rnx4~m_A_c"] = 0.1
	MonodAffinityConstantDict["MonodK~rnx5~m_B_c"] = 0.1
	MonodAffinityConstantDict["MonodK~rnx6~m_A_c"] = 0.1
	MonodAffinityConstantDict["MonodK~rnx7~m_C_c"] = 0.1

	# W_value in transcription control term, nomenclature: W_targetmRNA_actor
	W_value_dict = Dict{String, Float64}()
	W_value_dict["W~m_B_c~mRNA_E4_c"] = 0.1
	W_value_dict["W~m_A_c~mRNA_E2_c"] = 0.7
	W_value_dict["W~m_C_c~mRNA_E3_c"] = 0.2
	W_value_dict["W~m_A_e~mRNA_TA_c"] = 0.1
	W_value_dict["W~m_B_c~mRNA_TB_c"] = 1.0
	W_value_dict["W~m_A_c~mRNA_E1_c"] = 0.6
	W_value_dict["W~m_C_c~mRNA_TC_c"] = 1.1

	# disassociation constants in transfer function, nomenclature: KD_speciesName
	transferFunctionDisassociationConstantDict = Dict{String, Float64}()
	transferFunctionDisassociationConstantDict["KD~m_B_c"] = 1.0
	transferFunctionDisassociationConstantDict["KD~m_A_c"] = 1.0
	transferFunctionDisassociationConstantDict["KD~m_C_c"] = 1.0
	transferFunctionDisassociationConstantDict["KD~m_A_e"] = 1.0
	transferFunctionDisassociationConstantDict["KD~m_B_c"] = 1.0
	transferFunctionDisassociationConstantDict["KD~m_A_c"] = 1.0
	transferFunctionDisassociationConstantDict["KD~m_C_c"] = 1.0

	# transcription specific correction factor
	transcriptionSpecificCorrectionFactorDict = Dict{String, Float64}()
	transcriptionSpecificCorrectionFactorDict["mRNA_E4_c"] = 1.0  # mRNA_E4_c
	transcriptionSpecificCorrectionFactorDict["mRNA_E2_c"] = 1.0  # mRNA_E2_c
	transcriptionSpecificCorrectionFactorDict["mRNA_E3_c"] = 1.0  # mRNA_E3_c
	transcriptionSpecificCorrectionFactorDict["mRNA_TA_c"] = 1.0  # mRNA_TA_c
	transcriptionSpecificCorrectionFactorDict["mRNA_TB_c"] = 1.0  # mRNA_TB_c
	transcriptionSpecificCorrectionFactorDict["mRNA_E1_c"] = 1.0  # mRNA_E1_c
	transcriptionSpecificCorrectionFactorDict["mRNA_TC_c"] = 1.0  # mRNA_TC_c

	# background gene expression control term
	backgroundGeneExpressionControlTermDict = Dict{String, Float64}()
	backgroundGeneExpressionControlTermDict["mRNA_E4_c"] = 0.001  # mRNA_E4_c
	backgroundGeneExpressionControlTermDict["mRNA_E2_c"] = 0.001  # mRNA_E2_c
	backgroundGeneExpressionControlTermDict["mRNA_E3_c"] = 0.001  # mRNA_E3_c
	backgroundGeneExpressionControlTermDict["mRNA_TA_c"] = 0.001  # mRNA_TA_c
	backgroundGeneExpressionControlTermDict["mRNA_TB_c"] = 0.001  # mRNA_TB_c
	backgroundGeneExpressionControlTermDict["mRNA_E1_c"] = 0.001  # mRNA_E1_c
	backgroundGeneExpressionControlTermDict["mRNA_TC_c"] = 0.001  # mRNA_TC_c

	# translation specific correction factor
	translationSpecificCorrectionFactor = Dict{String, Float64}()
	translationSpecificCorrectionFactor["p_TA_c"] = 1.0  # p_TA_c
	translationSpecificCorrectionFactor["p_TC_c"] = 1.0  # p_TC_c
	translationSpecificCorrectionFactor["p_E2_c"] = 1.0  # p_E2_c
	translationSpecificCorrectionFactor["p_TB_c"] = 1.0  # p_TB_c
	translationSpecificCorrectionFactor["p_E1_c"] = 1.0  # p_E1_c
	translationSpecificCorrectionFactor["p_E4_c"] = 1.0  # p_E4_c
	translationSpecificCorrectionFactor["p_E3_c"] = 1.0  # p_E3_c

	# cooperativity number in transfer function
	cooperativity = 1

	# ------------------------------------------------------------------------------------------#
	# constants (from bionumbers)       units
	# ------------------------------------------------------------------------------------------#
	cell_diameter = 1.1                 # mum
	number_of_rnapII = 4600            	# copies/cells
	number_of_ribosome = 50000         	# copies/cells
	mRNA_half_life_TF = 0.083           # hrs
	protein_half_life = 70              # hrs
	doubling_time_cell = 0.33           # hrs
	max_translation_rate = 16.5         # aa/sec
	max_transcription_rate = 60.0       # nt/sec
	average_transcript_length = 1200   	# nt
	average_protein_length = 400       	# aa
	fraction_nucleus = 0.0             	# dimensionless
	av_number = 6.02e23                 # number/mol
	avg_gene_number = 2                 # number of copies of a gene
	polysome_number = 4									# number of ribsomoses per transcript
	# ------------------------------------------------------------------------------------------#
	#
	# ------------------------------------------------------------------------------------------#
	# Calculate constants using bionumber values
	# ------------------------------------------------------------------------------------------#
	# Calculate the volume (convert to L)
	V = ((1-fraction_nucleus)*(1/6)*(3.14159)*(cell_diameter)^3)*(1e-15)     # volume

	# Calculate the rnapII_concentration and ribosome_concentration
	rnapII_concentration = number_of_rnapII*(1/av_number)*(1/V)*1e9                   # nM
	ribosome_concentration = number_of_ribosome*(1/av_number)*(1/V)*1e9               # nM

	# degrdation rate constants -
	degradation_constant_mRNA = -(1/mRNA_half_life_TF)*log(0.5)                       # hr^-1
	degradation_constant_protein = -(1/protein_half_life)*log(0.5)                    # hr^-1

	# kcats for transcription and translation -
	kcat_transcription = max_transcription_rate*(3600/average_transcript_length)      # hr^-1
	kcat_translation = polysome_number*max_translation_rate*(3600/average_protein_length)             # hr^-1

	# Maximum specific growth rate -
	maximum_specific_growth_rate = (1/doubling_time_cell)*log(2)                      # hr^-1

	# What is the average gene concentration -
	avg_gene_concentration = avg_gene_number*(1/av_number)*(1/V)*1e9                  # nM

	# How fast do my cells die?
	death_rate_constant = 0.05*maximum_specific_growth_rate                            # hr^-1

	# Saturation constants for translation and trascription -
	saturation_transcription = 4600*(1/av_number)*(1/V)*1e9                           # nM
	saturation_translation = 150000*(1/av_number)*(1/V)*1e9                           # nM
	# -------------------------------------------------------------------------------------------#

	####################################
	# put all stuff in a Dictionary
	dataDictionary = Dict{String, Any}()
	dataDictionary["stoichiometric_matrix"] = stoichiometric_matrix
	dataDictionary["all_species_dict"] = all_species_dict
	dataDictionary["all_species_reversed_dict"] = all_species_reversed_dict
	dataDictionary["rnx_species_array"] = rnx_species_array
	dataDictionary["mRNA_species_array"] = mRNA_species_array
	dataDictionary["protein_species_array"] = protein_species_array
	dataDictionary["initial_condition"] = initial_condition
	dataDictionary["kcat_signaling"] = kcat_signaling
	dataDictionary["MonodAffinityConstantDict"] = MonodAffinityConstantDict
	dataDictionary["W_value_dict"] = W_value_dict
	dataDictionary["transferFunctionDisassociationConstantDict"] = transferFunctionDisassociationConstantDict
	dataDictionary["transcriptionSpecificCorrectionFactorDict"] = transcriptionSpecificCorrectionFactorDict
	dataDictionary["backgroundGeneExpressionControlTermDict"] = backgroundGeneExpressionControlTermDict
	dataDictionary["translationSpecificCorrectionFactor"] = translationSpecificCorrectionFactor
	dataDictionary["cooperativity"] = cooperativity
	dataDictionary["RNAPConcentration"] = rnapII_concentration
	dataDictionary["RIBOConcentration"] = ribosome_concentration
	dataDictionary["degradationConstantmRNA"] = degradation_constant_mRNA
	dataDictionary["degradationConstantProtein"] = degradation_constant_protein
	dataDictionary["kcatTranscription"] = kcat_transcription
	dataDictionary["kcatTranslation"] = kcat_translation
	dataDictionary["avgGeneConcentration"] = avg_gene_concentration
	dataDictionary["transcriptionSaturationConstant"] = saturation_transcription
	dataDictionary["translationSaturationConstant"] = saturation_translation
	dataDictionary["specificGrowthRate"] = maximum_specific_growth_rate - death_rate_constant

	return dataDictionary
end


text:
export unify
export UNIFICATION_FAILURE_LOGGING_LEVEL, logging_unification_failures
export unification_failure

"""
    UNIFICATION_FAILURE_LOGGING_LEVEL
See `@unification_failure`.
"""
UNIFICATION_FAILURE_LOGGING_LEVEL = false

function logging_unification_failures(f, level)
    global UNIFICATION_FAILURE_LOGGING_LEVEL
    was = UNIFICATION_FAILURE_LOGGING_LEVEL
    try
        UNIFICATION_FAILURE_LOGGING_LEVEL = level
        f()
    finally
        UNIFICATION_FAILURE_LOGGING_LEVEL = was
    end
end

"""
   @unification_failure(thing1, thing2)
Logs that `unify` failed for `thing` and `thing2`.
If the value of `UNIFICATION_FAILURE_LOGGING_LEVEL` is a
Logging.LogLevel` then it is the log level for this message,
otherwise no log entry is made.
"""
macro unification_failure(thing1, thing2, more...)
    args = [ :_file => __source__.file,
             :_line => __source__.line,
             :_module => __module__,
             more...]
    :(
        if UNIFICATION_FAILURE_LOGGING_LEVEL isa LogLevel
            #=
            # Using @logmsg gets the line number wrong.  I want the line
            # where @unification_failure appears, not the line in the
            # macro that logs the message.
            Logging.handle_message(
                Base.CoreLogging.current_logger_for_env(
                    UNIFICATION_FAILURE_LOGGING_LEVEL, :none, $m),
                UNIFICATION_FAILURE_LOGGING_LEVEL,
                "Unification failure",
                $m, :none, gensym("id"), $(source.line), "$(source.file)";
                thing1=$(esc(thing1)), thing2=$(esc(thing2)),
                more=$(esc(more)))
            =#
            @logmsg(UNIFICATION_FAILURE_LOGGING_LEVEL,
                    "Unification failure",
                    thing1=$(esc(thing1)), thing2=$(esc(thing2)),
                    $args...)
        end)
end

"""
    unify(continuation, expression1, expression2, bindings=EmptyBindings())

Unify the two expressions, calling `continuation` with each set of
bindings which satisfy the unification.
"""
function unify end

function unify(continuation, thing1, thing2)
    unify(continuation, thing1, thing2, EmptyBindings())
end

### Ignore

function unify(continuation, ::Ignore, ::Any, bindings::AbstractBindings)
  continuation(bindings)
end

function unify(continuation, ::Any, ::Ignore, bindings::AbstractBindings)
  continuation(bindings)
end

### Var

function unify(continuation, var::Var, other::Any, bindings::AbstractBindings)
    unify_var(continuation, var, other, bindings)
end

function unify(continuation, other::Any, var::Var, bindings::AbstractBindings)
    unify_var(continuation, var, other, bindings)
end

function unify(continuation, var1::Var, var2::Var, bindings::AbstractBindings)
    ##### Is this enough
    unify_var(continuation, var1, var2, bindings)
end

"""
`unify` with a `NoCirc` always succeeds unless the NoCirc is being unified
against its subject `var`.
"""
struct NoCirc <: AbstractVar
    var::AbstractVar
end

function unify(continuation, nc::NoCirc, ::Any, bindings::AbstractBindings)
    continuation(bindings)
end

function unify(continuation, nc::NoCirc, var::Var, bindings::AbstractBindings)
    if same(nc.var, var)
        return continuation(bindings)
    end
    continuation(bindings)
end

function unify(continuation, ::Any, nc::NoCirc, bindings)
    throw(ErrorException("$nc appeared in right hand side of unify."))
end

function unify_var(continuation, var::Var, other::Any, bindings::AbstractBindings)
    values, _ = lookupequiv(bindings, var)
    if length(values) > 1
        # Variable's value is already ambiguous.
        @unification_failure(var, other, values)
        return
    end
    if length(values) == 1
        return unify(continuation, first(values), other, bindings)
    end
    # To avoid a reference cycle, we must make sure `var` does not
    # appear anywhere in `other`.  It would be nice if there were only
    # one generic function that needed to code how to walk non-atomic
    # data.  We use `unify` for that.  We unify `other` against a NoCirc
    # that matches anything but the variable we are looking for.
    #=
    unify(NoCirc(var), other, bindings) do bindings
        ubind(continuation, var, other, bindings)
    end
    =#
    ubind(continuation, var, other, bindings)
end

### SubseqVar

### default

##### Is this method shadowed by the one on T, T where T?
function unify(continuation, thing1::Any, thing2::Any, bindings::AbstractBindings)
    if thing1 == thing2
        return continuation(bindings)
    end
    @unification_failure(thing1, thing2)
end

### Things that need to be == to unify

"""
    @unify_equal(typ, op)
Define a method on `unify` for type `typ` that only succeeds if the
two instances of `typ` satisfy `op`.
"""
macro unify_equal(typ, op)
    :(function unify(continuation, thing1::$typ, thing2::$typ, bindings::AbstractBindings)
          if $op(thing1, thing2)
              return continuation(bindings)
          end
          @unification_failure(thing1, thing2)
      end)
end

# Numbers of different types that might be equal:
@unify_equal(Number, ==)

@unify_equal(Symbol, ==)
@unify_equal(AbstractChar, ==)
@unify_equal(AbstractString, ==)

### Things that are finitely indexable:

struct UnifyIndexableElement
    thing
    index::Int

    UnifyIndexableElement(thing) = new(thing, firstindex(thing))
    UnifyIndexableElement(thing, index::Int) = new(thing, index)
end

function next(e::UnifyIndexableElement)
    UnifyIndexableElement(e.thing, e.index + 1)
end

function exhausted(e::UnifyIndexableElement)
    e.index > lastindex(e.thing)
end

function elt(e::UnifyIndexableElement)
    e.thing[e.index]
end

function unify_indexable(continuation, thing1, thing2,
                         bindings::AbstractBindings)
    unify_indexable(continuation,
                    UnifyIndexableElement(thing1),
                    UnifyIndexableElement(thing2),
                    bindings)
end

function unify_indexable(continuation,
                         e1::UnifyIndexableElement, e2::UnifyIndexableElement,
                         bindings::AbstractBindings)
    if exhausted(e1) && exhausted(e2)
        return continuation(bindings)
    end
    if exhausted(e1) || exhausted(e2)
        @unification_failure(e1, e2)
        return
    end
    if isa(elt(e1), SubseqVar) && isa(elt(e2), SubseqVar)
        ubind(elt(e1), elt(e2), bindings) do bindings
            unify_indexable(continuation, next(e1), next(e2), bindings)
        end
        return
    end
    if isa(elt(e1), SubseqVar)
        return unify_indexable(continuation, elt(e1), e1, e2, bindings)
    end
    if isa(elt(e2), SubseqVar)
        return unify_indexable(continuation, elt(e2), e2, e1, bindings)
    end
    unify(elt(e1), elt(e2), bindings) do bindings
        return unify_indexable(continuation, next(e1), next(e2), bindings)
    end
end

function unify_indexable(continuation,
                         var::SubseqVar, varUIE::UnifyIndexableElement,
                         e::UnifyIndexableElement,
                         bindings::AbstractBindings)
    # Try each length of subsequence of thing2
    for endindex in (e.index - 1):(lastindex(e.thing))
        ubind(var, view(e.thing, (e.index):endindex), bindings) do bindings
            unify_indexable(continuation,
                            next(varUIE),
                            UnifyIndexableElement(e.thing, endindex + 1),
                            bindings)
        end
    end
end

macro generate_unify_indexable_methods(types...)
    defs = []
    for t1 = types, t2 = types
        push!(defs,
              esc(:(function unify(continuation, thing1::$t1, thing2::$t2, bindings::AbstractBindings)
                        unify_indexable(continuation, thing1, thing2, bindings)
                    end)))
    end
    return :(begin $(defs...) end)
end

##### Maybe these types can be computed from methods(IndexStyle) that return IndexLinear.
@generate_unify_indexable_methods(AbstractVector, Tuple)

### Types with fields

# Hacky predicate dispatch mechanism.
abstract type UnificationStrategy end

function unify(continuation, thing1::T, thing2::T, bindings::AbstractBindings) where T
    # Short circuit for objects that are equal and structs with no fields
    if thing1 == thing2
        return continuation(bindings)
    end
    # Try each unification strategy
    for strategy in subtypes(UnificationStrategy)
        unify(continuation, strategy(),
              thing1, thing2, bindings::AbstractBindings)
    end
end

# I don't know how to tell if an object is implemented as a struct
# except to see if it has any fields.  Of course, a struct could have
# no fields, which would make it a singleton type.
struct UnifyFields <: UnificationStrategy end

function unify(continuation, strategy::UnifyFields, thing1, thing2, bindings::AbstractBindings)
    fields = fieldnames(typeof(thing1))
    if length(fields) == 0
        @unification_failure(thing1, thing2)
        return
    end
    function unify_fields(index, bindings)
        if index > length(fields)
            return continuation(bindings)
        end
        field = fields[index]
        unify(getfield(thing1, field),
              getfield(thing2, field),
              bindings) do bindings
                  return unify_fields(index + 1, bindings)
              end
    end
    unify_fields(1, bindings)
end


text:
@testset "random_fourier_features" begin
    rng = Random.MersenneTwister(54321)

    input_dims = 2
    x = rand(rng, 3, input_dims)
    x_rowvecs = RowVecs(x)
    x_colvecs = ColVecs(x')

    y = rand(rng, 4, input_dims)
    y_rowvecs = RowVecs(y)
    y_colvecs = ColVecs(y')

    num_feats = Int(1e7)

    function test_kernel_approx(rng, kernel, x, y, input_dims, num_feats)
        ϕ = sample_rff_basis(rng, kernel, input_dims, num_feats)
        km_true = kernelmatrix(kernel, x, y)
        km_approx = kernelmatrix(ϕ, x, y)
        @test km_true ≈ km_approx rtol = 1e-3 atol = 1e-3
    end

    @testset "SqExponentialKernel" begin
        rng = Random.MersenneTwister(54321)
        se_kernel = SqExponentialKernel()
        test_kernel_approx(rng, se_kernel, x_colvecs, y_colvecs, input_dims, num_feats)
    end
    @testset "ScaledKernel" begin
        rng = Random.MersenneTwister(54321)
        scaled_kernel = 2 * SqExponentialKernel()
        test_kernel_approx(rng, scaled_kernel, x_colvecs, y_colvecs, input_dims, num_feats)
    end
    @testset "ScaleTransform" begin
        rng = Random.MersenneTwister(54321)
        st_kernel = SqExponentialKernel() ∘ ScaleTransform(4.0)
        test_kernel_approx(rng, st_kernel, x_colvecs, y_colvecs, input_dims, num_feats)
    end
    @testset "ColVecs/RowVecs" begin
        rng = Random.MersenneTwister(54321)
        k = 3 * (SqExponentialKernel() ∘ ScaleTransform(3.0))
        ϕ = sample_rff_basis(rng, k, input_dims, 10)
        @test ϕ(x_colvecs) isa ColVecs
        @test ϕ(x_rowvecs) isa RowVecs
        @test ϕ(x_colvecs).X == ϕ(x_rowvecs).X'
        @test kernelmatrix(ϕ, x_colvecs) == kernelmatrix(ϕ, x_rowvecs)
        @test kernelmatrix(ϕ, x_colvecs, y_colvecs) == kernelmatrix(ϕ, x_rowvecs, y_colvecs)
        @test kernelmatrix(ϕ, x_colvecs, y_colvecs) == kernelmatrix(ϕ, x_colvecs, y_rowvecs)
    end
    @testset "rff_weight_space_approx" begin
        rng = Random.MersenneTwister(54321)
        k = 2 * (SqExponentialKernel() ∘ ScaleTransform(5.0))
        f = GP(k)
        rff_wsa = build_rff_weight_space_approx(rng, input_dims, num_feats)
        bfr = rff_wsa(f)
        @test cov(bfr(x_colvecs)) ≈ kernelmatrix(k, x_colvecs) rtol = 1e-3 atol = 1e-3

        f2 = GP(2.0, k)
        @test_throws ErrorException rff_wsa(f2)
    end
end


text:
module RetryRequest

import ..HTTP
import ..Layer, ..request
using ..Sockets
using ..IOExtras
using ..MessageRequest
using ..Messages
import ..@debug, ..DEBUG_LEVEL, ..sprintcompact

"""
    request(RetryLayer, ::URI, ::Request, body) -> HTTP.Response

Retry the request if it throws a recoverable exception.

`Base.retry` and `Base.ExponentialBackOff` implement a randomised exponentially
increasing delay is introduced between attempts to avoid exacerbating network
congestion.

Methods of `isrecoverable(e)` define which exception types lead to a retry.
e.g. `HTTP.IOError`, `Sockets.DNSError`, `Base.EOFError` and `HTTP.StatusError`
(if status is ``5xx`).
"""
abstract type RetryLayer{Next <: Layer} <: Layer end
export RetryLayer

function request(::Type{RetryLayer{Next}}, url, req, body;
                 retries::Int=4, retry_non_idempotent::Bool=false,
                 kw...) where Next

    retry_request = Base.retry(request,
        delays=ExponentialBackOff(n = retries),
        check=(s,ex)->begin
            retry = isrecoverable(ex, req, retry_non_idempotent)
            if retry
                @debug 1 "🔄  Retry $ex: $(sprintcompact(req))"
                reset!(req.response)
            else
                @debug 1 "🚷  No Retry: $(no_retry_reason(ex, req))"
            end
            return s, retry
        end)

    retry_request(Next, url, req, body; kw...)
end

isrecoverable(e) = false
isrecoverable(e::IOError) = true
isrecoverable(e::Sockets.DNSError) = true
isrecoverable(e::HTTP.StatusError) = e.status == 403 || # Forbidden
                                     e.status == 408 || # Timeout
                                     e.status >= 500    # Server Error

isrecoverable(e, req, retry_non_idempotent) =
    isrecoverable(e) &&
    !(req.body === body_was_streamed) &&
    !(req.response.body === body_was_streamed) &&
    (retry_non_idempotent || req.txcount == 0 || isidempotent(req))
    # "MUST NOT automatically retry a request with a non-idempotent method"
    # https://tools.ietf.org/html/rfc7230#section-6.3.1

function no_retry_reason(ex, req)
    buf = IOBuffer()
    show(IOContext(buf, :compact => true), req)
    print(buf, ", ",
        ex isa HTTP.StatusError ? "HTTP $(ex.status): " :
        !isrecoverable(ex) ?  "$ex not recoverable, " : "",
        (req.body === body_was_streamed) ? "request streamed, " : "",
        (req.response.body === body_was_streamed) ? "response streamed, " : "",
        !isidempotent(req) ? "$(req.method) non-idempotent" : "")
    return String(take!(buf))
end

end # module RetryRequest


text:
using GLFW, GLAbstraction, Reactive, ModernGL, Color,Compat, ImmutableArrays
include("../src/color.jl")
include("../src/types.jl")
include("../src/visualize/text.jl")
include("../src/share/text.jl")
include("../src/visualize/image.jl")
include("../src/visualize_interface.jl")

immutable Screen
    id::Symbol
    area
    parent::Screen
    children::Vector{Screen}
    inputs::Dict{Symbol, Any}
    renderlist::Vector{RenderObject}

    hidden::Signal{Bool}
    hasfocus::Signal{Bool}

    perspectivecam::PerspectiveCamera
    orthographiccam::OrthographicCamera
    nativewindow::GLFW.Window
    counter = 1
    function Screen(
    	area,
	    parent::Screen,
	    children::Vector{Screen},
	    inputs::Dict{Symbol, Any},
	    renderlist::Vector{RenderObject},

	    hidden::Signal{Bool},
	    hasfocus::Signal{Bool},

	    perspectivecam::PerspectiveCamera,
	    orthographiccam::OrthographicCamera,
	    nativewindow::GLFW.Window)
        new(symbol("display"*string(counter+=1)), area, parent, children, inputs, renderlist, hidden, hasfocus, perspectivecam, orthographiccam, nativewindow)
    end

    function Screen(
        area,
        children::Vector{Screen},
        inputs::Dict{Symbol, Any},
        renderlist::Vector{RenderObject},

        hidden::Signal{Bool},
        hasfocus::Signal{Bool},

        perspectivecam::PerspectiveCamera,
        orthographiccam::OrthographicCamera,
        nativewindow::GLFW.Window)
        parent = new()
        new(symbol("display"*string(counter+=1)), area, parent, children, inputs, renderlist, hidden, hasfocus, perspectivecam, orthographiccam, nativewindow)
    end
end

function Screen(
        parent::Screen;
        area 				      = parent.area,
        children::Vector{Screen}  = Screen[],
        inputs::Dict{Symbol, Any} = copy(parent.inputs),
        renderlist::Vector{RenderObject} = RenderObject[],

        hidden::Signal{Bool}   	  = parent.hidden,
        hasfocus::Signal{Bool}    = parent.hasfocus,

        nativewindow::GLFW.Window 	   = parent.nativewindow)

	insidescreen = lift(inputs[:mouseposition]) do mpos
		isinside(area.value, mpos...) && !any(children) do screen
			isinside(screen.area.value, mpos...)
		end
	end

	camera_input = merge(inputs, @compat(Dict(
		:mouseposition 	=> keepwhen(insidescreen, Vector2(0.0), inputs[:mouseposition]),
		:scroll_x 		=> keepwhen(insidescreen, 0, inputs[:scroll_x]),
		:scroll_y 		=> keepwhen(insidescreen, 0, inputs[:scroll_y]),
		:window_size 	=> lift(x->Vector4(x.x, x.y, x.w, x.h), area)
	)))
	ocamera      = OrthographicPixelCamera(camera_input)
	pcamera  	 = PerspectiveCamera(camera_input, Vec3(2), Vec3(0))
    screen = Screen(area, parent, children, inputs, renderlist, hidden, hasfocus, pcamera, ocamera, nativewindow)
	push!(parent.children, screen)
	screen
end
function GLAbstraction.isinside(x::Screen, position::Vector2)
	!any(screen->inside(screen.area.value, position...), x.children) && inside(x.area, position...)
end

function Screen(obj::RenderObject, parent::Screen)

	area 	 = boundingbox2D(obj)
	hidden   = Input(false)
	screen 	 = Screen(parent)
	mouse 	 = filter(inside, Input(Screen), parent.inputs[:mouseposition])

	hasfocus = lift(parent.inputs[:mouseposition], parent.inputs[:mousebuttonpressed], screen.area) do pos, buttons, area
		isinside(area, pos...) && !isempty(bottons)
	end
	buttons  = menubar(screen)
	push!(parent.children, screen)
	push!(screen.renderlist, buttons)
	push!(screen.renderlist, obj)
end

function Screen(style::Style{:Default}, parent=first(SCREEN_STACK))

	hidden   	= Input(true)
	screen 	 	= Screen(parent)
	mouse 	 	= filter(Input(Screen), parent.inputs[:mouseposition]) do screen, mpos
	end
	inputs 		= merge(parent.inputs, @compat(Dict(:mouseposition=>mouse)))
	opxcamera   = OrthographicPixelCamera(inputs)
	pcamera  	= PerspectiveCamera(inputs)
	hasfocus 	= lift(parent.inputs[:mouseposition], parent.inputs[:mousebuttonpressed], screen.area) do pos, buttons, area
		isinside(area, pos...) && !isempty(bottons)
	end
	screen 		= Screen(area, parent, children=Screen[], inputs, renderList, hidden, hasfocus, perspectivecam, orthographiccam)
	buttons     = menubar(screen, style)

	push!(parent.children, screen)
	push!(screen.renderlist, buttons)

end

function GLAbstraction.render(x::Screen)
    glViewport(x.area.value)
    render(x.renderlist)
    render(x.children)
end

GLFW.Init()
GLFW.WindowHint(GLFW.SAMPLES, 4)
w = 1200
h = 1000
window = GLFW.CreateWindow(w, h, "name")
GLFW.MakeContextCurrent(window)
GLFW.ShowWindow(window)
framebuffers 		= Input(Vector2{Int}(w, h))
window_size 		= Input(Vector4{Int}(0, 0, w, h))
inputs =  Dict(
	:insidewindow 					=> Input(false),
	:open 							=> Input(true),
	:hasfocus						=> Input(false),

	:window_size					=> window_size,
	:framebuffer_size 				=> framebuffers,
	:windowposition					=> Input(Vector2(0)),

	:unicodeinput					=> Input(Char[]),

	:buttonspressed					=> Input(IntSet()),
	:buttondown						=> Input(0),
	:buttonreleased					=> Input(0),

	:mousebuttonspressed			=> Input(IntSet()),
	:mousedown						=> Input(0),
	:mousereleased					=> Input(0),

	:mouseposition					=> Input(Vector2(0.0)),
	:mouseposition_glfw_coordinates	=> Input(Vector2(0.0)),

	:scroll_x						=> Input(0),
	:scroll_y						=> Input(0)
)
pcamera  = PerspectiveCamera(inputs, Vec3(2), Vec3(0))
pocamera = OrthographicPixelCamera(inputs)

ROOT_SCREEN = Screen(Input(Rectangle(0,0,1200,1000)), Screen[], inputs, RenderObject[], Input(false), inputs[:hasfocus], pcamera, pocamera, window)

function clear!(x::Vector{RenderObject})
    while !isempty(x)
        value = pop!(x)
        delete!(value)
    end
end
global gaggaaa = Texture("pic.jpg")

visualize_area = lift(ROOT_SCREEN.area) do x
Rectangle(div(x.w,7)*3, 0, div(x.w, 7)*3, x.h)
end
search_area = lift(visualize_area) do x
Rectangle(x.x, x.y, x.w, div(x.h,10))
end
visualize_screen = Screen(ROOT_SCREEN, area=visualize_area)

search_screen 	 = Screen(visualize_screen,  area=search_area)
w_height_search  = lift(search_screen.area) do x
    x.h
end
search_offset = lift(w_height_search) do x
    translationmatrix(Vec3(30,x-30,0))
end

search = visualize(Float32[1f0 for i=1:10, j=1:10], screen=search_screen)
push!(search_screen.renderlist, search)

@async for i = 1:100
	if i == 100
		value = eval(Main, :gaggaaa)
	    clear!(visualize_screen.renderlist)
	    obj     = visualize(value, screen=visualize_screen)
	    push!(visualize_screen.renderlist, obj)
	end
	sleep(0.01)
end

glClearColor(0,0,0,0)
function renderloop(screen)

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

	yield()
  	render(screen)

	GLFW.SwapBuffers(screen.nativewindow)
	GLFW.PollEvents()
end
while Romeo.ROOT_SCREEN.inputs[:open].value
    renderloop(Romeo.ROOT_SCREEN)
    sleep(0.0001)
end
GLFW.Terminate()


text:
module IsAvailable

@doc """ A function that gets the status of a package at the link https://juliapackages.com/p/ """ ->
function isavailable(name::String)
    url = "https://juliapackages.com/p/"*lowercase(name)
    try
        r = HTTP.get(url)
        if r.status == 200
            return true
        end
    catch
        return false
    end
end

@doc """ A function that gets the status of several packages in a list""" ->
function app(names::Array{String, 1})
    for name in names
        if isavailable(name)
            print(name*" is available. Do you wanna to install it? [type y(yes) or n(no)] ")
            answer = readline()
            validAnswers = ["yes", "no", "y", "n"]
            while lowercase(answer) ∉ validAnswers
                print("Please, type some valid option: [y] yes and [n] no. ")
                answer = readline()
            end
            if lowercase(answer) ∈ ["yes", "y"]
                try
                    Pkg.add(name)
                catch
                    print("Oops! I did it again. Check out the name of the package... It may be wrong")
                end
            elseif lowercase(answer) ∈ ["no", "n"]
                print("Okay... Bye!")
            end
        end
    end
end

end # module

max_stars_repo_name: Leticia-maria/isavailable

max_issues_repo_name: Leticia-maria/IsAvailable


text:
@testset "inchi" begin
    demomol = joinpath(dirname(@__FILE__), "..", "assets", "test", "demo.mol")
    txt = read(open(demomol), String)
    inchi1 = inchi(txt)
    ikey1 = inchikey(inchi1)
    demomol = joinpath(dirname(@__FILE__), "..", "assets", "test", "demo.mol")
    mol = sdftomol(demomol)
    inchi2 = inchi(mol)
    ikey2 = inchikey(inchi2)
    @test inchi1 == inchi2
    @test ikey1 == ikey2
end


text:
function get_derivatives(u::Vector{Float64},p::Vector{Float64})
    # derivatives: dF/d[bifurcation_param]
    dFdp::Vector{Float64} = zeros(V.NUM)

    dFdp[V.pERK] = p[C.kf]*(p[C.tERK] - u[V.pERK] - u[V.pERK_DUSP])

    return dFdp
end


text:
# Based partly on code in LegacyStrings that used to be part of Julia
# Licensed under MIT License, see LICENSE.md

# (Mostly written by Scott P. Jones in series of PRs contributed to the Julia project in 2015)

## Error messages for Unicode / UTF support

const UTF_ERR_SHORT =
  "invalid UTF-8 sequence starting at index <<1>> (0x<<2>>) missing one or more continuation bytes"
const UTF_ERR_CONT =
  "invalid UTF-8 sequence starting at index <<1>> (0x<<2>>) is not a continuation byte"
const UTF_ERR_LONG =
  "invalid UTF-8 sequence, overlong encoding starting at index <<1>> (0x<<2>>)"
const UTF_ERR_NOT_LEAD =
  "not a leading Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const UTF_ERR_NOT_TRAIL =
  "not a trailing Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const UTF_ERR_NOT_SURROGATE =
  "not a valid Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const UTF_ERR_MISSING_SURROGATE =
  "missing trailing Unicode surrogate code unit after index <<1>> (0x<<2>>)"
const UTF_ERR_INVALID =
  "invalid Unicode character starting at index <<1>> (0x<<2>> > 0x10ffff)"
const UTF_ERR_SURROGATE =
  "surrogate encoding not allowed in UTF-8 or UTF-32, at index <<1>> (0x<<2>>)"
const UTF_ERR_ODD_BYTES_16 =
  "UTF16String can't have odd number of bytes <<1>>"
const UTF_ERR_ODD_BYTES_32 =
  "UTF32String must have multiple of 4 bytes <<1>>"
const UTF_ERR_INVALID_ASCII =
  "invalid ASCII character at index <<1>> (0x<<2>> > 0x7f)"
const UTF_ERR_INVALID_LATIN1 =
  "invalid Latin1 character at index <<1>> (0x<<2>> > 0xff)"
const UTF_ERR_INVALID_CHAR =
  "invalid Unicode character (0x<<2>> > 0x10ffff)"
const UTF_ERR_INVALID_8 =
  "invalid UTF-8 data"
const UTF_ERR_INVALID_16 =
  "invalid UTF-16 data"
const UTF_ERR_INVALID_UCS2 =
  "invalid UCS-2 character (surrogate present)"
const UTF_ERR_INVALID_INDEX =
  "invalid character index <<1>>"

@static if isdefined(Base, :UnicodeError)
    Base.UnicodeError(msg) = UnicodeError(msg, 0%Int32, 0%UInt32)
else

struct UnicodeError <: Exception
    errmsg::AbstractString   ##< A UTF_ERR_ message
    errpos::Int32            ##< Position of invalid character
    errchr::UInt32           ##< Invalid character
    UnicodeError(msg, pos, chr) = new(msg, pos%Int32, chr%UInt32)
    UnicodeError(msg) = new(msg, 0%Int32, 0%UInt32)
end

_repmsg(msg, pos, chr) =
    replace(replace(msg, "<<1>>" => string(pos)), "<<2>>" =>  outhex(chr))
show(io::IO, exc::UnicodeError) =
    print(io, "UnicodeError: ", _repmsg(exc.errmsg, exc.errpos, exc.errchr))
end

const UTF_ERR_DECOMPOSE_COMPOSE = "only one of decompose or compose may be true"
const UTF_ERR_COMPAT_STRIPMARK  = "compat or stripmark true requires compose or decompose true"
const UTF_ERR_NL_CONVERSION     = "only one newline conversion may be specified"
const UTF_ERR_NORMALIZE         = " is not one of :NFC, :NFD, :NFKC, :NFKD"

@noinline boundserr(s, pos)      = throw(BoundsError(s, pos))
@noinline unierror(err)          = throw(UnicodeError(err))
@noinline unierror(err, pos, ch) = throw(UnicodeError(err, pos, ch))
@noinline unierror(err, v)       = unierror(string(":", v, err))
@noinline nulerr()               = unierror("cannot convert NULL to string")
@noinline neginderr(s, n)        = unierror("Index ($n) must be non negative")
@noinline codepoint_error(T, v)  = unierror(string("Invalid CodePoint: ", T, " 0x", outhex(v)))
@noinline argerror(startpos, endpos) =
    unierror(string("End position ", endpos, " is less than start position (", startpos, ")"))

@noinline ascii_err()    = throw(ArgumentError("Not a valid ASCII string"))
@noinline ncharerr(n)    = throw(ArgumentError(string("nchar (", n, ") must be not be negative")))
@noinline repeaterr(cnt) = throw(ArgumentError("repeat count $cnt must be >= 0"))

@static isdefined(Base, :string_index_err) && (const index_error = Base.string_index_err)

## Functions to check validity of UTF-8, UTF-16, and UTF-32 encoded strings,
#  and also to return information necessary to convert to other encodings

is_surrogate_lead(c::Unsigned) = ((c & ~0x003ff) == 0xd800)
is_surrogate_trail(c::Unsigned) = ((c & ~0x003ff) == 0xdc00)
is_surrogate_codeunit(c::Unsigned) = ((c & ~0x007ff) == 0xd800)
is_valid_continuation(c) = ((c & 0xc0) == 0x80)

## Return flags for check_string function

const UTF_LONG      =  1  ##< Long encodings are present
const UTF_LATIN1    =  2  ##< characters in range 0x80-0xFF present
const UTF_UNICODE2  =  4  ##< characters in range 0x100-0x7ff present
const UTF_UNICODE3  =  8  ##< characters in range 0x800-0xd7ff, 0xe000-0xffff
const UTF_UNICODE4  = 16  ##< non-BMP characters present
const UTF_SURROGATE = 32  ##< surrogate pairs present
const UTF_INVALID   = 64  ##< invalid sequences present

# Get a UTF-8 continuation byte, give error if invalid, return updated character value
@propagate_inbounds function check_continuation(dat, pos, ch, flag)
    @inbounds byt = get_codeunit(dat, pos)
    pos += 1
    if is_valid_continuation(byt)
        flag = false
    elseif !flag
        unierror(UTF_ERR_CONT, pos, byt)
    end
    (ch%UInt32 << 6) | (byt & 0x3f), pos, flag
end

"""
Validates and calculates number of characters in a UTF-8,UTF-16 or UTF-32 encoded vector/string

Warning: this function does not check the bounds of the start or end positions
Use `check_string` to make sure the bounds are checked

Input Arguments:

* `dat`    UTF-8 (`Vector{UInt8}`), UTF-16 (`Vector{UInt16}`) or UTF-32 (`Vector{UInt32}`, `AbstractString`) encoded string
* `pos`    start position
* `endpos` end position

Keyword Arguments:

* `accept_long_null`  = `false`  # Modified UTF-8 (`\\0` represented as `b\"\\xc0\\x80\"`)
* `accept_surrogates` = `false`  # `CESU-8`
* `accept_long_char`  = `false`  # Accept arbitrary long encodings
* `accept_invalids`   = `false`  # Accept invalid sequences (to be replaced on conversion)

Returns:

* (total characters, flags, 4-byte, 3-byte, 2-byte)

Throws:

* `UnicodeError`
"""
function unsafe_check_string end

_ret_check(totalchar, flags, invalids, latin1byte, num2byte, num3byte, num4byte) =
    (totalchar,
     ifelse(latin1byte == 0, 0, UTF_LATIN1) |
     ifelse(num2byte   == 0, 0, UTF_UNICODE2) |
     ifelse(num3byte   == 0, 0, UTF_UNICODE3) |
     ifelse(num4byte   == 0, 0, UTF_UNICODE4) |
     ifelse(invalids   == 0, 0, UTF_INVALID) | flags,
     num4byte, num3byte, num2byte, latin1byte, invalids)

function unsafe_check_string(dat::T, pos, endpos;
                             accept_long_null  = false,
                             accept_surrogates = false,
                             accept_long_char  = false,
                             accept_invalids   = false
                             ) where {T<:Union{AbstractVector{UInt8}, Ptr{UInt8}, String}}
    flags = 0%UInt
    totalchar = latin1byte = num2byte = num3byte = num4byte = invalids = 0
    @inbounds while pos <= endpos
        ch = get_codeunit(dat, pos)
        pos += 1
        totalchar += 1
        if ch > 0x7f
            # Check UTF-8 encoding
            if ch < 0xe0
                # 2-byte UTF-8 sequence (i.e. characters 0x80-0x7ff)
                if pos > endpos
                    accept_invalids || unierror(UTF_ERR_SHORT, pos, ch)
                    invalids += 1
                    break
                end
                ch, pos, flg = check_continuation(dat, pos, ch & 0x3f, accept_invalids)
                flg && (invalids += 1 ; continue)
                if ch > 0xff
                    num2byte += 1
                elseif ch > 0x7f
                    latin1byte += 1
                elseif accept_long_char
                    flags |= UTF_LONG
                elseif (ch == 0) && accept_long_null
                    flags |= UTF_LONG
                elseif accept_invalids
                    invalids += 1
                else
                    unierror(UTF_ERR_LONG, pos, ch)
                end
             elseif ch < 0xf0
                # 3-byte UTF-8 sequence (i.e. characters 0x800-0xffff)
                if pos + 1 > endpos
                    accept_invalids || unierror(UTF_ERR_SHORT, pos, ch)
                    invalids += 1
                    break
                end
                ch, pos, flg = check_continuation(dat, pos, ch & 0x0f, accept_invalids)
                flg && (invalids += 1 ; continue)
                ch, pos, flg = check_continuation(dat, pos, ch, accept_invalids)
                flg && (invalids += 1 ; continue)
                # check for surrogate pairs, make sure correct
                if is_surrogate_codeunit(ch)
                    if !is_surrogate_lead(ch)
                        accept_invalids || unierror(UTF_ERR_NOT_LEAD, pos-2, ch)
                        invalids += 1
                        continue
                    end
                    # next character *must* be a trailing surrogate character
                    if pos + 2 > endpos
                        accept_invalids || unierror(UTF_ERR_MISSING_SURROGATE, pos-2, ch)
                        invalids += 1
                        break
                    end
                    byt = get_codeunit(dat, pos)
                    pos += 1
                    if byt != 0xed
                        accept_invalids || unierror(UTF_ERR_NOT_TRAIL, pos, byt)
                        invalids += 1
                        continue
                    end
                    surr, pos, flg = check_continuation(dat, pos, 0x0000d, accept_invalids)
                    flg && (invalids += 1 ; continue)
                    surr, pos, flg = check_continuation(dat, pos, surr, accept_invalids)
                    flg && (invalids += 1 ; continue)
                    if !is_surrogate_trail(surr)
                        accept_invalids || unierror(UTF_ERR_NOT_TRAIL, pos-2, surr)
                        invalids += 1
                    elseif !accept_surrogates
                        accept_invalids || unierror(UTF_ERR_SURROGATE, pos-2, surr)
                        invalids += 1
                    else
                        flags |= UTF_SURROGATE
                        num4byte += 1
                    end
                elseif ch > 0x07ff
                    num3byte += 1
                elseif accept_long_char
                    flags |= UTF_LONG
                    num2byte += 1
                elseif accept_invalids
                    invalids += 1
                else
                    unierror(UTF_ERR_LONG, pos-2, ch)
                end
            elseif ch < 0xf5
                # 4-byte UTF-8 sequence (i.e. characters > 0xffff)
                if pos + 2 > endpos
                    accept_invalids || unierror(UTF_ERR_SHORT, pos, ch)
                    invalids += 1
                    break
                end
                ch, pos, flg = check_continuation(dat, pos, ch & 0x07, accept_invalids)
                flg && (invalids += 1 ; continue)
                ch, pos, flg = check_continuation(dat, pos, ch, accept_invalids)
                flg && (invalids += 1 ; continue)
                ch, pos, flg = check_continuation(dat, pos, ch, accept_invalids)
                flg && (invalids += 1 ; continue)
                if ch > 0x10ffff
                    accept_invalids || unierror(UTF_ERR_INVALID, pos-3, ch)
                    invalids += 1
                elseif ch > 0xffff
                    num4byte += 1
                elseif is_surrogate_codeunit(ch)
                    accept_invalids || unierror(UTF_ERR_SURROGATE, pos-3, ch)
                    invalids += 1
                elseif accept_long_char
                    # This is an overly long encoded character
                    flags |= UTF_LONG
                    if ch > 0x7ff
                        num3byte += 1
                    elseif ch > 0x7f
                        num2byte += 1
                    end
                elseif accept_invalids
                    invalids += 1
                else
                    unierror(UTF_ERR_LONG, pos-2, ch)
                end
            elseif accept_invalids
                invalids += 1
            else
                unierror(UTF_ERR_INVALID, pos, ch)
            end
        end
    end
    _ret_check(totalchar, flags, invalids, latin1byte, num2byte, num3byte, num4byte)
end

function unsafe_check_string(dat::Union{AbstractVector{T}, Ptr{T}}, pos, endpos;
                             accept_long_null  = false,
                             accept_surrogates = false,
                             accept_long_char  = false,
                             accept_invalids   = false) where {T<:Union{UInt16,UInt32}}
    flags = 0%UInt
    totalchar = latin1byte = num2byte = num3byte = num4byte = invalids = 0
    @inbounds while pos <= endpos
        ch = get_codeunit(dat, pos)%UInt32
        pos += 1
        totalchar += 1
        if ch > 0x7f
            if ch < 0x100
                latin1byte += 1
            elseif ch < 0x800
                num2byte += 1
            elseif ch > 0x0ffff
                if (ch > 0x10ffff)
                    accept_invalids || unierror(UTF_ERR_INVALID, pos, ch)
                    invalids += 1
                else
                    num4byte += 1
                end
            elseif !is_surrogate_codeunit(ch)
                num3byte += 1
            elseif is_surrogate_lead(ch)
                if pos > endpos
                    accept_invalids || unierror(UTF_ERR_MISSING_SURROGATE, pos, ch)
                    invalids += 1
                    break
                end
                # next character *must* be a trailing surrogate character
                ch = get_codeunit(dat, pos)
                pos += 1
                if !is_surrogate_trail(ch)
                    accept_invalids || unierror(UTF_ERR_NOT_TRAIL, pos, ch)
                    invalids += 1
                elseif typeof(dat) <: AbstractVector{UInt16} # fix this test!
                    num4byte += 1
                elseif accept_surrogates
                    flags |= UTF_SURROGATE
                    num4byte += 1
                elseif accept_invalids
                    invalids += 1
                else
                    unierror(UTF_ERR_SURROGATE, pos, ch)
                end
            elseif accept_invalids
                invalids += 1
            else
                unierror(UTF_ERR_NOT_LEAD, pos, ch)
            end
        end
    end
    _ret_check(totalchar, flags, invalids, latin1byte, num2byte, num3byte, num4byte)
end

function unsafe_check_string(str::T;
                             accept_long_null  = false,
                             accept_surrogates = false,
                             accept_long_char  = false,
                             accept_invalids   = false) where {T<:AbstractString}
    flags = 0%UInt
    totalchar = latin1byte = num2byte = num3byte = num4byte = invalids = 0
    pos = start(str)
    @inbounds while !done(str, pos)
        chr, nxt = next(str, pos)
        ch = chr%UInt32
        totalchar += 1
        if ch > 0x7f
            if ch < 0x100
                latin1byte += 1
            elseif ch < 0x800
                num2byte += 1
            elseif ch > 0x0ffff
                if (ch > 0x10ffff)
                    accept_invalids || unierror(UTF_ERR_INVALID, pos, ch)
                    invalids += 1
                else
                    num4byte += 1
                end
            elseif !is_surrogate_codeunit(ch)
                num3byte += 1
            elseif is_surrogate_lead(ch)
                if done(str, nxt)
                    accept_invalids || unierror(UTF_ERR_MISSING_SURROGATE, pos, ch)
                    invalids += 1
                    break
                end
                # next character *must* be a trailing surrogate character
                chr, nxt = next(str, nxt)
                if !is_surrogate_trail(chr)
                    accept_invalids || unierror(UTF_ERR_NOT_TRAIL, pos, chr)
                    invalids += 1
                elseif accept_surrogates
                    flags |= UTF_SURROGATE
                    num4byte += 1
                elseif accept_invalids
                    invalids += 1
                else
                    unierror(UTF_ERR_SURROGATE, pos, ch)
                end
            elseif accept_invalids
                invalids += 1
            else
                unierror(UTF_ERR_NOT_LEAD, pos, ch)
            end
        end
        pos = nxt
    end
    _ret_check(totalchar, flags, invalids, latin1byte, num2byte, num3byte, num4byte)
end

@inline function skipascii(beg::Ptr{UInt8}, fin::Ptr{UInt8})
    align = reinterpret(UInt, beg)
    pnt = reinterpret(Ptr{UInt64}, align & ~CHUNKMSK)
    v = unsafe_load(pnt)
    (align &= CHUNKMSK) != 0 && (v &= ~_mask_bytes(align))
    while (pnt += CHUNKSZ) < fin
        # find first byte that is not ASCII
        (v &= hi_mask) == 0 || return Int(pnt - CHUNKSZ - beg + (trailing_zeros(v)>>>3))
        v = unsafe_load(pnt)
    end
    Int(pnt - CHUNKSZ - beg +
        (trailing_zeros((pnt == fin ? v : (v | ~_mask_bytes(fin - pnt - CHUNKSZ))) & hi_mask)>>>3))
end

function fast_check_string(beg::Ptr{UInt8}, len)
    pnt = beg
    fin = pnt + len
    flags = 0%UInt
    asciichar = latin1byte = num2byte = num3byte = num4byte = 0
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch < 0x80
            cnt = skipascii(pnt, fin)
            asciichar += cnt
            pnt += cnt
            continue
        # Check UTF-8 encoding
        elseif ch < 0xe0
            # 2-byte UTF-8 sequence (i.e. characters 0x80-0x7ff)
            (pnt += 1) < fin || unierror(UTF_ERR_SHORT, pnt - beg, ch)
            (ch > 0xc1 && checkcont(pnt)) || unierror(UTF_ERR_INVALID, pnt - beg, ch)
            ch > 0xc3 ? (num2byte += 1) : (latin1byte += 1)
        elseif ch < 0xf0
            # 3-byte UTF-8 sequence (i.e. characters 0x800-0xffff)
            (pnt += 2) < fin || unierror(UTF_ERR_SHORT, pnt - beg - 1, ch)
            b2 = get_codeunit(pnt - 1)
            (is_valid_continuation(b2) && checkcont(pnt)) ||
                unierror(UTF_ERR_INVALID, pnt - beg - 1, ch)
            if ch == 0xe0 # Might be overlong
                b2 < 0xa0 && unierror(UTF_ERR_LONG, pnt - beg - 1, get_utf8_3byte(pnt, ch))
            elseif ch == 0xed # Might be surrogate pair
                b2 > 0x9f && unierror(UTF_ERR_SURROGATE, pnt - beg - 1, get_utf8_3byte(pnt, ch))
            end
            num3byte += 1
        elseif ch < 0xf5
            # 4-byte UTF-8 sequence (i.e. characters > 0xffff)
            (pnt += 3) < fin || unierror(UTF_ERR_SHORT, pnt - beg - 2, ch)
            b2 = get_codeunit(pnt - 2)
            #println(ch,", ",b2,", ",get_codeunit(pnt-1),", ",get_codeunit(pnt))
            (is_valid_continuation(b2) && checkcont(pnt-1) && checkcont(pnt)) ||
                unierror(UTF_ERR_INVALID, pnt - beg - 2, ch)
            if ch == 0xf0
                b2 < 0x90 && unierror(UTF_ERR_LONG, pnt - beg - 2, get_utf8_4byte(pnt, ch))
            elseif ch == 0xf4
                b2 > 0x8f && unierror(UTF_ERR_INVALID, pnt - beg - 2, get_utf8_4byte(pnt, ch))
            end
            num4byte += 1
        else
            unierror(UTF_ERR_INVALID, pnt - beg + 1, ch)
        end
        pnt += 1
    end
    _ret_check(asciichar+latin1byte+num2byte+num3byte+num4byte,
               flags, 0, latin1byte, num2byte, num3byte, num4byte)
end

const _bmp_mask = 0xd800_d800_d800_d800
@inline _mask_allsurr(v)  = xor((v | v<<1 | v<<2 | v<<3 | v<<4) & _hi_bit_16, _hi_bit_16)

@inline _get_bmp_mask(v::UInt64) = _mask_allsurr(xor(v, _bmp_mask))

const msk_ascii_16 = 0xff80ff80ff80ff80
const msk_latin_16 = 0xff00ff00ff00ff00
const msk_num2b_16 = 0xf800f800f800f800

# v -> 1111 1yyy l000 0000 1111 1yyy l000 0000 1111 1yyy l000 0000 1111 1yyy l000 0000

@inline function countmask(v, cnta, cntl, cnt2, cnt3)
    # First mask off ASCII bits: 0xff80ff80ff80ff80, count all 0 words (0-4)
    # Then mask off Latin1 bits: 0xff00ff00ff00ff00, count all 0 words (0-4)
    # Then mask off num2byte:    0xf800f800f800f800, count all 0 words (0-4)
    va = v & msk_ascii_16
    vl = v & msk_latin_16
    v2 = v & msk_num2b_16
    println("v=$v, 0x$(hex(va)), 0x$(hex(vl)), 0x$(hex(v2)), $cnta, $cntl, $cnt2, $cnt3")
    cnta, cntl, cnt2, cnt3
end

@inline function skipbmp(beg::Ptr{UInt16}, fin, cnta, cntl, cnt2, cnt3)
    align = reinterpret(UInt, beg)
    pnt = reinterpret(Ptr{UInt64}, align & ~CHUNKMSK)
    v = unsafe_load(pnt)
    (align &= CHUNKMSK) != 0 && (v &= ~_mask_bytes(align))
    while (pnt += CHUNKSZ) < fin
        # If vm is 0, it means there are non-BMP characters present
        if (vm = _get_bmp_mask(v)) == 0
            # Get how many characters are before the non-BMP character
            cnt = (leading_zeros(vm)>>>3)&~1
            v &= _mask_bytes(cnt)
            return pnt - CHUNKSZ + cnt, countmask(v, cnta, cntl, cnt2, cnt3)...
        end
        cnta, cntl, cnt2, cnt3 = countmask(v, cnta, cntl, cnt2, cnt3)
        v = unsafe_load(pnt)
    end
    v = (pnt == fin ? v : (v & _mask_bytes(fin - pnt - CHUNKSZ)))
    if (vm = _get_bmp_mask(v)) == 0
        # Get how many characters are before the non-BMP character
        cnt = (leading_zeros(vm)>>>3)&~1
        v &= _mask_bytes(cnt)
        fin = pnt - CHUNKSZ + cnt
    end
    fin, countmask(v, cnta, cntl, cnt2, cnt3)...
end

# Check for valid UTF-16
function fast_check_string(beg::Ptr{UInt16}, len)
    pnt = beg
    fin = bytoff(pnt, len)
    flags = 0%UInt
    asciichar = latin1byte = num2byte = num3byte = num4byte = 0
    while pnt < fin
        ch = get_codeunit(pnt)
        #=
        if !is_surrogate_codepoint(ch)
            pnt, asciichar, latin1byte, num2byte, num3byte =
                skipbmp(pnt, fin, asciichar, latin1byte, num2byte, num3byte)
        =#
        if ch <= 0x7f
            asciichar += 1
        elseif ch <= 0xff
            latin1byte += 1
        elseif ch < 0x7ff
            num2byte += 1
        elseif !is_surrogate_codeunit(ch)
            num3byte += 1
        elseif !is_surrogate_lead(ch)
            unierror(UTF_ERR_NOT_LEAD, Int((pnt - beg)>>>1), ch)
        elseif (pnt += 2) >= fin
            unierror(UTF_ERR_SHORT, Int((pnt - beg - 2)>>>1), ch)
        else
            c2 = get_codeunit(pnt)
            is_surrogate_trail(c2) || unierror(UTF_ERR_NOT_TRAIL, Int((pnt - beg)>>>1), c2)
            pnt += 2
            num4byte += 1
        end
        pnt += 2
    end
    _ret_check(asciichar+latin1byte+num2byte+num3byte+num4byte,
               flags, 0, latin1byte, num2byte, num3byte, num4byte)
end

# Check for valid UTF-32
function fast_check_string(beg::Ptr{UInt32}, len)
    pnt = beg
    fin = bytoff(pnt, len)
    flags = 0%UInt
    asciichar = latin1byte = num2byte = num3byte = num4byte = 0
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch <= 0x7f
            asciichar += 1
        elseif ch <= 0xff
            latin1byte += 1
        elseif ch <= 0x7ff
            num2byte += 1
        elseif ch <= 0xd7ff
            num3byte += 1
        elseif ch <= 0xdfff
            unierror(UTF_ERR_SURROGATE, (pnt-beg+4)>>>2, ch)
        elseif ch <= 0xffff
            num3byte += 1
        elseif ch <= 0x10ffff
            num4byte += 1
        else
            unierror(UTF_ERR_INVALID, (pnt-beg+4)>>>2, ch)
        end
        pnt += 4
    end
    _ret_check(asciichar+latin1byte+num2byte+num3byte+num4byte,
               flags, 0, latin1byte, num2byte, num3byte, num4byte)
end

"""
Calculate the total number of characters, as well as number of
latin1, 2-byte, 3-byte, and 4-byte sequences in a validated UTF-8 string
"""
function count_chars(::Type{UTF8Str}, ::Type{S}, dat::Union{AbstractVector{S}, Ptr{S}},
                     pos, len) where {S<:CodeUnitTypes}
    totalchar = latin1byte = num2byte = num3byte = num4byte = 0
    @inbounds while pos <= len
        ch = get_codeunit(dat, pos)
        totalchar += 1
        if ch < 0x80 # ASCII characters
            pos += 1
        elseif ch < 0xc4 # 2-byte Latin 1 characters (0x80-0xff)
            latin1byte += 1
            pos += 2
        elseif ch < 0xe0 # 2-byte BMP sequence (i.e. characters 0x100-0x7ff)
            pos += 2
            num2byte += 1
        elseif ch < 0xf0 # 3-byte BMP sequence (0x800-0xffff)
            pos += 3
            num3byte += 1
        else # 4-byte non-BMP sequence (0x10000 - 0x10ffff)
            pos += 4
            num4byte += 1
        end
    end
    _ret_check(totalchar, 0%UInt, 0, latin1byte, num2byte, num3byte, num4byte)
end

"""
Calculate the total number of characters, as well as number of
latin1, 2-byte, 3-byte, and 4-byte sequences in a validated UTF-16, UCS2, or UTF-32 string
"""
function count_chars(::Type{T}, ::Type{S}, dat::Union{AbstractVector{S}, Ptr{S}},
                     pos, endpos) where {S<:CodeUnitTypes, T<:Str}
    totalchar = latin1byte = num2byte = num3byte = num4byte = 0
    @inbounds while pos <= endpos
        ch = get_codeunit(dat, pos)%UInt32
        pos += 1
        totalchar += 1
        if ch <= 0x7f # ASCII characters
        elseif ch <= 0xff # 2-byte Latin 1 characters (0x80-0xff)
            latin1byte += 1
        elseif ch <= 0x7ff # 2-byte BMP sequence (i.e. characters 0x100-0x7ff)
            num2byte += 1
        elseif T == UTF16Str
            if is_surrogate_lead(ch)
                pos += 1
                num4byte += 1
            else
                num3byte += 1
            end
        elseif ch <= 0xffff # 3-byte BMP sequence (0x800-0xffff)
            num3byte += 1
        else # 4-byte non-BMP sequence (0x10000 - 0x10ffff)
            num4byte += 1
        end
    end
    _ret_check(totalchar, 0%UInt, 0, latin1byte, num2byte, num3byte, num4byte)
end

count_chars(T, dat, len) = count_chars(T, codeunit(T), dat, 1, len)

"""
Calculate the total number of bytes > 0x7f
"""
function count_latin(len, pnt::Ptr{UInt8})
    # Todo: optimize this to work on chunks when pnt is aligned
    cnt = 0
    fin = pnt + len
    while pnt < fin
        cnt += (get_codeunit(pnt) > 0x7f)
        pnt += 1
    end
    cnt
end

"""
Validates and calculates number of characters in a UTF-8,UTF-16 or UTF-32 encoded vector/string

This function checks the bounds of the start and end positions
Use `unsafe_check_string` to avoid that overhead if the bounds have already been checked

Input Arguments:

* `dat`    UTF-8 (`Vector{UInt8}`), UTF-16 (`Vector{UInt16}`) or UTF-32 (`Vector{UInt32}`, `AbstractString`) encoded string

Optional Input Arguments:

* `startpos` start position (defaults to 1)
* `endpos`   end position   (defaults to `lastindex(dat)`)

Keyword Arguments:

* `accept_long_null`  = `false`  # Modified UTF-8 (`\\0` represented as `b\"\\xc0\\x80\"`)
* `accept_surrogates` = `false`  # `CESU-8`
* `accept_long_char`  = `false`  # Accept arbitrary long encodings
* `accept_invalids`   = `false`  # Accept invalid sequences (to be replaced on conversion)

Returns:

* (total characters, flags, 4-byte, 3-byte, 2-byte)

Throws:

* `UnicodeError`
"""
function check_string end

# No need to check bounds if using defaults
check_string(dat; kwargs...) = unsafe_check_string(dat, 1, lastindex(dat); kwargs...)

# Make sure that beginning and end positions are bounds checked
function check_string(dat, startpos, endpos = lastindex(dat); kwargs...)
    @boundscheck checkbounds(dat, startpos)
    @boundscheck checkbounds(dat, endpos)
    endpos < startpos && argerror(startpos, endpos)
    unsafe_check_string(dat, startpos, endpos; kwargs...)
end

byte_string_classify(data) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), data, length(data))
byte_string_classify(data::Vector{UInt8}) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), data, length(data))
byte_string_classify(s::Str{<:Byte_CSEs}) = byte_string_classify(s.data)
    # 0: neither valid ASCII nor UTF-8
    # 1: valid ASCII
    # 2: valid UTF-8

is_valid(::Type{ASCIIStr},  s::Vector{UInt8}) = is_ascii(s)
is_valid(::Type{UTF8Str},   s::Vector{UInt8}) = byte_string_classify(s) != 0
is_valid(::Type{<:Str{LatinCSE}},  s::Vector{UInt8}) = true
# This should be optimized, stop at first character > 0x7f
is_valid(::Type{<:Str{_LatinCSE}}, s::Vector{UInt8}) = !is_ascii(s)

is_valid(::Type{UniStr}, s::String) = is_unicode(s)
is_valid(::Type{<:Str{C}}, s::String) where {C<:Union{UTF8CSE,UTF16CSE,UTF32CSE}} = is_unicode(s)
is_valid(::Type{<:Str{ASCIICSE}}, s::String) = is_ascii(s)
is_valid(::Type{<:Str{UCS2CSE}}, s::String)  = is_bmp(s)
is_valid(::Type{<:Str{LatinCSE}}, s::String) = is_latin(s)

function _copysub(pnt::Ptr{T}, len) where {T<:CodeUnitTypes}
    buf, out = _allocate(eltype(pnt), len)
    _memcpy(out, pnt, len)
    buf
end
_copysub(str::String)    = str
_copysub(str::Str)       = str.data
_copysub(str::SubString) = @preserve str _copysub(pointer(str), ncodeunits(str))
_copysub(vec::Vector)    = @preserve vec _copysub(pointer(vec), length(vec))

function _cvtsize(::Type{T}, dat, len) where {T <: CodeUnitTypes}
    buf, pnt = _allocate(T, len)
    @inbounds for i = 1:len
        set_codeunit!(pnt, get_codeunit(dat, i))
        pnt += sizeof(T)
    end
    buf
end

function _cvtsize(::Type{T}, pnt::Ptr{S}, len) where {S<:CodeUnitTypes,T<:CodeUnitTypes}
    buf, out = _allocate(T, len)
    fin = bytoff(pnt, len)
    #println("buf=$buf, out=$out, pnt=$pnt, fin=$fin, len=$len")
    while pnt < fin
        set_codeunit!(out, get_codeunit(pnt))
        out += sizeof(T)
        pnt += sizeof(S)
    end
    buf
end

(*)(s1::Union{C1, S1}, ss::Union{C2, S2}...) where {C1<:Chr,C2<:Chr,S1<:Str,S2<:Str} =
    string(s1, ss...)

thisind(str::MaybeSub{<:Str}, i::Integer) = thisind(str, Int(i))

function filter(fun, str::MaybeSub{T}) where {C<:CSE,T<:Str{C}}
    out = get_iobuffer(sizeof(str))
    @inbounds for ch in codepoints(str)
        fun(ch) && _write(C, out, ch)
    end
    Str{C}(String(take!(out)))
end

# Todo: These should be optimized based on the traits, and return internal substrings, once
# I've implemented those

first(str::Str, n::Integer) = str[1:min(end, nextind(str, 0, n))]
last(str::Str, n::Integer)  = str[max(1, prevind(str, ncodeunits(str)+1, n)):end]

const Chrs = @static V6_COMPAT ? Union{Char,AbstractChar} : Chr

function repeat(ch::CP, cnt::Integer) where {CP <: Chrs}
    C = codepoint_cse(CP)
    cnt > 1 && return Str(C, _repeat(CodePointStyle(C), C, codepoint(ch), cnt))
    cnt < 0 && repeaterr(cnt)
    cnt == 0 ? empty_str(C) : _convert(C, codepoint(ch))
end

(^)(ch::CP, cnt::Integer) where {CP <: Chrs} = repeat(ch, cnt)

# low level mem support functions

const (WidChr,OthChr) = @static sizeof(Cwchar_t) == 4 ? (UInt32,UInt16) : (UInt16,UInt32)

_fwd_memchr(ptr::Ptr{UInt8}, byt::UInt8, len::Integer) =
    ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), ptr, byt, len)

_fwd_memchr(ptr::Ptr{WidChr}, wchr::WidChr, len::Integer) =
    ccall(:wmemchr, Ptr{WidChr}, (Ptr{WidChr}, Int32, Csize_t), ptr, wchr, len)

_fwd_memchr(beg::Ptr{OthChr}, wchr::OthChr, len::Integer) =
    _fwd_memchr(beg, ch, bytoff(beg, len))

_fwd_memchr(ptr::Ptr{UInt8}, byt::UInt8, fin::Ptr{UInt8}) =
    ptr < fin ? _fwd_memchr(ptr, byt, fin - ptr) : C_NULL

_fwd_memchr(ptr::Ptr{WidChr}, wchr::WidChr, fin::Ptr{WidChr}) =
    ptr < fin ? _fwd_memchr(ptr, wchr, chroff(fin - ptr)) : C_NULL

function _fwd_memchr(pnt::Ptr{T}, wchr::T, fin::Ptr{T}) where {T<:OthChr}
    while pnt < fin
        get_codeunit(pnt) == ch && return pnt
        pnt += sizeof(T)
    end
    C_NULL
end

_rev_memchr(ptr::Ptr{UInt8}, byt::UInt8, len::Integer) =
    ccall(:memrchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), ptr, byt, len)

function _rev_memchr(beg::Ptr{T}, ch::T, len::Integer) where {T<:Union{UInt16,UInt32}}
    pnt = bytoff(beg, len)
    while (pnt -= sizeof(T)) >= beg
        get_codeunit(pnt) == ch && return pnt
    end
    C_NULL
end

_memcmp(a::Ptr{UInt8}, b::Ptr{UInt8}, len) =
    ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt), a, b, len)
_memcmp(a::Ptr{WidChr}, b::Ptr{WidChr}, len) =
    ccall(:wmemcmp, Int32, (Ptr{WidChr}, Ptr{WidChr}, UInt), a, b, len)

function _memcmp(apnt::Ptr{OthChr}, bpnt::Ptr{OthChr}, len)
    fin = bytoff(apnt, len)
    while apnt < fin
        (c1 = get_codeunit(apnt)) == (c2 = get_codeunit(bpnt)) || return _cmp(c1, c2)
        apnt += sizeof(OthChr)
        bpnt += sizeof(OthChr)
    end
    0
end

# These should probably be handled by traits, or dispatched by getting the codeunit type for each
_memcmp(a::Union{String, Str{<:Byte_CSEs}},
        b::Union{String, Str{<:Byte_CSEs}, SubString{String}, SubString{<:Str{<:Byte_CSEs}}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::SubString{<:Union{String, Str{<:Byte_CSEs}}}, b::Union{String, Str{<:Byte_CSEs}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::SubString{<:Union{String, Str{<:Byte_CSEs}}},
        b::SubString{<:Union{String, Str{<:Byte_CSEs}}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)

_memcmp(a::Str{<:Word_CSEs}, b::MaybeSub{<:Str{<:Word_CSEs}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::Str{<:Quad_CSEs}, b::MaybeSub{<:Str{<:Quad_CSEs}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::SubString{<:Str{<:Word_CSEs}}, b::Str{<:Word_CSEs}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::SubString{<:Str{<:Quad_CSEs}}, b::Str{<:Quad_CSEs}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::SubString{<:Str{<:Word_CSEs}}, b::SubString{<:Str{Word_CSEs}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)
_memcmp(a::SubString{<:Str{<:Quad_CSEs}}, b::SubString{<:Str{Quad_CSEs}}, siz) =
    _memcmp(pointer(a), pointer(b), siz)

_memcpy(dst::Ptr{UInt8}, src::Ptr, siz) =
    ccall(:memcpy, Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), dst, src, siz)
_memcpy(a::Ptr{WidChr}, b::Ptr{WidChr}, len) =
    ccall(:wmemcpy, Ptr{WidChr}, (Ptr{WidChr}, Ptr{WidChr}, UInt), a, b, len)
_memcpy(a::Ptr{OthChr}, b::Ptr{OthChr}, len) =
    ccall(:memcpy, Ptr{OthChr}, (Ptr{OthChr}, Ptr{OthChr}, UInt), a, b, bytoff(OthChr, len))

@inline _memset(pnt::Ptr{UInt8}, ch::UInt8, cnt) =
    ccall(:memset, Ptr{UInt8}, (Ptr{UInt8}, Cint, Csize_t), pnt, ch, cnt)
@inline _memset(pnt::Ptr{WidChr}, ch::WidChr, cnt) =
    ccall(:wmemset, Ptr{WidChr}, (Ptr{WidChr}, Cuint, Csize_t), pnt, ch, cnt)

@inline function _memset(pnt::Ptr{OthChr}, ch::OthChr, cnt)
    fin = bytoff(pnt, cnt)
    while pnt < fin
        set_codeunit!(pnt, ch)
        pnt += sizeof(OthChr)
    end
end

# Little-endian output here
@inline get_utf8_16(ch) =
    (ch >>> 6) | ((ch & 0x3f)%UInt16<<8) | 0x80c0
@inline get_utf8_32(ch) =
    (ch & 0xc0000 >>> 18) | (ch & 0x3f000 >>> 4) |
    (ch & 0xfc0 << 10) | (ch & 0x3f)<<24 | 0x808080f0
@inline get_utf16_32(ch) =
    (0xd7c0 + (ch >>> 10))%UInt16 << 6 | (0xdc00 + (ch & 0x3ff))%UInt32

@inline function _repeat_chr(::Type{T}, ch, cnt) where {T<:CodeUnitTypes}
    #println("_repeat_chr($T, $ch, $cnt)")
    buf, pnt = _allocate(T, cnt)
    _memset(pnt, ch%T, cnt)
    buf
end

@inline function _repeat_3(ch, cnt)
    buf, pnt = _allocate(UInt8, cnt*3)
    b1, b2, b3 = get_utf8_3(ch)
    fin = pnt + cnt*3
    while pnt < fin
        set_codeunit!(pnt,     b1)
        set_codeunit!(pnt + 1, b2)
        set_codeunit!(pnt + 2, b3)
        pnt += 3
    end
    buf
end

_repeat(::CodeUnitSingle, ::Type{C}, ch::T, cnt) where {T,C<:CSE} =
    _repeat_chr(basetype(T), ch, cnt)

function _repeat(::CodeUnitMulti, ::Type{UTF8CSE}, ch, cnt)
    if ch <= 0x7f
        _repeat_chr(UInt8, ch, cnt)
    elseif ch <= 0x7ff
        _repeat_chr(UInt16, get_utf8_16(ch), cnt)
    elseif ch <= 0xffff
        _repeat_3(ch, cnt)
    else
        _repeat_chr(UInt32, get_utf8_32(ch), cnt)
    end
end

_repeat(::CodeUnitMulti, ::Type{UTF16CSE}, ch, cnt) =
    ch <= 0xffff ? _repeat_chr(UInt16, ch, cnt) : _repeat_chr(UInt32, get_utf16_32(ch), cnt)

function repeat(str::T, cnt::Integer) where {C<:CSE,T<:Str{C}}
    cnt < 2 && return cnt == 1 ? str : (cnt == 0 ? empty_str(C) : repeaterr(cnt))
    CU = codeunit(T)
    @preserve str begin
        len = ncodeunits(str)
        totlen = len * cnt
        buf, out = _allocate(CU, totlen)
        if len == 1 # common case: repeating a single codeunit string
            _memset(out, get_codeunit(pointer(str)), cnt)
        else
            pnt = pointer(str)
            fin = bytoff(out, totlen)
            siz = bytoff(CU, len)
            while out < fin
                _memcpy(out, pnt, len)
                out += siz
            end
        end
    end
    Str(C, buf)
end
(^)(str::T, cnt::Integer) where {T<:Str} = repeat(str, cnt)

# Definitions for C compatible strings, that don't allow embedded
# '\0', and which are terminated by a '\0'

containsnul(str::Str{<:Byte_CSEs}) = containsnul(unsafe_convert(Ptr{Cchar}, str), sizeof(str))

# Check 4 characters at a time
function containsnul(str::Str{<:Word_CSEs})
    (siz = sizeof(str)) == 0 && return true
    @preserve str begin
        pnt, fin = _calcpnt(str, siz)
        while (pnt += CHUNKSZ) <= fin
            ((v = unsafe_load(pnt))%UInt16 == 0 || (v>>>16)%UInt16 == 0 ||
             (v>>>32)%UInt16 == 0 || (v>>>48) == 0) && return true
        end
        pnt - CHUNKSZ != fin &&
            ((v = (unsafe_load(pnt) | ~_mask_bytes(siz)))%UInt16 == 0 ||
             (v>>>16)%UInt16 == 0 || (v>>>32)%UInt16 == 0)
    end
end

function containsnul(str::Str{<:Quad_CSEs})
    (siz = sizeof(str)) == 0 && return true
    @preserve str begin
        pnt, fin = _calcpnt(str, siz)
        while (pnt += CHUNKSZ) <= fin
            ((v = unsafe_load(pnt))%UInt32 == 0 || (v>>>32) == 0) && return true
        end
        pnt - CHUNKSZ != fin && unsafe_load(reinterpret(Ptr{UInt32}, pnt)) == 0x00000
    end
end


text:
using Clustering
using DataFrames
using BioFeatures
using Distances

function clusters_by_size_order_permutation(cluster_object)
    cluster_assignments =  assignments(cluster_object)
    cluster_counts      =  counts(cluster_object)
    println("CLUSTER COUNTS: ", cluster_counts)

    cluster_size = {i => cluster_counts[i] for i in cluster_assignments}
    cluster_assignments_size_labels = Int64[]

    for c in cluster_assignments
         c_size = cluster_size[c]
         push!( cluster_assignments_size_labels, c_size )
    end
    #sort_order = sortperm(cluster_assignments_size_labels)
  println("finished cluster_assignments_size_labels")
  return cluster_assignments_size_labels
end
#=
function clusters_by_size_order_permutation(cluster_assignments)
    clusters = unique(cluster_assignments)
    println(clusters)
    num_clusters = length(clusters)

    # obtain counts in each cluster label
    h = hist(cluster_assignments, [i for i=0:num_clusters])
    println(h)
    h = h[2]
    println("Got these clusters $h")
    # now make dictionary mapping the current labels (1 to n )
    # to a label that is the size of how many are in the cluster
    cluster_size = {i => h[i] for i = 1:num_clusters}
    println("cluster size: $cluster_size")
    cluster_assignments_size_labels = Int64[]
    println("do for loop btw: cluster_assignments_size_labels")
    for c in cluster_assignments
         c_size = cluster_size[c]
         push!( cluster_assignments_size_labels, c_size )
    end
    # println("pushed on to $cluster_assignments_size_labels")
  #sort_order = sortperm(cluster_assignments_size_labels)
  return cluster_assignments_size_labels
end
function order_vector(vector,cluster_assignments)
    #sort_order = sortperm(cluster_assignments)
    cluster_assignments_size_labels = clusters_by_size_order_permutation(cluster_assignments)
    sort_order = sortperm(cluster_assignments_size_labels)
    sort!(cluster_assignments_size_labels)
    return ( vector[sort_order], cluster_assignments_size_labels)
end
=#
function order_vector(vector,cluster_object)
    #sort_order = sortperm(cluster_assignments)
    cluster_assignments_size_labels = clusters_by_size_order_permutation(cluster_object)
    sort_order = sortperm(cluster_assignments_size_labels)
    sort!(cluster_assignments_size_labels)
    return ( vector[sort_order], cluster_assignments_size_labels)
end

function order_matrix(tiles_matrix,cluster_object)
    #cluster_assignments =  assignments(cluster_object)
    #cluster_counts      =  counts(cluster_object)
    #cluster_size = {i => cluster_counts[i] for i in cluster_assignments}
    #cluster_assignments_size_labels = Int64[]
    #println("do for loop btw: cluster_assignments_size_labels")
    #for c in cluster_assignments
    #     c_size = cluster_size[c]
    #     push!( cluster_assignments_size_labels, c_size )
    #end
    #sort_order = sortperm(cluster_assignments_size_labels)
   cluster_assignments_size_labels = clusters_by_size_order_permutation(cluster_object)
   sort_order = sortperm(cluster_assignments_size_labels)
   return tiles_matrix[ sort_order,:]
end

#=
function order_matrix( tiles_matrix, cluster_assignments)

    # we want each row to be together with its other cluster members
    # AND we want those groups to be sorted by the number that are in each cluster
    # This means changing the number of the cluster (e.g. 1,2,3 and 4) to the amount
    # that are in the cluster (e.g. 13,445,34 and 67) and sorting on that .

    clusters = unique(cluster_assignments)
    println("CLUSTERS: $clusters")
    num_clusters = length(clusters)

    # obtain counts in each cluster label
    #h = hist(cluster_assignments, [i for i=0:num_clusters])
    h = Dict()
    for cl in cluster_assignments
        if haskey(h,cl)
            h[cl]+=1
        else
           h[cl] = 1
       end
    end
    println("COUNTS PER CLUSTER: $h")
    #h = h[2]
    println("Got these clusters $h")
    # now make dictionary mapping the current labels (1 to n )
    # to a label that is the size of how many are in the cluster
    cluster_size = {i => h[i] for i in cluster_assignments}
    println("cluster size: $cluster_size")
    cluster_assignments_size_labels = Int64[]
    println("do for loop btw: cluster_assignments_size_labels")
    for c in cluster_assignments
         c_size = cluster_size[c]
         push!( cluster_assignments_size_labels, c_size )
    end
    # println("pushed on to $cluster_assignments_size_labels")
    sort_order = sortperm(cluster_assignments_size_labels)
    return tiles_matrix[ sort_order,:] #this doesn't work
end
=#

function cluster_tiles(tiles_matrix;method="kmeans", num_clusters_range=5:5, kmeans_repeats=1,  num_clusters=5,eps=5,minpts=100)

   cluster_assignments = nothing
   tiles_matrix_transpose = tiles_matrix'
   sil_vals = nothing
   if (method == "kmeans")
      #cluster_object = kmeans_cluster_tiles(tiles_matrix_transpose,num_clusters)
      (cluster_object, sil_vals) = kmeans_cluster_tiles(tiles_matrix_transpose,kmeans_repeats, num_clusters_range)
   elseif( method == "affprop"  )
      cluster_object = aff_prop_cluster_tiles(tiles_matrix)
   elseif(method == "dbscan")
     cluster_object = dbscan_hammings(tiles_matrix_transpose,eps,minpts)
     cluster_object.assignments =cluster_object.assignments + 1
     cluster_counts = sum(cluster_object.counts)
     count_label_one = length(cluster_object.assignments) - cluster_counts
     unshift!(cluster_object.counts,count_label_one)
   end

   cluster_assignments = assignments(cluster_object)
   println("CLUSTER ASSIGNMENTS: ", cluster_assignments)

   #tiles_matrix_ordered_by_cluster = order_matrix(tiles_matrix_transpose, cluster_assignments)
   tiles_matrix_ordered_by_cluster = order_matrix(tiles_matrix, cluster_object)

   return (tiles_matrix_ordered_by_cluster, cluster_object,sil_vals)
end

function kmeans_cluster_tiles(tiles_matrix,num_clusters)
    println("kmeans ready")
    tiles_matrix_FLOAT = convert(Array{Float64,2},tiles_matrix)
    kmeans_cluster_object = kmeans(tiles_matrix_FLOAT, num_clusters)
    #cluster_assignments = kmeans_cluster_object.assignments
    #return cluster_assignments
    return kmeans_cluster_object
end

# THIS DOES LOTS MORE THAN THE ABOVE
function kmeans_cluster_tiles(tiles_matrix, num_repetitions, cluster_range)
   println("kmeans (ADVANCED) ready")
    tiles_matrix_FLOAT = convert(Array{Float64,2},tiles_matrix)
  current_mean_sil_val = 0
  current_kmeans_cluster_object = nothing
  #sil_dict = Dict{Int64, Array{Float64}}()
  k_array = Int64[]
  mean_sil_array=Float64[]
  rep_array = Int64[]
  cost_array=Float64[]
  for k in cluster_range
      for r =1:num_repetitions
           println("cluster k $k, repetition $r")
           kmeans_cluster_object = kmeans(tiles_matrix_FLOAT, k)
           distance_matrix = pairwise(SqEuclidean(),tiles_matrix)
           sil_val = silhouettes(kmeans_cluster_object,distance_matrix)
           mean_sil_val = mean(sil_val)
           if mean_sil_val > current_mean_sil_val
               current_kmeans_cluster_object = kmeans_cluster_object
               current_mean_sil_val = mean_sil_val
           end
           println("mean sil val: ", mean_sil_val)
           push!(k_array,k)
           push!(mean_sil_array, mean_sil_val)
           push!(rep_array, r)
           push!(cost_array,kmeans_cluster_object.totalcost)
           println("mean cost: ", kmeans_cluster_object.totalcost)
           println("current sil val: ", current_mean_sil_val )
      end
  end
  sil_df = DataFrame(k=k_array,r=rep_array,ave_sil=mean_sil_array,cost=cost_array)
  return (current_kmeans_cluster_object,sil_df)
end

function dbscan_hammings(tiles_matrix, eps=5,minpts=10)
    println("dbscan ready")
      distance_matrix = pairwise(Hamming(),tiles_matrix)
    ds = dbscan(distance_matrix,eps,minpts)
    return ds #.assignments
end

function aff_prop_cluster_tiles(tiles_matrix, distance=SqEuclidean())
    println("affinity propogation ready")
    #similarity_matrix = -pairwise(Hamming(),tiles_matrix)
    similarity_matrix = -pairwise(distance,tiles_matrix)
    af=affinityprop(similarity_matrix)
    return af
end

function cluster_tiles_plot()
end

function cluster_tiles_barplot()
end


text:
"
General Purpose stuff for tax simulation models
"
module TBComponents

using Printf
using IterableTables
using IteratorInterfaceExtensions
using TableTraits

# inequality stuff
export OutputDict, OutputDictArray
export DEFAULT_ATKINSON_ES, DEFAULT_ENTROPIES, DEFAULT_FGT_ALPHAS
export makegini, makepoverty, makeinequality, binify, adddecomposedtheil

# reweighting
export DistanceFunctionType
export doreweighting, dochisquarereweighting

# general tax routines

export WEEKS_PER_YEAR, weeklyise, annualise, delete_thresholds_up_to
export TaxResult, RateBands, IncomesDict
export calctaxdue, calc_indirect, IndirResult,*,times

# budget constraint stuff
export Point2DG, Point2D, BudgetConstraint, BCSettings, DEFAULT_SETTINGS
export makebc, pointstoarray, annotate_bc

# equivalence scales

export Equivalence_Scale_Type, oxford, modified_oecd, eq_square_root, mcclements, eq_per_capita
export Equiv_Person_Type, eq_head, eq_spouse_of_head, eq_other_adult, eq_dependent_child
export EQ_Person
export get_equivalence_scale

include( "common_types.jl" )
include( "poverty_inequality.jl" )
include( "reweighter.jl" )
include( "taxcalcs.jl" )
include( "piecewise_linear_generator.jl" )
include( "equivalence_scales.jl" )

end # module


text:
"""
    helpers.jl -- Feel free to modify anything you like in this file!

Provides the problem structure and some useful functions. In particular, students may be
interested in `count` - to help with writing their algorithm, and `main` - to test their
methods. Note that `starter_code/localtest.jl` provides an identical test to that found
in the autograder. `main`, on the other hand, returns other metrics that may be of
interest during development.
"""

# Statistics is the only allowed external package (Random is part of the standard library)
# If you would like to use any other (standard) packages, make sure to do so in `project1.jl`
# *NOT HERE*.
using Random
using Statistics

# A global counter that keeps track of how many times each function has been called.
# It may seem like a clever hack to edit this dictionary as part of your optimize
# method to get infinite function evaluations, but beware...
# you'll regret it when it goes through the autograder
const COUNTERS = Dict{String, Int}()

"""
    @counted

A function defined with this macro increments the global counter `COUNTERS`
each time it's called.

Example:
    @counted f(x) = 2x  #each time `f(x)` is called, we also have `COUNTERS["f"] += 1`
"""
macro counted(f)
    name = f.args[1].args[1]
    name_str = String(name)
    body = f.args[2]
    update_counter = quote
        if !haskey(COUNTERS, $name_str)
            COUNTERS[$name_str] = 0
        end
        COUNTERS[$name_str] += 1
    end
    insert!(body.args, 1, update_counter)
    return f
end

# simple.jl defines the 3 simple problems. It's included down here rather than at the
# top because it relies on @counted and COUNTERS to track the evaluation counts.
include("simple.jl")

const PROBS = Dict("simple1" => (f=rosenbrock, g=rosenbrock_gradient, x0=rosenbrock_init, n=20),
                   "simple2" => (f=himmelblau, g=himmelblau_gradient, x0=himmelblau_init, n=40),
                   "simple3" => (f=powell, g=powell_gradient, x0=powell_init, n=100))

"""
    count(f::Function)
    count(f::Function, g::Function)

Check how many times the function f has been called, or calculate f + 2g.
"""
Base.count(f::Function)              = get(COUNTERS, string(nameof(f)), 0)
Base.count(f::Function, g::Function) = count(f) + 2*count(g)

"""
    get_score(f, g, x_star_hat, n)

The score is computed as `f(x⋆)` using the potential optimum `x_star_hat`.
If `count(f, g) > n`, the score in necessarily `Inf` (didn't adhere to constraints).
Also returns the number of evaluations for use in `main` and `localtest`.
"""
function get_score(f, g, x_star_hat, n)
    num_evals = count(f, g)

    score = num_evals <= n ? f(x_star_hat) : Inf

    return num_evals, score
end

"""
    main(probname, repeat, opt_func)

Evaluates a problem given by `probname` `repeat` times using `opt_func`
as the optimization (pass in your `optimize`).

## Arguments:
    - `probname`: Name of optimization problem (e.g. "simple1")
    - `repeat`: Number of Monte Carlo evaluations
    - `opt_func`: Optimization algorithm
## Returns:
    - (`mean_score`, `max_evals`)
"""
function main(probname::String, repeat::Int, opt_func, seed = 42)
    scores = zeros(repeat)
    nevals = zeros(Int, repeat)

    prob = PROBS[probname]

    # Repeat the optimization with a different initialization
    for i in 1:repeat
        empty!(COUNTERS) # fresh eval-count each time
        Random.seed!(seed + i)
        x_star_hat = opt_func(prob.f, prob.g, prob.x0(), prob.n, probname)
        nevals[i], scores[i] = get_score(prob.f, prob.g, x_star_hat, prob.n)
    end

    return scores, nevals
end


text:
import Base.exp

abstract type Function end

mutable struct Variable
    data::Array
    grad::Union{Array,Nothing}
    creator::Union{Function,Nothing}
    Variable(data) = new(data, nothing, nothing)
end
function backward(v::Variable)::Nothing
    if isnothing(v.grad)
        v.grad = ones(size(v.data))
    end
    funcs::Array{Function} = [v.creator]
    while !isempty(funcs)
        f::Function = pop!(funcs) # No need to check for nothingness at the first time
        x::Variable, y::Variable = f.input, f.output
        x.grad = backward(f, y.grad)
        if !isnothing(x.creator)
            push!(funcs, x.creator)
        end
    end
end

function _function(f::Function, forward::Core.Function, inputs::Array{Variable})::Array{Variable}
    xs::Array{Array} = [x.data for x in inputs]
    ys::Tuple{Vararg{Array}} = forward(xs)
    outputs = [Variable(y) for y in ys]
    for output in outputs
        output.creator = f
    end
    f.inputs = inputs
    f.outputs = outputs
    return outputs
end

mutable struct Add <: Function
    inputs::Union{Array{Variable},Nothing}
    outputs::Union{Array{Variable},Nothing}
    Add() = new()
end
(f::Add)(inputs::Array{Variable})::Array{Variable} = _function(f, xs -> (xs[1] + xs[2],), inputs)