text:
using ArgParse
using Parameters
using Printf
using Base: Float64, Float32, Int, Int32
import Base.read

const DefaultInterations = 8
const DefaultNPoses = 65536
const RefNPoses = 65536
const DefaultWGSize = 64
const DefaultPPWI = 4

const Zero = 0.0f0
const Quarter = 0.25f0
const Half = 0.5f0
const One = 1.0f0
const Two = 2.0f0
const Four = 4.0f0
const Cnstnt = 45.0f0

const HbtypeF = 70
const HbtypeE = 69
const Hardness = 38.0f0
const Npnpdist = 5.5f0
const Nppdist = 1.0f0

const Float32Max = typemax(Float32)

struct Vec3f32
  x::Float32
  y::Float32
  z::Float32
end

struct Vec4f32
  x::Float32
  y::Float32
  z::Float32
  w::Float32
end

struct Atom
  x::Float32
  y::Float32
  z::Float32
  type::Int32
end

function read(io::IO, ::Type{Atom})
  Atom(read(io, Float32), read(io, Float32), read(io, Float32), read(io, Int32))
end

struct FFParams
  hbtype::Int32
  radius::Float32
  hphb::Float32
  elsc::Float32
end

function read(io::IO, ::Type{FFParams})
  FFParams(read(io, Int32), read(io, Float32), read(io, Float32), read(io, Float32))
end

@with_kw mutable struct Params
  device::String = "1"
  list::Bool = false
  numposes::UInt = DefaultNPoses
  iterations::UInt = DefaultInterations
  wgsize::UInt = DefaultWGSize
  ppwi::UInt = DefaultPPWI
  deck::String = "../data/bm1"
  enzyme::Bool = false
  verify::Bool = false
end

struct Deck
  protein::AbstractArray{Atom}
  ligand::AbstractArray{Atom}
  forcefield::AbstractArray{FFParams}
  poses::AbstractArray{Float32,2}
end

const DeviceWithRepr = Tuple{Any,String,Any}

function parse_options(given::Params)
  s = ArgParseSettings()
  @add_arg_table s begin
    "--list", "-l"
    help = "List available devices"
    action = :store_true
    "--device", "-d"
    help = "Select device at DEVICE or the device name substring if DEVICE is not an integer, NOTE: Julia is 1-indexed"
    arg_type = String
    default = given.device
    "--numposes", "-n"
    help = "Compute energies for N poses"
    arg_type = Int
    default = convert(Int, given.numposes)
    "--iterations", "-i"
    help = "Repeat kernel I times"
    arg_type = Int
    default = convert(Int, given.iterations)
    "--wgsize", "-w"
    help = "Run with work-group size WGSIZE"
    arg_type = Int
    default = convert(Int, given.wgsize)
    "--ppwi", "-p"
    help = "Compute PPWI poses per work-item"
    arg_type = Int
    default = convert(Int, given.ppwi)
    "--deck"
    help = "Use the DECK directory as input deck"
    arg_type = String
    default = given.deck
    "--enzyme"
    help = "Use the enzyme"
    arg_type = Bool
    default = false
    "--verify"
    help = "Verify derivatives"
    arg_type = Bool
    default = false
  end
  args = parse_args(s)
  # surely there's a better way than doing this:
  for (arg, val) in args
    setproperty!(given, Symbol(arg), val)
  end
end

function read_structs(path::String, ::Type{T})::Vector{T} where {T}
  io = open(path, "r")
  size = filesize(io)
  [read(io, T) for _ = 1:size÷sizeof(T)]
end

function print_timings(params::Params, deck::Deck, elapsedSeconds::Float64, wgsize::UInt)

  # Average time per iteration
  averageIterationSeconds = elapsedSeconds / params.iterations

  natlig = length(deck.ligand)
  natpro = length(deck.protein)

  # Compute FLOP/s
  ops_per_wg = (wgsize * 27 + natlig * (2 + wgsize * 18 + natpro * (10 + wgsize * 30)) + wgsize)

  total_ops = ops_per_wg * (params.numposes / wgsize)
  flops = total_ops / averageIterationSeconds
  gflops = flops / 1.0e9

  total_finsts = 25.0 * natpro * natlig * params.numposes
  finsts = total_finsts / averageIterationSeconds
  gfinsts = finsts / 1e9

  interactions = params.numposes * natlig * natpro
  interactions_per_sec = interactions / averageIterationSeconds

  # Print stats
  @printf "- Kernel time:    %.03f ms\n" (elapsedSeconds * 1000.0)
  @printf "- Average time:   %.03f ms\n" (averageIterationSeconds * 1000.0)
  @printf "- Interactions/s: %.03f billion\n" interactions_per_sec / 1e9
  @printf "- GFLOP/s:        %.03f\n" gflops
  @printf "- GFInst/s:       %.03f\n" gfinsts
end

function main()

  params::Params = Params()
  parse_options(params)

  ds = devices()

  if params.list
    for (i, (_, name, type)) in enumerate(ds)
      println("[$i] ($type) $name")
    end
    exit(0)
  end

  deviceIndex = try
    index = parse(Int, params.device)
    if index < 1 || index > length(ds)
      error("Device $(index) out of range (1..$(length(ds))), NOTE: Julia is 1-indexed")
    end
    index
  catch
    index = findfirst(x -> occursin(params.device, "($(x[3])) $(x[2])"), ds)
    if index === nothing
      error("No device match the substring `$(params.device)`, see --list for available devices")
    end
    index
  end

  poses = permutedims( # reshape is column order so we flip it back again
    reshape(
      read_structs("$(params.deck)/poses.in", Float32),  # read poses as one long array
      (params.numposes, 6), # reshape it to 6 slices of numposes sized rows
    ),
  )

  deck = Deck(
    read_structs("$(params.deck)/protein.in", Atom),
    read_structs("$(params.deck)/ligand.in", Atom),
    read_structs("$(params.deck)/forcefield.in", FFParams),
    poses,
  )

  println("Poses     : ", size(deck.poses)[2])
  println("Iterations: ", params.iterations)
  println("Ligands   : ", length(deck.ligand))
  println("Protein   : ", length(deck.protein))
  println("Forcefield: ", length(deck.forcefield))
  println("Deck      : ", params.deck)
  println("WGsize    : ", params.wgsize)
  println("PPWI      : ", params.ppwi)

  set_zero_subnormals(true)
  # GC.enable(false)

  (energies, rumtimeSeconds, ppwi) = run(params, deck, ds[deviceIndex])

  # GC.enable(true)
  set_zero_subnormals(false)

  print_timings(params, deck, rumtimeSeconds, ppwi)

  output = open("energies.out", "w+")
  println("\nEnergies:")
  for i = 1:size(deck.poses)[2]
    @printf(output, "%7.2f\n", energies[i])
    if (i < 16)
      @printf("%7.2f\n", energies[i])
    end
  end

  ref_energies = open("$(params.deck)/ref_energies.out", "r")
  n_ref_poses = size(deck.poses)[2]
  if n_ref_poses > RefNPoses
    println("Only validating the first $(RefNPoses) poses")
    n_ref_poses = RefNPoses
  end

  maxdiff = 0.0
  for i = 1:n_ref_poses
    line = readline(ref_energies)
    if line == ""
      error("ran out of ref energies lines to verify")
    end
    e = parse(Float32, line)
    if (abs(e) < 1.0 && abs(energies[i]) < 1.0)
      continue
    end
    diff = abs(e - energies[i]) / e
    if (diff > 0.01)
      println("![$(i)] $(e) $(energies[i]) ~ $(diff)")
      maxdiff = diff
    end
  end

  @printf "Largest difference was %.03f%%.\n\n" 100 * maxdiff # Expect numbers to be accurate to 2 decimal places
  close(ref_energies)

end

max_stars_repo_name: wsmoses/Enzyme-BUDE

max_issues_repo_name: wsmoses/Enzyme-BUDE

max_forks_repo_name: wsmoses/Enzyme-BUDE


text:
using UpwindDifferences
using Test

using SparseArrays
using LinearAlgebra
using DelimitedFiles
using ProgressMeter

function evaluateparameters(ps)
    θ, ρ, r, y, dt = ps[1], ps[2], ps[3], ps[4], ps[5]
    xs = LinRange(0.01, 10.0, 1000)

    reward(x, c) = c ^ (1.0 - θ) / (1.0 - θ)
    drift(x, c) = r * x + y - c
    policy(x, dv) = dv ^ (-1.0 / θ)
    zerodrift(x) = r * x + y

    R, GF, GB = zeros(1000), zeros(1000), zeros(1000)
    UR = UpwindResult(R, GF, GB)
    bs = (zeros(1000-1), zeros(1000-1))
    dv = zeros(1000-1)

    function backwardsiterate!(v0, v1)
        upwind!(UR, v1, xs, reward, drift, policy, zerodrift)
        A = sparse(Tridiagonal(-UR.GB[2:end] ./ step(xs), (UR.GB .- UR.GF) ./ step(xs), UR.GF[1:end-1] ./ step(xs)))
        v0 .= ((ρ + 1.0 / dt) * I - A) \ (UR.R .+ v1 ./ dt)
        return v0
    end

    vinit = map(x -> reward(x, zerodrift(x)) ./ ρ, xs)
    vs = (copy(vinit), copy(vinit))
    tol, maxiter = 1e-10, 10000
    err = tol
    for iter in 1:maxiter
        backwardsiterate!(vs[2], vs[1])
        err = maximum(abs.(vs[2] .- vs[1]))
        if err < tol
            return (value = vs[2], iter = iter, err = err, termination = :tolerance)
        else
            vs[1] .= vs[2]
        end
    end
    return (value = vs[1], iter = maxiter, err = err, termination = :iterations)
end

@testset "UpwindDifferences.jl" begin

    allparameters = readdlm("testdata/hjb_parameters.tsv")
    allvalues = readdlm("testdata/hjb_values.tsv")
    @showprogress for i in axes(allparameters, 1)
        ps = allparameters[i, :]
        vnew = evaluateparameters(ps).value
        vold = allvalues[i, :]
        @test maximum(abs.(vnew .- vold)) < 1e-8
    end
end

ext: jl


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

_rad2sec{T<:AbstractFloat}(rad::T) = rad2deg(rad)*3600.0

"""
    rad2sec(rad) -> seconds

### Purpose ###

Convert from radians to seconds.

### Argument ###

* `rad`: number of radians.  It can be either a scalar or a vector.

### Output ###

The number of seconds corresponding to `rad`.  If `rad` is an array, an array of
the same length is returned.

### Example ###

``` julia
rad2sec(1)
# => 206264.80624709636
```

### Notes ###

Use `sec2rad` to convert seconds to radians.
"""
rad2sec(rad::Real) = _rad2sec(float(rad))

@vectorize_1arg Real rad2sec


text:
"""
std_setup : functions called during standard set-up of optimization. Includes CAMEOS_main() which does most of the work
			computing the tensor with appropriate values. Also functions to deal with input/output presentation of all sequences
			and handling mismatches between HMM/MRF models of sequences.
"""

module std_setup

using Logging, StatsBase, Distributions, Unicode

islower(s) = all(c->islowercase(c) | isspace(c), s);

function debug(x)
	println("DEBUG: $x")
end

function err(x)
	println("ERR: $x")
end

include("types.jl")
#include("math.jl")
include("bio_seq.jl")
include("hmm_logic.jl")
include("dyn_frwd.jl")
include("mrf.jl")
#include("lookup.jl")

function CAMEOS_main(x, x_prot, y, y_prot, frame, extra_out = false)
	start_time = time()

	X_len, X_hmm_state, X_hmm_insert, X_hmm_match = create_HMM(x)
	X_hmm = types.wholeHMM(X_hmm_state, X_hmm_insert, X_hmm_match)

	Y_len, Y_hmm_state, Y_hmm_insert, Y_hmm_match = create_HMM(y)
	Y_hmm = types.wholeHMM(Y_hmm_state, Y_hmm_insert, Y_hmm_match)

	X_len += 2
	Y_len += 2 #for the stop and since we're looking one back...

	mark_name = x_prot
	m_pd, m_nd = load_hmm_seq(mark_name, X_hmm, X_len)
	mark_prot = m_pd[mark_name]
	mark_nuc = m_nd[mark_name]

	deg_name = y_prot
	d_pd, d_nd = load_hmm_seq(deg_name, Y_hmm, Y_len)
	deg_prot = d_pd[deg_name]
	deg_nuc = d_nd[deg_name]

	min_len = min(Y_len, X_len) - 1
	x_init = follow_prot_give_log_2d(Y_hmm, deg_prot, Y_len - min_len)
	y_init = follow_prot_give_log_2d(X_hmm, mark_prot, X_len - min_len)

	all_results = fill_in_tensor(X_hmm, Y_hmm, X_len, Y_len, x_init, y_init, mark_prot, deg_prot, min_len, frame)
	if extra_out
		return all_results, X_hmm, Y_hmm, X_len, Y_len, mark_prot, deg_prot, mark_nuc, deg_nuc
	else
		return all_results
	end
end

function align_consensus(prot_seq, hmm_file)
	hmm_con = get_hmm_consensus(hmm_file)

	alignment = read(pipeline(`printf ">pro_se\n$prot_seq\n"`, `hmmalign --outformat=A2M $hmm_file - `), String)
	alignment = join(split(alignment, '\n')[2:end], "")
	alignment = fix_hmmalign(alignment, hmm_con)

	return alignment
end

function i_j_from_path(path)
	first_elem = path[end]
	last_elem = path[1]
	cur_i = first_elem[1]; cur_j = first_elem[2]
	end_i = last_elem[1]; end_j = last_elem[2]
	return cur_i, cur_j, end_i, end_j
