text:
module Json

import JSON3, HTTP, Reexport

Reexport.@reexport using Genie
Reexport.@reexport using Genie.Renderer

module JSONParser

import JSON3

const parse = JSON3.read

const json = JSON3.write

# ugly but necessary
JSON3.StructTypes.StructType(::T) where {T<:DataType} = JSON3.StructTypes.Struct()

end

using .JSONParser

const JSON = JSONParser
const JSON_FILE_EXT = ".json.jl"
const JSONString = String

export JSONString, json

function render(viewfile::Genie.Renderer.FilePath; context::Module = @__MODULE__, vars...) :: Function
  Genie.Renderer.registervars(; context = context, vars...)

  path = viewfile |> string
  partial = true

  f_name = Genie.Renderer.function_name(string(path, partial)) |> Symbol
  mod_name = Genie.Renderer.m_name(string(path, partial)) * ".jl"
  f_path = joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name)
  f_stale = Genie.Renderer.build_is_stale(path, f_path)

  if f_stale || ! isdefined(context, f_name)
    Genie.Renderer.build_module(
      Genie.Renderer.Html.to_julia(read(path, String), nothing, partial = partial, f_name = f_name),
      path,
      mod_name,
      output_path = false
    )

    Base.include(context, f_path)
  end

  () -> try
    getfield(context, f_name)() |> first
  catch ex
    if isa(ex, MethodError) && string(ex.f) == string(f_name)
      Base.invokelatest(getfield(context, f_name)) |> first
    else
      rethrow(ex)
    end
  end |> JSONParser.json
end

function render(data::Any; forceparse::Bool = false, context::Module = @__MODULE__) :: Function
  () -> JSONParser.json(data)
end

function Genie.Renderer.render(::Type{MIME"application/json"}, datafile::Genie.Renderer.FilePath; context::Module = @__MODULE__, vars...) :: Genie.Renderer.WebRenderable
  Genie.Renderer.WebRenderable(render(datafile; context = context, vars...), :json)
end

function Genie.Renderer.render(::Type{MIME"application/json"}, data::String; context::Module = @__MODULE__, vars...) :: Genie.Renderer.WebRenderable
  Genie.Renderer.WebRenderable(render(data; context = context, vars...), :json)
end

function Genie.Renderer.render(::Type{MIME"application/json"}, data::Any; context::Module = @__MODULE__, vars...) :: Genie.Renderer.WebRenderable
  Genie.Renderer.WebRenderable(render(data), :json)
end

### json API

function json(resource::Genie.Renderer.ResourcePath, action::Genie.Renderer.ResourcePath; context::Module = @__MODULE__,
              status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), vars...) :: Genie.Renderer.HTTP.Response
  json(Genie.Renderer.Path(joinpath(Genie.config.path_resources, string(resource), Renderer.VIEWS_FOLDER, string(action) * JSON_FILE_EXT));
        context = context, status = status, headers = headers, vars...)
end

function json(datafile::Genie.Renderer.FilePath; context::Module = @__MODULE__,
              status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), vars...) :: Genie.Renderer.HTTP.Response
  Genie.Renderer.WebRenderable(Genie.Renderer.render(MIME"application/json", datafile; context = context, vars...), :json, status, headers) |> Genie.Renderer.respond
end

function json(data::String; context::Module = @__MODULE__,
              status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), vars...) :: Genie.Renderer.HTTP.Response
  Genie.Renderer.WebRenderable(Genie.Renderer.render(MIME"application/json", data; context = context, vars...), :json, status, headers) |> Genie.Renderer.respond
end

function json(data::Any; status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders()) :: Genie.Renderer.HTTP.Response
  Genie.Renderer.WebRenderable(Genie.Renderer.render(MIME"application/json", data), :json, status, headers) |> Genie.Renderer.respond
end

### === ###
### EXCEPTIONS ###

function Genie.Router.error(error_message::String, ::Type{MIME"application/json"}, ::Val{500}; error_info::String = "") :: HTTP.Response
  json(Dict("error" => "500 Internal Error - $error_message", "info" => error_info), status = 500)
end

function Genie.Router.error(error_message::String, ::Type{MIME"application/json"}, ::Val{404}; error_info::String = "") :: HTTP.Response
  json(Dict("error" => "404 Not Found - $error_message", "info" => error_info), status = 404)
end

function Genie.Router.error(error_code::Int, error_message::String, ::Type{MIME"application/json"}; error_info::String = "") :: HTTP.Response
  json(Dict("error" => "$error_code Error - $error_message", "info" => error_info), status = error_code)
end

end


text:
# Consensus and fine splitting functions

"""
    denoise(seqs; <keyword arguments> )

                    rough_radius::Float64 = 0.01,
                    fine_radius::Float64 = 1.0,
                    reassign::Bool = true,
                    clust_multiplier::Float64 = 1.0,
                    min_devs::Int64 = 6,
                    initial_devs::Int64 = 20,
                    user_min_clust::Int64 = 5,
                    k::Int64 = 6,
                    distfunc = euclidean,
                    center = mean,
                    verbose::Int64 = 1,
                    cycle_lim = 20,
                    triangle::Bool = false,
                    degap_param::Bool = false,
                    Polish::Bool = true,
                    p_value::Float64 = 0.01

#Arguments
- `seqs`: array containing the sequences to be operated on
- `rough_radius::Float64`: distance between SNPs used in fine clustering
- `fine_radius::Float64`: distance between SNPs used in rough clustering
- `reassign::Bool = true`: indicates whether or not to reassign the rough singletons
back to consensus sequences
- `clust_multiplier::Float64 = 1.0`: minimum non-error cluster size
- `min_devs::Int64 = 6`: number of high standard deviation kmer to fine cluster on
- `initial_devs::Int64 = 20`: initial amount of high standard deviation kmer to test
- `user_min_clust::Int64 = 5`: minimum size of a cluster, competes with clust_multiplier
- `k::Int64 = 6`: length of kmer sequences calculated
- `distfunc = euclidean`: type of distance function calculated on the kmer sequences
- `center = mean`: sets the center of a centroid to the mean of the centroid
- `verbose::Int64 = 1`: when true, the status of the program will be printed to 'stdout'
- `cycle_lim = 20`: maximum number of cycles allowed in clustering call
- `triangle::Bool = false`: controls whether or not the triangle inequality is used in calculations
- `degap_param::Bool = false`: controls whether degap is used in consensus calculations
- `Polish::Bool = false`: if true, draws a consensus again from singletons after reassignment
- `p_value::Float64 = 0.01`: statistical threshold for splitting decisions

Performs clustering and sequence consensus.

Order of process:
- Rough clustering
- Fine clustering
- Consensus of clusters

Recommended `radius` for rough clustering: 0.01
Recommended `radius` for fine clustering: 1.0

"""
function denoise(seqs; rough_radius::Float64=0.01, fine_radius::Float64=1.0, reassign::Bool=true,
                clust_multiplier::Float64=1.0, min_devs::Int64=6, initial_devs::Int64=20, user_min_clust::Int64=5,
                k::Int64=6, distfunc=euclidean, center=mean, verbose::Int64=1, cycle_lim=3, triangle::Bool=false,
                degap_param::Bool=false, Polish::Bool=true, p_value::Float64=0.01)

    #transforms input sequences into vectors of kmer frequencies
    kmer_vecs = [kmer_count(seq, k) for seq in seqs];

    fine_indices =  kmer_split_clustering(kmer_vecs, seqs, rough_radius=rough_radius, fine_radius=fine_radius,
                        reassign=reassign, clust_multiplier=clust_multiplier, min_devs=min_devs, initial_devs=initial_devs,
                        user_min_clust=user_min_clust, k=k, distfunc=distfunc, center=center, verbose=verbose,
                        cycle_lim=cycle_lim, triangle=triangle, p_value=p_value)

    return kmer_split_consensus(kmer_vecs, seqs, fine_indices, verbose=verbose, distfunc=distfunc, reassign=reassign,
                            k=k, degap_param=degap_param, user_min_clust=user_min_clust, Polish=Polish)
end

"""
    kmer_split_clustering(kmer_vecs, seqs; rough_radius::Float64 = 0.01, fine_radius::Float64 = 1.0,  reassign::Bool = true,
                        clust_multiplier::Float64 = 1.0, min_devs::Int64 = 6, initial_devs::Int64 = 20,
                        user_min_clust::Int64 = 5, k::Int64 = 6, distfunc = euclidean, center = mean,
                        verbose::Int64 = 1, cycle_lim = 20, triangle::Bool = false, p_value::Float64 = 0.01)

Returns the fine clustering result from the given sequences.

"""
function kmer_split_clustering(kmer_vecs, seqs; rough_radius::Float64=0.01, fine_radius::Float64=1.0, reassign::Bool=true,
                        clust_multiplier::Float64=1.0, min_devs::Int64=6, initial_devs::Int64=20, user_min_clust::Int64=5,
                        k::Int64=6, distfunc=euclidean, center=mean, verbose::Int64=1, cycle_lim=20, triangle::Bool=true,
                        p_value::Float64=0.01)

    #clusters the indices of input sequences into groups originating from the same templates
    return fine_clustering(kmer_vecs, rough_radius=rough_radius, fine_radius=fine_radius, clust_multiplier=clust_multiplier,
                           min_devs=min_devs, initial_devs=initial_devs, user_min_clust=user_min_clust, distfunc=distfunc,
                           center=center, verbose=verbose, cycle_lim=cycle_lim, triangle=triangle, p_value=p_value)
end

"""
    kmer_split_consensus(kmer_vecs, seqs, fine_indices; verbose::Int64=1, distfunc=euclidean, reassign::Bool=true,
                       k::Int64=6, degap_param::Bool=false, user_min_clust::Int64=5, Polish::Bool=true)

Returns the consensus of the sequences after clustering.

"""
function kmer_split_consensus(kmer_vecs, seqs, fine_indices; verbose::Int64=1, distfunc=euclidean, reassign::Bool=true,
                            k::Int64=6, degap_param::Bool=false, user_min_clust::Int64=5, Polish::Bool=true)
    if verbose > 0
        println("Finding consensus sequences...")
    end
    #convert the clusters of indices into clusters of sequences
    seq_clusters = [seqs[fine_indices[ind]] for ind in 1:length(fine_indices)]

    #determine consensus sequences that spawned each of the sequence clusters
    consensus = [consensus_seq(clust, degap_param = degap_param) for clust in seq_clusters]
    consensus, sizes = FAD_clean(consensus, [length(indices) for indices in fine_indices])
    #adds back singletons from a preliminary clustering step in fine_clustering in order to obtain more
    #accurate frequencies
    if reassign
        if verbose > 0
            println("Reassigning singleton inputs to consensus sequences to improve accuracy of frequencies")
        end
        ave_dist_arr = Array{Float64, 1}(length(seqs))
        consensus_vecs =  [kmer_count(seq, k) for seq in consensus];
        num_C = length(consensus_vecs)
        num_vecs = length(kmer_vecs)
        C_arr = consensus_vecs
        cluster_indices = [[] for i in 1:num_C]
        for i in 1:num_vecs
            curr_input = kmer_vecs[i]
            dists = [distfunc(C_arr[j], curr_input) for j in 1:num_C]
            min_index = indmin(dists)
            push!(cluster_indices[min_index], i)
        end
        sizes = [length(cluster) for cluster in cluster_indices]
        arr = find(user_min_clust.<sizes) #consensus seqs with zero inputs disappear!
        if verbose > 0
            println("After reassigning, there are $(length(arr)) templates")
        end
        if Polish
            cluster_indices=cluster_indices[arr]
            sizes=sizes[arr]
            seq_clusters = [seqs[cluster_indices[ind]] for ind in 1:length(cluster_indices)]
            consensus = [consensus_seq(clust, degap_param = degap_param) for clust in seq_clusters]
            return consensus, sizes, cluster_indices
        else
            return consensus[arr], sizes[arr], cluster_indices[arr]
        end
    else
        return consensus, sizes, fine_indices
    end
end

function FAD_clean(seqs, sizes; alpha = 0.01, neigh_thresh = 1.0,method = 2,err_rate=0.02, phreds = nothing)
    expected_zero_errors=1.0
    if !(method in [1,2])
        warn("Please pick a method in 1,2,3")
        return [],[]
    end
    #Method list:
    #1: No statistical considerations. All smaller relatives get eaten.
    #2: Use Poisson test to decide to add not-too-small relatives back in.