end

#This is logic to format the output of the path through the sequence.
function format_output_map(cur_seq, frame, deg_nuc, mark_nuc, deg_map, mark_map, cur_i, cur_j, end_i, end_j, deg_len, mark_len)
	deg_trns = 1 #Translation start for essential gene.
	deg_trne = 1 #Translation end for essential gene.
	mark_trns = 1 #Same as above.
	mark_trne = 1 #Same...

	if frame == "p1"
		final_mark = ""; final_deg = ""; final_seq = ""
		if mark_map[cur_j] > 1
			final_mark *= lowercase(mark_nuc[1:((mark_map[cur_j - 1]) * 3) - 1]);
			final_seq *= lowercase(mark_nuc[1:((mark_map[cur_j - 1]) * 3) - 1])
			final_deg *= "A" * uppercase(cur_seq[2:end]);
			final_seq *= "A" * uppercase(cur_seq[2:end]);
			final_mark *= "A" * uppercase(cur_seq[2:end]) #we always start with A.
			mark_trns = 1
			mark_trne = length(final_seq) - 1
			deg_trns = length(mark_nuc[1:((mark_map[cur_j - 1]) * 3) - 1]) + 1 #Not sure about the +1...
			deg_trne = length(final_seq) - 1
		elseif deg_map[cur_i] > 1
			final_deg *= lowercase(deg_nuc[1:((deg_map[cur_i - 1]) * 3) - 2] * "At")
			final_seq *= lowercase(deg_nuc[1:((deg_map[cur_i - 1]) * 3) - 2] * "At")
			final_deg *= uppercase(cur_seq);
			final_seq *= uppercase(cur_seq);
			final_mark *= uppercase(cur_seq[2:end])
			deg_trns = 1
			deg_trne = length(final_seq) - 1
			mark_trns = length(deg_nuc[1:((deg_map[cur_i - 1]) * 3)]) + 2 #- 1 + the At stuff...
			mark_trne = length(final_seq) - 1
		else
			#I guess it's possible that both sequences are at 1.
			#I think then we just don't have anything in prefix.
			final_seq *= uppercase(cur_seq)
			final_deg *= uppercase(cur_seq)
			final_mark *= uppercase(cur_seq[2:end])

			deg_trns = 1
			mark_trns = 2
			deg_trne = length(final_seq) - 1
			mark_trne = length(final_seq) - 1
		end
		if (end_j in keys(mark_map)) && mark_map[end_j] < maximum(keys(mark_map)) # div(mark_len, 3)
			final_mark *= lowercase(mark_nuc[(mark_map[end_j] * 3) - 2 : end]);
			final_seq *= lowercase(mark_nuc[(mark_map[end_j] * 3) - 2 : end])
			mark_trne = length(final_seq)
		elseif (end_i in keys(deg_map)) && deg_map[end_i] < maximum(keys(deg_map)) #div(deg_len, 3)
			final_deg *= lowercase(deg_nuc[(deg_map[end_i] * 3) - 1 : end]);
			final_seq *= lowercase(deg_nuc[(deg_map[end_i] * 3) - 1 : end])
			deg_trne = length(final_seq)
		end
	elseif frame == "p2"
		final_mark = ""; final_deg = ""; final_seq = "";
		if mark_map[cur_j] > 1
			final_mark *= lowercase(mark_nuc[1:((mark_map[cur_j - 1]) * 3) - 1])
			final_seq *= lowercase(mark_nuc[1:((mark_map[cur_j - 1]) * 3) - 1])
			final_mark *= "A" * uppercase(cur_seq[2:end])
			final_seq *= "A" * uppercase(cur_seq[2:end])
			final_deg *= "A" * uppercase(cur_seq[2:end])
			mark_trns = 1
			mark_trne = length(final_seq) - 1
			deg_trns = length(mark_nuc[1:((mark_map[cur_j - 1]) * 3) - 1]) + 1
			deg_trne = length(final_seq) - 1
		elseif deg_map[cur_i] > 1
			final_deg *= lowercase(deg_nuc[1:((deg_map[cur_i - 1]) * 3)])
			final_seq *= lowercase(deg_nuc[1:((deg_map[cur_i - 1]) * 3)])
			final_deg *= "A" * uppercase(cur_seq[2:end])
			final_seq *= "A" * uppercase(cur_seq[2:end])
			final_mark *= uppercase(cur_seq[2:end])
			deg_trns = 1
			deg_trne = length(final_seq) - 1
			mark_trns = length(deg_nuc[1:((deg_map[cur_i - 1]) * 3)]) + 2
			mark_trne = length(final_seq) - 1
		end
		if (end_j in keys(mark_map)) && mark_map[end_j] < div(mark_len, 3)
			final_mark *= lowercase(mark_nuc[(mark_map[end_j] * 3) - 2:end])
			final_seq *= lowercase(mark_nuc[(mark_map[end_j] * 3) - 2:end])
			mark_trne = length(final_seq)
		elseif (end_i in keys(deg_map)) && deg_map[end_i] < div(deg_len, 3)
			final_deg *= lowercase(deg_nuc[(deg_map[end_i] * 3) - 1:end])
			final_seq *= lowercase(deg_nuc[(deg_map[end_i] * 3) - 1:end])
			deg_trne = length(final_seq)
		end
	end
	return final_deg, final_mark, final_seq, deg_trns, deg_trne, mark_trns, mark_trne
end

function NST_get_explicit_mapping(wt_seq, hmm_file)
	#Returns a mapping from HMM nodes to position in a protein sequence.
	prot_seq = align_consensus(wt_seq, hmm_file)
	mapping = Dict{Int64, Int64}()
	cur_prot = 0
	cur_hmm = 0
	for cha in prot_seq
		if islower(cha)
			cur_prot += 1
		elseif cha == '-'
			cur_hmm += 1
		else
			cur_prot += 1
			cur_hmm += 1
		end
		mapping[cur_hmm] = cur_prot
	end
	mapping[cur_hmm + 1] = cur_prot + 1 #for stop codon
	return mapping
end

function NST_full_general_main(x_prot, x_hmm, y_prot, y_hmm, frame, num_samples, prob_thresh, soft_temp, soft_temp_start, abs_bad, fin_x_nuc="", fin_x_prot="", fin_y_nuc="", fin_y_prot="")
	gen_bbt, X_hmm_obj, Y_hmm_obj, X_len, Y_len, x_prot_seq, y_prot_seq, mark_nuc, deg_nuc = CAMEOS_main(x_hmm, x_prot, y_hmm, y_prot, frame, true)
	println("CAMEOS tensor built")
	cumulative_score = 0.0
	failures = 0
	success_count = 0
	min_len = min(Y_len, X_len) - 1
	deg_map = Dict{Int64, Int64}()
	mark_map = Dict{Int64, Int64}()
	if fin_x_prot != ""
		deg_map = NST_get_explicit_mapping(fin_x_prot, x_hmm)
	else
		deg_map = Dict{Int64, Int64}(i => i for i in 1:X_len)
	end
	if fin_y_prot != ""
		mark_map = NST_get_explicit_mapping(fin_y_prot, y_hmm)
	else
		mark_map = Dict{Int64, Int64}(i => i for i in 1:Y_len)
	end

	all_samples = types.SampleNucs[]
	all_results = Any[]
	#best_result is deterministic run of CAMEOS.
	X_tail = tail_prot_give_log_2d(X_hmm_obj, x_prot_seq, X_len)
	Y_tail = tail_prot_give_log_2d(Y_hmm_obj, y_prot_seq, Y_len)

	best_result = sample_start_point(X_hmm_obj, x_prot_seq, X_len, X_tail, Y_hmm_obj, y_prot_seq, Y_len, Y_tail, 2, 1.0, 1.0, 1e-10, gen_bbt)
	for r in best_result
		push!(all_results, r)
		break #just add one.
	end

	results = sample_start_point(X_hmm_obj, x_prot_seq, X_len, X_tail, Y_hmm_obj, y_prot_seq, Y_len, Y_tail, num_samples, prob_thresh, soft_temp, soft_temp_start, gen_bbt)
	for r in results
		push!(all_results, r)
	end

	for res in all_results
		full_seq, path, cumulative_score, cur_score = res
		cur_i, cur_j, end_i, end_j = i_j_from_path(path)

		final_mark, final_deg, final_seq, mark_trns, mark_trne, deg_trns, deg_trne = format_output_map(full_seq, frame, fin_x_nuc, fin_y_nuc, deg_map, mark_map, cur_i, cur_j, end_i, end_j, X_len * 3, Y_len * 3)

		push!(all_samples, types.SampleNucs(final_seq, deg_trns, deg_trne, mark_trns, mark_trne, path))
	end

	return all_samples
end

function NST_full_general_main_range(x_prot, x_hmm, y_prot, y_hmm, frame, num_samples, prob_thresh, soft_temp, soft_temp_start, abs_bad, x_range, y_range, fin_x_nuc="", fin_x_prot="", fin_y_nuc="", fin_y_prot="")
	#This function only goes through a single range.
	gen_bbt, X_hmm_obj, Y_hmm_obj, X_len, Y_len, x_prot_seq, y_prot_seq, mark_nuc, deg_nuc = CAMEOS_main(x_hmm, x_prot, y_hmm, y_prot, frame, true)
	println("CAMEOS tensor built")
	cumulative_score = 0.0
	failures = 0
	success_count = 0
	min_len = min(Y_len, X_len) - 1
	deg_map = Dict{Int64, Int64}()
	mark_map = Dict{Int64, Int64}()
	if fin_x_prot != ""
		deg_map = NST_get_explicit_mapping(fin_x_prot, x_hmm)
	else
		deg_map = Dict{Int64, Int64}(i => i for i in 1:X_len)
	end
	if fin_y_prot != ""
		mark_map = NST_get_explicit_mapping(fin_y_prot, y_hmm)
	else
		mark_map = Dict{Int64, Int64}(i => i for i in 1:Y_len)
	end

	all_samples = types.SampleNucs[]
	all_results = Any[]

	results = sample_start_point_on_range(X_hmm_obj, x_prot_seq, X_len, Y_hmm_obj, y_prot_seq, Y_len, num_samples, x_range, y_range, prob_thresh, soft_temp, soft_temp_start, gen_bbt)
	for r in results
		push!(all_results, r)
	end

	for res in all_results #[1:5]
		full_seq, path, cumulative_score, cur_score = res
		cur_i, cur_j, end_i, end_j = i_j_from_path(path)
		final_mark, final_deg, final_seq, mark_trns, mark_trne, deg_trns, deg_trne = format_output_map(full_seq, frame, fin_x_nuc, fin_y_nuc, deg_map, mark_map, cur_i, cur_j, end_i, end_j, X_len * 3, Y_len * 3) #ength(fin_x_nuc), length(fin_y_nuc))
		push!(all_samples, types.SampleNucs(final_seq, deg_trns, deg_trne, mark_trns, mark_trne, path))
	end

	return all_samples #, new_all_samples
end

function produce_vectors(trns, trne, skip_item, insert_item)
	tallied_skip = Dict{Int64, Int64}()
	tallied_insert = Dict{Int64, Int64}()

	#Being cautious b/c countmap on an empty list gives an error.
	if length(skip_item) > 0
		tallied_skip = countmap(skip_item) #this looks like: {210 => 4, 228 => 2} kind of thing..
	end
	if length(insert_item) > 0
		tallied_insert = countmap(insert_item) #similarly looks like {195 => 3, 204 => 2}
	end

	ts_keys = sort(collect(keys(tallied_skip)))
	ti_keys = sort(collect(keys(tallied_insert)))

	grem_aa = 0
	bp_counter = 0
	skip_counter = 0
	skip_covered = Int64[]
	insert_covered = Int64[] #trn_info is keep_track[name] in main body

	skip_vec = Int64[]
	insert_vec = Tuple{Int64, Int64}[]
	insert_counter = 0
	bp = trns - 1

	while (bp - insert_counter) <= trne
		bp += 1
		if skip_counter > 0
			push!(skip_vec, bp - insert_counter)
			skip_counter -= 1
    		continue
		end

		if bp_counter % 3 == 0
			grem_aa = div(bp_counter, 3)
		end #So at bp_counter=0, we are at cur_aa = 1. Great.

		if (grem_aa in ts_keys) && !(grem_aa in skip_covered) #then this is going to be skipped
			push!(skip_vec, bp - insert_counter)
			push!(skip_covered, grem_aa) #so we just do this once.
			skip_counter = (3 * tallied_skip[grem_aa]) - 1 #-1 because we're currently adding our base.
			bp_counter += 0 #because we're finishing off the codon.
			continue #don't want to evaluate anything else. Don't want bp_counter to increase
		end

		if (grem_aa in ti_keys) && !(grem_aa in insert_covered) #then we're going to add one or more "..." here
			push!(insert_covered, grem_aa)
			push!(insert_vec, (bp - 1 - insert_counter, tallied_insert[grem_aa]))
			insert_counter += (3 * tallied_insert[grem_aa])
		end

		bp_counter += 1
	end

	return skip_vec, insert_vec
end