#    counts = countmap(seqs);
    seq_keys = seqs;
seq_freqs_all = reverse(sort([(sizes[k],seq_keys[k],kmer_count(seq_keys[k],6)) for k in 1:length(seq_keys)]));
    seq_freqs = seq_freqs_all[[k[1]>=2 for k in seq_freqs_all]];

    current_set = [seq_freqs[1]];
    for i in 2:length(seq_freqs)
        dists = [corrected_kmer_dist_full(k[3],seq_freqs[i][3]) for k in current_set]
        neighbour_set = dists.<=neigh_thresh

        #if no neighbours within radius, add to current_set
        if sum(neighbour_set)==0
            push!(current_set,seq_freqs[i])
        else
            if method == 2
                inds = collect(1:length(current_set))[neighbour_set]
                best = indmax(current_set[inds])
                p_val = (1-cdf(Poisson((current_set[inds][best][1]/expected_zero_errors)*err_rate),seq_freqs[i][1]))*length(seq_freqs[i][2])
                if p_val < alpha
                    push!(current_set,seq_freqs[i])
                end
            end
        end
    end

    frequencies = zeros(length(current_set))
    for s in seq_freqs_all
        dists = [corrected_kmer_dist_full(k[3],s[3],k=6) for k in current_set]
        frequencies[indmin(dists)] += s[1]
    end
    return [k[2] for k in current_set],frequencies

end

"""
    get_ave_dists(fine_indices, original_seqs, consensus_seqs; distfunc=euclidean, center=median, k::Int64=6)

Finds the average distance using the declared distance function
"""
function get_ave_dists(fine_indices, original_seqs, consensus_seqs; distfunc=euclidean, center=median, k::Int64=6)
    kmer_vecs = [kmer_count(seq, k) for seq in original_seqs];
    num_seqs = length(consensus_seqs)
    ave_dist_arr = Array{Float64, 1}(num_seqs)
    consensus_vecs = [kmer_count(seq, k) for seq in consensus_seqs];
    for ind in 1:num_seqs
        arr = consensus_vecs[ind]
        clust = fine_indices[ind]
        ave_dist_arr[ind] = center([distfunc(arr, kmer_vecs[index]) for index in clust])
    end
    return ave_dist_arr
end

"""
    denoise_pipeline(infile::String; <Keyword Arguments> )

                         prefix = "",
                         outfile = "",
                         ref_file = "",
                         error_filter = 0.01,
                         homopoly_filter = true,
                         minlength = 20,
                         maxlength = 1000000,
                         trim_seqs = true,
                         align_ref_frames = true,
                         write_bad_seqs = false,
                         write_aminoacids = true,
                         minsize = 5,
                         cycle_lim = 10,
                         merge_differing_homopolymers = false,
                         cluster_filedump_path = "",
                         post_filedump = false,
                         rough_radius::Float64=0.01,
                         fine_radius::Float64=1.0,
                         reassign::Bool=true,
                         clust_multiplier::Float64=1.0,
                         min_devs::Int64=6,
                         initial_devs::Int64=20,
                         user_min_clust::Int64=5,
                         k::Int64=6,
                         distfunc=euclidean,
                         center=mean,
                         verbose::Int64=1,
                         triangle::Bool=false,
                         degap_param::Bool=false,
                         Polish::Bool=true,
                         p_value::Float64=0.01)

# Arguments
- `prefix = ""`: prefix for names of sequences in output file (ie in case of multiple time points, etc.).
- `outfile = ""`: default is "[infile name without extension]-clusters.fasta".
- `ref_file = ""`: panel of reference seqs -- if provided, used to orient and trim nearest sequences in input file.
- `error_filter = 0.01`: if the input is a ".fastq" file, any sequence with a site with error greater than this value is removed before clustering.
- `homopoly_filter = true`: if true, remove sequences with erroneous run-on homopolymer regions (or trim the sequences if this region is on an end) via hidden markov model inference.
- `minlength = 20`: removes sequences shorter than this prior to clustering.
- `maxlength = 1000000`: removes sequences longer than this prior to clustering.
- `trim_seqs = true`: if a reference panel is provided (`ref_file`) and this is true, locally trims consensus seqs to nearest reference sequence
- `align_ref_frames = true`: corrects codon reading frames of inferred consensus sequences.
- `write_bad_seqs = false`: write sequences with uncorrectable reading frames to file
- `write_aminoacids = true`: write amino acid translations of consensus sequences to file
- `minsize = 5`: throw away clusters/consensuses with fewer than this many reads.
- `cycle_lim = 10` : TODO
- `merge_differing_homopolymers = false`: aggregates clusters with consensus seqs that only differ by a single gap in a homopolymer region.
- `cluster_filedump_path = ""`: if nonempty, writes original sequences to output files determined by resulting cluster in given directory. This string may also include a file name prefix.
- `post_filedump`: writes inferred consensus seqs to file after clustering, before aligning and trimming.
- `rough_radius::Float64`: distance between SNPs used in fine clustering
- `fine_radius::Float64`: distance between SNPs used in rough clustering
- `reassign::Bool = true`: indicates whether or not to reassign the rough singletons
back to consensus sequences
- `clust_multiplier::Float64 = 1.0`: minimum non-error cluster size
- `min_devs::Int64 = 6`: number of high standard deviation kmer to fine cluster on
- `initial_devs::Int64 = 20`: initial amount of high standard deviation kmer to test
- `user_min_clust::Int64 = 5`: minimum size of a cluster, competes with clust_multiplier
- `k::Int64 = 6`: length of kmer sequences calculated
- `distfunc = euclidean`: type of distance function calculated on the kmer sequences
- `center = mean`: sets the center of a centroid to the mean of the centroid
- `verbose::Int64 = 1`: when true, the status of the program will be printed to 'stdout'
- `cycle_lim = 20`: maximum number of cycles allowed in clustering call
- `triangle::Bool = false`: controls whether or not the triangle inequality is used in calculations
- `degap_param::Bool = false`: controls whether degap is used in consensus calculations
- `Polish::Bool = false`: if true, draws a consensus again from singletons after reassignment
- `p_value::Float64 = 0.01`: statistical threshold for splitting decisions

A clustering pipeline for sequence files.
Order of processes:
- Error filter (if .fastq)
- Length filter
- Homopolymer filter
- Length filter again
- Orient (either to reference panel or inferred coarse cluster centroid)
- Rough clustering
- Fine clustering
- Consensus of clusters

Recommended `radius` for rough clustering: 0.01
Recommended `radius` for fine clustering: 1.0

`infile` and `ref_file` must be either .fasta or .fastq.
If file is .fastq, then removes sequences with greater than `error_filter`
error.

"""
function denoise_pipeline(infile::String;
                         prefix = "", outfile = "", ref_file = "",
                         error_filter = 0.01,
                         homopoly_filter = true,
                         minlength = 20,
                         maxlength = 1000000,
                         trim_seqs = true,
                         align_ref_frames = true,
                         write_bad_seqs = false,
                         write_aminoacids = true,
                         minsize = 5,
                         cycle_lim = 10,
                         merge_differing_homopolymers = false,
                         cluster_filedump_path = "",
                         post_filedump = false,
                         rough_radius::Float64=0.01,
                         fine_radius::Float64=1.0,
                         reassign::Bool=true,
                         clust_multiplier::Float64=1.0,
                         min_devs::Int64=6,
                         initial_devs::Int64=20,
                         user_min_clust::Int64=5,
                         k::Int64=6,
                         distfunc=euclidean,
                         center=mean,
                         verbose::Int64=1,
                         triangle::Bool=false,
                         degap_param::Bool=false,
                         Polish::Bool=true,
                         p_value::Float64=0.01)
    seqs = []
    if infile[end-5:end] == ".fasta"
        seqs, names = read_fasta_with_names_in_other_order(infile, seqtype=String)
        phreds = nothing
    elseif infile[end-5:end] == ".fastq"
        seqs, phreds, names = read_fastq(infile, seqtype=String)
        println("read $(length(seqs)) seqs")
        # maybe put this somewhere else?
        seqs, phreds, names = quality_filter(seqs, phreds, names, errorRate=error_filter)
        println("kept $(length(seqs)) seqs")
    else
        error("Unsupported file type: $(infile[end-5:end])")
    end
    if verbose > 0
        println(length(seqs), " total sequences read from ", infile)
    end

    if verbose > 0
        println("Filtering...")
    end

    minlength = max(minlength, 1)
    @time @sync seqs, phreds, names = length_filter(Array{String, 1}(seqs), phreds, names, minlength, maxlength)
    # TODO: write to specific filtered file option
    if verbose > 0
        println(length(seqs), " total sequences after length filter")
    end

    if homopoly_filter
        @time @sync seqs, phreds, names = homopolymer_filter(Array{String, 1}(seqs), phreds, names)
        # TODO: write to specific filtered file option
        if verbose > 0
            println(length(seqs), " total sequences after homopolymer filter")
        end

        @time @sync seqs, phreds, names = length_filter(Array{String, 1}(seqs), phreds, names, minlength, maxlength)
        # TODO: write to specific filtered file option
        if verbose > 0
            println(length(seqs), " total sequences after length filter")
        end
    end

    if length(ref_file) > 0
        @time @sync seqs, phreds, names = orient_to_refs_file(Array{String, 1}(seqs), phreds, names, ref_file)
    else
        # orient to coarse cluster centroid
        @time @sync seqs, phreds, names = orient_strands(Array{String, 1}(seqs), phreds, names, get_coarse_centroid(seqs))
    end
    # TODO: write to specific filtered file option
    if verbose > 0
        println(length(seqs), " total sequences after orienting")
    end

    if verbose > 0
        println("Clustering and Consensusing...")
    end

    #clusts, sizes, fine_indices = denoise(seqs, triangle=true, cycle_lim=cycle_lim, user_min_clust=10)
    clusts, sizes, fine_indices = denoise(seqs, rough_radius=rough_radius, fine_radius=fine_radius, reassign=reassign,
                clust_multiplier=clust_multiplier, min_devs=min_devs, initial_devs=initial_devs, user_min_clust=user_min_clust,
                k=k, distfunc=distfunc, center=center, verbose=1, cycle_lim=cycle_lim, triangle=triangle,
                degap_param=degap_param, Polish=Polish, p_value=p_value)

    if length(clusts) == 0
        warn("No consensus sequences found; try lowering minimum cluster size")
    end

    if length(prefix) > 0
        prefix = prefix * "_"
    end
    if outfile == ""
        outfile = infile[1:end-6] * "-clusters.fasta"
    elseif outfile[end-5:end] != ".fasta"
        println("Output filetype must be fasta, writing to default filename")
        outfile = infile[1:end-6] * "-clusters.fasta"
    end

    if verbose > 0
        println("Trimming and aligning...")
    end
    if trim_seqs && length(ref_file) > 0
        if post_filedump
            names = [prefix * "seq" * string(i) * "_" * string(sizes[i]) for i in 1:length(sizes)]
            write_fasta(outfile[1:end-6] * "-untrimmed.fasta", clusts, names=names)
        end
        clusts = trim_to_refs(Array{String, 1}(clusts), ref_file)
        tmp = merge_consensuses([[clusts[i], sizes[i]] for i in 1:length(sizes)])
        clusts, sizes = [t[1] for t in tmp], [t[2] for t in tmp]
    end

    if align_ref_frames
        if post_filedump
            names = [prefix * "seq" * string(i) * "_" * string(sizes[i]) for i in 1:length(sizes)]
            write_fasta(outfile[1:end-6] * "-uncorrected.fasta", clusts, names=names)
        end
        goodclsts, badclsts = align_reference_frames([clusts, sizes])
        clusts, sizes = goodclsts
        badclusts, badsizes = badclsts
    end

    names = [prefix * "seq" * string(i) * "_" * string(sizes[i]) for i in 1:length(sizes)]
    write_fasta(outfile, clusts, names=names)
    if write_aminoacids
        write_fasta(outfile[1:end-6] * "-aminoacids.fasta", [translate_to_aa(s[1:end-end%3]) for s in clusts], names=names)
    end
    if align_ref_frames && write_bad_seqs
        badnames = ["seq" * string(i) * "_" * string(badsizes[i]) for i in 1:length(badsizes)]
        write_fasta(outfile[1:end-6] * "-badseqs.fasta", badclusts, names=names)
    end

    if verbose > 0
        numseqs = length(sizes) > 0 ? sum(sizes) : 0
        medsize = length(sizes) > 0 ? median(sizes) : 0
        println("\nConsensus sequences written to ", outfile)
        println("Total sequences included: ", numseqs)
        println("Median cluster size: ", medsize)
        println("Number of clusters/consensus seqs: ", length(sizes))
    end
end

simple_clusterpipeline = denoise_pipeline

#==
    OLD PIPELINE
==#

"""
    cluster_and_consensus(seqs::Array{String}, radii; <keyword arguments>)

Cluster sequences, compute the consensus of each, and locally refine consensus sequences to corresponding reads.
Returns consensus sequences and cluster sizes for each.
`radii` may be a single float64 or an array of float64.

- `k = 6`: size of kmers computed in kmer vector for approximate dirichlet process clustering.
- `distfunc = corrected_kmer_dist(k)'`: distance function used in clustering.
- `center = mean`: method of computing cluster centroids.
- `shift = 7`: window size use to locally refine reference sequences.
- `minsize::Int = 10`: minimum size of cluster to return.
- `maxsize::Int = 100`: maximum number of sequences used to refine consensus sequence of cluster.
- `cycle_lims = fill(20, length(radii))'`: number of iterations of clustering.
- `filedump = false`: write original sequences to files separated by cluster.
- `filedump_path = ""`: path to write clusters of seqs.
- `names = nothing`: names of sequences to write to file. If not given, names are "seq1", etc.
- `phreds = nothing`: phreds of sequences to write to file. If given, writes ".fastq", else writes ".fasta".
- `verbose::Bool = false`: progress updates printed to std_out.
"""
function cluster_and_consensus(seqs::Array{String}, radii;
                                    k = 6,
                                    distfunc = corrected_kmer_dist(k),
                                    center = mean,
                                    shift = 7,
                                    minsize::Int = 10,
                                    maxsize::Int = 100,
                                    cycle_lims = fill(20, length(radii)),
                                    filedump = false,
                                    filedump_path = "",
                                    names = nothing,
                                    phreds = nothing,
                                    verbose::Bool = false)

    if verbose
        println("Computing kmer vectors...")
    end
    kmer_vecs = [kmer_count(seq, k) for seq in seqs]
    if verbose
        println("Clustering...")
    end
    μs, sizes, cluster_indices, centroid_indices = dp_centers(kmer_vecs, radii,
                                                              distfunc=distfunc,
                                                              center=center,
                                                              cycle_lims=cycle_lims,
                                                              verbose=verbose)
    filt = sizes .>= minsize
    final_sizes = sizes[filt]
    if verbose
        println("Keeping $(sum(filt)) clusters")
    end
    filt_indices = cluster_indices[filt]
    consensus_seqs = Array{String}(length(filt_indices))
    for i in 1:length(filt_indices)
        if verbose
            println("Computing consensus for cluster $i with size $(sizes[filt][i])")
        end

        consensus_seqs[i] = refine_ref(seqs[centroid_indices[filt][i]],
                                       sample(seqs[filt_indices[i]],
                                              min(maxsize, length(filt_indices[i])),
                                              replace=false),
                                       shift=shift)
        if filedump
            # write fasta for each split cluster - includes < minsize clusters
            ss = seqs[filt_indices[i]]
            nms = names==nothing ? ["seq$(iii)" for iii in 1:length(ss)] : names[filt_indices[i]]
            pds = phreds==nothing ? nothing : phreds[filt_indices[i]]
            if pds != nothing
                write_fastq(filedump_path * "clust-$(i).fastq", ss, pds, names=nms)
            else
                write_fasta(filedump_path * "clust-$(i).fasta", ss, names=nms)
            end
        end
    end
    # TODO: merge identical consensuses? (cant use merge_consensuses)
    return consensus_seqs, final_sizes
end

"""
    fine_cluster_and_consensus(seqs::Array{String}, radii; <keyword arguments>).

                                    polylen = 3,
                                    minsize::Int = 10,
                                    maxsize::Int = 100,
                                    cycle_lims = fill(20, length(radii)),
                                    merge_differing_homopolymers = false,
                                    filedump = false,
                                    filedump_path = "",
                                    names = nothing,
                                    phreds = nothing,
                                    verbose::Bool = false)

Cluster sequences, compute the consensus of each, and split into sub-clusters on above-threshold-discrepancies.
Returns consensus sequences and cluster sizes for each sub-cluster.
`radii` may be a single float64 or an array of float64.

- `k = 6`: size of kmers computed in kmer vector for approximate dirichlet process clustering.
- `distfunc = corrected_kmer_dist(k)'`: distance function used in clustering.
- `center = mean`: method of computing cluster centroids.
- `shift = 7`: window size use to locally refine reference sequences.
- `polylen`: min length of homopolymer region to be determined as a homopolymer region (see `peak_split`).
- `minsize::Int = 10`: minimum size of cluster to return.
- `cycle_lims = fill(20, length(radii))'`: number of iterations of clustering.
- `filedump = false`: write original sequences to files separated by cluster.
- `filedump_path = ""`: path to write clusters of seqs.
- `merge_differing_homopolymers = false`: if inferred consensus seqs only differ by a single gap in a homopolymer region, merge their clusters. (Note: this occurs after filedump)
- `names = nothing`: names of sequences to write to file. If not given, names are "seq1", etc.
- `phreds = nothing`: phreds of sequences to write to file. If given, writes ".fastq", else writes ".fasta".
- `verbose::Bool = false`: progress updates printed to std_out.
"""
function fine_cluster_and_consensus(seqs::Array{String}, radii;
                                    k = 6,
                                    distfunc = corrected_kmer_dist(k),
                                    center = mean,
                                    shift = 7,
                                    polylen = 3,
                                    minsize::Int = 10,
                                    cycle_lims = fill(20, length(radii)),
                                    merge_differing_homopolymers = false,
                                    filedump = false,
                                    filedump_path = "",
                                    names = nothing,
                                    phreds = nothing,
                                    verbose::Bool = false)

    if verbose
        println("Computing kmer vectors...")
    end
    kmer_vecs = [kmer_count(seq, k) for seq in seqs]
    if verbose
        println("Clustering...")
    end
    μs, sizes, cluster_indices, centroid_indices = dp_centers(kmer_vecs, radii,
                                                              distfunc=distfunc,
                                                              center=center,
                                                              cycle_lims=cycle_lims,
                                                              verbose=verbose)
    filt = sizes .>= minsize
    final_sizes = sizes[filt]
    if verbose
        println("Keeping $(sum(filt)) clusters")
    end
    filt_indices = cluster_indices[filt]

    consensus_seqs = []
    for i in 1:length(filt_indices)
        if verbose
            println("Computing consensus for cluster $i with size $(sizes[filt][i])")
        end
        # tmp array with only this cluster
        tmp_consensuses = []
        reads = []
        nms = names==nothing ? ["seq$(iii)" for iii in 1:length(filt_indices[i])] : names[filt_indices[i]]
        pds = phreds==nothing ? nothing : phreds[filt_indices[i]]

        recursive_split(seqs[filt_indices[i]],
                        tmp_consensuses,
                        centroid=seqs[centroid_indices[filt][i]],
                        shift=shift,
                        minreads=minsize,
                        polylen=polylen,
                        k=k,
                        reads=reads,
                        phreds=pds,
                        names=nms)
        if filedump
            for j in 1:length(reads)
                # check if phreds exist
                if reads[j][2] != nothing
                    write_fastq(filedump_path * "-clust-$(i)-split-$(j).fastq", reads[j][1], reads[j][2], names=reads[j][3])
                else
                    write_fasta(filedump_path * "-clust-$(i)-split-$(j).fasta", reads[j][1], names=reads[j][3])
                end
            end
        end
        # add these new consensuses to total list
        append!(consensus_seqs, tmp_consensuses)
    end
    if merge_differing_homopolymers
        consensus_seqs = merge_diff_homopolymers(consensus_seqs, polylen=polylen)
    else
        consensus_seqs = merge_consensuses(consensus_seqs)
    end
    filter!(clust->clust[2] >= minsize, consensus_seqs)
    return [c[1] for c in consensus_seqs], [c[2] for c in consensus_seqs]
end

#-----Maybe move these to different files------

"""
    trim_to_refs(seqs::Array{String, 1}, refspath::String)

Trims each sequence in `seqs` to nearest reference in `refspath` file after local
alignment. Distance to references determined by amino acid similarity (kmer dot product).
`refspath` may be either .fasta or .fastq.
"""
function trim_to_refs(seqs::Array{String, 1}, refspath::String)
    if length(refspath) < 5
        return seqs
    end
    if refspath[end-5:end] == ".fasta"
        refs = read_fasta_with_names_in_other_order(refspath, seqtype=String)
    elseif refspath[end-5:end] == ".fastq"
        refs = read_fastq(refspath, seqtype=String)
    else
        error("Unsupported file type: $(infile[end-5:end])")
    end
    refs = refs[1]
    return trim_to_refs(seqs, refs)
end

"""
    trim_to_refs(seqs::Array{String, 1}, refs::Array{String, 1}; kmersize = 6)

Trims each sequence in `seqs` to nearest reference in `refs` array after local
alignment. Distance determined by amino acid similarity (kmer dot product).
"""
function trim_to_refs(seqs::Array{String, 1}, refs::Array{String, 1}; kmersize = 6)
    if length(refs) == 0
        println("Not trimming: no reference sequences")
        return seqs
    end
    ref_kmers = [sparse_aa_kmer_count(rf, kmersize) for rf in refs]
    trimmed = seqs[:]
    function trim(sq::String)
        best_match = 0
        best_match_i = 0
        vec = sparse_aa_kmer_count(sq, kmersize)
        # compare to every reference, update min dists
        matches = [dot(vec, rfk) for rfk in ref_kmers]
        best_match_i = indmax(matches)
        if best_match_i > 0
            ali = loc_kmer_seeded_align(refs[best_match_i], sq)
            return degap(ali[2])
        else
            # this probably happens because it does not share any kmers with any references
            # fix by shortening kmersize?
            # make flag to throw these away?
            return sq
        end
    end
    trimmed = map(trim, trimmed)
    return trimmed
end

ext: jl


text:
# The basics:
# 1) If the type of the extracted data is not defined, then it is String.
# 2) If not provided, the custom parser is `identity` for String, and
#    `parse` for any type `<: Number`.
# 3) If the `parse` throws, the extractor does not try to be smart.
# 4) In the case `extraction` returns NoDataFound, and the default is
#    NoDefault, then throw an error; otherwise return the default.
# 5) If no default is given, then NoDefault is used.

abstract struct AbstractExtractor{T} end
struct GenericExtractor{T} <: AbstractExtractor{T}
	extract
	default :: Union{T, NoDefault}
end
struct RegexExtractor{T} <: AbstractExtractor{T}
	regex :: Regex
	default :: Union{T, NoDefault}
end
function default(e :: GenericExtractor{T}) :: Union{T, NoDefault}
	return e.default
end
function default(e :: RegexExtractor{T}) :: Union{T, NoDefault}
	return e.default
end
function raw_extractor(e, ::Type{T}, default :: Union{T, NoDefault})

end


text:
import Generate

fhir_versions = [:R4]

for fhir_version in fhir_versions
    url = "https://www.hl7.org/fhir/$(fhir_version)/definitions.json.zip"
    json_definitions = Generate.download_fhir_json_schema(url);
    schema_string = json_definitions["fhir.schema.json.zip"]["fhir.schema.json"];
    output_file = joinpath(dirname(dirname(@__DIR__)), "src", "autogenerated", "autogenerated-$(fhir_version).jl")
    Generate.output_fhir_types(; schema_string = schema_string, output_file = output_file)
end


text:
# ---
# title: 1576. Replace All ?'s to Avoid Consecutive Repeating Characters
# id: problem1576
# author: Indigo
# date: 2021-01-27
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/replace-all-s-to-avoid-consecutive-repeating-characters/description/>
# hidden: true
# ---
#
# Given a string `s` containing only lower case English letters and the '?'
# character, convert **all** the '?' characters into lower case letters such
# that the final string does not contain any **consecutive repeating
# **characters. You **cannot** modify the non '?' characters.
#
# It is **guaranteed** that there are no consecutive repeating characters in the
# given string **except** for '?'.
#
# Return the final string after all the conversions (possibly zero) have been
# made. If there is more than one solution, return any of them. It can be shown
# that an answer is always possible with the given constraints.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "?zs"
#     Output: "azs"
#     **Explanation** : There are 25 solutions for this problem. From "azs" to "yzs", all are valid. Only "z" is an invalid modification as the string will consist of consecutive repeating characters in "zzs".
#
# **Example 2:**
#
#
#
#     Input: s = "ubv?w"
#     Output: "ubvaw"
#     **Explanation** : There are 24 solutions for this problem. Only "v" and "w" are invalid modifications as the strings will consist of consecutive repeating characters in "ubvvw" and "ubvww".
#
#
# **Example 3:**
#
#
#
#     Input: s = "j?qg??b"
#     Output: "jaqgacb"
#
#
# **Example 4:**
#
#
#
#     Input: s = "??yw?ipkj?"
#     Output: "acywaipkja"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 100`
#   * `s` contains only lower case English letters and `'?'`.
#
#
## @lc code=start
using LeetCode