function try_run_on_SampleNucs(sample_prots, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem)
	#Same code as run_on_SampleNucs with Exception Handling to avoid cases where problems with protein translation lead to a bug.
	#Bugs are rare and tend to occur with nonsense proteins that confuse alignment via hmmalign and result in very weird sequences.
	#We don't care about these anyway, because they'll be low-scoring. We are fine throwing them away.
	mark_grem_prot = mrf.get_gremlin_consensus(mark_grem, mark_name)
	deg_grem_prot = mrf.get_gremlin_consensus(deg_grem, deg_name)
	sample_count = 1
	mark_sub_prots = Dict{Int64, AbstractString}()
	deg_sub_prots = Dict{Int64, AbstractString}()
	keep_track = Dict{Int64, Any}()
	titles = Int64[]
	deg_prot_trans = ""
	mark_prot_trans = ""
	for sample_nuc in sample_prots
		sample_name = sample_count
		try
			deg_prot_trans = bio_seq.translate(sample_nuc.final_seq[sample_nuc.deg_trns:sample_nuc.deg_trne])
			mark_prot_trans = bio_seq.translate(sample_nuc.final_seq[sample_nuc.mark_trns:sample_nuc.mark_trne])
		catch Exception
			println("Exception hit")
			continue
		end

		#If we have no problem translating...
		mark_sub_prots[sample_name] = mark_prot_trans
		deg_sub_prots[sample_name] = deg_prot_trans
		keep_track[sample_name] = sample_nuc
		push!(titles, sample_name)
		sample_count += 1
	end
	deg_skip, deg_insert = redo_alignment(align_grem(deg_grem_prot, deg_hmm), sample_prots, false)
	mark_skip, mark_insert = redo_alignment(align_grem(mark_grem_prot, mark_hmm), sample_prots, true)
	ex_sam = types.ExtendedSampleNucs[]
	for name in titles
		trn_info = keep_track[name]
		try
			deg_skip_vec, deg_insert_vec = produce_vectors(trn_info.deg_trns, trn_info.deg_trne, deg_skip[name], deg_insert[name])
			mark_skip_vec, mark_insert_vec = produce_vectors(trn_info.mark_trns, trn_info.mark_trne, mark_skip[name], mark_insert[name])
			new_sn = types.ExtendedSampleNucs(trn_info.final_seq, trn_info.deg_trns, trn_info.deg_trne, "+", trn_info.mark_trns, trn_info.mark_trne, "+", trn_info.path, deg_skip_vec, deg_insert_vec, mark_skip_vec, mark_insert_vec)
			push!(ex_sam, new_sn)
		catch ne_err
			#I don't care about an exception, just won't use that sample.
		end
	end
	return ex_sam, mark_grem_prot, deg_grem_prot
end

#This is a helper function to deal with insertions and other sequence features that we "censor" meaning that
#model isn't aware that we're making slight modifications to sequence with understanding that it won't affect downstream outputs.
function interpret_censorship(nucs, first_skip, first_insert, first_tran, second_skip, second_insert, second_tran)
	global_skip = union(first_skip, second_skip)

	first_map_local_skip = Int64[]
	for mapped_x in first_skip
		push!(first_map_local_skip, mapped_x)
	end
	first_map_global_skip = Int64[]
	for mapped_x in global_skip
		push!(first_map_global_skip, mapped_x)
	end
	mapped_first_insert = Tuple[]
	for (mapped_x, y) in first_insert
		push!(mapped_first_insert, (mapped_x, y))
	end

	if first_tran[3] == "+"
		first_nuc = types.MRF_nuc(nucs[first_tran[1]:first_tran[2]], first_map_global_skip, first_map_local_skip, mapped_first_insert)
	else
		first_nuc = types.MRF_nuc(bio_seq.rev_comp(nucs[first_tran[2]:first_tran[1]]), first_map_global_skip, first_map_local_skip, mapped_first_insert)
	end

	aa_count = 0
	first_map = Dict{Int64, Int64}()

	for c in first_tran[1]:first_tran[2]
		if !(c in first_map_local_skip)
			aa_count += 1
			first_map[c] = div(aa_count-1, 3) + 1
		end
	end

	second_map_local_skip = Int64[]
	for mapped_x in second_skip
		push!(second_map_local_skip, mapped_x)
	end
	second_map_global_skip = Int64[]
	for mapped_x in global_skip
		push!(second_map_global_skip, mapped_x)
	end
	mapped_second_insert = Tuple[]
	for (mapped_x, y) in second_insert
		push!(mapped_second_insert, (mapped_x, y))
	end

	if second_tran[3] == "+"
		second_nuc = types.MRF_nuc(nucs[second_tran[1]:second_tran[2]], second_map_global_skip, second_map_local_skip, mapped_second_insert)
	elseif second_tran[3] == "-"
		second_nuc = types.MRF_nuc(bio_seq.rev_comp(nucs[second_tran[2]:second_tran[1]]), second_map_global_skip, second_map_local_skip, mapped_second_insert)
	end

	second_map = Dict{Int64, Int64}()
	aa_count = 0
	for c in second_tran[1]:second_tran[2]
		if !(c in second_map_local_skip)
			aa_count += 1
			second_map[c] = div(aa_count-1, 3) + 1
		end
	end

	return first_nuc, first_map, second_nuc, second_map
end

function hmm_score_helper(deg_in_read)
	true_lines = Any[]
	for line in deg_in_read
		if line[1] == '#'
			continue
		else
			push!(true_lines, split(line, " ", keepempty=false))
		end
	end
	deg_hmm_scores = Dict{AbstractString, Float64}()
	for x in true_lines
		prot_id = x[3]
		score = parse(Float64, x[6])
		deg_hmm_scores[prot_id] = score
	end
	return deg_hmm_scores
end

function full_set_up(out_path, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem, pop_size, num_samples, bad_threshold, rand_barcode, frame="p1")
	my_prots = NST_full_general_main(mark_name, mark_hmm, deg_name, deg_hmm, frame, pop_size * 50, 0.90, 1.0, 1.0, 1200, gen_hmm_trace(mark_name, mark_hmm)..., gen_hmm_trace(deg_name, deg_hmm)...)

	println("Evaluating HMM seeds")
	@debug("We've generated samples.")

	sample_count = 1
	mark_sub_prots = Dict{Int64, String}()
	deg_sub_prots = Dict{Int64, String}()
	keep_track = Dict{Int64, Any}()
	titles = Int64[]
	for sample_nuc in my_prots
		sample_name = sample_count
		deg_prot_trans = bio_seq.translate(sample_nuc.final_seq[sample_nuc.deg_trns:sample_nuc.deg_trne], 1)
		mark_prot_trans = bio_seq.translate(sample_nuc.final_seq[sample_nuc.mark_trns:sample_nuc.mark_trne], 1)
		mark_sub_prots[sample_name] = mark_prot_trans
		deg_sub_prots[sample_name] = deg_prot_trans
		keep_track[sample_name] = sample_nuc
		push!(titles, sample_name)
		sample_count += 1
	end

	deg_lines = String[]
	mark_lines = String[]
	for sam in titles
		if length(strip(deg_sub_prots[sam])) > 1 && length(strip(mark_sub_prots[sam])) > 1 #hmmscan doesn't like empty lines.
			push!(deg_lines, ">sample_$(sam)\n$(deg_sub_prots[sam])\n")
			push!(mark_lines, ">sample_$(sam)\n$(mark_sub_prots[sam])\n")
		end
	end

	deg_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/deg_out_$(rand_barcode).fa", "w")
	write(deg_out_file, "$(join(deg_lines))")
	close(deg_out_file)

	mark_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/mark_out_$(rand_barcode).fa", "w")
	write(mark_out_file, "$(join(mark_lines))")
	close(mark_out_file)

	deg_hmm_run = read(pipeline(`hmmscan --cpu 1 --tblout=$out_path/$(mark_name)_$(deg_name)_$frame/deg_$(rand_barcode).tbl $deg_hmm $out_path/$(mark_name)_$(deg_name)_$frame/deg_out_$(rand_barcode).fa `), String)
	mark_hmm_run= read(pipeline(`hmmscan --cpu 1 --tblout=$out_path/$(mark_name)_$(deg_name)_$frame/mark_$(rand_barcode).tbl $mark_hmm $out_path/$(mark_name)_$(deg_name)_$frame/mark_out_$(rand_barcode).fa `), String)

	#Okay now read the results and see which samples are top quality.
	deg_in_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/deg_$(rand_barcode).tbl")
	deg_in_read = readlines(deg_in_file)
	close(deg_in_file)
	deg_hmm_scores = hmm_score_helper(deg_in_read)

	mark_in_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/mark_$(rand_barcode).tbl")
	mark_in_read = readlines(mark_in_file)
	close(mark_in_file)
	mark_hmm_scores = hmm_score_helper(mark_in_read)

	#Okay, now we just sum up the contributions of each score in deg/mark for total.
	#We will then choose the top-N scoring samples for further work.

	total_sam = Dict{Int64, Float64}()
	for sam in titles
		sam_key = "sample_$sam"
		if sam_key in keys(mark_hmm_scores) && sam_key in keys(deg_hmm_scores)
			total_score = mark_hmm_scores[sam_key] + deg_hmm_scores[sam_key]
			total_sam[sam] = total_score
		end
	end

	total_sam_num = length(collect(keys(total_sam)))

	safety_margin = round(Int64, pop_size * 1.1) #10% more than pop_size in case of random failures.
	#We sort the dictionary and only keep some of the prots.
	sorted_by_score = sort(collect(total_sam), by= x -> -x[2]) #negative so that top is first.
	elite_candidates = [sorted_by_score[i][1] for i in 1:5]

	#top_quart = sorted_by_score[1:div(total_sam_num, 4)]
	#top_half = sorted_by_score[1:div(total_sam_num, 2)]
	top_third = sorted_by_score[1:div(total_sam_num, 3)]

	top_key_pairs = sample(top_third, safety_margin, replace=false) #clearly safety_margin can never be bigger than top_half/top_quart.
	top_candidates = [top_key_pairs[i][1] for i in 1:safety_margin]

	for elite in elite_candidates
		push!(top_candidates, elite)
	end

	top_candidates = reverse(top_candidates) #to get "elites" in first.

	my_top_prots = types.SampleNucs[]
	for topper in top_candidates
		push!(my_top_prots, my_prots[topper]) #should be an int.
	end

	deg_lines = String[]
	mark_lines = String[]

	#Let's output the top samples just to be able to query its diversity.
	for iiiii in top_candidates
		push!(deg_lines, ">sample_$(iiiii) ($(total_sam[iiiii]))\n$(deg_sub_prots[iiiii])\n")
		push!(mark_lines, ">sample_$(iiiii) ($(total_sam[iiiii]))\n$(mark_sub_prots[iiiii])\n")
	end

	deg_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/deg_select_$(rand_barcode).fa", "w")
	write(deg_out_file, "$(join(deg_lines))")
	close(deg_out_file)

	mark_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/mark_select_$(rand_barcode).fa", "w")
	write(mark_out_file, "$(join(mark_lines))")
	close(mark_out_file)

	@debug("From new_sample_traceback, we get $(length(my_top_prots)) samples with which to work.")
	my_ex_prots, mark_grem_prot, deg_grem_prot = try_run_on_SampleNucs(my_top_prots, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem)
	@debug("We set up with $(length(my_ex_prots)) proteins from auto_align_grem.")

	mark_grem_model = mrf.init_model(mark_grem)
	deg_grem_model = mrf.init_model(deg_grem)

	chromosomes = types.Chromosome[]

	iss = Any[]
	deg_prot_seqs = String[]
	mark_prot_seqs = String[]

	successes = 0
	for ex in my_ex_prots
		if successes >= pop_size #this may happen if our sampling from CAMEOS is out of step with our desired population size.
			break
		end
		#try

		deg_nuc, deg_map, mark_nuc, mark_map = interpret_censorship(ex.final_seq, ex.deg_skip, ex.deg_insert, (ex.deg_trns, ex.deg_trne, ex.deg_d), ex.mark_skip, ex.mark_insert, (ex.mark_trns, ex.mark_trne, ex.mark_d))
		deg_trans = uppercase(bio_seq.translate_constrained_maybe_map(deg_nuc, ex.deg_trns, do_map = false))
		mark_trans = uppercase(bio_seq.translate_constrained_maybe_map(mark_nuc, ex.mark_trns, do_map = false))

		if length(strip(deg_trans, '*')) == deg_grem_model.nNodes && length(strip(mark_trans, '*')) == mark_grem_model.nNodes
			push!(deg_prot_seqs, deg_trans)
			push!(mark_prot_seqs, mark_trans)
			push!(iss, (ex.path, ex.final_seq, deg_nuc, deg_map, ex.deg_trns, ex.deg_trne, ex.deg_d, ex.deg_skip, ex.deg_insert, mark_nuc, mark_map, ex.mark_trns, ex.mark_trne, ex.mark_d, ex.mark_skip, ex.mark_insert))
			#I just want to write all these proteins and assess their likelihoods and stuff as a big batch.
			successes += 1
		end
	end

	@debug("We had $successes sampling successes...")

	deg_prot_mat = zeros(Int64, successes, deg_grem_model.nNodes * 21)
	mark_prot_mat = zeros(Int64, successes, mark_grem_model.nNodes * 21)
	for i in 1:successes
		deg_prot_mat[i, 1:end] = bio_seq.convert_protein(deg_prot_seqs[i][1:deg_grem_model.nNodes])
		mark_prot_mat[i, 1:end] = bio_seq.convert_protein(mark_prot_seqs[i][1:mark_grem_model.nNodes])
	end

	deg_probs, deg_ull, deg_pv_w1, deg_pv_w2 = mrf.cpu_assess(deg_prot_mat, deg_grem_model.w1, deg_grem_model.w2, deg_grem_model.nNodes, successes)
	mark_probs, mark_ull, mark_pv_w1, mark_pv_w2 = mrf.cpu_assess(mark_prot_mat, mark_grem_model.w1, mark_grem_model.w2, mark_grem_model.nNodes, successes)

	for i in 1:successes
		new_chr = types.Chromosome(iss[i]..., deg_probs[i], mark_probs[i])
		push!(chromosomes, new_chr)
	end
	return mark_grem_model, deg_grem_model, chromosomes, mark_grem_prot, deg_grem_prot
end

function full_sample_set_up(out_path, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem, pop_size, num_samples, bad_threshold, rand_barcode, mark_range, deg_range, frame="p1")
	my_sampled_prots = NST_full_general_main_range(mark_name, mark_hmm, deg_name, deg_hmm, frame, pop_size * 50, 0.90, 1.0, 1.0, 1200, mark_range, deg_range, gen_hmm_trace(mark_name, mark_hmm)..., gen_hmm_trace(deg_name, deg_hmm)...)

	sample_count = 1
	mark_sub_prots = Dict{String, String}()
	deg_sub_prots = Dict{String, String}()
	keep_track = Dict{String, Any}()
	titles = String[]
	for sample_nuc in my_sampled_prots
		sample_name = "$sample_count"
		deg_prot_trans = bio_seq.translate(sample_nuc.final_seq[sample_nuc.deg_trns: sample_nuc.deg_trne], 1)
		mark_prot_trans = bio_seq.translate(sample_nuc.final_seq[sample_nuc.mark_trns: sample_nuc.mark_trne], 1)
		mark_sub_prots[sample_name] = mark_prot_trans
		deg_sub_prots[sample_name] = deg_prot_trans
		keep_track[sample_name] = sample_nuc
		push!(titles, sample_name)
		sample_count += 1
	end

	deg_lines = String[]
	mark_lines = String[]
	for sam in titles
		if length(strip(deg_sub_prots[sam])) > 1 && length(strip(mark_sub_prots[sam])) > 1 #hmmscan doesn't like empty lines.
			push!(deg_lines, ">sample_$(sam)\n$(deg_sub_prots[sam])\n")
			push!(mark_lines, ">sample_$(sam)\n$(mark_sub_prots[sam])\n")
		end
	end

	deg_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/deg_out_$(rand_barcode).fa", "w")
	write(deg_out_file, "$(join(deg_lines))")
	close(deg_out_file)

	mark_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/mark_out_$(rand_barcode).fa", "w")
	write(mark_out_file, "$(join(mark_lines))")
	close(mark_out_file)

	deg_hmm_run = read(pipeline(`hmmscan --cpu 1 --tblout=$out_path/$(mark_name)_$(deg_name)_$frame/deg_$(rand_barcode).tbl $deg_hmm $out_path/$(mark_name)_$(deg_name)_$frame/deg_out_$(rand_barcode).fa `), String)
	mark_hmm_run= read(pipeline(`hmmscan --cpu 1 --tblout=$out_path/$(mark_name)_$(deg_name)_$frame/mark_$(rand_barcode).tbl $mark_hmm $out_path/$(mark_name)_$(deg_name)_$frame/mark_out_$(rand_barcode).fa `), String)

	#Okay now read the results and see which samples are top quality.
	deg_in_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/deg_$(rand_barcode).tbl")
	deg_in_read = readlines(deg_in_file)
	close(deg_in_file)
	deg_hmm_scores = hmm_score_helper(deg_in_read)

	mark_in_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/mark_$(rand_barcode).tbl")
	mark_in_read = readlines(mark_in_file)
	close(mark_in_file)
	mark_hmm_scores = hmm_score_helper(mark_in_read)

	total_sam = Dict{String, Float64}()
	for sam in titles
		sam_key = "sample_$sam"
		if sam_key in keys(mark_hmm_scores) && sam_key in keys(deg_hmm_scores)
			total_score = mark_hmm_scores[sam_key] + deg_hmm_scores[sam_key]
			total_sam[sam] = total_score
		end
	end

	total_sam_num = length(collect(keys(total_sam)))

	safety_margin = round(Int64, pop_size * 1.1) #10% more than pop_size in case of random failures.
	#Okay, well now let's sort the dictionary and only keep some of the prots!
	sorted_by_score = sort(collect(total_sam), by= x -> -x[2]) #negative so that top is first.

	num_from_elite = 5

	num_from_general_sample = safety_margin
	num_from_range_sample = 0

	top_third = sorted_by_score[1:div(total_sam_num, 3)]
	top_key_pairs = sample(top_third, num_from_general_sample, replace=false) #clearly safety_margin can never be bigger than top_half/top_quart.
	top_candidates = [top_key_pairs[i][1] for i in 1:num_from_general_sample]

	elite_candidates = [sorted_by_score[i][1] for i in 1:num_from_elite]

	for elite in elite_candidates
		push!(top_candidates, elite)
	end

	top_candidates = reverse(top_candidates) #to get elites in first.

	my_top_prots = types.SampleNucs[]
	for topper in top_candidates
		push!(my_top_prots, keep_track[topper]) #should be an int.
	end

	deg_lines = String[]
	mark_lines = String[]

	#Let's output the top samples just to be able to query its diversity.
	for iiiii in top_candidates
		push!(deg_lines, ">sample_$(iiiii) ($(total_sam[iiiii]))\n$(deg_sub_prots[iiiii])\n")
		push!(mark_lines, ">sample_$(iiiii) ($(total_sam[iiiii]))\n$(mark_sub_prots[iiiii])\n")
	end

	deg_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/deg_select_$(rand_barcode).fa", "w")
	write(deg_out_file, "$(join(deg_lines))")
	close(deg_out_file)

	mark_out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/mark_select_$(rand_barcode).fa", "w")
	write(mark_out_file, "$(join(mark_lines))")
	close(mark_out_file)

	@debug("From new_sample_traceback, we get $(length(my_top_prots)) samples with which to work.")
	my_ex_prots, mark_grem_prot, deg_grem_prot = try_run_on_SampleNucs(my_top_prots, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem)
	@debug("We set up with $(length(my_ex_prots)) proteins from auto_align_grem.")

	mark_grem_model = mrf.init_model(mark_grem)
	deg_grem_model = mrf.init_model(deg_grem) #MAKE SURE THIS IS FLOAT32.

	chromosomes = types.Chromosome[]

	iss = Any[]
	deg_prot_seqs = String[]
	mark_prot_seqs = String[]

	successes = 0
	for ex in my_ex_prots
		if successes >= pop_size #this may happen if our sampling from CAMEOS is out of step with our desired population size.
			break
		end
		deg_nuc, deg_map, mark_nuc, mark_map = interpret_censorship(ex.final_seq, ex.deg_skip, ex.deg_insert, (ex.deg_trns, ex.deg_trne, ex.deg_d), ex.mark_skip, ex.mark_insert, (ex.mark_trns, ex.mark_trne, ex.mark_d))
		deg_trans = uppercase(bio_seq.translate_constrained_maybe_map(deg_nuc, ex.deg_trns, do_map = false))
		mark_trans = uppercase(bio_seq.translate_constrained_maybe_map(mark_nuc, ex.mark_trns, do_map = false))

		if length(strip(deg_trans, '*')) == deg_grem_model.nNodes && length(strip(mark_trans, '*')) == mark_grem_model.nNodes
			push!(deg_prot_seqs, deg_trans)
			push!(mark_prot_seqs, mark_trans)
			push!(iss, (ex.path, ex.final_seq, deg_nuc, deg_map, ex.deg_trns, ex.deg_trne, ex.deg_d, ex.deg_skip, ex.deg_insert, mark_nuc, mark_map, ex.mark_trns, ex.mark_trne, ex.mark_d, ex.mark_skip, ex.mark_insert))
			#I just want to write all these proteins and assess their likelihoods and stuff as a big batch.
			successes += 1
		end
	end

	@debug("We had $successes sampling successes...")

	deg_prot_mat = zeros(Int64, successes, deg_grem_model.nNodes * 21)
	mark_prot_mat = zeros(Int64, successes, mark_grem_model.nNodes * 21)
	for i in 1:successes
		deg_prot_mat[i, 1:end] = bio_seq.convert_protein(deg_prot_seqs[i][1:deg_grem_model.nNodes])
		mark_prot_mat[i, 1:end] = bio_seq.convert_protein(mark_prot_seqs[i][1:mark_grem_model.nNodes])
	end

	deg_probs, deg_ull, deg_pv_w1, deg_pv_w2 = mrf.cpu_assess(deg_prot_mat, deg_grem_model.w1, deg_grem_model.w2, deg_grem_model.nNodes, successes)
	mark_probs, mark_ull, mark_pv_w1, mark_pv_w2 = mrf.cpu_assess(mark_prot_mat, mark_grem_model.w1, mark_grem_model.w2, mark_grem_model.nNodes, successes)

	for i in 1:successes
		#I can just splat the tuple with "..."
		new_chr = types.Chromosome(iss[i]..., deg_probs[i], mark_probs[i])
		push!(chromosomes, new_chr)
	end
	return mark_grem_model, deg_grem_model, chromosomes, mark_grem_prot, deg_grem_prot
end

function short_set_up(mark_name, deg_name, mark_grem, deg_grem, mark_hmm, deg_hmm, pop_size)
	mark_grem_model = mrf.init_model(mark_grem)
	deg_grem_model = mrf.init_model(deg_grem)

	seqs = bio_seq.load_fasta("gremlin_consensus.fa")
	return mark_grem_model, deg_grem_model
end

function get_explicit_mapping(grem_seq, hmm_file)
	prot_seq = align_consensus(grem_seq, hmm_file)
	#prot_seq = join(split(align, "\n")[2:end])
	mapping = Dict{Int64, Int64}()
	cur_grem = 0
	cur_hmm = 0
	for cha in prot_seq
		if islower(cha)
			cur_grem += 1
		elseif cha == '-'
			cur_hmm += 1
		else
			cur_grem += 1
			cur_hmm += 1
		end
		mapping[cur_hmm] = cur_grem
	end
	mapping[cur_hmm + 1] = cur_grem #for stop codon, but we just map to last character.
	mapping[cur_hmm + 2] = cur_grem
	return mapping
end

end


text:
"""
```
imgb = binarize(Balanced(), img)
binarize!(Balanced(), img)
```

By balancing both side of the images histogram the threshold is iterativly
updated to balance both by removing a bin untill sides untill converge to a value.

# Arguments

The function arguments are described in more detail below.

##  `img`

An `AbstractGray` image which is to be binarized into background (`Gray(0)`)
and foreground (`Gray(1)`) pixels.

## Reference

"BI-LEVEL IMAGE THRESHOLDING - A Fast Method", Proceedings of the First International Conference on Bio-inspired Systems and Signal Processing, 2008. Available: [10.5220/0001064300700076](https://doi.org/10.5220/0001064300700076)
"""
function binarize(algorithm::Balanced,  img::AbstractArray{T,2}) where T <: AbstractGray
  binarize!(algorithm, copy(img))
end

"""
```
t = binarize!(Balanced(), img)
```
Same as [`binarize`](@ref binarize(::Balanced, ::AbstractArray{<:AbstractGray,2})) except that it modifies the image that was passed as an argument.
"""
function binarize!(algorithm::Balanced,  img::AbstractArray{T,2}) where T <: AbstractGray
  edges, counts = build_histogram(img,  256)
  t = find_threshold(HistogramThresholding.Balanced(), counts[1:end], edges)
  map!(img,img) do val
    val < t ? Gray(0) : Gray(1)
  end
end


text:
using PyPlot
using Statistics

Statistics.mean(x::Real,y::Real) = (x+y)/2

x = [0.5,1.0,2.0]
xint = [mean(x[i],x[i+1]) for i in 1:length(x)-1]

fig,(ax1,ax2,ax3,ax4,ax5) = subplots(5,1,figsize=(5,5),sharex=true)
axs = [ax1,ax2,ax3,ax4,ax5]