function modify_string(s::String)
    char_vec = collect(s)
    for i in 1:length(char_vec)
        if char_vec[i] == '?'
            char_vec[i] = 'a'
            (i > 1 && char_vec[i - 1] == 'a') && (char_vec[i] += 1)
            (i < length(char_vec) && char_vec[i + 1] == char_vec[i]) && (char_vec[i] += 1)
            (i > 1 && char_vec[i - 1] == char_vec[i]) && (char_vec[i] += 1)
        end
    end
    return join(char_vec)
end
## @lc code=end


text:
@testset "Import MAT model" begin
    filepath = joinpath("data", "agora-model.mat")

    download_data_file(
        "https://www.vmh.life/files/reconstructions/AGORA/1.03/reconstructions/mat/Mycoplasma_hominis_ATCC_23114.mat",
        filepath,
        "03362073aa917f0691a0c896948f6e8eebe47f02dcbe0c3f00275fa87396e220",
    )

    cp = load_model(CoreModel, filepath)
    @test cp isa CoreModel
    @test size(cp.S) == (475, 496)

    cp = load_model(CoreModel, joinpath("data", "toyModel1.mat"))
    @test size(cp.S) == (6, 7)

    cp = load_model(CoreModel, joinpath("data", "toyModel2.mat"))
    @test size(cp.S) == (6, 7)

    cp = load_model(CoreModel, joinpath("data", "toyModel3.mat"))
    @test size(cp.S) == (9, 12)
end

@testset "Save MAT model" begin
    filepath = joinpath("data", "toyModel1.mat")
    loaded = load_model(CoreModel, filepath)
    save_model(loaded, "test_model.mat")
    wrote = load_model(CoreModel, "test_model.mat")
    @test wrote isa CoreModel
    @test isequal(wrote, loaded)
end


text:
module JuliaMBD

export ConstantBlock, GainBlock
export Line, Var
export connect
export _toexpr

include("_block.jl")

end

max_stars_repo_name: JuliaReliab/JuliaMBD

max_issues_repo_name: JuliaReliab/JuliaMBD

max_forks_repo_name: JuliaReliab/JuliaMBD


text:
# How a regular module may be written

module WebCrawler

using Dates

# public interface
export Target
export add_site!, crawl_sites!, current_sites, reset_crawler!

# target sites
Base.@kwdef mutable struct Target
    url::String
    finished::Bool = false
    finish_time::Union{DateTime,Nothing} = nothing
end

# tracks current set of targets
const sites = Target[]

function add_site!(site::Target)
    push!(sites, site)
end

function crawl_sites!()
    for s in sites
        index_site!(s)
    end
end

function current_sites()
    copy(sites)
end

function index_site!(site::Target)
    site.finished = true
    site.finish_time = now()
    println("Site $(site.url) crawled.")
end

function reset_crawler!()
    empty!(sites)
end

end # module

using Main.WebCrawler
add_site!(Target(url = "http://cnn.com"));
add_site!(Target(url = "http://yahoo.com"));
current_sites()

crawl_sites!()
current_sites()

#=
julia> using Main.WebCrawler

julia> add_site!(Target(url = "http://cnn.com"));

julia> add_site!(Target(url = "http://yahoo.com"));

julia> current_sites()
2-element Array{Target,1}:
 Target("http://cnn.com", false, nothing)
 Target("http://yahoo.com", false, nothing)

julia> crawl_sites!()
Site http://cnn.com crawled.
Site http://yahoo.com crawled.

julia> current_sites()
2-element Array{Target,1}:
 Target("http://cnn.com", true, 2019-09-23T22:17:51.222)
 Target("http://yahoo.com", true, 2019-09-23T22:17:51.226)

=#

# ---------------------------------------------------------------------------
# Using closure to hide non-public fields/functions

module WebCrawler

using Dates

# public interface
export Target
export add_site!, crawl_sites!, current_sites, reset_crawler!

# target sites
Base.@kwdef mutable struct Target
    url::String
    finished::Bool = false
    finish_time::Union{DateTime,Nothing} = nothing
end

let sites = Target[]

    global function add_site!(site::Target)
        push!(sites, site)
    end

    global function crawl_sites!()
        for s in sites
            index_site!(s)
        end
    end

    global function current_sites()
        copy(sites)
    end

    function index_site!(site::Target)
        site.finished = true
        site.finish_time = now()
        println("Site $(site.url) crawled.")
    end

    global function reset_crawler!()
        empty!(sites)
    end

end # let

end # module

# test access
using Main.WebCrawler
add_site!(Target(url = "http://cnn.com"));
add_site!(Target(url = "http://yahoo.com"));
current_sites()

crawl_sites!()
current_sites()

#=
julia> using Main.WebCrawler

julia> add_site!(Target(url = "http://cnn.com"));

julia> add_site!(Target(url = "http://yahoo.com"));

julia> current_sites()
2-element Array{Target,1}:
 Target("http://cnn.com", false, nothing)
 Target("http://yahoo.com", false, nothing)

julia> crawl_sites!()
Site http://cnn.com crawled.
Site http://yahoo.com crawled.

julia> current_sites()
2-element Array{Target,1}:
 Target("http://cnn.com", true, 2019-09-29T11:02:51.034)
 Target("http://yahoo.com", true, 2019-09-29T11:02:51.037)
=#

# test hidden variable/functions
Main.WebCrawler.index_site!()
Main.WebCrawler.sites

#=
julia> Main.WebCrawler.index_site!()
ERROR: UndefVarError: index_site! not defined

julia> Main.WebCrawler.sites
ERROR: UndefVarError: sites not defined
=#


text:
map_str = ""
passports = []
open("day04/day04.in") do io
    global map_str, passports
    for line in readlines(io)
        line = line * " "
        if cmp(line, " ") == 0
            push!(passports, map_str)
            map_str = ""
        else
            map_str = map_str * line
        end
    end
end

pass_map = Dict(
    "byr" => 2, "iyr" => 3, "eyr" => 5, "hgt" => 7, "hcl" => 11, "ecl" => 13,
    "pid" => 17, "cid" => 19
    )
pass_valid = Int64(0)
for passport in passports
    global pass_valid
    checkmult = Int64(1)
    validity = true
    for item in split(passport, ' ')
        if occursin(":", item)
            item_type, item_value = split(item, ":")
            if length(item) > 2
                checkmult *= pass_map[item_type]
            end
            if item_type == "byr"
                item_value = parse(Int64, item_value)
                validity &= 1920 <= item_value <= 2002
            elseif item_type == "iyr"
                item_value = parse(Int64, item_value)
                validity &= 2010 <= item_value <= 2020
            elseif item_type == "eyr"
                item_value = parse(Int64, item_value)
                validity &= 2020 <= item_value <= 2030
            elseif item_type == "hgt"
                if occursin("in", item_value)
                    item_value = split(item_value, "in")[1]
                    item_value = parse(Int64, item_value)
                    validity &= 59 <= item_value <= 76
                elseif occursin("cm", item_value)
                    item_value = split(item_value, "cm")[1]
                    item_value = parse(Int64, item_value)
                    validity &= 150 <= item_value <= 193
                else
                    validity = false
                end
            elseif item_type == "hcl"
                validity &= item_value[1] == '#' && length(item_value) == 7
                for ch in item_value[2:end]
                    validity &= '0' <= ch <= '9' || 'a' <= ch <= 'f'
                end
            elseif item_type == "ecl"
                validity &= item_value in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
            elseif item_type == "pid"
                validity &= length(item_value) == 9
                for ch in item_value[1:end]
                    validity &= '0' <= ch <= '9'
                end
            end
        end
    end
    if (checkmult == 9699690 || checkmult == 510510 ) && validity
        pass_valid += 1
    end
end
println(pass_valid)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
mutable struct AbstractParametersHSPARK{IT,DT,TT,D,S,R,P,VT,FT,UT,GT,ϕT,tabType} <: Parameters{DT,TT}
    f_v::VT
    f_f::FT
    f_u::UT
    f_g::GT
    f_ϕ::ϕT

    Δt::TT

    tab::tabType

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function AbstractParametersHSPARK{IT,DT,D,S,R,P}(f_v::VT, f_f::FT, f_u::UT, f_g::GT, f_ϕ::ϕT, Δt::TT, tableau::tabType) where {IT,DT,TT,D,S,R,P,VT,FT,UT,GT,ϕT,tabType}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new{IT,DT,TT,D,S,R,P,VT,FT,UT,GT,ϕT,tabType}(f_v, f_f, f_u, f_g, f_ϕ, Δt, tableau, zero(TT), q, p, λ)
    end
end