# BITSHAVE
for i in 1:length(x)
    ax1.scatter(x[i],0,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax1.plot([0,x[1]],[0,0],"C5")
ax1.plot([x[1],x[2]],[0,0],"C0")
ax1.plot([x[2],x[3]],[0,0],"C1")
ax1.plot([x[3],3],[0,0],"C2")
ax1.set_yticks([])
ax1.set_title("Bitshave (round-to-zero)",loc="left")

ax1.plot([1.2,1.2],[-0.2,0.2],"k--",lw=0.5)
ax1.text(1.22,-0.2,L"x is round to $x_2$")

# BITSET
for i in 1:length(x)-1
    ax2.scatter(x[i+1]-0.02,0,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax2.plot([0,x[1]],[0,0],"C5")
ax2.plot([x[1],x[2]],[0,0],"C0")
ax2.plot([x[2],x[3]],[0,0],"C1")
ax2.plot([x[3],3],[0,0],"C2")
ax2.set_yticks([])
ax2.set_title("Bitset (round-away-from-zero)",loc="left")

ax2.plot([1.2,1.2],[-0.2,0.2],"k--",lw=0.5)
ax2.text(1.22,-0.2,L"x is round to $x_3^-$")

# HALFSHAVE
for i in 1:length(xint)
    ax3.scatter(xint[i],0,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax3.plot([0,x[1]],[0,0],"C5")
ax3.plot([x[1],x[2]],[0,0],"C0")
ax3.plot([x[2],x[3]],[0,0],"C1")
ax3.plot([x[3],3],[0,0],"C2")
ax3.set_yticks([])
ax3.set_title("Halfshave",loc="left")

ax3.plot([1.2,1.2],[-0.2,0.2],"k--",lw=0.5)
ax3.text(1.22,-0.2,L"x is round to $x_2^*$")

# ROUND NEAREST
for i in 1:length(x)
    ax4.scatter(x[i],0,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax4.plot([0,xint[1]],[0,0],"C0")
ax4.plot([xint[1],xint[2]],[0,0],"C1")
ax4.plot([xint[2],3],[0,0],"C2")
ax4.set_yticks([])
ax4.set_title("Round-to-nearest",loc="left")

ax4.plot([1.2,1.2],[-0.2,0.2],"k--",lw=0.5)
ax4.text(1.22,-0.2,L"x is round to $x_2$")

# STOCHASTIC ROUNDING
ax5.plot([-0.5,0,x[1]],[0,1,0],"C5")
ax5.plot([0,x[1],x[2]],[0,1,0],"C0")
ax5.plot([x[1],x[2],x[3]],[0,1,0],"C1")
ax5.plot([x[2],x[3],4],[0,1,0],"C2")
ax5.plot([x[3],4,5],[0,1,0],"C4")

for i in 1:length(x)
    ax5.scatter(x[i],1,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax5.plot([1.2,1.2],[0,1],"k--",lw=0.5)
ax5.text(1.18,-0.3,L"x = ")
ax5.text(1.31,-0.15,L"$x_2$ at 80%")
ax5.text(1.31,-0.45,L"$x_3$ at 20% chance.")

ax5.set_xlim(x[1]-0.2,x[3]+0.2)
ax5.set_yticks([0,1])
ax5.set_ylim([-0.1,1.1])
ax5.set_yticklabels([L"0\%",L"100\%"])
ax5.set_xticks(x)
ax5.set_xticklabels([L"$x_1$",L"$x_2$",L"$x_3$"])
ax5.set_title("Stochastic rounding",loc="left")

for (iax,ax) in enumerate(axs)
    ax.set_title(string(Char(iax+96)),loc="right",fontweight="bold")
    ax.set_frame_on(false)
end

tight_layout()
savefig("figs/schematic_all.pdf")
savefig("figs/schematic_all.png",dpi=200)


text:
using NNlib: conv, crosscor, ∇conv_filter, ∇conv_data, ∇maxpool, maxpool, depthwiseconv, ∇depthwiseconv_filter, ∇depthwiseconv_data

@testset "conv2d" begin
    x = reshape(Float64[1:20;], 5, 4, 1, 1)
    w = reshape(Float64[1:4;], 2, 2, 1, 1)

    @test dropdims(conv(x, w), dims = (3,4)) == [
        29 79 129;
        39 89 139;
        49 99 149;
        59 109 159.]

    @test dropdims(conv(view(x, :, :, :, :), w), dims = (3,4)) == [
        29 79 129;
        39 89 139;
        49 99 149;
        59 109 159.]

    @test dropdims(crosscor(x, w), dims = (3,4)) == [
	 51  101  151;
	 61  111  161;
 	 71  121  171;
	 81  131  181.]

    @test dropdims(conv(Float32.(x), Float32.(w)), dims=(3,4)) == Float32.([
        29 79 129;
        39 89 139;
        49 99 149;
        59 109 159.])

    @test dropdims(conv(x, w; stride=2), dims = (3,4)) == [
        29 129;
        49 149.]

    @test dropdims(conv(Float32.(x), Float32.(w); stride=2), dims = (3,4)) == Float32.([
        29 129;
        49 149.])

    @test dropdims(conv(x, w; pad=1), dims = (3,4)) == [
        1.0   9.0   29.0   49.0   48.0;
        4.0  29.0   79.0  129.0  115.0;
        7.0  39.0   89.0  139.0  122.0;
        10.0  49.0   99.0  149.0  129.0;
        13.0  59.0  109.0  159.0  136.0;
        10.0  40.0   70.0  100.0   80.0
    ]

    @test dropdims(conv(Float32.(x), Float32.(w); pad=1), dims = (3,4)) == Float32.([
        1.0   9.0   29.0   49.0   48.0;
        4.0  29.0   79.0  129.0  115.0;
        7.0  39.0   89.0  139.0  122.0;
        10.0  49.0   99.0  149.0  129.0;
        13.0  59.0  109.0  159.0  136.0;
        10.0  40.0   70.0  100.0   80.0
    ])

    @test dropdims(conv(x, w; dilation=2), dims = (3,4)) == [
        48 98;
        58 108;
        68 118.]

	# NaN tests for dilation forward pass

	ys = []
	for idx in 1:1000
    	push!(ys, conv(x, w; dilation=2))
	end
	@test !any([any(isnan.(ys[idx])) for idx in 1:1000])

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed convolution code won't change often)

    @test size(∇conv_filter(reshape(rand(4,3), 4, 3, 1, 1), x)) == size(w)
    @test size(∇conv_data(reshape(rand(4,3), 4, 3, 1, 1), w)) == size(x)

    # Test that stride/pad work backward as well
    y = conv(x, w; stride=2, pad=1, dilation=2)
    @test size(y) == (3, 2, 1, 1)
    @test size(∇conv_filter(y, x; size=size(w), stride=2, pad=1, dilation=2)) == size(w)
    @test size(∇conv_data(y, w; size=size(x), stride=2, pad=1, dilation=2)) == size(x)

	# NaN tests for dilation backward pass: filters
	dy = randn(size(ys[1]))
	dws = []
	for idx in 1:1000
	    push!(dws, ∇conv_filter(dy, x; size=size(w), dilation=2))
	end

	# NaN tests for dilation backward pass: input
	dxs = []
	for idx in 1:1000
	    push!(dxs, ∇conv_data(dy, w; size=size(x), dilation=2))
	end

	@test !any([any(isnan.(dws[idx])) for idx in 1:1000])
	@test !any([any(isnan.(dxs[idx])) for idx in 1:1000])

end

@testset "depthwiseconv2d" begin
    x = reshape(Float64[1:18;], 3, 3, 2, 1)
    w = reshape(Float64[1:16;], 2, 2, 2, 2)

    @test depthwiseconv(x, w)[:] == [23.0, 33.0, 53.0, 63.0, 71.0, 97.0, 149.0, 175.0, 497.0, 539.0, 623.0, 665.0, 689.0, 747.0, 863.0, 921.0]

    @test depthwiseconv(x, w, stride = 2, pad = 1)[:] == [1.0, 7.0, 19.0, 63.0, 5.0, 27.0, 63.0, 175.0, 90.0, 218.0, 287.0, 665.0, 130.0, 310.0, 403.0, 921.0]

    @test depthwiseconv(x, w, stride = 2)[:] == [23.0, 71.0, 497.0, 689.0]

    @test depthwiseconv(x, w, pad = 1)[:] == [1.0, 4.0, 7.0, 6.0, 7.0, 23.0, 33.0, 24.0, 19.0, 53.0, 63.0, 42.0, 21.0, 52.0, 59.0, 36.0, 5.0, 16.0, 27.0, 18.0, 27.0, 71.0, 97.0, 60.0, 63.0, 149.0, 175.0, 102.0, 49.0, 112.0, 127.0, 72.0, 90.0, 199.0, 218.0, 120.0, 227.0, 497.0, 539.0, 294.0, 287.0, 623.0, 665.0, 360.0, 176.0, 379.0, 402.0, 216.0, 130.0, 283.0, 310.0, 168.0, 319.0, 689.0, 747.0, 402.0, 403.0, 863.0, 921.0, 492.0, 240.0, 511.0, 542.0, 288.0]

    # the correctness of the gradients are being verified by calling
    # the corresponding counvolution gradients

    dy = reshape(Float64[1:16;], 2,2,4,1)
    local z = ∇depthwiseconv_data(dy,x,w)
    for i in 1:2
        X = copy(x[:,:,i:i,:]);
        W = copy(permutedims(w[:,:,:,i:i],[1,2,4,3]));
        DY = copy(dy[:,:,2i-1:2i,:]);
        res = ∇conv_data(DY,W;size=size(X))
        @test dropdims(z[:,:,i:i,:], dims=(3,4)) == dropdims(res, dims=(3,4))
    end

    z = ∇depthwiseconv_filter(dy, x, w)
    for i in 1:2
        X = copy(x[:,:,i:i,:]);
        W = copy(permutedims(w[:,:,:,i:i],[1,2,4,3]))
        DY = copy(dy[:,:,2i-1:2i,:])
        res = ∇conv_filter(DY,X; size=size(W))
        @test dropdims(z[:,:,:,i:i]; dims=(4)) == dropdims(res; dims=(3))
    end

    @test size(∇depthwiseconv_filter(rand(2,2,4,1), x, w)) == size(w)
    @test size(∇depthwiseconv_data(rand(2,2,4,1), x, w)) == size(x)

    # Test for the stride/pad for backward pass
    y = depthwiseconv(x,w,stride=2,pad=1)
    @test size(y) == (2,2,4,1)
    @test size(∇depthwiseconv_filter(rand(Float64, size(y)), x, w, stride=2, pad=1)) == size(w)
    @test size(∇depthwiseconv_data(rand(Float64, size(y)), x, w, stride=2, pad=1)) == size(x)
end

@testset "maxpool2d" begin

    x = reshape(Float64[1:20;], 5, 4, 1, 1)

    @test dropdims(maxpool(x, (2,2)), dims = (3,4)) == [7 17; 9 19]
    @test dropdims(maxpool(x, (2,2); stride=(2,2)), dims = (3,4)) == [7 17; 9 19]
    @test dropdims(maxpool(x, (2,2); pad=(1,1)), dims = (3,4)) == [
        1.0  11.0  16.0;
        3.0  13.0  18.0;
        5.0  15.0  20.0;
    ]

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed maxpooling code won't change often)

    y = maxpool(x, (2,2))
    dy = reshape(rand(2,2), 2, 2, 1, 1)
    @test size(∇maxpool(dy, y, x, (2,2))) == size(x)

end

@testset "conv3d" begin

    x = reshape(Float64[1:60;], 5, 4, 3, 1, 1)
    w = reshape(Float64[1:8;], 2, 2, 2, 1, 1)
    res = zeros(4,3,2)
    res[:, :, 1] = [
        322.0  502.0  682.0;
        358.0  538.0  718.0;
        394.0  574.0  754.0;
        430.0  610.0  790.0]
    res[:, :, 2] = [
        1042.0  1222.0  1402.0;
        1078.0  1258.0  1438.0;
        1114.0  1294.0  1474.0;
        1150.0  1330.0  1510.0]
    @test dropdims(conv(x, w), dims = (4,5)) == res

    @test dropdims(conv(Float32.(x), Float32.(w)), dims = (4,5)) == Float32.(res)

    @test dropdims(conv(x, w; stride=2), dims = (3,4,5)) == [
        322.0 682.0;
        394.0 754.0]

    @test dropdims(conv(Float32.(x), Float32.(w); stride=2), dims = (3,4,5)) == Float32.([
        322.0 682.0;
        394.0 754.0])

    res = zeros(6,5,4)
    res[:, :, 1] = [
        1.0   9.0   29.0   49.0   48.0;
        4.0  29.0   79.0  129.0  115.0;
        7.0  39.0   89.0  139.0  122.0;
        10.0  49.0   99.0  149.0  129.0;
        13.0  59.0  109.0  159.0  136.0;
        10.0  40.0   70.0  100.0   80.0]
    res[:, :, 2] = [
        26.0  126.0  206.0  286.0  220.0;
        80.0  322.0  502.0  682.0  502.0;
        94.0  358.0  538.0  718.0  524.0;
        108.0  394.0  574.0  754.0  546.0;
        122.0  430.0  610.0  790.0  568.0;
        80.0  260.0  360.0  460.0  320.0]
    res[:, :, 3] = [
        146.0   446.0   526.0   606.0   420.0;
        360.0  1042.0  1222.0  1402.0   942.0;
        374.0  1078.0  1258.0  1438.0   964.0;
        388.0  1114.0  1294.0  1474.0   986.0;
        402.0  1150.0  1330.0  1510.0  1008.0;
        240.0   660.0   760.0   860.0   560.0]
    res[:, :, 4] = [
        205.0   517.0   577.0   637.0  392.0;
        456.0  1133.0  1263.0  1393.0  847.0;
        467.0  1159.0  1289.0  1419.0  862.0;
        478.0  1185.0  1315.0  1445.0  877.0;
        489.0  1211.0  1341.0  1471.0  892.0;
        270.0   660.0   730.0   800.0  480.0]
    @test dropdims(conv(x, w; pad=1), dims = (4,5)) == res

    @test dropdims(conv(Float32.(x), Float32.(w); pad=1), dims = (4,5)) == Float32.(res)

    @test dropdims(conv(x, w; dilation=2), dims = (3,4,5)) == [
        608 788;
        644 824;
        680 860.
    ]

    @test dropdims(conv(Float32.(x), Float32.(w); dilation=2), dims = (3,4,5)) == Float32.([
        608 788;
        644 824;
        680 860.
    ])
	# NaN tests for dilation forward pass

	ys = []
	for idx in 1:1000
    	push!(ys, conv(x, w; dilation=2))
	end
	@test !any([any(isnan.(ys[idx])) for idx in 1:1000])

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed convolution code won't change often)

    @test size(∇conv_filter(reshape(rand(4,3,2), 4, 3, 2, 1, 1), x; size=size(w))) == size(w)
    @test size(∇conv_data(reshape(rand(4,3,2), 4, 3, 2, 1, 1), w; size=size(x))) == size(x)

	# NaN tests for dilation backward pass: filters
	dy = randn(size(ys[1]))
	dws = []
	for idx in 1:1000
	    push!(dws, ∇conv_filter(dy, x; size=size(w), dilation=2))
	end

	# NaN tests for dilation backward pass: input
	dxs = []
	for idx in 1:1000
	    push!(dxs, ∇conv_data(dy, w; size=size(x), dilation=2))
	end

	@test !any([any(isnan.(dws[idx])) for idx in 1:1000])
	@test !any([any(isnan.(dxs[idx])) for idx in 1:1000])

end

@testset "maxpool3d" begin

    x = reshape(Float64[1:60;], 5, 4, 3, 1, 1)

    @test dropdims(maxpool(x, (2,2,2)), dims = (3,4,5)) == [27 37; 29 39.]
    @test dropdims(maxpool(x, (2,2,2); stride=(2,2,2)), dims = (3,4,5)) == [27 37; 29 39.]
    res = zeros(3,3,2)
    res[:, :, 1] = [
        1.0  11.0  16.0;
        3.0  13.0  18.0;
        5.0  15.0  20.0]
    res[:, :, 2] = [
        41.0  51.0  56.0;
        43.0  53.0  58.0;
        45.0  55.0  60.0]
    @test dropdims(maxpool(x, (2,2,2), pad=(1,1,1)), dims = (4,5)) == res

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed maxpooling code won't change often)

    y = maxpool(x, (2,2,2))
    dy = reshape(rand(2,2), 2, 2, 1, 1, 1)
    @test size(∇maxpool(dy, y, x, (2,2,2))) == size(x)

end


text:
export iauRz
"""
Rotate an r-matrix about the z-axis.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   psi    double          angle (radians)

Given and returned:
   r      double[3][3]    r-matrix, rotated

Notes:

   1. Calling this function with positive psi incorporates in the
   supplied r-matrix r an additional rotation, about the z-axis,
   anticlockwise as seen looking towards the origin from positive z.

   2. The additional rotation can be represented by this matrix:

      (  + cos(psi)   + sin(psi)     0  )
      (                                 )
      (  - sin(psi)   + cos(psi)     0  )
      (                                 )
      (       0            0         1  )

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauRz(double psi, double r[3][3])
function iauRz(psi::Real, r::AbstractMatrix{<:Real})

   # Transpose matrix upfront
   r = Array{Float64, 2}(r') # Transpose input up front

   ccall((:iauRz, libsofa_c), Cvoid,
           (Cdouble, Ptr{Cdouble}),
           psi, pointer(r))

   return convert(Array{Float64, 2}, r')
end


text:
# adapted from Winston.jl

using Cairo

# function to draw marker symbols

const symbol_funcs = Dict(
    "asterisk" => (c, x, y, r) -> (
        move_to(c, x, y+r);
        line_to(c, x, y-r);
        move_to(c, x+0.866r, y-0.5r);
        line_to(c, x-0.866r, y+0.5r);
        move_to(c, x+0.866r, y+0.5r);
        line_to(c, x-0.866r, y-0.5r)
    ),
    "circle" => (c, x, y, r) -> (
        new_sub_path(c);
        circle(c, x, y, r)
    ),
    "cross" => (c, x, y, r) -> (
        move_to(c, x+r, y+r);
        line_to(c, x-r, y-r);
        move_to(c, x+r, y-r);
        line_to(c, x-r, y+r)
    ),
    "diamond" => (c, x, y, r) -> (
        move_to(c, x, y+r);
        line_to(c, x+r, y);
        line_to(c, x, y-r);
        line_to(c, x-r, y);
        close_path(c)
    ),
    "dot" => (c, x, y, r) -> (
        new_sub_path(c);
        rectangle(c, x-r, y-r, 2r, 2r)
    ),
    "plus" => (c, x, y, r) -> (
        move_to(c, x+r, y);
        line_to(c, x-r, y);
        move_to(c, x, y+r);
        line_to(c, x, y-r)
    ),
    "square" => (c, x, y, r) -> (
        new_sub_path(c);
        rectangle(c, x-0.866r, y-0.866r, 1.732r, 1.732r)
    ),
    "triangle" => (c, x, y, r) -> (
        move_to(c, x, y+r);
        line_to(c, x+0.866r, y-0.5r);
        line_to(c, x-0.866r, y-0.5r);
        close_path(c)
    ),
    "down-triangle" => (c, x, y, r) -> (
        move_to(c, x, y-r);
        line_to(c, x+0.866r, y+0.5r);
        line_to(c, x-0.866r, y+0.5r);
        close_path(c)
    ),
    "right-triangle" => (c, x, y, r) -> (
        move_to(c, x+r, y);
        line_to(c, x-0.5r, y+0.866r);
        line_to(c, x-0.5r, y-0.866r);
        close_path(c)
    ),
    "left-triangle" => (c, x, y, r) -> (
        move_to(c, x-r, y);
        line_to(c, x+0.5r, y+0.866r);
        line_to(c, x+0.5r, y-0.866r);
        close_path(c)
    ),
)

# function symbols(cc::Cairo.CairoContext, x, y;
#                  kind = "circle", size = 1, line_width = 0.1,
#                  line_rgba = (0,0,0,1), fill_rgba = nothing)

#     symbol_func = get(symbol_funcs, kind, "circle")

#     Cairo.save(cc)
#     Cairo.set_line_width(cc, line_width)
#     Cairo.set_dash(cc, Float64[])
#     Cairo.new_path(cc)
#     for i = 1:min(length(x),length(y))
#         symbol_func(cc, x[i], y[i], size)
#     end
#     if fill_rgba != nothing
#         Cairo.set_source_rgba(cc, fill_rgba...)
#         Cairo.fill_preserve(cc)
#     end
#     Cairo.set_source_rgba(cc, line_rgba...)
#     Cairo.stroke(cc)
#     Cairo.restore(cc)
# end

# function curve(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray;
#                line_width = 1,
#                line_rgba = (0,0,0,1))

#     n = min(length(x), length(y))
#     n > 0 || return

#     Cairo.save(cc)
#     Cairo.set_line_width(cc, line_width)
#     Cairo.set_source_rgba(cc, line_rgba...)
#     new_path(cc)

#     lo = 1
#     while lo < n
#         while lo <= n && !(isfinite(x[lo]) && isfinite(y[lo]))
#             lo += 1
#         end

#         hi = lo + 1
#         while hi <= n &&  (isfinite(x[hi]) && isfinite(y[hi]))
#             hi += 1
#         end
#         hi -= 1

#         if lo < hi
#             Cairo.move_to(cc, x[lo], y[lo])
#             for i = (lo+1):hi
#                 Cairo.line_to(cc, x[i], y[i])
#                 if i < hi && (i & 127) == 0
#                     Cairo.stroke(cc)
#                     Cairo.move_to(cc, x[i], y[i])
#                 end
#             end
#             Cairo.stroke(cc)
#         end

#         lo = hi + 1
#     end
#     Cairo.restore(cc)
# end

# function bars(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray;
#               line_width = 1,
#               line_rgba = (0,0,0,1))

#     Cairo.save(cc)
#     Cairo.set_line_width(cc, line_width)
#     Cairo.set_source_rgba(cc, line_rgba...)
#     new_path(cc)

#     for i = 1:min(length(x),length(y))
#         Cairo.move_to(cc, x[i], 0)
#         Cairo.line_to(cc, x[i], y[i])
#     end

#     Cairo.stroke(cc)
#     Cairo.restore(cc)
# end

# ---- implement layer geometries

function points(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray;
                alpha = 1.0, color = "black", fill = "red", shape = "circle", size = 2, linewidth = 0.1)

    @assert length(x) == length(y) "x and y must be of same length: $(length(x)) != $(length(y))"
    #info("points: x = $x, y = $y")
    alphas = along(alpha, x)
    colors = along(color, x)
    fills = along(fill, x)
    shapes = along(shape, x)
    sizes = along(size, x)
    #info("fills = $fills")
    Cairo.save(cc)
    Cairo.set_line_width(cc, linewidth)
    Cairo.set_dash(cc, Float64[])
    Cairo.new_path(cc)

    for i = 1:length(x)
        symbol_funcs[shapes[i]](cc, x[i], y[i], sizes[i])
        if fills[i] != nothing
            Cairo.set_source_rgba(cc, rgba(fills[i], alphas[i])...)
            Cairo.fill_preserve(cc)
        end
        Cairo.set_source_rgba(cc, rgba(colors[i], alphas[i])...)
        Cairo.stroke(cc)
    end
    Cairo.restore(cc)

end

function lines(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray;
               alpha = 1.0, color = "black", linetype = "dash", size = 0.1)

    @assert length(x) == length(y) "x and y must be of same length: $(length(x)) != $(length(y))"
    n = length(x)

    Cairo.save(cc)
    Cairo.set_line_width(cc, size)
    Cairo.set_source_rgba(cc, rgba(color, alpha)...)
    Cairo.set_line_type(cc, linetype)

    new_path(cc)

    lo = 1
    while lo < n
        while lo <= n && !(isfinite(x[lo]) && isfinite(y[lo]))
            lo += 1
        end

        hi = lo + 1
        while hi <= n &&  (isfinite(x[hi]) && isfinite(y[hi]))
            hi += 1
        end
        hi -= 1

        if lo < hi
            Cairo.move_to(cc, x[lo], y[lo])
            for i = (lo+1):hi
                Cairo.line_to(cc, x[i], y[i])
                if i < hi && (i & 127) == 0
                    Cairo.stroke(cc)
                    Cairo.move_to(cc, x[i], y[i])
                end
            end
            Cairo.stroke(cc)
        end

        lo = hi + 1
    end
    Cairo.restore(cc)
end

function cols(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray, basey;
              alpha = 1.0, color = "black", linetype = "dash", size = 0.1)

    @assert length(x) == length(y) "x and y must be of same length: $(length(x)) != $(length(y))"
    alphas = along(alpha, x)
    colors = along(color, x)
    linetypes = along(linetype, x)
    sizes = along(size, x)

    Cairo.save(cc)
    new_path(cc)

    for i = 1:length(x)
        Cairo.set_source_rgba(cc, rgba(colors[i], alphas[i])...)
        Cairo.set_line_width(cc, sizes[i])
        Cairo.set_line_type(cc, linetypes[i])
        Cairo.move_to(cc, x[i], basey)
        Cairo.line_to(cc, x[i], y[i])
        Cairo.stroke(cc)
    end

    Cairo.restore(cc)

end

function along(x, y::AbstractArray)
    if ~(typeof(x) <: AbstractArray)
        return [x for i in 1:length(y)]
    else
        @assert length(x) == length(y) "arrays of incompatible length"
        x
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Icestorm"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/YosysHQ/icestorm.git", "c495861c19bd0976c88d4964f912abe76f3901c3")
]

dependencies = Dependency[
]

# Bash recipe for building across all platforms
script = raw"""
cd icestorm
make PREFIX=${prefix} ICEPROG=0 -j${nproc}
make install PREFIX=${prefix} ICEPROG=0
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> !Sys.iswindows(p), supported_platforms(;experimental=true))
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("iceunpack", :iceunpack),
    ExecutableProduct("icebox_diff", :icebox_diff),
    ExecutableProduct("icebox_stat", :icebox_stat),
    ExecutableProduct("icebox_colbuf", :icebox_colbuf),
    ExecutableProduct("icebox_maps", :icebox_maps),
    ExecutableProduct("icebox_explain", :icebox_explain),
    ExecutableProduct("icebox_chipdb", :icebox_chipdb),
    ExecutableProduct("icepll", :icepll),
    ExecutableProduct("icebox_html", :icebox_html),
    ExecutableProduct("icemulti", :icemulti),
    ExecutableProduct("icebox_asc2hlc", :icebox_asc2hlc),
    ExecutableProduct("icebox_hlc2asc", :icebox_hlc2asc),
    ExecutableProduct("icebox_vlog", :icebox_vlog),
    ExecutableProduct("icebram", :icebram),
    ExecutableProduct("icepack", :icepack),
    ExecutableProduct("icetime", :icetime)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _radec(ra::T, dec::T, hours::Bool) where {T<:AbstractFloat}
    # Compute right ascension.
    if hours
        ra_hr, ra_min, ra_sec = sixty(mod(ra, 24))
    else
        ra_hr, ra_min, ra_sec = sixty(mod(ra, 360) / 15)
    end
    # Compute declination.
    dec_deg, dec_min, dec_sec = sixty(dec)
    return ra_hr, ra_min, ra_sec, dec_deg, dec_min, dec_sec
end

"""
    radec(ra::Real, dec::Real[, hours=true]) -> ra_hours, ra_minutes, ra_seconds, dec_degrees, dec_minutes, dec_seconds

### Purpose ###

Convert right ascension and declination from decimal to sexagesimal units.

### Explanation ###

The conversion is to sexagesimal hours for right ascension, and sexagesimal
degrees for declination.

### Arguments ###

* `ra`: decimal right ascension, scalar or array.  It is expressed in degrees,
  unless the optional keyword `hours` is set to `true`.
* `dec`: declination in decimal degrees, scalar or array, same number of
  elements as `ra`.
* `hours` (optional boolean keyword): if `false` (the default), `ra` is assumed
  to be given in degrees, otherwise `ra` is assumed to be expressed in hours.

### Output ###

A 6-tuple of `AbstractFloat`:

    (ra_hours, ra_minutes, ra_seconds, dec_degrees, dec_minutes, dec_seconds)

If `ra` and `dec` are arrays, also each element of the output 6-tuple are arrays
of the same dimension.

### Example ###

Position of Sirius in the sky is (ra, dec) = (6.7525, -16.7161), with right
ascension expressed in hours.  Its sexagesimal representation is given by

```jldoctest
julia> using AstroLib

julia> radec(6.7525, -16.7161, hours=true)
(6.0, 45.0, 9.0, -16.0, 42.0, 57.9600000000064)
```
"""
radec(ra::Real, dec::Real; hours::Bool=false) =
    _radec(promote(float(ra), float(dec))..., hours)

function radec(ra::AbstractArray{R}, dec::AbstractArray{D};
               hours::Bool=false) where {R<:Real, D<:Real}
    @assert length(ra) == length(dec)
    typera = float(R)
    ra_hr   = similar(ra, typera)
    ra_min  = similar(ra, typera)
    ra_sec  = similar(ra, typera)
    dec_deg = similar(ra, typera)
    dec_min = similar(ra, typera)
    dec_sec = similar(ra, typera)
    for i in eachindex(ra)
        ra_hr[i], ra_min[i], ra_sec[i], dec_deg[i], dec_min[i], dec_sec[i] =
            radec(ra[i], dec[i], hours=hours)
    end
    return ra_hr, ra_min, ra_sec, dec_deg, dec_min, dec_sec
end


text:
using NeXLParticle
using Test
using Pkg
using Pkg.Artifacts

# This is the path to the Artifacts.toml we will manipulate
artifacts_toml = joinpath(@__DIR__, "Artifacts.toml")

# Query the `Artifacts.toml` file for the hash bound to the name "iris"
# (returns `nothing` if no such binding exists)
zep_hash = artifact_hash("zeptest", artifacts_toml)

# If the name was not bound, or the hash it was bound to does not exist, create it!
if isnothing(zep_hash) || !artifact_exists(zep_hash)
    # create_artifact() returns the content-hash of the artifact directory once we're finished creating it
    zep_hash = create_artifact() do artifact_dir
        println("Downloading test data to $artifact_dir")
        tarball = joinpath(artifact_dir, "zeptest.tar.gz")
        download("https://drive.google.com/uc?export=download&id=1TZh4zbw2VY6QFlTfTpiPIG1U5uWepXRh",tarball)
        Pkg.probe_platform_engines!()
        Pkg.unpack(tarball, artifact_dir, verbose=true)
        rm(tarball)
    end
    # Now bind that hash within our `Artifacts.toml`.  `force = true` means that if it already exists,
    # just overwrite with the new content-hash.  Unless the source files change, we do not expect
    # the content hash to change, so this should not cause unnecessary version control churn.
    bind_artifact!(artifacts_toml, "zeptest", zep_hash)
end

# Get the path of the iris dataset, either newly created or previously generated.
# this should be something like `~/.julia/artifacts/dbd04e28be047a54fbe9bf67e934be5b5e0d357a`
zeptest = artifact_path(zep_hash)

@testset "ZepTest" begin
    zep = Zeppelin(joinpath(zeptest,"test.hdz"))
    @test eachparticle(zep)==1:250
    hdz = header(zep)
    @test hdz["MAG0"]=="347 78 250 98.24 5.112"
    @test hdz["INSTRUMENT"]=="NIST's TESCAN MIRA3 in 217 F101"
    es = elms(zep)
    @test length(es)==32
    @test n"Sb" in es
    @test beamenergy(zep)==20.0e3
    @test probecurrent(zep)==0.961832
    @test magdata(zep,0)["Fields"]==78.0
    @test magdata(zep,0)["Area"]==5.112
    @test length(classes(zep))==52
    @test maxparticle(zep,1:10)[123] == 42.0
    relms = ( n"Ag", n"Al", n"Ba", n"Bi", n"Br", n"C", n"Ca", n"Ce", n"Cl", n"Co", n"Cr", n"Cu", n"F", n"Fe", n"K", n"Mg", n"Mn", n"Na", n"Nd", n"Ni",
        n"O", n"P", n"Pb", n"S", n"Sb", n"Si", n"Sn", n"Sr", n"Ti", n"V", n"W", n"Zn", n"Zr" )
    refs = Dict{Element,Spectrum}( elm => loadspectrum(joinpath(zeptest,"Standards","$(elm.symbol) std.msa")) for elm in relms)
    det = matching(refs[n"Fe"],132.0)
    res = quantify(zep, det, refs, withUncertainty=true)
    @test isapprox(res[6,:FE], 94.87, atol=0.01)
    @test isapprox(res[26,:CA], 23.48, atol=0.01)
    @test isapprox(res[169, :BA], 39.37, atol=0.01)
    @test res[160,:FIRSTELM]==n"Si"
    @test res[172,:SECONDELM]==n"Ti"
    @test res[179,:THIRDELM]==n"Zn"
    @test ismissing(res[177,:FOURTHELM])
    @test res[195,:FOURTHELM]==n"K"
    NeXLParticle.writeZep(res, res.headerfile)
    res2 = Zeppelin(res.headerfile)
    @test isapprox(res2[6,:FE], 94.87, atol=0.01)
    @test isapprox(res2[26,:CA], 23.48, atol=0.01)
    @test isapprox(res2[169, :BA], 39.37, atol=0.01)
    @test res2[160,:FIRSTELM]==n"Si"
    @test res2[172,:SECONDELM]==n"Ti"
    @test res2[179,:THIRDELM]==n"Zn"
    @test ismissing(res2[177,:FOURTHELM])
    @test res2[195,:FOURTHELM]==n"K"
#@time res=NeXLParticle.quantify(zep, det, refs, withUncertainty=false);
end


text:
struct KryslCuthillMcKeePermuter
    sortbydeg::Bool
    assumesymmetry::Bool
    checksymmetry::Bool
end

KryslCuthillMcKeePermuter(; sortbydeg=false, assumesymmetry=false, checksymmetry=true) = KryslCuthillMcKeePermuter(sortbydeg, assumesymmetry, checksymmetry)

function permute_stripe(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, method::KryslCuthillMcKeePermuter; kwargs...) where {Tv, Ti}
    prm = SymRCM.symrcm(SparseMatrixCSC(A), sortbydeg=method.sortbydeg)
    return DomainPermutation(prm)
end

function permute_stripe(A::SparseMatrixCSC, method::KryslCuthillMcKeePermuter; kwargs...)
    return permute_plaid(A, method)[end]
end

function permute_plaid(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, method::KryslCuthillMcKeePermuter; kwargs...) where {Tv, Ti}
    prm = SymRCM.symrcm(SparseMatrixCSC(A), sortbydeg=method.sortbydeg)
    return (DomainPermutation(prm), DomainPermutation(prm))
end

function permute_plaid(A::SparseMatrixCSC, method::KryslCuthillMcKeePermuter; adj_A = nothing, kwargs...)
    if method.assumesymmetry || (method.checksymmetry && issymmetric(pattern(A)))
        prm = SymRCM.symrcm(A, sortbydeg=method.sortbydeg)
        return (DomainPermutation(prm), DomainPermutation(prm))
    else
        m, n = size(A)
        if adj_A === nothing
            adj_A = adjointpattern(A)
        end
        B = hvcat((2, 2), spzeros(m, m), A, adj_A, spzeros(n, n))
        prm = SymRCM.symrcm(B, sortbydeg=method.sortbydeg)
        iprm = Int[]
        jprm = Int[]
        for v in prm
            if v > m
                push!(jprm, v - m)
            else
                push!(iprm, v)
            end
        end
        return (DomainPermutation(iprm), DomainPermutation(jprm))
    end
end


text:
# cuda16 # Array,Array->Array (Multi dimensional broadcast up to some dimension with loop unrolling)
fp = open("cuda16.cu","w")
#using Knet: broadcast_ops

function cuda16src(f, j=f, ex="$f(xi,yi)")
  sprint() do s
    # it can handle arrays with 3 to 5 dimensions
    for dim_count=3:5
      for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
          print(s,"__global__ void _$(F)_16_$(dim_count)($T *x,$T *y, $T *z,")
          # place input variables
          for counter=0:dim_count-1
            print(s,"int stridex_$counter,")
          end
          for counter=0:dim_count-1
            print(s,"int stridey_$counter,")
          end
          for counter=0:dim_count-1
            print(s,"int stridez_$counter,")
          end
          # z is global index calculated from block and tread id
          print(s,

              """int N_z) {
                  int index_z = threadIdx.x + (blockIdx.x * blockDim.x);
                  int index_x,index_y;
              """)

              print(s,
              """
                  int coords[$(dim_count)];

                  while (index_z < N_z) {
                      int temp_index = index_z;
              """)
              for counter=dim_count-1:-1:0
                print(s,"\n\tcoords[$counter] = temp_index / stridez_$counter;")
                print(s,"\n\ttemp_index = temp_index % stridez_$counter;")
              end
              print(s,
              """\n
                      index_x =0;
                      index_y = 0;
              """)
              for counter=0:dim_count-1
                print(s,"\n\tindex_x+= stridex_$counter*coords[$counter];")
                print(s,"\n\tindex_y+= stridey_$counter*coords[$counter];")
              end
              print(s,
              """\n
                      $T xi = x[index_x];
                      $T yi = y[index_y];
                      z[index_z]=$ex;
                      //z[index_z]=index_z;
                      index_z+=blockDim.x * gridDim.x;

                  }
              }

              extern "C" {
                $DLLEXPORT void $(F)_16_$(dim_count)($T *x,$T *y,$T *z,""")
              for counter=0:dim_count-1
                print(s,"int stridex_$counter,")
              end
              for counter=0:dim_count-1
                print(s,"int stridey_$counter,")
              end
              for counter=0:dim_count-1
                print(s,"int stridez_$counter,")
              end

              print(s,
                """int Nz) {

                  _$(F)_16_$(dim_count)<<<256,256>>>(x,y,z,""")
              for counter=0:dim_count-1
                print(s,"stridex_$counter,")
              end
              for counter=0:dim_count-1
                print(s,"stridey_$counter,")
              end
              for counter=0:dim_count-1
                print(s,"stridez_$counter,")
              end
              print(s,
                    """Nz);
                }
              }
              """)
      end
    end
  end
end

for a in broadcast_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda16src(a...))
end

close(fp)


text:
## Banded LU decomposition

# This file is based on a part of Julia. License is MIT: https://julialang.org/license

####################
# Banded LU Factorization #
# This is like LU but factors does not actually
# permute the rows.
####################
struct BandedLU{T,S<:AbstractMatrix{T}} <: Factorization{T}
    factors::S
    ipiv::Vector{BlasInt}
    info::BlasInt

    function BandedLU{T,S}(factors, ipiv, info) where {T,S<:AbstractMatrix{T}}
        require_one_based_indexing(factors)
        new{T,S}(factors, ipiv, info)
    end
end
function BandedLU(factors::AbstractMatrix{T}, ipiv::Vector{BlasInt}, info::BlasInt) where {T}
    BandedLU{T,typeof(factors)}(factors, ipiv, info)
end
function BandedLU{T}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}, info::Integer) where {T}
    BandedLU(convert(AbstractMatrix{T}, factors),
       convert(Vector{BlasInt}, ipiv),
       BlasInt(info))
end

# iteration for destructuring into components
Base.iterate(S::BandedLU) = (S.L, Val(:U))
Base.iterate(S::BandedLU, ::Val{:U}) = (S.U, Val(:p))
Base.iterate(S::BandedLU, ::Val{:p}) = (S.p, Val(:done))
Base.iterate(S::BandedLU, ::Val{:done}) = nothing

adjoint(F::BandedLU) = Adjoint(F)
transpose(F::BandedLU) = Transpose(F)

lu(S::BandedLU) = S

function BandedLU{T}(F::BandedLU) where T
    M = convert(AbstractMatrix{T}, F.factors)
    BandedLU{T,typeof(M)}(M, F.ipiv, F.info)
end
BandedLU{T,S}(F::BandedLU) where {T,S} = BandedLU{T,S}(convert(S, F.factors), F.ipiv, F.info)
Factorization{T}(F::BandedLU{T}) where {T} = F
Factorization{T}(F::BandedLU) where {T} = BandedLU{T}(F)

copy(A::BandedLU{T,S}) where {T,S} = BandedLU{T,S}(copy(A.factors), copy(A.ipiv), A.info)

size(A::BandedLU)    = size(getfield(A, :factors))
size(A::BandedLU, i) = size(getfield(A, :factors), i)

Base.propertynames(F::BandedLU, private::Bool=false) =
    (:L, :U, :p, :P, (private ? fieldnames(typeof(F)) : ())...)

LinearAlgebra.issuccess(F::BandedLU) = F.info == 0

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::BandedLU)
    if issuccess(F)
        summary(io, F); println(io)
        println(io, "L factor:")
        show(io, mime, F.L)
        println(io, "\nU factor:")
        show(io, mime, F.U)
    else
        print(io, "Failed factorization of type $(typeof(F))")
    end
end

/(B::AbstractMatrix, A::BandedLU) = copy(transpose(transpose(A) \ transpose(B)))

# Conversions
AbstractMatrix(F::BandedLU) = (F.L * F.U)[invperm(F.p),:]
AbstractArray(F::BandedLU) = AbstractMatrix(F)
Matrix(F::BandedLU) = Array(AbstractArray(F))
Array(F::BandedLU) = Matrix(F)

##
function lu!(A::BandedMatrix{T}, pivot::Union{Val{false}, Val{true}} = Val(true);
             check::Bool = true) where T<:BlasFloat
    if pivot === Val(false)
        return banded_lufact!(A, pivot; check = check)
    end
    m= size(A,1)
    l,u = bandwidths(A) # l of the bands are ignored and overwritten
    _, ipiv = LAPACK.gbtrf!(l, u-l, m, bandeddata(A))
    return BandedLU{T,typeof(A)}(A, ipiv, zero(BlasInt))
end

_lu!(::AbstractBandedLayout, axes, A, pivot::Union{Val{false}, Val{true}} = Val(true); check::Bool = true) =
    banded_lufact!(A, pivot; check = check)

function _lu(::AbstractBandedLayout, axes, A, pivot::Union{Val{false}, Val{true}} = Val(true); check::Bool = true) where {T<:Real}
    l,u = bandwidths(A)
    lu!(BandedMatrix{float(eltype(A))}(A,(l,l+u)), pivot; check = check)