function update_params!(params::AbstractParametersHSPARK, sol::AtomicSolutionPDAE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
    params.λ .= sol.λ
end

equation(int::AbstractIntegratorHSPARK) = int.equation
timestep(int::AbstractIntegratorHSPARK) = int.params.Δt
tableau(int::AbstractIntegratorHSPARK) = int.tableau
pstages(int::AbstractIntegratorHSPARK) = int.tableau.r

function initialize!(int::AbstractIntegratorHSPARK, sol::AtomicSolutionPDAE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::AbstractIntegratorHSPARK, sol::AtomicSolutionPDAE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in eachdim(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in eachdim(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = 0
        end
    end

    if isdefined(tableau(int), :λ) && tableau(int).λ.c[1] == 0
        for k in eachdim(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(k-1)+3] = sol.λ[k]
        end
    end
end


text:
using PooksoftAlphaVantageDataStore

# ok, so lets build by user object -
my_current_dir = pwd()   # where am I?
path_to_config_file = my_current_dir*"/test/configuration/Configuration.json"

# build the api user model -
user_model = build_api_user_model(path_to_config_file) |> checkresult

# make an API call -
stock_symbol = "aapl"
data_type = :json
outputsize = :compact
df = execute_sts_adjusted_daily_api_call(user_model, stock_symbol; data_type = data_type, outputsize = outputsize, logger=nothing) |> checkresult


text:
using FileIO;
using DelimitedFiles;
using ScikitLearn;
@sk_import svm: SVC
using Random
using Random:seed!

function trainClassSVM(trainingInputs::Array{Float64,2}, trainingTargets::Array{Bool,2},
    testInputs::Array{Float64,2}, testTargets::Array{Bool,2};
    kernel::String="rbf", kernelDegree::Int64=3, kernelGamma::Int64=2, C::Int64=1)

    # Comprobamos que el numero de filas (numero de patrones) coincide tanto en entrenamiento como en test
    @assert(size(trainingInputs,1)==size(trainingTargets,1));
    @assert(size(testInputs,1)==size(testTargets,1));
    # Comprobamos que el numero de columnas coincide en los grupos de entrenamiento y test
    @assert(size(trainingInputs,2)==size(testInputs,2));
    @assert(size(trainingTargets,2)==size(testTargets,2));

    # Creamos nuestra máquina de soporte vectorial
    model = SVC(kernel=kernel, degree=kernelDegree, gamma=kernelGamma, C=C);
    # Entrenamos el modelo con el conjunto de entrenamiento
    model = fit!(model, trainingInputs, trainingTargets);
    # Pasamos el conjunto de test
    testOutputs = predict(model, testInputs);

    return model;
end;


text:
module NotSoSimpleMultigrid

# dependencies
using SimpleMultigrid, StaticArrays, LinearAlgebra, SparseArrays, Statistics, Printf

import SimpleMultigrid: stencil2mat, TransferKind, R₁, P₁, Grid, MultigridCycle, V, W, F, Smoother, zero_x, MultigridIterable, μ_cycle!, F_cycle!, smooth!, residu, factor_twos, coarse_grid_correction!

# export statements
export laplace2d, elliptic2d, anisotropic2d, rotated_anisotropic2d, laplace3d, elliptic3d, anisotropic3d
export Injection, FullWeighting, Cubic
export GaussSeidel, Jacobi
export coarsen
export MultigridIterable, V, W, F, V_cycle, W_cycle, F_cycle, \

# include source files
include("compose_matrix.jl")
include("grid_transfer_operators.jl")
include("grids.jl")
include("multigrid.jl")

end # module


text:
module V4

using MSTM_jll: mpiexec, mstm
using DataFrames
using DocStringExtensions
using Printf
using Scanf
using UUIDs
using ...STMMRunner
export run_mstm

abstract type MSTMOutput end

struct NearField
    spheres::DataFrame
    field::DataFrame
end

struct MSTMFixedOutput <: MSTMOutput
    q_ext_unpolarized::Float64
    q_abs_unpolarized::Float64
    q_scat_unpolarized::Float64
    q_ext_parallel::Float64
    q_abs_parallel::Float64
    q_scat_parallel::Float64
    q_ext_perpendicular::Float64
    q_abs_perpendicular::Float64
    q_scat_perpendicular::Float64
    q_scat_down_unpolarized::Float64
    q_scat_up_unpolarized::Float64
    q_scat_down_parallel::Float64
    q_scat_up_parallel::Float64
    q_scat_down_perpendicular::Float64
    q_scat_up_perpendicular::Float64
    q_ext_down_unpolarized::Float64
    q_ext_up_unpolarized::Float64
    q_ext_down_parallel::Float64
    q_ext_up_parallel::Float64
    q_ext_down_perpendicular::Float64
    q_ext_up_perpendicular::Float64
    q_scat_waveguide_unpolarized::Float64
    q_scat_waveguide_parallel::Float64
    q_scat_waveguide_perpendicular::Float64
    scattering_matrix::DataFrame
    near_field::Union{NearField,Nothing}
end

struct MSTMLatticeOutput <: MSTMOutput
    ref_unpolarized::Float64
    abs_unpolarized::Float64
    scat_unpolarized::Float64
    ref_parallel::Float64
    abs_parallel::Float64
    scat_parallel::Float64
    ref_perpendicular::Float64
    abs_perpendicular::Float64
    scat_perpendicular::Float64
    scattering_matrix::DataFrame
end

struct MSTMRandomOutput <: MSTMOutput
    q_ext::Float64
    q_abs::Float64
    q_scat::Float64
    scattering_matrix::DataFrame
end

struct MSTMMonteCarloOutput <: MSTMOutput
    q_ext::Float64
    q_abs::Float64
    q_scat::Float64
    q_scat_down_unpolarized::Float64
    q_scat_up_unpolarized::Float64
    q_scat_down_parallel::Float64
    q_scat_up_parallel::Float64
    q_scat_down_perpendicular::Float64
    q_scat_up_perpendicular::Float64
    scattering_matrix_azimuthal_average::DataFrame
    scattering_matrix_diffuse::DataFrame
    scattering_matrix_expansion_coefficients::DataFrame
end

"""
Use the given configuration to run MSTM v4.

$(SIGNATURES)

- If `keep = true`, the working directory will not be removed after the run.
"""
function run_mstm(cfg::STMMConfig; keep::Bool = false)
    current_dir = pwd()

    id = string(uuid1())
    if cfg.working_directory == ""
        dir = joinpath(pwd(), "mstm_tmp_$id")
    else
        dir = cfg.working_directory
    end
    results = nothing
    proc = nothing

    try
        if !ispath(dir)
            @debug "[Run MSTM] Creating workding directory"
            mkdir(dir)
        end

        @debug "[Run MSTM] Entering working directory"
        cd(dir)

        @debug "[Run MSTM] Writing MSTM input file"
        write_input(cfg)

        # We cannot simply do `$(mpiexec()) ... $(mstm())` here due to the limit of Cmd interpolation.
        # See https://github.com/JuliaLang/julia/issues/39282 for more details.
        @debug "[Run MSTM] Running MSTM..."
        proc = open(
            setenv(`$(mpiexec().exec[1]) -n $(cfg.number_processors) $(mstm().exec[1])`, mstm().env),
            stdout;
            write = true
        )
        wait(proc)

        @debug "[Run MSTM] Collecting MSTM output"
        results = collect_output(cfg)
    catch e
        @error "[Run MSTM] $e"
        @error "[Run MSTM] Runtime error occurred"
    finally
        @debug "[Run MSTM] Leaving working directory"
        cd(current_dir)

        if !isnothing(proc)
            @debug "[Run MSTM] Killing MSTM process"
            kill(proc)
            sleep(0.2)
            if process_running(proc)
                @warn "Failed to kill, please retry manually."
            end
        end

        if !keep && occursin(String(id), dir)
            @assert ispath(dir)
            @debug "[Run MSTM] Removing temporary directory"
            rm(dir; force = true, recursive = true)
        end

        return results
    end
end

function format_sphere(sphere::Sphere)::String
    radius_and_center = @sprintf "%.4e %.4e %.4e %.4e " sphere.x sphere.y sphere.z sphere.r

    if !isempty(sphere.t_matrix)
        return radius_and_center * sphere.t_matrix
    else
        mₗ = sphere.m / (1.0 - sphere.β * sphere.m)
        mᵣ = sphere.m / (1.0 + sphere.β * sphere.m)

        return radius_and_center *
               (@sprintf "(%.4e,%.4e) (%.4e,%.4e)" mₗ.re mₗ.im mᵣ.re mᵣ.im)
    end
end

function format_spheres(cfg::STMMConfig)
    return """
number_spheres
$(length(cfg.spheres))
sphere_data
$(join(map(format_sphere, cfg.spheres),"\n"))
end_of_sphere_data"""
end

function format_fft_options(cfg::STMMConfig)
    if !cfg.use_fft_translation
        return "fft_translation_option\nfalse"
    else
        return """
        fft_translation_option
        true
        min_fft_nsphere
        $(length(cfg.spheres))
        cell_volume_fraction
        $(cfg.cell_volume_fraction)
        node_order
        $(cfg.node_order)"""
    end
end

# FIXME Optically active layers are not supported in MSTM v4 yet.
# See https://github.com/dmckwski/MSTM/issues/6
function format_layers(cfg::STMMConfig)::String
    N = length(cfg.layers)

    return """
number_plane_boundaries
$(N - 1)
layer_thickness
$(N >= 3 ? join(map(x -> x.d, cfg.layers[2:N-1]), ",") : "")
layer_ref_index
$(join(map(x -> (@sprintf "(%.4e,%.4e)" x.m.re x.m.im), cfg.layers), ","))"""
end

function format_orientation(cfg::STMMConfig)::String
    if cfg.orientation == RandomOrientation
        if cfg.use_monte_carlo_integration
            return "incidence_average\ntrue\nnumber_incident_directions\n$(cfg.number_incident_directions)"
        else
            return "random_orientation\ntrue"
        end
    else
        return "incident_alpha_deg\n$(cfg.α)\nincident_beta_deg\n$(cfg.β)\nazimuthal_average\n$(cfg.azimuthal_average)"
    end
end

function format_near_field(cfg::STMMConfig)::String
    if cfg.calculate_near_field && cfg.orientation == FixedOrientation
        return """
    calculate_near_field
    true
    near_field_output_file
    $(cfg.near_field_output_file)
    near_field_minimum_border
    $(cfg.near_field_x_min),$(cfg.near_field_y_min),$(cfg.near_field_z_min)
    near_field_maximum_border
    $(cfg.near_field_x_max),$(cfg.near_field_y_max),$(cfg.near_field_z_max)
    near_field_step_size
    $(cfg.near_field_step_size)
    near_field_calculation_model
    $(cfg.near_field_calculation_model)
    store_surface_vector
    $(cfg.store_surface_vector)""" * (cfg.store_surface_vector ? """
    \nnear_field_expansion_spacing
    $(cfg.near_field_expansion_spacing)
    near_field_expansion_order
    $(cfg.near_field_expansion_order)""" : "")
    else
        return "calculate_near_field\nfalse"
    end
end

function format_lattice(cfg::STMMConfig)
    if !cfg.periodic
        return "periodic_lattice\nfalse"
    else
        if cfg.orientation == RandomOrientation
            throw(ErrorException("Lattice mode is incompatible with random operation!"))
        end

        x_min = minimum(s.x - s.r for s in cfg.spheres)
        x_max = maximum(s.x + s.r for s in cfg.spheres)
        y_min = minimum(s.y - s.r for s in cfg.spheres)
        y_max = maximum(s.y + s.r for s in cfg.spheres)
        dx = x_max - x_min
        dy = y_max - y_min
        cx = max(dx, cfg.cell_size[1])
        cy = max(dy, cfg.cell_size[2])
        return "periodic_lattice\ntrue\ncell_width\n$cx,$cy"
    end
end

# We set `print_sphere_data` to `false` since we do not need them.
function write_input(cfg::STMMConfig)
    inp = """
    $(format_spheres(cfg))
    $(format_layers(cfg))
    output_file
    $(cfg.output_file)
    run_file
    $(cfg.run_print_file)
    print_sphere_data
    false
    mie_epsilon
    $(cfg.mie_epsilon)
    translation_epsilon
    $(cfg.translation_epsilon)
    solution_epsilon
    $(cfg.solution_epsilon)
    t_matrix_convergence_epsilon
    $(cfg.t_matrix_epsilon)
    max_iterations
    $(cfg.max_iterations)
    max_t_matrix_order
    $(cfg.max_t_matrix_order)
    store_translation_matrix
    $(cfg.store_translation_matrix)
    $(format_fft_options(cfg))
    $(format_orientation(cfg))
    $(format_near_field(cfg))
    incident_frame
    $(cfg.frame == IncidentFrame)
    normalize_s11
    $(cfg.normalize_scattering_matrix)
    gaussian_beam_constant
    $(cfg.gaussian_beam_constant)
    gaussian_beam_focal_point
    $(@sprintf("(%.4e,%.4e,%.4e)", cfg.gaussian_beam_focal_point...))
    calculate_scattering_matrix
    $(cfg.calculate_scattering_matrix)
    scattering_map_model
    $(cfg.scattering_map_model)
    scattering_map_increment
    $(cfg.Δθ)
    scattering_map_dimension
    $(cfg.scattering_map_dimension)
    $(format_lattice(cfg))
    end_of_options
    """

    open("mstm.inp", "w") do io
        write(io, inp)
    end
end

function read_ints(s::AbstractString)::Vector{Int}
    return map(x -> parse(Int, x), filter(!isempty, split(s, r"\s+")))
end

function read_floats(s::AbstractString)::Vector{Float64}
    return map(x -> parse(Float64, x), filter(!isempty, split(s, r"\s+")))
end

function collect_output(cfg::STMMConfig)::MSTMOutput
    @assert isfile(cfg.output_file)

    out = split(read(open(cfg.output_file), String), "\n")
    i = 1

    if cfg.orientation == FixedOrientation
        q = cfg.periodic ? zeros(9) : zeros(24)
    else
        q = cfg.use_monte_carlo_integration ? zeros(9) : zeros(3)
    end

    θ = Float64[]
    kx = Float64[]
    ky = Float64[]
    s11 = Float64[]
    s12 = Float64[]
    s13 = Float64[]
    s14 = Float64[]
    s21 = Float64[]
    s22 = Float64[]
    s23 = Float64[]
    s24 = Float64[]
    s31 = Float64[]
    s32 = Float64[]
    s33 = Float64[]
    s34 = Float64[]
    s41 = Float64[]
    s42 = Float64[]
    s43 = Float64[]
    s44 = Float64[]
    type = String[]

    while !occursin("calculation results for run", out[i])
        i += 1
    end

    # Skip header
    i += cfg.orientation == FixedOrientation ? 4 : 6

    if cfg.orientation == FixedOrientation
        if cfg.periodic
            # Lattice mode
            q[1:9] = collect(@scanf out[i+1] "%e%e%e%e%e%e%e%e%e" zeros(9)...)[2:end]

            i += 6
        elseif length(cfg.layers) == 1
            # One layer

            # Read total efficiencies
            q[1:9] = read_floats(out[i+1])
            i += 2

            # Read hemispherical scattering efficiencies
            q[10:15] = read_floats(out[i+1])
            i += 3
        else
            # Multiple layers

            # Read total efficiencies
            q[1:9] = read_floats(out[i+1])
            i += 2

            # Read down/up extinction efficiencies
            q[16:21] = read_floats(out[i+1])
            i += 2

            # Read hemispherical scattering efficiencies
            q[10:15] = read_floats(out[i+1])
            i += 2

            # Read waveguide scattering efficiencies
            q[22:24] = read_floats(out[i+1])
            i += 3
        end

        current_type = "total"
        if occursin("reflection", out[i])
            current_type = "reflection"
            i += 4
        elseif occursin("kx", out[i])
            current_type = "backward"
            i += 1
        else
            i += 3
        end

        # Read scattering matrix
        while !isempty(out[i])
            if occursin("transmission", out[i])
                current_type = "transmission"
                i += 2
                continue
            elseif occursin("forward", out[i])
                current_type = "forward"
                i += 2
                if cfg.periodic
                    i += 2
                end
                continue
            end

            v = read_floats(out[i])
            i += 1

            push!(type, current_type)

            # TODO: need to check whether lattice output can support θ instead kx-ky convention
            if cfg.scattering_map_model == 0 && !cfg.periodic
                push!(θ, v[1])
            else
                push!(kx, v[1])
                push!(ky, v[2])
            end

            push!(s44, v[end])
            push!(s43, v[end-1])
            push!(s42, v[end-2])
            push!(s41, v[end-3])
            push!(s34, v[end-4])
            push!(s33, v[end-5])
            push!(s32, v[end-6])
            push!(s31, v[end-7])
            push!(s24, v[end-8])
            push!(s23, v[end-9])
            push!(s22, v[end-10])
            push!(s21, v[end-11])
            push!(s14, v[end-12])
            push!(s13, v[end-13])
            push!(s12, v[end-14])
            push!(s11, v[end-15])
        end

        if cfg.scattering_map_model == 0 && !cfg.periodic
            scattering_matrix = DataFrame(;
                θ,
                s11,
                s12,
                s13,
                s14,
                s21,
                s22,
                s23,
                s24,
                s31,
                s32,
                s33,
                s34,
                s41,
                s42,
                s43,
                s44,
                type
            )
        else
            scattering_matrix = DataFrame(;
                kx,
                ky,
                s11,
                s12,
                s13,
                s14,
                s21,
                s22,
                s23,
                s24,
                s31,
                s32,
                s33,
                s34,
                s41,
                s42,
                s43,
                s44,
                type
            )
        end

        if cfg.calculate_near_field && isfile(cfg.near_field_output_file)
            nf = split(read(open(cfg.near_field_output_file), String), "\n")

            Ns = read_ints(nf[3])[1]
            sx = Float64[]
            sy = Float64[]
            sz = Float64[]
            r = Float64[]
            for j = 1:Ns
                v = read_floats(nf[j+3])
                push!(sx, v[1])
                push!(sy, v[2])
                push!(sz, v[3])
                push!(r, v[4])
            end
            spheres = DataFrame(; x = sx, y = sy, z = sz, r)

            Nb = read_ints(nf[Ns+4])[1]
            Nx, Ny, Nz = read_ints(nf[Ns+Nb+7])

            # Read near field data
            x = Float64[]
            y = Float64[]
            z = Float64[]
            Ex₌ = ComplexF64[]
            Ex⊥ = ComplexF64[]
            Ey₌ = ComplexF64[]
            Ey⊥ = ComplexF64[]
            Ez₌ = ComplexF64[]
            Ez⊥ = ComplexF64[]
            Hx₌ = ComplexF64[]
            Hx⊥ = ComplexF64[]
            Hy₌ = ComplexF64[]
            Hy⊥ = ComplexF64[]
            Hz₌ = ComplexF64[]
            Hz⊥ = ComplexF64[]
            for j = 1:Nx*Ny*Nz
                v = read_floats(nf[j+7+Ns+Nb])
                push!(x, v[1])
                push!(y, v[2])
                push!(z, v[3])
                push!(Ex₌, complex(v[4], v[5]))
                push!(Ey₌, complex(v[6], v[7]))
                push!(Ez₌, complex(v[8], v[9]))
                push!(Hx₌, complex(v[10], v[11]))
                push!(Hy₌, complex(v[12], v[13]))
                push!(Hz₌, complex(v[14], v[15]))
                push!(Ex⊥, complex(v[16], v[17]))
                push!(Ey⊥, complex(v[18], v[19]))
                push!(Ez⊥, complex(v[20], v[21]))
                push!(Hx⊥, complex(v[22], v[23]))
                push!(Hy⊥, complex(v[24], v[25]))
                push!(Hz⊥, complex(v[26], v[27]))
            end

            field = DataFrame(;
                x,
                y,
                z,
                Ex₌,
                Ex⊥,
                Ey₌,
                Ey⊥,
                Ez₌,
                Ez⊥,
                Hx₌,
                Hx⊥,
                Hy₌,
                Hy⊥,
                Hz₌,
                Hz⊥
            )

            near_field = NearField(spheres, field)
        else
            near_field = nothing
        end

        return cfg.periodic ? MSTMLatticeOutput(q..., scattering_matrix) : MSTMFixedOutput(q..., scattering_matrix, near_field)
    elseif !cfg.use_monte_carlo_integration
        # Random orientation, analytical average

        # Read efficiencies
        q[1:3] = read_floats(out[i+1])
        i += 5

        # Read scattering matrix
        Nθ = Int(round((cfg.θ_max - cfg.θ_min) / cfg.Δθ)) + 1
        for j = 1:Nθ
            v = read_floats(out[i+j])

            push!(θ, v[1])
            push!(s11, v[2])
            push!(s12, v[3])
            push!(s22, v[4])
            push!(s33, v[5])
            push!(s34, v[6])
            push!(s44, v[7])
        end
        scattering_matrix = DataFrame(; θ, s11, s12, s22, s33, s34, s44)

        return MSTMRandomOutput(q..., scattering_matrix)
    else
        # Random orientation, Monte Carlo

        # Read efficiencies
        q[1:3] = read_floats(out[i+1])
        i += 2

        q[4:9] = read_floats(out[i+1])
        i += 4

        # Read scattering matrix
        Nθ, _ = read_ints(out[i])
        i += 1

        for j = 1:Nθ
            v = read_floats(out[i+j])
            push!(θ, v[1])
            push!(s11, v[2])
            push!(s12, v[3])
            push!(s22, v[4])
            push!(s33, v[5])
            push!(s34, v[6])
            push!(s44, v[7])
        end

        scattering_matrix_azimuthal_average = DataFrame(; θ, s11, s12, s22, s33, s34, s44)

        for v in [θ, s11, s12, s22, s33, s34, s44]
            resize!(v, 0)
        end

        i += Nθ + 1
        @assert occursin("diffuse", out[i])
        i += 1

        for j = 1:Nθ
            v = read_floats(out[i+j])
            push!(θ, v[1])
            push!(s11, v[2])
            push!(s12, v[3])
            push!(s22, v[4])
            push!(s33, v[5])
            push!(s34, v[6])
            push!(s44, v[7])
        end

        scattering_matrix_diffuse = DataFrame(; θ, s11, s12, s22, s33, s34, s44)

        i += Nθ + 1
        @assert occursin(
            "azimuthal averaged scattering matrix expansion coefficients",
            out[i],
        )
        i += 2

        n = Int[]
        a11 = Float64[]
        a44 = Float64[]
        a12 = Float64[]
        a34 = Float64[]
        a22p = Float64[]
        a22m = Float64[]

        while !isempty(out[i])
            v = read_floats(out[i])

            push!(n, v[1])
            push!(a11, v[2])
            push!(a44, v[3])
            push!(a12, v[4])
            push!(a34, v[5])
            push!(a22p, v[6])
            push!(a22m, v[7])

            i += 1
        end

        scattering_matrix_expansion_coefficients =
            DataFrame(; n, a11, a44, a12, a34, a22p, a22m)

        return MSTMMonteCarloOutput(
            q...,
            scattering_matrix_azimuthal_average,
            scattering_matrix_diffuse,
            scattering_matrix_expansion_coefficients,
        )
    end
end

end


text:
using Random
Random.seed!(209)

@testset "0D grids" begin
  @test Space((1,)).graph == Agents.Graph(1)
end

@testset "1D grids" begin
  a = Space((5,1))
  ae = collect(Agents.LightGraphs.edges(a.graph))

  b = Space((5,1), periodic=true)
  be = collect(Agents.LightGraphs.edges(b.graph))

  @test ae == [Agents.LightGraphs.Edge(1,2),Agents.LightGraphs.Edge(2,3), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(4,5)]
  @test be == [Agents.LightGraphs.Edge(1,2), Agents.LightGraphs.Edge(1,5), Agents.LightGraphs.Edge(2,3), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(4,5)]
end

@testset "2D grids" begin
  @test Space((2,1)).graph == Agents.Graph(2,1)

  a = Space((2,3))
  b = Space((2,3), periodic=true) # 2D grid

  @test a.dimensions == (2,3)
  @test b.dimensions == (2,3)
  @test length(a.agent_positions) == 6
  @test length(b.agent_positions) == 6

  @test Agents.nv(a) == 6
  @test Agents.ne(a) == 7
  @test Agents.nv(b) == 6
  @test Agents.ne(b) == 9

  ae = collect(Agents.LightGraphs.edges(a.graph))
  be = collect(Agents.LightGraphs.edges(b.graph))

  @test ae == [Agents.LightGraphs.Edge(1,2), Agents.LightGraphs.Edge(1,3), Agents.LightGraphs.Edge(2,4), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(3,5), Agents.LightGraphs.Edge(4,6), Agents.LightGraphs.Edge(5,6)]
  @test be == [Agents.LightGraphs.Edge(1,2), Agents.LightGraphs.Edge(1,3), Agents.LightGraphs.Edge(1,5), Agents.LightGraphs.Edge(2,4), Agents.LightGraphs.Edge(2,6), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(3,5), Agents.LightGraphs.Edge(4,6), Agents.LightGraphs.Edge(5,6)]
end

@testset "2D Moore" begin
  a = Space((3, 3), moore=true)
  @test Agents.nv(a) == 9
  @test Agents.ne(a) == 20

  a = Space((3, 4), moore=true)
  @test Agents.nv(a) == 12
  @test Agents.ne(a) == 29

  b = Space((3, 2), moore=true, periodic=true)
  @test Agents.nv(b) == 6
  @test Agents.ne(b) == 15

  b = Space((3, 3), moore=true, periodic=true)
  @test Agents.nv(b) == 9
  @test Agents.ne(b) == 36
end

@testset "3D grid" begin
  g1 = Space((2,3,2))
  g2 = Space((2,3,3))
  g3 = Space((2,3,2), periodic=true)
  g4 = Space((2,3,3), periodic=true)
  @test Agents.ne(g1) == 20
  @test Agents.ne(g2) == 33
  @test Agents.ne(g3) == 24
  @test Agents.ne(g4) == 39
end

@testset "grid coord/vertex conversions" begin
  @test coord2vertex((2,2,1), (3,4,1)) == 5
  @test coord2vertex((2,2), (3,4)) == 5
  @test coord2vertex((2,2,1), (3,4,1)) == 5
  @test coord2vertex((2,2), (3,4)) == 5
  @test vertex2coord(5, (3,4,1)) == (2,2,1)
  @test vertex2coord(5, (3,4)) == (2,2)

  @test coord2vertex((2,2,1), (2,3,3)) == 4
  @test coord2vertex((2,2,2), (2,3,3)) == 10
  @test coord2vertex((2,2,3), (2,3,3)) == 16
  @test coord2vertex(((1,3,1)), (2,3,3)) == 5
  @test coord2vertex((1,3,2), (2,3,3)) == 11
  @test coord2vertex((1,3,3), (2,3,3)) == 17

  @test vertex2coord(5, (2,3,3)) == (1,3,1)
  @test vertex2coord(7, (2,3,3)) == (1,1,2)
  @test vertex2coord(13, (2,3,3)) == (1,1,3)
  @test vertex2coord(13, (2,3,3)) == (1,1,3)
  @test vertex2coord(15, (2,3,3)) == (1,2,3)
  @test vertex2coord(18, (2,3,3)) == (2,3,3)
  @test vertex2coord(18, (2,3,3)) == (2,3,3)
end

mutable struct Agent1 <: AbstractAgent
  id::Int
  pos::Tuple{Int,Int}
end

mutable struct Agent2 <: AbstractAgent
  id::Int
  pos::Tuple{Int,Int}
  p::Int
end

@testset "Agent-Space interactions" begin

  model = model_initiation(f=0.1, d=0.8, p=0.1, griddims=(20, 20), seed=2)  # forest fire model

  agent = model.agents[1]
  move_agent!(agent, (3,4), model)  # node number 63
  @test agent.pos == (3,4)
  @test agent.id in model.space.agent_positions[63]

  agent = model.agents[2]
  move_agent!(agent, 83, model)  # pos (3,5)
  @test agent.pos == (3,5)
  @test agent.id in model.space.agent_positions[83]

  new_pos = move_agent!(agent, model)
  @test agent.id in model.space.agent_positions[coord2vertex(new_pos, model)]

  add_agent!(agent, (2,9), model)
  @test agent.pos == (2,9)
  @test agent.id in model.space.agent_positions[coord2vertex((2,9), model)]
  @test agent.id in model.space.agent_positions[coord2vertex(new_pos, model)]

  model1 = ABM(Agent1, Space((3,3)))
  add_agent!(1, model1)
  @test model1.agents[1].pos == (1, 1)
  add_agent!((2,1), model1)
  @test model1.agents[2].pos == (2, 1)

  model2 = ABM(Agent2, Space((3,3)))
  add_agent!(1, model2, 3)
  @test model2.agents[1].pos == (1,1)
  add_agent!((2,1), model2, 2)
  @test model2.agents[2].pos == (2,1)

  @test agent.id in get_node_contents(agent, model)

  ii = model.agents[length(model.agents)]
  @test id2agent(ii.id, model) == model.agents[ii.id]

  agent = model.agents[1]
  agent_pos = coord2vertex(agent.pos, model)
  kill_agent!(agent, model)
  @test_throws KeyError id2agent(1, model)
  @test !in(1, Agents.agent_positions(model)[agent_pos])
end


text:
@testset "Replication-related" begin
  println(testname("Replication-related tests"))
  print("  [  ] Streaming changes ")
  count = 0
  maxch = 5
  for ch in @task changes_streaming(db; limit=maxch)
    count += 1
  end
  @test count == maxch + 1 # In stream mode, last item is the CouchDB "last_seq" so need to add 1.
  println("\r  [OK] Streaming changes")

  print("  [  ] Static changes ")
  data = changes(db; limit=maxch, conflicts=true, include_docs=true, attachments=true, att_encoding_info=true)
  @test maxch == length(data["results"]) # In static mode, "last_seq" is a key in the dict.
  println("\r  [OK] Static changes")

  print("  [  ] Filtered changes ")
  data2 = changes(db; doc_ids=[data["results"][1]["id"], data["results"][2]["id"], data["results"][3]["id"]])
  @test length(data2["results"]) == 3
  println("\r  [OK] Filtered changes")

  print("  [  ] Streaming changes, filtered ")
  count = 0
  for ch in @task changes_streaming(db; doc_ids=[data["results"][1]["id"], data["results"][2]["id"], data["results"][3]["id"]])
    count += 1
  end
  @test count > 0
  println("\r  [OK] Streaming changes, filtered")

  print("  [  ] revs_diff ")
  fakerev = "2-1f0e2f0d841ba6b7e3d735b870ebeb8c"
  fakerevs = Dict(data["results"][1]["id"] => [data["results"][1]["changes"][1]["rev"], fakerev])
  diff = revs_diff(db; data=fakerevs)
  @test haskey(diff, data["results"][1]["id"])
  @test diff[data["results"][1]["id"]]["missing"][1] == fakerev
  println("\r  [OK] revs_diff")

  print("  [  ] bulk_get (note: needs CouchDB2 or Cloudant DBNext) ")
  fetchdata = [
    Dict{String, String}("id" => data["results"][1]["id"], "rev" => data["results"][1]["changes"][1]["rev"]),
  ]
  response = bulk_get(db; data=fetchdata)
  @test length(response["results"]) == 1
  println("\r  [OK] bulk_get (note: needs CouchDB2 or Cloudant DBNext)")
end


text:
#  this function takes in a string and a char
#  or sub_string and removes all instances
#  of char from the passed_string

function remove_char(passed_string, search_character_string)

aaa=split(passed_string, search_character_string)
temp="x"
  for i in 1:size(aaa, 1)
    if length(aaa[i])!=0
      temp=string(temp,aaa[i])
    end
  end

  return_string = temp[2:end]
  return return_string

end


text:
"""
    SARSALambdaSolver

SARSA-λ implementation for tabular MDPs, assign credits using eligibility traces

Parameters:
- `exp_policy::ExplorationPolicy`:
    Exploration policy to select the actions
- `n_episodes::Int64`:
    Number of episodes to train the Q table
    default: `100`
- `max_episode_length::Int64`:
    Maximum number of steps before the episode is ended
    default: `100`
- `learning_rate::Float64`:
    Learning rate
    defaul: `0.001`
- `lambda::Float64`:
    Exponential decay parameter for the eligibility traces
    default: `0.5`
- `eval_every::Int64`:
    Frequency at which to evaluate the trained policy
    default: `10`
- `n_eval_traj::Int64`:
    Number of episodes to evaluate the policy
- `rng::AbstractRNG`: random number generator
- `verbose::Bool`:
    print information during training
    default: `true`
"""
@with_kw mutable struct SARSALambdaSolver{E<:ExplorationPolicy} <: Solver
   n_episodes::Int64 = 100
   max_episode_length::Int64 = 100
   learning_rate::Float64 = 0.001
   exploration_policy::E
   Q_vals::Union{Nothing, Matrix{Float64}} = nothing
   eligibility::Union{Nothing, Matrix{Float64}} = nothing
   lambda::Float64 = 0.5
   eval_every::Int64 = 10
   n_eval_traj::Int64 = 20
   rng::AbstractRNG = Random.GLOBAL_RNG
   verbose::Bool = true
end

function solve(solver::SARSALambdaSolver, mdp::MDP)
    rng = solver.rng
    if solver.Q_vals === nothing
        Q = zeros(length(states(mdp)), length(actions(mdp)))
    else
        Q = solver.Q_vals
    end
    if solver.eligibility === nothing
        ecounts = zeros(length(states(mdp)), length(actions(mdp)))
    else
        ecounts = solver.eligibility
    end
    exploration_policy = solver.exploration_policy
    on_policy = ValuePolicy(mdp, Q)
    sim = RolloutSimulator(rng=rng, max_steps=solver.max_episode_length)
    k = 0
    for i = 1:solver.n_episodes
        s = rand(rng, initialstate(mdp))
        a = action(exploration_policy, on_policy, k, s)
        t = 0
        while !isterminal(mdp, s) && t < solver.max_episode_length
            sp, r = @gen(:sp, :r)(mdp, s, a, rng)
            k += 1
            ap = action(exploration_policy, on_policy, k, sp)
            si = stateindex(mdp, s)
            ai = actionindex(mdp, a)
            spi = stateindex(mdp, sp)
            api = actionindex(mdp, ap)
            delta = r + discount(mdp) * Q[spi,api] - Q[si,ai]
            ecounts[si,ai] += 1
            for es in states(mdp)
                for ea in actions(mdp)
                    esi, eai = stateindex(mdp, es), actionindex(mdp, ea)
                    Q[esi,eai] += solver.learning_rate * delta * ecounts[esi,eai]
                    ecounts[esi,eai] *= discount(mdp) * solver.lambda
                end
            end
            s, a = sp, ap
            t += 1
        end
        if i % solver.eval_every == 0
            r_tot = 0.0
            for traj in 1:solver.n_eval_traj
                r_tot += simulate(sim, mdp, on_policy, rand(rng, initialstate(mdp)))
            end
            solver.verbose ? println("On Iteration $i, Returns: $(r_tot/solver.n_eval_traj)") : nothing
        end
    end
    return on_policy
end

POMDPLinter.@POMDP_require solve(solver::SARSALambdaSolver, problem::Union{MDP,POMDP}) begin
    P = typeof(problem)
    S = statetype(P)
    A = actiontype(P)
    @req initialstate(::P, ::AbstractRNG)
    @req gen(::DDNOut{(:sp, :r)}, ::P, ::S, ::A, ::AbstractRNG)
    @req stateindex(::P, ::S)
    @req states(::P)
    ss = states(problem)
    @req length(::typeof(ss))
    @req actions(::P)
    as = actions(problem)
    @req length(::typeof(as))
    @req actionindex(::P, ::A)
    @req discount(::P)
end


text:
#=============================================
    MAKERGAULx.jl

    different formula for each cycle (dfc)

=============================================#

## write functions to fit MAKERGAUL model here, which will be called in `mod_bl_q` in "amp.jl"

import JuMP: Model, @variable, @constraint, @NLconstraint, @NLobjective,
    solve, getvalue, getobjectivevalue

## bounds
const fb_B_MULTIPLE     = 1.9
const d0_LB             = 1e-14
const d0_UB             = Inf_T ## used: 0.1 (Bultmann 2013)
const eu0_inh_LB        = 0.0001
const eu0_UB_MULTIPLE   = 10
const inh_UB_MULTIPLE   = 10

## start values
const eu0_START             = 7e3  ## used: eu0_inh_LB, 0.01, 1, 50
const MAKERGAUL_d0_START    = 1    ## used: 0, 1e-14 (change_d3), 0.01, 1 (change_d2), 50 (change_d1)
const inh_START             = 4e-6 ## used: eu0_inh_LB (Infeasible for flat line), 0, 0.05, 1, (Invalid_Number_Detected for flat line), 10 (Infeasible for flat line), 50 `:Optimal` when `max_of_idx == 1`, "Invalid_Number_Detected" for the rest: eu0_START = 0.01, inh_START = 0; eu0_START = 50, inh_START = 0.05; eu0_START = 50, inh_START = 1; eu0_START = 50, inh_START = 50;

## function definitions

function pred_from_nm1(
    ::Union{Type{Val{MAKERGAUL3}}, Type{Val{MAKERGAUL4}}},
    eu_nm1      ::Real,
    d_nm1       ::Real,
    inh         ::Real
)
    eu_n = eu_nm1 / (1 + inh * d_nm1)
    d_n = d_nm1 + d_nm1 * eu_n / (eu_n + d_nm1)
    return [eu_n d_n]
end ## pred_from_nm1()

function pred_from_cycs( ## 0.7to1.2e-5 sec for 40 cycles on PC
    ::Type{Val{MAKERGAUL3}},
    cycs2fit    ::AbstractVector,
    fb          ::Real,
    eu0         ::Real,
    d0          ::Real,
    inh         ::Real
)
    max_cyc = maximum(cycs)
    pred_ary_eu_d_w0 = [eu0 d0]
    i = 1
    while i <= max_cyc
        pred_ary_eu_d_w0 = vcat(
            pred_ary_eu_d_w0,
            pred_from_nm1(MAKERGAUL3(), pred_ary_eu_d_w0[i,:]..., inh)
        )
        i += 1
    end
    fs_w0 = pred_ary_eu_d_w0[:, 2] + fb
    return fs_w0[2:end][map(Int_T, cycs)]
end ## pred_from_cycs(::Type{Val{MAKERGAUL3}})

function pred_from_cycs(
    ::Type{Val{MAKERGAUL4}},
    cycs        ::AbstractVector,
    fb          ::Real,
    bl_k        ::Real,
    eu0         ::Real,
    d0          ::Real,
    inh         ::Real
)
    max_cyc = maximum(cycs)
    pred_ary_eu_d_w0 = [eu0 d0]
    i = 1
    while i <= max_cyc
        pred_ary_eu_d_w0 = vcat(
            pred_ary_eu_d_w0,
            pred_from_nm1(MAKERGAUL4(), pred_ary_eu_d_w0[i,:]..., inh)
        )
        i += 1
    end
    return fb + bl_k * cycs .+ pred_ary_eu_d_w0[:, 2][2:end][map(Int_T, cycs)]
end ## pred_from_cycs(:MAKERGAUL4)

function fit(
    ::Type{Val{MAKERGAUL3}},
    cycs        ::AbstractVector, ## continous integers or not
    obs_fluos   ::AbstractVector,
    wts         ::AbstractVector =ones(length(obs_fluos));
    kwargs_Model... ## argument for `JuMP.Model`
)
    ## find maximum observed fluorescence
    max_of, max_of_idx = findmax(obs_fluos)
    idc2fit = 1:max_of_idx
    obs2fit = obs_fluos[idc2fit]
    cycs2fit = 1:cycs[max_of_idx]

    jmp_model = Model(;kwargs_Model...)

    min_of = minimum(obs_fluos)
    of_diff = max_of - min_of
    fb_b_abs = fb_B_MULTIPLE * of_diff

    @variable(jmp_model,
        # min_of - fb_b_abs <= fb <= min_of + fb_b_abs,
        fb,
        start=min_of)
    @variable(jmp_model, # eu0
        # eu0_inh_LB <= eu0 <= eu0_UB_MULTIPLE * of_diff,
        eu0 >= 0,
        start=eu0_START)
    # @variable(jmp_model, d0_LB <= d0 <= d0_UB, start=MAKERGAUL_d0_START) ## change_a1
    @variable(jmp_model, d0 >= 0, start=MAKERGAUL_d0_START) ## change_a2
    @variable(jmp_model, ## inh
        # eu0_inh_LB <= inh <= inh_UB_MULTIPLE * of_diff,
        inh >= 0,
        start=inh_START)
    @variable(jmp_model, f[cycs2fit])
    # @variable(jmp_model, eu[cycs2fit]) ## change_b1
    @variable(jmp_model, eu[cycs2fit] >= 0) ## change_b2
    # @variable(jmp_model, d[cycs2fit]) ## change_c1
    @variable(jmp_model, d[cycs2fit] >= 0) ## change_c2

    @constraint(jmp_model, f_constr[cyc in cycs2fit], f[cyc] == fb + d[cyc])
    # change_e1: with division
    @NLconstraint(jmp_model, eu_constr_01, eu[1] == eu0 / (1 + inh * d0))
    @NLconstraint(jmp_model, d_constr_01, d[1] == d0 + d0 * eu[1] / (eu[1] + d0))
    @NLconstraint(jmp_model, eu_constr_2p[cyc in cycs2fit[2:end]], eu[cyc] == eu[cyc-1] / (1 + inh * d[cyc-1]))
    # @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] == d[cyc-1] + d[cyc-1] * eu[cyc] / (eu[cyc] + d[cyc-1])) # "Invalid_Number_Detected"
    @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] == d[cyc-1] + 1 / (1 / d[cyc-1] + 1 / eu[cyc]))
    #
    ## change_e2: get rid of division by multiplying both sides by denominator
    # @NLconstraint(jmp_model, eu_constr_01, eu[1] * (1 + inh * d0) == eu0)
    # @NLconstraint(jmp_model, d_constr_01, d[1] * (eu[1] + d0) == d0 * (d0 + 2 * eu[1]))
    # @NLconstraint(jmp_model, eu_constr_2p[cyc in cycs2fit[2:end]], eu[cyc] * (1 + inh * d[cyc-1]) == eu[cyc-1])
    # # @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] == d[cyc-1] + d[cyc-1] * eu[cyc] / (eu[cyc] + d[cyc-1])) # "Invalid_Number_Detected"
    # @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] * (eu[cyc] + d[cyc-1]) == d[cyc-1] * (d[cyc-1] + 2 * eu[cyc]))

    ## ssre by mean of predicted
    # len_cycs = length(cycs)
    # @NLobjective(
    #     jmp_model, Min,
    #     sum( ((f[cycs[idx]] - obs_fluos[idx]) / (sum(f[cycs[idx]] / len_cycs for idx in idc2fit) - obs_fluos[idx])) ^ 2 for idx in idc2fit )
    # )
    ## ssre by mean of observed
    # mean_of = mean(obs2fit)
    # @NLobjective(
    #     jmp_model, Min,
    #     sum( ((f[cycs[idx]] - obs_fluos[idx]) / (mean_of - obs_fluos[idx])) ^ 2 for idx in idc2fit )
    # )
    ## sse # the only one with optimal so far
    @NLobjective(jmp_model, Min, sum((f[cycs[idx]] - obs_fluos[idx]) ^ 2 for idx in idc2fit)) ## can solve to optimal when ssre can't

    # @NLobjective(
    #     jmp_model, Min,
    #     @eval $(parse(get_mak2_rsq_str(obs_fluos[1:max_d2_idx]))) ## `OutOfMemoryError()`
    # )

    status = solve(jmp_model)
    coef_strs = ["fb", "eu0", "d0", "inh"]
    coefs = map(getvalue, [fb, eu0, d0, inh])
    obj_val = getobjectivevalue(jmp_model)

    return MAKERGAUL3Fit(
        max_of_idx,
        coef_strs,
        coefs,
        status,
        obj_val,
        jmp_model,
        # init_coefs
    )
end ## fit(::Type{Val{MAKERGAUL3}})

function fit(
    ::Type{Val{MAKERGAUL4}},
    cycs        ::AbstractVector, ## continous integers or not
    obs_fluos   ::AbstractVector,
    wts         ::AbstractVector=ones(length(obs_fluos));
    kwargs_Model... ## argument for `JuMP.Model`
)
    ## find maximum observed fluorescence
    max_of, max_of_idx = findmax(obs_fluos)
    idc2fit = 1:max_of_idx
    obs2fit = obs_fluos[idc2fit]
    cycs2fit = 1:cycs[max_of_idx]

    ## fit a linear model to the estimated baseline portion of the curve
    d2_vec = finite_diff(cycs, obs_fluos; nu=2)
    max_d2_idx = indmax(d2_vec)
    idc2fit_4bl = 1:(max(1, max_d2_idx - 1))
    fb_start, bl_k_start = linreg(cycs[idc2fit_4bl], obs_fluos[idc2fit_4bl])

    jmp_model = Model(;kwargs_Model...)

    min_of = minimum(obs_fluos)
    of_diff = max_of - min_of
    fb_b_abs = fb_B_MULTIPLE * of_diff

    @variable(jmp_model,
        # min_of - fb_b_abs <= fb <= min_of + fb_b_abs,
        fb,
        start=fb_start)
    @variable(jmp_model,
        bl_k,
        start=bl_k_start)
    @variable(jmp_model, ## eu0
        # eu0_inh_LB <= eu0 <= eu0_UB_MULTIPLE * of_diff,
        eu0 >= 0,
        start=eu0_START)
    # @variable(jmp_model, d0_LB <= d0 <= d0_UB, start=MAKERGAUL_d0_START) ## change_a1
    @variable(jmp_model, d0 >= 0, start=MAKERGAUL_d0_START) ## change_a2
    @variable(jmp_model, ## inh
        # eu0_inh_LB <= inh <= inh_UB_MULTIPLE * of_diff,
        inh >= 0,
        start=inh_START)
    @variable(jmp_model, f[cycs2fit])
    # @variable(jmp_model, eu[cycs2fit])    ## change_b1
    @variable(jmp_model, eu[cycs2fit] >= 0) ## change_b2
    # @variable(jmp_model, d[cycs2fit])     ## change_c1
    @variable(jmp_model, d[cycs2fit] >= 0)  ## change_c2

    @constraint(jmp_model,
        f_constr[cyc in cycs2fit],
        f[cyc] == fb + bl_k * cyc + d[cyc])
    ## change_e1: with division
    @NLconstraint(jmp_model, eu_constr_01, eu[1] == eu0 / (1 + inh * d0))
    @NLconstraint(jmp_model, d_constr_01, d[1] == d0 + d0 * eu[1] / (eu[1] + d0))
    @NLconstraint(jmp_model, eu_constr_2p[cyc in cycs2fit[2:end]], eu[cyc] == eu[cyc-1] / (1 + inh * d[cyc-1]))
    # @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] == d[cyc-1] + d[cyc-1] * eu[cyc] / (eu[cyc] + d[cyc-1])) # "Invalid_Number_Detected"
    @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] == d[cyc-1] + 1 / (1 / d[cyc-1] + 1 / eu[cyc]))
    #
    ## change_e2: get rid of division by multiplying both sides by denominator
    # @NLconstraint(jmp_model, eu_constr_01, eu[1] * (1 + inh * d0) == eu0)
    # @NLconstraint(jmp_model, d_constr_01, d[1] * (eu[1] + d0) == d0 * (d0 + 2 * eu[1]))
    # @NLconstraint(jmp_model, eu_constr_2p[cyc in cycs2fit[2:end]], eu[cyc] * (1 + inh * d[cyc-1]) == eu[cyc-1])
    # # @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] == d[cyc-1] + d[cyc-1] * eu[cyc] / (eu[cyc] + d[cyc-1])) # "Invalid_Number_Detected"
    # @NLconstraint(jmp_model, d_constr_2p[cyc in cycs2fit[2:end]], d[cyc] * (eu[cyc] + d[cyc-1]) == d[cyc-1] * (d[cyc-1] + 2 * eu[cyc]))

    ## ssre by mean of predicted
    # len_cycs = length(cycs)
    # @NLobjective(
    #     jmp_model, Min,
    #     sum( ((f[cycs[idx]] - obs_fluos[idx]) / (sum(f[cycs[idx]] / len_cycs for idx in idc2fit) - obs_fluos[idx])) ^ 2 for idx in idc2fit )
    # )
    ## ssre by mean of observed
    # mean_of = mean(obs2fit)
    # @NLobjective(
    #     jmp_model, Min,
    #     sum( ((f[cycs[idx]] - obs_fluos[idx]) / (mean_of - obs_fluos[idx])) ^ 2 for idx in idc2fit )
    # )
    ## sse ## the only one with optimal so far
    @NLobjective(jmp_model, Min, sum((f[cycs[idx]] - obs_fluos[idx]) ^ 2 for idx in idc2fit)) ## can solve to optimal when ssre can't

    # @NLobjective(
    #     jmp_model, Min,
    #     @eval $(parse(get_mak2_rsq_str(obs_fluos[1:max_d2_idx]))) ## `OutOfMemoryError()`
    # )

    status = solve(jmp_model)
    coef_strs = ["fb", "bl_k", "eu0", "d0", "inh"]
    coefs = map(getvalue, [fb, bl_k, eu0, d0, inh])
    obj_val = getobjectivevalue(jmp_model)

    return MAKERGAUL4Fit(
        max_of_idx,
        fb_start,
        bl_k_start,
        coef_strs,
        coefs,
        status,
        obj_val,
        jmp_model,
        # init_coefs
    )