end

# Jesus christ someone loves to write "StridedMatrix" for no reason in Base!
function getproperty(F::BandedLU{T}, d::Symbol) where T
    m, n = size(F)
    if d == :L
        # not clear how to get it from F.factors so we
        # form it in the most insane way possible
        Ai = F\Matrix(I,n,m)
        tril!((F.P/Ai)/F.U)
    elseif d == :U
        return triu!(getfield(F, :factors)[1:min(m,n), 1:n])
    elseif d == :p
        return ipiv2perm(getfield(F, :ipiv), m)
    elseif d == :P
        return Matrix{T}(I, m, m)[:,invperm(F.p)]
    else
        getfield(F, d)
    end
end

# function banded_lufact!(A::AbstractMatrix{T}, ::Val{Pivot} = Val(true);
#                          check::Bool = true) where {T,Pivot}
#     m, n = size(A)
#     l, u = bandwidths(A)
#     minmn = min(m,n)
#     info = 0
#     ipiv = Vector{BlasInt}(undef, minmn)
#     @inbounds begin
#         for k = 1:minmn
#             # find index max
#             kp = k
#             if Pivot
#                 amax = abs(zero(T))
#                 for i = k:min(k+l,m)
#                     absi = abs(A[i,k])
#                     if absi > amax
#                         kp = i
#                         amax = absi
#                     end
#                 end
#             end
#             ipiv[k] = kp
#             if !iszero(A[kp,k])
#                 if k != kp
#                     # Interchange
#                     for i = rowrange(A,k)
#                         tmp = A[k,i]
#                         A[k,i] = A[kp,i]
#                         A[kp,i] = tmp
#                     end
#                 end
#                 # Scale first column
#                 Akkinv = inv(A[k,k])
#                 for i = k+1:min(k+l,m)
#                     A[i,k] *= Akkinv
#                 end
#             elseif info == 0
#                 info = k
#             end
#             # Update the rest
#             for j = k+1:min(k+u,n)
#                 for i = k+1:min(k+l,m)
#                     A[i,j] -= A[i,k]*A[k,j]
#                 end
#             end
#         end
#     end
#     check && checknonsingular(info)
#     return BandedLU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
# end

_apply_ipiv_rows!(A::BandedLU, B::AbstractVecOrMat) = _ipiv_rows!(A, 1 : length(A.ipiv), B)
_apply_inverse_ipiv_rows!(A::BandedLU, B::AbstractVecOrMat) = _ipiv_rows!(A, length(A.ipiv) : -1 : 1, B)

function _ipiv_rows!(A::BandedLU, order::OrdinalRange, B::AbstractVecOrMat)
    for i = order
        if i != A.ipiv[i]
            _swap_rows!(B, i, A.ipiv[i])
        end
    end
    B
end

function _swap_rows!(B::AbstractVector, i::Integer, j::Integer)
    B[i], B[j] = B[j], B[i]
    B
end

function _swap_rows!(B::AbstractMatrix, i::Integer, j::Integer)
    for col = 1 : size(B, 2)
        B[i,col], B[j,col] = B[j,col], B[i,col]
    end
    B
end


text:
struct Production{S<:System}
    type::Type{S}
    args
end
Base.iterate(p::Production) = (p, nothing)
Base.iterate(p::Production, ::Nothing) = nothing
Base.eltype(::Type{Production}) = Production

mutable struct Produce{P,V,S} <: State{P}
    name::Symbol # used in recurisve collecting in collect()
    value::V
    productions::Vector{Production{<:S}}
end

Produce(; _name, _type, _...) = begin
    T = _type
    #TODO: maybe now we can get rid of V
    if T <: System
        P = T
        V = Union{T,Nothing}
        v = nothing
        S = typefor(P)
    elseif T <: Vector{<:System}
        P = T
        V = T
        v = V[]
        S = typefor(eltype(P))
    end
    Produce{P,V,S}(_name, v, Production{S}[])
end

constructortags(::Val{:Produce}) = ()

produce(s::Type{<:System}; args...) = Production(typefor(s), args)
produce(::Nothing; args...) = nothing
unittype(s::Produce) = nothing

Base.getindex(s::Produce{<:System}, i::Int) = i == 1 ? s.value : throw(BoundsError(s, i))
Base.length(s::Produce{<:System}) = isnothing(s.value) ? 0 : 1
Base.iterate(s::Produce{<:System}) = isempty(s) ? nothing : (s.value, nothing)
Base.iterate(s::Produce{<:System}, ::Nothing) = nothing
Base.eltype(::Type{Produce{S}}) where {S<:System} = S

Base.getindex(s::Produce{<:Vector}, i::Int) = getindex(s.value, i)
Base.getindex(s::Produce{<:Vector}, ::Nothing) = s
Base.length(s::Produce{<:Vector}) = length(s.value)
Base.iterate(s::Produce{<:Vector}, i=1) = i > length(s) ? nothing : (s[i], i+1)
Base.eltype(::Type{<:Produce{Vector{S}}}) where {S<:System} = S

Base.isempty(s::Produce) = length(s) == 0

priority(::Type{<:Produce}) = PrePriority()

produce!(s::Produce, ::Nothing) = nothing
produce!(s::Produce, p::Production) = push!(s.productions, p)
produce!(s::Produce, P::Vector) = produce!.(Ref(s), P)

produce!(s::Produce, p::System) = (s.value = p)
produce!(s::Produce{V}, p::System) where {V<:Vector} = push!(s.value, p)

export produce

updatetags!(d, ::Val{:Produce}; _...) = begin
    #HACK: fragile type string check
    d[:single] = !occursin("Vector{", string(d[:_type]))
end

genvartype(v::VarInfo, ::Val{:Produce}; V, _...) = begin
    if istag(v, :single)
        @q Produce{$V,Union{$V,Nothing},typefor($V)}
    else
        @q Produce{$V,$V,typefor(eltype($V))}
    end
end

gendefault(v::VarInfo, ::Val{:Produce}) = nothing

genupdate(v::VarInfo, ::Val{:Produce}, ::PreStage; kw...) = begin
    @gensym s a P c p b
    @q let $s = $(symstate(v)),
           $a = $C.value($s),
           $P = $s.productions,
           $c = context
        for $p in $P
            $b = $p.type(; context=$c, $p.args...)
            $C.produce!($s, $b)
        end
        empty!($P)
        $C.update!($a, $C.PreStage())
    end
end

# Produce referenced in args expected to be raw state, not extracted by value(), for querying
genupdate(v::VarInfo, ::Val{:Produce}, ::PreStep; kw...) = symstate(v)

genupdate(v::VarInfo, ::Val{:Produce}, ::MainStep; kw...) = begin
    @gensym s a
    @q let $s = $(symstate(v)),
           $a = $C.value($s)
        $C.update!($a)
        $s
    end
end

genupdate(v::VarInfo, ::Val{:Produce}, ::PostStep; kw...) = begin
    @gensym s P
    q = if istag(v, :single)
        @q let $s = $(symstate(v))
            if isempty($s)
                let $P = $(genbody(v))
                    $C.produce!($s, $P)
                end
            end
        end
    else
        @q let $s = $(symstate(v)),
               $P = $(genbody(v))
            $C.produce!($s, $P)
        end
    end
    w = gettag(v, :when)
    isnothing(w) ? q : @q if $C.value($w); $q end
end

genupdate(v::VarInfo, ::Val{:Produce}, ::PostStage; kw...) = begin
    @gensym s a
    @q let $s = $(symstate(v)),
           $a = $C.value($s)
        $C.update!($a, $C.PostStage())
    end
end


text:
abstract type AbstractRating{T <: Number}
end

struct PredictRating{T <: Number} <: AbstractRating{T}
    value::Real
    target::T
    PredictRating(x::Real, preference::Preference{T}) where T <: Number = new{T}(correct(x, preference), round(x, preference))
end

struct Rating{T <: Number} <: AbstractRating{T}
    value::T
    history::Vector{Persa.Rating{T}}
    Rating(x::T, preference::Preference{T}) where T <: Number = new{T}(correct(x, preference), Vector{Persa.Rating{T}}())
end

struct MissingRating{T <: Number} <: AbstractRating{T}
end

function Base.getindex(rating::Persa.Rating, i::Int)
    if i < 1 || (i - 1) > length(rating.history)
        throw(BoundsError(rating, i))
    elseif i == 1
        return rating
    end

    return rating.history[i-1]
end

MissingRating() = MissingRating{Number}()

function rerating(r1::Persa.Rating{T}, r2::Persa.Rating{T}) where T
    push!(r1.history, r2)
    return r1
end

value(rating::Rating) = rating.value
value(rating::PredictRating) = rating.value
value(rating::MissingRating) = missing

Base.string(x::Rating) = string(value(x))
Base.print(io::IO, x::Rating) = print(io, string(x))
Base.show(io::IO, x::Rating) = print(io, "Rating: ", x)

Base.string(x::PredictRating) = string(value(x), " (", x.target ,")")
Base.print(io::IO, x::PredictRating) = print(io, string(x))
Base.show(io::IO, x::PredictRating) = print(io, "Rating: ", x)

Base.string(x::MissingRating) = string(value(x))
Base.print(io::IO, x::MissingRating) = print(io, string(x))
Base.show(io::IO, x::MissingRating) = print(io, "Rating: ", x)

Base.isnan(rating::AbstractRating{T}) where T <: Number = false
Base.isnan(rating::MissingRating{T}) where T <: Number = true

Base.ismissing(rating::MissingRating{T}) where T <: Number = true
Base.ismissing(rating::AbstractRating{T}) where T <: Number = false

Base.zero(::Type{AbstractRating{T}}) where T <: Number = MissingRating{T}()

function convert(values::Array{T}, preference::Preference{T})::Array{AbstractRating{T}} where T <: Number
    [Rating(value, preference) for value in values]
end

function convert(values::Array{T})::Array{AbstractRating{T}} where T <: Number
    convert(values, Preference(values))
end

#Missing Operations
missingOperationsError(operation) = "Operation \'$operation\' is invalid with missing values"
Base.:(+)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("-"))
Base.:(-)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("-"))
Base.:(*)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("*"))
Base.:(/)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("/"))

Base.:(+)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("-"))
Base.:(-)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("-"))
Base.:(*)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("*"))
Base.:(/)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("/"))

Base.:(+)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("+"))
Base.:(-)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("-"))
Base.:(*)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("*"))
Base.:(/)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("/"))

Base.:(+)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("+"))
Base.:(-)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("-"))
Base.:(*)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("*"))
Base.:(/)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("/"))

Base.:(==)(r1::MissingRating, r2::MissingRating) where {T <: Number} = error(missingOperationsError("=="))

Base.:(==)(r1::MissingRating, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("=="))
Base.:(!=)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("!="))
Base.:(>)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError(">"))
Base.:(<)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("<"))
Base.:(>=)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError(">="))
Base.:(<=)(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = error(missingOperationsError("<="))

Base.:(==)(r1::AbstractRating{T}, r2::MissingRating) where {T <: Number} = error(missingOperationsError("=="))
Base.:(!=)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("!="))
Base.:(>)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError(">"))
Base.:(<)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("<"))
Base.:(>=)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError(">="))
Base.:(<=)(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = error(missingOperationsError("<="))

Base.:(==)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("=="))
Base.:(!=)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("!="))
Base.:(>)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError(">"))
Base.:(<)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("<"))
Base.:(>=)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError(">="))
Base.:(<=)(r1::MissingRating{T}, x::Number) where {T <: Number} = error(missingOperationsError("<="))

Base.:(==)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("=="))
Base.:(!=)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("!="))
Base.:(>)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError(">"))
Base.:(<)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("<"))
Base.:(>=)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError(">="))
Base.:(<=)(x::Number, r1::MissingRating{T}) where {T <: Number} = error(missingOperationsError("<="))

Base.:(+)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) + value(r2)
Base.:(-)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) - value(r2)
Base.:(*)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) * value(r2)
Base.:(/)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) / value(r2)

Base.:(+)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) + x
Base.:(-)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) - x
Base.:(*)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) * x
Base.:(/)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) / x

Base.:(+)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x + value(r1)
Base.:(-)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x - value(r1)
Base.:(*)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x * value(r1)
Base.:(/)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x / value(r1)

Base.:(==)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) == value(r2)
Base.:(!=)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) != value(r2)
Base.:(>=)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) >= value(r2)
Base.:(<=)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) <= value(r2)
Base.:(>)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) > value(r2)
Base.:(<)(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) < value(r2)

Base.:(==)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) == x
Base.:(!=)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) != x
Base.:(>=)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) >= x
Base.:(<=)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) <= x
Base.:(>)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) > x
Base.:(<)(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) < x

Base.:(==)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x == value(r1)
Base.:(!=)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x != value(r1)
Base.:(>=)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x >= value(r1)
Base.:(<=)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x <= value(r1)
Base.:(>)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x > value(r1)
Base.:(<)(x::Number, r1::AbstractRating{T}) where {T <: Number} = x < value(r1)