end ## fit(::Type{Val{MAKERGAUL4}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
struct Program
    context::AbstractContext
    vertexarray::AbstractArray
    uniforms::Dict{Symbol, Any}
    vertex_source::String
    fragment_source::String
end

struct InstancedProgram
    program::Program
    per_instance::AbstractArray
end

name_type_iter(x) = (s = Tables.schema(x); zip(s.names, s.types))

function getter_function(io, T, t_str, name, plot)
    println(io, t_str, " get_$(name)(){return $name;}")
end

function getter_function(io, ::AbstractSampler, t_str, name, plot)
end

function InstancedProgram(
        context::AbstractContext,
        vertshader, fragshader,
        instance::AbstractVector,
        per_instance::AbstractVector;
        uniforms...
    )
    instance_attributes = sprint() do io
        println(io, "\n// Per instance attributes: ")
        for (name, T) in name_type_iter(per_instance)
            t_str = type_string(context, T)
            println(io, "in ", t_str, " $name;")
            getter_function(io, T, t_str, name, uniforms)
        end
        println(io)
    end
    p = Program(
        context,
        instance_attributes * vertshader,
        fragshader,
        instance; uniforms...
    )
    InstancedProgram(p, per_instance)
end

function Program(
        context::AbstractContext,
        vertshader, fragshader,
        instance::AbstractVector;
        uniforms...
    )
    c_uniforms = Dict{Symbol, Any}()
    uniform_block = sprint() do io
        println(io, "\n// Uniforms: ")
        for (name, v) in uniforms
            vc = convert_uniform(context, v)
            t_str = type_string(context, vc)
            println(io, "uniform ", t_str, " $name;")
            getter_function(io, vc, t_str, name, uniforms)
            c_uniforms[name] = vc
        end
        println(io)
    end
    src = sprint() do io
        println(io, "// Instance inputs: ")
        for (name, T) in name_type_iter(instance)
            t_str = type_string(context, T)
            println(io, "in ", t_str, " $name;")
            getter_function(io, T, t_str, name, uniforms)
        end
        println(io, uniform_block)
        println(io)
        println(io, vertshader)
    end
    precision = """
    #version 300 es
    precision mediump int;
    precision mediump float;
    precision mediump sampler2D;
    precision mediump sampler3D;

    """
    Program(
        context, instance, c_uniforms,
        precision * src,
        precision * """
        out vec4 fragment_color;
        """ * uniform_block * fragshader
    )
end