text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AzureReachabilityReportLocation <: SwaggerModel
    country::Any # spec type: Union{ Nothing, String } # spec name: country
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    city::Any # spec type: Union{ Nothing, String } # spec name: city

    function AzureReachabilityReportLocation(;country=nothing, state=nothing, city=nothing)
        o = new()
        validate_property(AzureReachabilityReportLocation, Symbol("country"), country)
        setfield!(o, Symbol("country"), country)
        validate_property(AzureReachabilityReportLocation, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(AzureReachabilityReportLocation, Symbol("city"), city)
        setfield!(o, Symbol("city"), city)
        o
    end
end # type AzureReachabilityReportLocation

const _property_map_AzureReachabilityReportLocation = Dict{Symbol,Symbol}(Symbol("country")=>Symbol("country"), Symbol("state")=>Symbol("state"), Symbol("city")=>Symbol("city"))
const _property_types_AzureReachabilityReportLocation = Dict{Symbol,String}(Symbol("country")=>"String", Symbol("state")=>"String", Symbol("city")=>"String")
Base.propertynames(::Type{ AzureReachabilityReportLocation }) = collect(keys(_property_map_AzureReachabilityReportLocation))
Swagger.property_type(::Type{ AzureReachabilityReportLocation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AzureReachabilityReportLocation[name]))}
Swagger.field_name(::Type{ AzureReachabilityReportLocation }, property_name::Symbol) =  _property_map_AzureReachabilityReportLocation[property_name]

function check_required(o::AzureReachabilityReportLocation)
    (getproperty(o, Symbol("country")) === nothing) && (return false)
    true
end

function validate_property(::Type{ AzureReachabilityReportLocation }, name::Symbol, val)
end


text:
# windows users are required to specify full path to the exe file unless they configured PATH variable
TESSERACT_PATH = "tesseract"
image_path = joinpath("data", "book-covers", "test", "007184418X.jpg")

text_str = readstring(`$TESSERACT_PATH $image_path stdout`)
text_array = readlines(`$TESSERACT_PATH $image_path stdout`)

println(text_str)
println(text_array)

ext: jl


text:
module PlotlyJSTest

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using PlotlyJS
typealias M PlotlyJS

tests = length(ARGS) > 0 ? ARGS : ["traces", "api", "styles"]

for fn in tests
    include(string(fn, ".jl"))
end

end


text:
module ContStressDemo

using TopOpt, LinearAlgebra, StatsFuns
# using Makie
# using TopOpt.TopOptProblems.Visualization: visualize

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 4.0 # filter radius
xmin = 0.0001 # minimum density
problem_size = (160, 40)
x0 = fill(1.0, prod(problem_size)) # initial design
p = 4.0 # penalty
stress_threshold = 1.0 # maximum von Mises stress
compliance_threshold = 800 # maximum compliance

problem = PointLoadCantilever(Val{:Linear}, problem_size, (1.0, 1.0), E, v, f)
#problem = HalfMBB(Val{:Linear}, problem_size, (1.0, 1.0), E, v, f)

solver = FEASolver(Direct, problem; xmin=xmin)

cheqfilter = DensityFilter(solver; rmin=rmin)
stress = TopOpt.MicroVonMisesStress(solver)
comp = TopOpt.Compliance(problem, solver)

function obj(x)
    # minimize volume
    return sum(cheqfilter(x)) / length(x)
end
function constr1(x)
    # L-20 norm aggregation to approximate maximum function
    return norm(stress(cheqfilter(x)), 20) - stress_threshold
end
function constr2(x)
    # compliance upper-bound
    return comp(cheqfilter(x)) - compliance_threshold
end

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr1)
Nonconvex.add_ineq_constraint!(m, constr2)

options = MMAOptions(; maxiter=1000, tol=Tolerance(; kkt=1e-4, f=1e-4))
TopOpt.setpenalty!(solver, p)
@time r = Nonconvex.optimize(
    m, MMA87(; dualoptimizer=ConjugateGradient()), x0; options=options
)

@show obj(r.minimizer)
@show constr1(r.minimizer)
@show constr2(r.minimizer)
@show maximum(stress(cheqfilter(r.minimizer)))
topology = cheqfilter(r.minimizer)
# fig = visualize(problem, solver.u;
#     topology = topology, default_exagg_scale=0.0, scale_range=10.0)
# Makie.display(fig)

end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# all other glasses should follow the format below, new glasses must be added to OTHER_GLASSES and OTHER_GLASS_NAMES where the index in the array matches the numeric part of the GlassID

module CARGILLE
using ..GlassCat: Glass, GlassID, OTHER

const OG0608 = Glass(GlassID(OTHER, 1), -2, 1.4451400, 0.0043176, -1.80659e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 800.0, -1.0, 0, -1.0, [(0.32, 0.03, 10.0), (0.365, 0.16, 100.0), (0.4047, 0.40, 100.0), (0.480, 0.71, 100.0), (0.4861, 0.72, 100.0), (0.5461, 0.80, 100.0), (0.5893, 0.90, 100.0), (0.6328, 0.92, 100.0), (0.6439, 0.95, 100.0), (0.6563, 0.96, 100.0), (0.6943, 0.99, 100.0), (0.840, 0.99, 100.0), (0.10648, 0.74, 100.0), (0.1300, 0.39, 100.0), (0.1550, 0.16, 100.0)], 1.457518, -1.0, -1.0, 0, 57.18978, 0, 0.878, -1)
const OG0607 = Glass(GlassID(OTHER, 2), -2, 1.44503, 0.0044096, -2.85878e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [(0.32, 0.15, 10.0), (0.365, 0.12, 100.0), (0.4047, 0.42, 100.0), (0.480, 0.78, 100.0), (0.4861, 0.79, 100.0), (0.5461, 0.86, 100.0), (0.5893, 0.90, 100.0), (0.6328, 0.92, 100.0), (0.6439, 0.90, 100.0), (0.6563, 0.92, 100.0), (0.6943, 0.98, 100.0), (0.840, 0.99, 100.0), (0.10648, 0.61, 100.0), (0.1300, 0.39, 100.0), (0.1550, 0.11, 100.0)], 1.457587, -1.0, -1.0, 0, 57.19833, 0, 0.878, -1)
const OG081160 = Glass(GlassID(OTHER, 3), -2, 1.49614, 0.00692199, -8.07052e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.000885983052189022, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.014, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [(0.32, 0.04, 100.0), (0.365, 0.13, 100.0), (0.4047, 0.26, 100.0), (0.480, 0.48, 100.0), (0.4861, 0.49, 100.0), (0.5461, 0.60, 100.0), (0.5893, 0.68, 100.0), (0.6328, 0.71, 100.0), (0.6439, 0.73, 100.0), (0.6563, 0.74, 100.0), (0.6943, 0.76, 100.0), (0.840, 0.83, 100.0), (0.10648, 0.86, 100.0), (0.1300, 0.89, 100.0), (0.1550, 0.90, 100.0)], 1.515549, -1.0, -1.0, 0, 36.82493, 0, 1.11, -1)

end

const OTHER_GLASSES = [CARGILLE.OG0607, CARGILLE.OG0608, CARGILLE.OG081160]
const OTHER_GLASS_NAMES = ["CARGILLE.OG0607", "CARGILLE.OG0608", "CARGILLE.OG081160"]


text:
module Observations

using ..OceanTurbulenceParameterEstimation
using ..OceanTurbulenceParameterEstimation.Grids: OneDimensionalEnsembleGrid

using Oceananigans
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Fields
using JLD2

abstract type AbstractTimeSeries end

const AbstractTimeSeriesBatch = Vector{<:AbstractTimeSeries}

export # truth_data.jl
       OneDimensionalTimeSeries,
       OneDimensionalTimeSeriesBatch,
       column_ensemble_interior,

       # set_field.jl
       set!

include("file_wrangling.jl")
include("one_dimensional_time_series.jl")
include("one_dimensional_time_series_batch.jl")
include("set_field.jl")

end #module


text:
push!(LOAD_PATH,"../src/") # load module from local directory

using TheDancer

N, Edisorder, Jdisorder, modelJ, B, dipolestrength, r = init!(50, 0.0, 0.0, 0.001, 298, 7.5)

function main()

    SCFcycles=200
    Unitarycycles=500

    for dampening in [0.07,0.025,0.05]
        SCFthenUnitary(dampening, SCFcycles, Unitarycycles, PNG=false) # PNG=true for .pngs for movie making
    end
end

main() # Party like it's C99!


text:
@blueprint struct Documenter
    assets::Vector{String}
    # logo should be user input from interface
    make_jl::String # make.jl template
    index_md::String # index.md template
    devbranch::Maybe{String} # if nothing use default branch
end

function compile(bp::Documenter, ctx::Context)
end


text:
export turbconv_sources

turbconv_sources(m::TurbulenceConvectionModel) = ()


text:
"""
Description

      The colors() method either lists available colors or plots a specified color(s).

Usage

      colors() # shows color definitions that are available
      colors("goldenrod") # plots the color
      colors("goldenrod", "khaki") # plots both colors

Method Signature(s)

      colors()
      colors(c::ASCIIString)
      colors(c1::ASCIIString, c2::ASCIIString)

Details

      The colors function has two multiple-dispatch definitions. Without an arugment, the method returns
      a list of available colors. When a string is passed, that defined color is plotted. When two strings
      are passed, both colors are platted adjacent to one another. The definitions are derived from the
      Colors.jl package's src/names_data.jl file. This list includes both X11 and SVG definitions
      and when a collision occurs, the SVG definition takes precedence.

References

      http://www.html.am/html-codes/color/color-code-chart.cfm
      https://github.com/JuliaGraphics/Colors.jl
"""
list    = include(Pkg.dir("Colors/src/names_data.jl"))

palette = ASCIIString[]

for (k,v) in list
    push!(palette, k)
end

function colors()
    println()
    print_with_color(:yellow, "The following colors are available")
    println()
    for p in palette
        print_with_color(:blue, "$p  ")
    end
    println()
end

function colors(c::ASCIIString)
    x   = linspace(0,1)
    y   = zeros(x)
    abv = FillAbove(x, y, color=c)
    blw = FillBelow(x, y, color=c)
    p   = Winston.plot()
    add(p, abv, blw)
end

function colors(c1::ASCIIString, c2::ASCIIString)
    x   = linspace(0,1)
    y   = zeros(x)
    abv = FillAbove(x, y, color=c1)
    blw = FillBelow(x, y, color=c2)
    p   = Winston.plot()
    add(p, abv, blw)
end


text:
const days_per_year = 365.2422
const global_base_date = Date(2000,1,1)
const global_base_date_as_day = convert(Dates.Day, global_base_date)

"""
    years_between(until::Date, from::Date)
    years_between(until::Dates.Day, from::Dates.Day)

Returns the number of years between two dates. For the purposes of this calculation
there are 365.2422 days in a year.
"""
function years_between(until::Date, from::Date)
    return (Dates.days(until) -Dates.days(from))/ days_per_year
end
function years_between(until::Dates.Day, from::Dates.Day)
    return (convert(Int, until)-convert(Int, from))/ days_per_year
end

"""
    years_between(a::Date, b::Date)
    years_between(a::Dates.Day, b::Dates.Day)

Returns the number of years that have elapsed since 1-Jan-2000. For the purposes of this calculation
there are 365.2422 days in a year.
"""
function years_from_global_base(a::Date)
    return years_between(a, global_base_date)
end

function years_from_global_base(a::Dates.Day)
    return years_between(a, global_base_date_as_day)
end

"""
    Period length is designed to convert TimePeriod objects to a float in a consistent way to years_from_global_base
"""
function period_length(a::Dates.DatePeriod, base::Date = global_base_date)
    return years_between(base+a, base)
end


text:
using ParticleDA, MPI

include(joinpath(@__DIR__, "model", "model.jl"))
using .Model

# Initialise MPI
MPI.Init()

# Get the number or ranks, so that we can set a number of particle as an integer
# multiple of them.
my_size = MPI.Comm_size(MPI.COMM_WORLD)

params = Dict(
    "filter" => Dict(
        "nprt" => my_size,
        "enable_timers" => true,
        "verbose"=> false,
        "n_time_step" => 5,
        "output_filename" => "warmup.h5",
    ),
    "model" => Dict(
        "llw2d" => Dict(
            "nx" => 20,
            "ny" => 20,
            "nobs" => 4,
            "padding" => 0,
        ),
    ),
)

# Warmup
run_particle_filter(Model.init, params, BootstrapFilter())
# Flush a newline
println()

# Run the command
params["filter"]["output_filename"] = "bootstrap_filter.h5"
params["filter"]["verbose"] = true
params["filter"]["nprt"] = 2 * my_size
params["model"]["llw2d"]["nobs"] = 36
run_particle_filter(Model.init, params, BootstrapFilter())
# Flush a newline
println()

# Run the command
params["filter"]["output_filename"] = "optimal_filter.h5"
run_particle_filter(Model.init, params, OptimalFilter())
# Flush a newline
println()


text:
#
# JSON client with Julia
#
# author: Atsushi Sakai
#

using HTTP
using JSON

function get_config()
    ip = "localhost" #default
    if haskey(ENV,"SERVER_IP")
        ip = ENV["SERVER_IP"]
    end

	port = "8000" #default
    if haskey(ENV,"SERVER_PORT")
        port = ENV["SERVER_PORT"]
    end

    println("request to ip:", ip)
    println("request to port:", port)

    return ip, port
end

function main()
	req = Dict("Name" => "julia json request")
	println("request:")
	println(req)

    ip, port = get_config()

    url ="http://"*ip*":"*port
  	r = HTTP.request("POST", url,
               ["Content-Type" => "application/json"],
               JSON.json(req))

    println("response:")
  	println(JSON.parse(String(r.body)))

    println(PROGRAM_FILE," Done!!")
end

main()

max_stars_repo_name: AtsushiSakai/JSONServerAndClient

max_issues_repo_name: AtsushiSakai/JSONServerAndClient

max_forks_repo_name: AtsushiSakai/JSONServerAndClient


text:
##########
##########
# Functions to perform perturbations on the food webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

# set up =======================================================================
using DrWatson
@quickactivate "trait-based-rewiring"

# logger set up
using Logging: global_logger
using TerminalLoggers: TerminalLogger
using DifferentialEquations
global_logger(TerminalLogger())

include(srcdir("00_species_list_creation.jl"))
include(srcdir("01_eqn_component_functions.jl"))
include(srcdir("04_perturbation_sims_functions.jl"))

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# pulse perturbations ==========================================================

# little line to set the network id (this happens from the shell line)
network_id = 1

species_list = readdlm(datadir("./intermediate-files/species-lists-pre-perturbation/species_list_0001.csv"), ',', Any)
network = create_interaction_matrix(species_list)

# set vector callback
extinction_callback =
    VectorContinuousCallback(extinction_threshold,
                            extinction_affect!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))
steady_callback = TerminateSteadyState(1e-8,
                                        1e-6,
                                        DiffEqCallbacks.allDerivPass,
                                        min_t = 1205)

# set callback set
save_positions = (true, true)
cb_inc = DiscreteCallback(inc_mort_cond, affect_inc!, save_positions=save_positions)
save_positions = (false, true)
cb_dec = DiscreteCallback(dec_mort_cond, affect_dec!, save_positions=save_positions)
callback_set = CallbackSet(extinction_callback,
                            steady_callback, #keep this out until fixed
                            cb_inc,
                            cb_dec)

#  allocate x0
u0 = Float64[]
for i in 1:size(species_list, 1)
    append!(u0, species_list[i, biomass_pos])
end

# set perturbation tspan
tspan = (1000.0, 2000.0)

# set constant stop values
const stop1 = [1100.0]
const stop2 = [1200.0]
const tstop = [1100.;1200.]

# loop through options of the c value to see what works
# c_list = [80 100 110 120 130 140 150 160 170 180 190 200 210 220]

for i in 1:length(c_list)

    c = c_list[i]
    x = 100
    f = 1
    z = c/(x*(1/f))

    # pick set of random species to perturb
    all = Array(1:size(species_list, 1))
    perturb_sp = rand(Array(1:size(species_list, 1)),
                        Int(round(size(species_list, 1)*(1/2)))) #changing to perturb half of them
    non_perturb = setdiff(Set(all), Set(perturb_sp))

    # set p
    # the amount (fourth element) has to be initially set to one to keep the
    # "amount" equal to one until the affect_inc! function comes in
    p = [species_list, network, 0.3, 1.0, perturb_sp, non_perturb]

    prob_pulse = ODEProblem(biomass_perturb_random,
                                        u0,
                                        tspan,
                                        p,
                                        tstops=tstop)
    @time sol_prob_pulse = solve(prob_pulse,
                                        Tsit5(),
                                        callback = callback_set,
                                        reltol = 1e-8,
                                        #saveat = save_times,
                                        progress = true,
                                        progress_steps = 1)
    writedlm(datadir("testing-files/c-values/solution_ext-1e04$c.csv"), sol_prob_pulse, ",")
    plot(sol_prob_pulse)
    # write out complete version of the final c states
    num_ext = sol_prob_pulse.u[end]
    num_ext = string("number of extinct species for ", c, " is ", length(num_ext[num_ext .== 0]))
    open(datadir("testing-files/c-values/zeros_ext-1e04$c.txt"), "w") do file
        write(file, num_ext)
    end
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
function make_all_charts(parent_directory::String, corpus::CorpusUtils.Corpus, save_images::Bool; show_images::Bool=false)
	if save_images Plots.plotlyjs() elseif show_images Plots.plotly() end
	all_tokens = sum(map(d -> sum(values(d.value)), corpus.documents))
	info("all tokens: $all_tokens")
	all_likelihoods = Plots.plot(title="ELBOs", size=(1024, 768), ylabel="lower bound on marginal likelihood", xlabel="# documents seen", guidefont=font(16, "sans-serif"), titlefont=font(18, "sans-serif"), tickfont=font(14, "sans-serif"))
	all_test = Plots.plot(title="Predictive likelihood", ylabel="predictive likelihood", xlabel="# documents seen", size=(1024, 768),  guidefont=font(16, "sans-serif"), titlefont=font(18, "sans-serif"), tickfont=font(14, "sans-serif"))
	efn = ""
	ppx_values = Dict{String, Float64}()
	for dir in readdir(parent_directory)
		if isdir(joinpath(parent_directory, dir))
			kernel_params = match(r"\[(.*)\]", dir)[1]
			kernel_name = replace(dir, r"Kernel(Real)?\[.*\]", "")
			for alpha_dir in readdir(joinpath(parent_directory, dir))
				cur_alpha_dir = joinpath(parent_directory, dir, alpha_dir)
				alpha = replace(alpha_dir,"alpha", "")
				if isdir(cur_alpha_dir)
					for model_param_var in readdir(cur_alpha_dir)
						cur_model_param_path = joinpath(cur_alpha_dir, model_param_var)
						if isdir(cur_model_param_path)
							for num_topics in readdir(cur_model_param_path)
								num_topics_path = joinpath(cur_model_param_path, num_topics)
								if isdir(num_topics_path)
									K = parse(Int64, replace(num_topics, "_topics", ""))
									inducing_point_model = false
									if isfile(joinpath(num_topics_path, "inducing_times.dat")) inducing_point_model = true end
							    	println("processing $num_topics_path (is inducing point model: $inducing_point_model)")
							    	# plot_label = "$kernel_name($kernel_params), ⍺=$alpha, K=$K"
							    	my_kernel_name = replace(kernel_name, r"Brownian", "Wiener")
							    	my_kernel_name = replace(my_kernel_name, r"OrnsteinUhlenbeck", "Ornstein Uhlenbeck")
							    	plot_label = "$my_kernel_name kernel"
							    	(c, e) = tab_split(joinpath(num_topics_path, "elbo"))
							    	Plots.plot!(all_likelihoods, c, e, label=plot_label, linewidth=2, legendfont=font(16, "sans-serif"))
							    	(tc, tppx) = tab_split(joinpath(num_topics_path, "test"))
							    	ppx_values[plot_label] = tppx[end]
							    	Plots.plot!(all_test, tc, tppx, label=plot_label, linewidth=2, legendfont=font(16, "sans-serif"))
							    	make_topic_charts(num_topics_path, corpus, save_images, show_images)
							    end
							end
						end
					end
				end
			end
	    end
	end
	display(ppx_values)
	ppx_file = open(joinpath(parent_directory, "perplexities"), "w")
	for (key, val) in ppx_values
		println("$key: $val")
		@printf(ppx_file, "%s\t%f\n", key, val)
	end
	close(ppx_file)
	if save_images
		Plots.pdf(all_likelihoods, joinpath(parent_directory, "elbos"))
		Plots.pdf(all_test, joinpath(parent_directory, "tests"))
		Plots.plotly()
	else
		display(all_likelihoods)
		display(all_test)
	end
end

function make_topic_charts(path::String, corpus::CorpusUtils.Corpus, save_images::Bool, show_images::Bool)
	if save_images
		if isdir(joinpath(path, "word_trajectories"))
			return
		end
		Plots.plotlyjs()
	elseif show_images
		Plots.plotly()
	end
	all_colors = Colors.distinguishable_colors(length(corpus.lexicon))
	all_times = readlines(joinpath(path, "times.dat"))

	ticks = DynamicData.get_unique_timestamps(corpus.documents)

	if length(ticks) < 20
		tick_pos = collect(1:length(ticks))
	else
		tick_pos = collect(1:10:length(ticks))
	end
	tick_labels = String[]

	juliandays = map(Dates.datetime2julian, ticks)
	time_range = maximum(juliandays) - minimum(juliandays)
	println("time range: $time_range")
	if time_range > 21*365
		tick_labels = map(string, map(Dates.year, ticks))
	elseif time_range > 2*365
		tick_labels = map(ts -> Dates.format(ts, "yyyy-mm"), ticks)
	else
		tick_labels = map(ts -> Dates.format(ts, "yyyy-mm-dd"), ticks)
	end
	ticks = tuple(tick_pos, tick_labels[tick_pos])

	topic_probabilities = get_topic_distributions(path, save_images)
	K = length(topic_probabilities)
	for k in 1:K
		(T,V) = size(topic_probabilities[k])
		word_means = vec(mean(topic_probabilities[k],1))
		words_to_plot = sortperm(word_means, rev=true)

		p = Plots.plot(size=(1024, 768), xlabel="Time", ylabel="Probability in topic", xticks=ticks, grid=false, tickfont=font(14, "sans-serif"), guidefont=font(16, "sans-serif"), rotation=90, legend=:legend, color_palette=all_colors[words_to_plot])
		for w in words_to_plot
			for t in 1:T
				if isnan(topic_probabilities[k][t,w]) || isinf(topic_probabilities[k][t,w])
					println("warning, word prob is nan or inf")
					topic_probabilities[k][t,w] = 0.
				end
			end
			f_w = topic_probabilities[k][:,w]
			Plots.plot!(p, collect(1:T), f_w, label=corpus.lexicon.words[w].value, linewidth=4, legendfont=font(20, "sans-serif"))
		end

		if save_images
			img_path = joinpath(path, "word_trajectories")
			mkpath(img_path)
			Plots.pdf(p, joinpath(img_path, "topic$k"))
		elseif show_images
			display(p)
		end
	end
	if save_images Plots.plotly() end
end


text:
include("../src/coach.jl")

using Test

# @test length(read_noise()) > 0
# @test length(read_noise()) == 84671552

# @test length(alter_noise(read_noise())) > 0
# @test length(alter_noise(read_noise())) < 84671552

# @test create_multiplier(Int8(0)) == 1.0
# @test create_multiplier(Int8(127)) > 1.12
# @test create_multiplier(Int8(-127)) < 0.88

noise_array = alter_noise(read_noise())

workouts = generate_workout_sequence(8.0,
                          42.2,
                          noise_array,
                          4,
                          20,
                          200,
                          10,
                          200,
                          20,
                          600,
                          Dates.Date(2018,9,4),
                          true,
                          0.1,
                          Dates.Sunday)

for i in workouts
    println(i)
end


text:
@deprecate numericallyapunov(args...; kwargs...) lyapunov_from_data(args...; kwargs...)
@deprecate grassberger(args...; kwargs...) grassberger_dim(args...; kwargs...)


text:
module Assignment04

export normalizeDNA,
       basecomposition,
       gc_content,
       kmercount

include("assignment.jl")

end


text:
@testset "tests" begin

@test parse_one("", Epsilon()) == []
@test parse_one("", Insert("foo")) == ["foo"]
@test parse_one("", Drop(Insert("foo"))) == []
@test_throws ParserException parse_one("x", Equal("a"))
@test parse_one("a", Equal("a")) == ["a"]
@test parse_one("aa", Equal("a")) == ["a"]
@test_throws ParserException parse_one("a", Repeat(Equal("a"), 2, 2))
@test parse_one("aa", Repeat(Equal("a"), 2, 2)) == ["a", "a"]
@test parse_one("aa", Repeat(Equal("a"), 1, 2)) == ["a", "a"]
@test parse_one("", Repeat(Equal("a"), 0, 0)) == []
@test_throws ParserException parse_one("a", Repeat(Equal("a"), 2, 2; greedy=false); debug=true)
@test parse_one("aa", Repeat(Equal("a"), 2, 2; greedy=false)) == ["a", "a"]
@test parse_one("aa", Repeat(Equal("a"), 1, 2; greedy=false)) == ["a"]
@test parse_one("", Repeat(Equal("a"), 0, 0; greedy=false)) == []

@test parse_one("ab", Series(Pattern(r"a"), Dot(); flatten=false)) == Any[["a"], ['b']]
@test parse_one("ab", Series(Pattern(r"a"), Dot())) == Any["a", 'b']
@test parse_one("ab", Seq(Pattern(r"a"), Dot())) == ["a", 'b']
@test parse_one("abc", Pattern("(.)(.)(.)", 1, 3)+Eos()) == ["a", "c"]
@test parse_one("true", p"([Tt][Rr][Uu][Ee])|([Ff][Aa][Ll][Ss][Ee])"+Eos()) == ["true"]
@test parse_one("abc", Seq(Equal("a"))) == ["a"]
@test parse_one("abc", Seq(Equal("a"), Equal("b"))) == ["a", "b"]
@test parse_one("abc", Seq(p"."[1:2], Equal("c"))) == ["a", "b", "c"]
@test parse_one("abc", Seq(p"."[1:2], Equal("b"))) == ["a", "b"]
@test parse_one("abc", Seq!(p"."[1:2], Equal("c"))) == ["a", "b", "c"]
@test_throws ParserException  parse_one("abc", Seq!(p"."[1:2], Equal("b"))) == ["a", "b"]
@test parse_one("abc", Seq(p"."[1:2], p"."[1:2])) == ["a", "b", "c"]
@test parse_one("abc", Seq(p"."[1:2,:&], p"."[1:2])) == Any[["a"], ["b"], "c"]
@test parse_one("abc", Seq(p"."[1:2,:&,:?], p"."[1:2])) == Any[["a"], "b", "c"]
@test_throws ErrorException parse_one("abc", Seq(p"."[1:2,:&,:?,:x], p"."[1:2]))
@test parse_one("abc", Seq(p"."[1:2], p"."[1:2], Equal("c"))) == ["a", "b", "c"]
@test parse_one("ab", p"." + e"b") == ["a", "b"]
@test parse_one("abc", p"." + e"b" + e"c") == ["a", "b", "c"]
@test parse_one("abc", p"." + E"b" + e"c") == ["a", "c"]
@test parse_one("b", Alt(e"a", e"b", e"c")) == ["b"]
@test parse_one("b", Alt!(e"a", e"b", e"c")) == ["b"]
@test collect(parse_all("b", Trace(Alt(Epsilon(), Repeat(e"b", 0, 1))))) == Array[[], ["b"], []]

@test collect(parse_all("b", Alt(Epsilon(), Repeat(e"b", 0, 1; greedy=false)))) == Array[[], [], ["b"]]
@test parse_one("abc", p"." + (e"b" | e"c")) == ["a", "b"]
@test length(collect(parse_all("abc", p"."[0:3]))) == 4
@test length(collect(parse_all("abc", p"."[1:2]))) == 2
@test parse_one("abc", p"."[3] > tuple) == [("a", "b", "c")]
@test parse_one("abc", p"."[3] > vcat) == Any[Any["a", "b", "c"]]
@test_throws ParserException parse_one("abc", And(Equal("a"), Lookahead(Equal("c")), Equal("b")))
@test parse_one("abc", And(Equal("a"), Not(Lookahead(Equal("c"))), Equal("b"))) == Any[["a"], [], ["b"]]
@test parse_one("1.2", PFloat64()) == [1.2]
m1 = Delayed()
m1.matcher = Seq(Dot(), Opt(m1))
@test parse_one("abc", m1) == ['a', 'b', 'c']
@test collect(parse_all("abc", Repeat(Fail(); flatten=false))) == Any[[]]
@test collect(parse_all("abc", Repeat(Fail(); flatten=false, greedy=false))) == Any[[]]
@test parse_one("12c", Lookahead(p"\d") + PInt()) == [12]
@test parse_one("12c", Lookahead(p"\d") + PInt() + Dot()) == [12, 'c']
@test_throws ParserException parse_one("12c", Not(Lookahead(p"\d")) + PInt() + Dot())
@test collect(parse_all("123abc", Seq!(p"\d"[0:end], p"[a-z]"[0:end]))) == Any[Any["1", "2", "3", "a", "b", "c"]]
@test parse_one("€", p"."; debug=true) == ["€"]

# check that repeat is exactly the same as regexp

for i in 1:10
    for greedy in (true, false)
        lo = rand(0:3)
        hi = lo + rand(0:2)
        r = Regex("a{$lo,$hi}" * (greedy ? "" : "?"))
        n = rand(0:4)
        s = repeat("a", n)
        m = match(r, s)
        if m === nothing
            @test_throws ParserException parse_one(s, Repeat(Equal("a"), lo, hi; greedy=greedy))
        else
            @test length(m.match) == length(parse_one(s, Repeat(Equal("a"), lo, hi; greedy=greedy)))
        end
    end
end

# backtracking in sequences of repeats

@test parse_one("ab", Seq(Equal("a"), Equal("b"))) == ["a", "b"]
@test parse_one("abc", Dot() + Dot() + Dot()) == ['a', 'b', 'c']

for backtrack in (true, false)

    @test map(x -> [length(x[1]), length(x[2])],
              collect(parse_all("aaa",
                                Seq((Repeat(Equal("a"), 0, 3; backtrack=backtrack) > tuple),
                                    (Repeat(Equal("a"), 0, 3; backtrack=backtrack) > tuple))))) ==
    Array[[3,0],
          [2,1],[2,0],
          [1,2],[1,1],[1,0],
          [0,3],[0,2],[0,1],[0,0]]

    @test map(x -> [length(x[1]), length(x[2])],
              collect(parse_all("aaa",
                                Seq((Repeat(Equal("a"), 0, 3; backtrack=backtrack, greedy=false) > tuple),
                                    (Repeat(Equal("a"), 0, 3; backtrack=backtrack, greedy=false) > tuple))))) ==
    Array[[0,0],[0,1],[0,2],[0,3],
          [1,0],[1,1],[1,2],
          [2,0],[2,1],
          [3,0]]
end

end


text:
# Original example:
# http://gnuplot.sourceforge.net/demo/hidden2.html

using Gnuplot

x = LinRange(-10, 10, 25)
y = LinRange(-10, 10, 25)

@gsp    "set xyplane at 0"
@gsp :- "unset key" :-
@gsp :- "set palette rgbformulae 31,-11,32"
@gsp :- "set style fill solid 0.5"
@gsp :- "set cbrange [-1:1]"
@gsp :- title="Mixing pm3d surfaces with hidden-line plots"
@gsp :- "set hidden3d front"

f = [sin(-sqrt((x+5)^2+(y-7)^2)*0.5) for x in x, y in y]
z = [x*x-y*y for x in x, y in y]
@gsp :- x y f "w pm3d" x y z "w l lc rgb 'black'"


text:
function Base.show(io::IO, u::LinkedEigenfunction)
    println(io, "$(length(u.us)) linked $(eltype(u.us))")
    if !haskey(io, :compact) || !io[:compact]
        recur_io = IOContext(io, :compact => true)
        for v in u.us
            print(io, "- ")
            show(recur_io, v)
        end
        print(io, "number of set coefficients: $(length(coefficients(u)))")
    end
end

coefficients(u::LinkedEigenfunction) = coefficients(first(u.us))

function set_domain!(u::LinkedEigenfunction, domain::AbstractDomain)
    foreach(v -> set_domain!(v, domain), u.us)
    return u
end

function recompute!(u::LinkedEigenfunction)
    foreach(recompute!, u.us)
    return u
end

function set_eigenfunction!(
    u::LinkedEigenfunction{<:AbstractPlanarEigenfunction},
    coefficients::Vector,
)
    foreach(v -> set_eigenfunction!(v, coefficients), u.us)
    return u
end

function (u::LinkedEigenfunction)(
    xy::AbstractVector{T},
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    res = [zero(first(xy)) for _ in eachindex(ks)]

    for i in eachindex(u.us)
        if !(boundary ∈ u.excluded_boundaries[i])
            res .+=
                u.extra_coefficients[i] .*
                u.us[i](xy, λ, ks, boundary = boundary, notransform = notransform)
        end
    end

    return res
end


text:
using Luxor
using Test
using Random

Random.seed!(42)

e1 = PathLine(Point(100, 45))
e2 = PathMove(Point(0, 5))
e3 = PathClose()
e4 = PathCurve(Point(0, 45), Point(50, 10), Point(100, -15))

function path_test(fname)
    Drawing(500, 500, fname)
    background("grey30")
    origin()
    fontsize(40)
    sethue("gold")

    p = Path([e1, e2, e3, e4])
    drawpath(p, action=:stroke)

    pts = star(O, 200, 6, 0.5, vertices=true)
    p = Path(pts)
    drawpath(p, action=:stroke)

    textpath("hello", O, :stroke)
    p = storepath()
    drawpath(p, action=:stroke)

    textpath("hello", O + (0, 100), action=:stroke)
    p = storepath()
    drawpath(p, :stroke)

    textpath("hello", O + (0, 150), action=:fill)
    p = storepath()
    drawpath(p, action=:fill)

    sethue("cyan")
    scale(0.9)
    drawpath(polytopath(pts), action=:stroke)

    sethue("magenta")
    scale(0.8)
    drawpath(polytopath(pts), :stroke)

    bp = BezierPath([BezierPathSegment(
        rand(BoundingBox() ),
        rand(BoundingBox() ),
        rand(BoundingBox() ),
        rand(BoundingBox() ),
    ) for i in 1:10])

    sethue("white")
    setlinejoin("round")
    setline(15)
    drawbezierpath(bp, :stroke)

    setline(5)
    sethue("black")
    drawpath(bezierpathtopath(bp), :stroke)

    @test finish() == true
    println("...finished test: output in $(fname)")
end

path_test("path-test.png")


text:
using IterTools, DataStructures

function update_WF_params(wms::Array{Ty,1}, α::Array{Ty,1}, Λ::Array{Array{Int64,1},1}, y::Array{Int64,2}; debug = false) where Ty<:Number
    #y is a matrix of dimension J*K, with K the dimension of the process
    # and J the number of observations
    # Julia is in row major, so the first index indicates the row (index of observation)
    # and the second the column (index of the dimension) (as in matrix mathematical notation)
    @assert length(wms) == size(Λ, 1)

    nJ = sum(y, dims = 2) |> vec#sum_j=1^K n_ij
    nK = sum(y, dims = 1) |> vec#sum_i=1^J n_ij
    sy = sum(y)
    J = size(y,1)
    sα = sum(α)

    first_term = sum(SpecialFunctions.logfactorial.(nJ) - sum(SpecialFunctions.logfactorial.(y), dims = 2))

    function lpga(m::Array{Int64,1})
        sm = sum(m)
        second_term = lgamma_local(sα + sm)
        third_term = sum(lgamma_local.(α + m + nK))
        fourth_term = -lgamma_local(sα + sm + sy)
        fifth_term = -sum(lgamma_local.(α + m))
        return first_term + second_term + third_term + fourth_term + fifth_term
    end

    filter_ = wms .== 0

    lwms_hat = log.(wms) .+ map(lpga, Λ)

    wms_hat = exp.(lwms_hat)
    wms_hat[filter_] .= 0
    wms_hat = wms_hat |> normalise

    if debug&&any(isnan, wms_hat)
        println("NAs in update step")
        println("wms")
        println(wms)
        println("pga")
        println(map(pga, Λ))
        println("wms_hat")
        println(wms_hat)
    end

    return [m .+ nK for m in Λ], wms_hat
end

function get_next_WF_filtering_distribution_and_loglikelihood_precomputed(current_Λ, current_wms, current_time, next_time, α, sα, next_y, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff)
    predicted_Λ, predicted_wms = predict_WF_params_precomputed(current_wms, sα, current_Λ, next_time-current_time, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff)

    μν_prime_im1 = logsumexp(log(predicted_wms[k]) + logμπh_WF(α, predicted_Λ[k], vec(next_y)) for k in eachindex(predicted_Λ))# for the time being, because one only deal with single observations

    filtered_Λ, filtered_wms = update_WF_params(predicted_wms, α, predicted_Λ, next_y)

    return filtered_Λ, filtered_wms, μν_prime_im1
end

function predict_WF_params_precomputed(wms::Array{Ty,1}, sα::Ty, Λ::Array{Array{Int64,1},1}, t::Ty, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff; wm = 1) where {Ty<:Number}

    res = Accumulator{Array{Int64,1}, Float64}()

    for k in eachindex(Λ)
        res = merge(res, WF_prediction_for_one_m_precomputed(Λ[k], sα, t, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff; wm = wms[k]))
    end

    ks = keys(res) |> collect

    return ks, [res[k] for k in ks]

end

function WF_prediction_for_one_m_precomputed(m::Array{Int64,1}, sα::Ty, t::Ty, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff; wm = 1) where {Ty<:Number}
    # gm = map(x -> 0:x, m) |> vec |> x -> Iterators.product(x...)
    gm = indices_of_tree_below(m)

    function fun_n(n)
        i = m.-n
        # println(i)
        si = sum(i)
        sm = sum(m)
        return wm*(logpmmi_precomputed(i, m, sm, si, t, sα, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff) |> exp)
    end

    Dict( collect(n) => fun_n(n) for n in gm ) |> Accumulator

end

function filter_WF_adaptive_precomputation_ar(α, data, do_the_pruning::Function; silence = false, return_precomputed_terms = false, trim0 = false)
    # println("filter_WF_mem2")

    @assert length(α) == length(data[collect(keys(data))[1]])
    Δt = assert_constant_time_step_and_compute_it(data)

    smmax = values(data) |> sum |> sum
    log_ν_ar = Array{Float64}(undef, smmax, smmax)
    log_Cmmi_ar = Array{Float64}(undef, smmax, smmax)
    log_binomial_coeff_ar_offset = Array{Float64}(undef,    smmax+1, smmax+1)

    sα = sum(α)
    times = keys(data) |> collect |> sort
    Λ_of_t = Dict()
    wms_of_t = Dict()

    filtered_Λ, filtered_wms = update_WF_params([1.], α, [repeat([0], inner = length(α))], data[times[1]])
    Λ_pruned, wms_pruned = do_the_pruning(filtered_Λ, filtered_wms)

    if trim0
        Λ_pruned, wms_pruned = Λ_pruned[wms_pruned.>0], wms_pruned[wms_pruned.>0]
    end

    Λ_of_t[times[1]] = filtered_Λ
    wms_of_t[times[1]] = filtered_wms
    new_sm_max = maximum(sum.(Λ_pruned))
    precompute_next_terms_ar!(0, new_sm_max, log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset, sα, Δt)
    sm_max_so_far = new_sm_max

    for k in 1:(length(times)-1)
        if (!silence)
            println("Step index: $k")
            println("Number of components: $(length(filtered_Λ))")
        end
        last_sm_max = maximum(sum.(Λ_pruned))
        new_sm_max = last_sm_max + sum(data[times[k+1]])

        if sm_max_so_far < new_sm_max
            precompute_next_terms_ar!(sm_max_so_far, new_sm_max, log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset, sα, Δt)
            sm_max_so_far = max(sm_max_so_far,new_sm_max)
        end

        filtered_Λ, filtered_wms = get_next_filtering_distribution_precomputed(Λ_pruned, wms_pruned, times[k], times[k+1], α, sα, data[times[k+1]], log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset)
        Λ_pruned, wms_pruned = do_the_pruning(filtered_Λ, filtered_wms)
        Λ_of_t[times[k+1]] = filtered_Λ
        wms_of_t[times[k+1]] = filtered_wms
    end

    if return_precomputed_terms
        return Λ_of_t, wms_of_t, log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset, sm_max_so_far
    else
        return Λ_of_t, wms_of_t
    end

end

function get_next_filtering_distribution_precomputed(current_Λ, current_wms, current_time, next_time, α, sα, next_y, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff)
    predicted_Λ, predicted_wms = predict_WF_params_precomputed(current_wms, sα, current_Λ, next_time-current_time, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff)
    filtered_Λ, filtered_wms = update_WF_params(predicted_wms, α, predicted_Λ, next_y)

    return filtered_Λ, filtered_wms
end

filter_WF = filter_WF_adaptive_precomputation_ar

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
export simulate!
export propagate_init!

#---------------constructor first time step-------------------------------------
function propagate_init!(s::Simulation) #check seems okay #everything good.
    compute_eigen_H!(s)
    H_derivatives!(s)
    generate_therm_surf!(s)

    #set initial wavefunction
    @views s.ψ .= convert(Array{ComplexF64, 2}, s.Γ[:, s.surfp])

    # KE = sum(0.5 ./m_spread .* s.v.^2)
    # PE = s.hp[1] + sum(s.Γ[s.surfp])

    #calculate initial forces
    update_forces!(s)
    #fix position of the last layer of atoms
    s.F[399:530, :] .= 0.0
    s.v[399:530, :] .= 0.0

    #convert to Forces
    s.F .= -s.F
    if logopt == 1
        storage_init!(s)
    end
end

#stores data from initialization into first columns of storage arrays
@inline function storage_init!(s::Simulation)
    get_temperature!(s, 1)
    s.storage_K[1] = 0.5 * sum(s.v.^2 .* m_spread)
    s.storage_P[1] = s.hp[1] + sum(s.λ[s.surfp])
    rtemp = norm(s.Δ_no)
    @views s.storage_e[1, 1] = 0.5 * (mass_arr[1] + mass_arr[2]) *
     sum(((mass_arr[1].*s.v[1, :] + mass_arr[2].*s.v[2, :])/(mass_arr[1] + mass_arr[2])).^2)
    @views T_vib = 0.50 * μ * sum(((s.v[1, :] .- s.v[2, :]).*s.Δ_no/rtemp).^2)
    U_vib = F_n *(1.0 - exp(-δ_n * (rtemp -  r_0_NO)))^2
    E_vib = U_vib + T_vib
    s.storage_e[2, 1] = E_vib
    @views T_tot = 0.5 * mass_arr[1] * sum(s.v[1, :].^2) + 0.5*mass_arr[2]*sum(s.v[2, :].^2)
    @views s.storage_e[3, 1] = T_tot - T_vib - s.storage_e[1, 1]
    s.storage_e[4, 1] = sum(s.λ[s.surfp]) - sum(s.λ[s.surfpinit])       #das hier macht keinen Sinn
    @views s.storage_xno[1:3, 1] = s.x[1, :]
    @views s.storage_xno[4:6, 1] = s.x[2, :]
    @views s.storage_vno[1:3, 1] = s.v[1, :]
    @views s.storage_vno[4:6, 1] = s.v[2, :]#check seems okay
end

@inline function compute_eigen_H!(s::Simulation)
    s.H .= Symmetric(h0 .+ vm .* s.hp[3] ./ sqrt_de)
    s.H[1, 1] = s.H[1, 1] + s.hp[2] - s.hp[1]
    s.λ .= eigvals(s.H)
    s.Γ .= eigvecs(s.H)#check seems okay

end

@inbounds @inline function H_derivatives!(s::Simulation) #seems ok, fixed error
    for i in 1:Ms+1
        s.dhdea[:, i] .= view(s.Γ, 1, :) .* s.Γ[1, i]
    end
    temp = vm * s.Γ  #needs preallocation
    mul!(s.dhdv, transpose(s.Γ), temp)
    s.dhdv .= s.dhdv ./ sqrt_de
end

@inline @inbounds function generate_therm_surf!(s::Simulation) #check seems okay
    #generate thermalized populations for electronic states
    for j in 1:Ne
        s.surfp[j] = j
    end

    for j in 1:Ne
        s.occnum[s.surfp[j]] = 1
    end

    k = 1
    for j in 1:Ms+1
        if s.occnum[j] == 0
            s.surfh[k] = j
            k = k + 1
        end
    end

    for t in 1:Ne
        for k in 1:Ne*(Ms + 1 - Ne)
            hoprand = rand(Float64, 3)
            #hoprand = collect([0.2, 0.4, 0.6])
            jp = Int(floor(hoprand[1]*Ne)) + 1
            jh = Int(floor(hoprand[2] * (Ms + 1 - Ne))) + 1
            rtemp = exp(-(s.λ[s.surfh[jh]] - s.dhp_neutral[s.surfp[jp]])) / (kb * tsurf)
            if hoprand[3] < rtemp
                m = s.surfh[jh]
                s.surfh[jh] = s.surfp[jp]
                s.surfp[jp] = m
            end
        end
    end

    s.surfpinit .= deepcopy(s.surfp)
end

@inline @inbounds function update_forces!(s::Simulation) #check seems okay
    for j in 1:Ne
        s.F .= s.F .+ (s.dhdea[s.surfp[j], s.surfp[j]] .* (s.dhp_ion .- s.dhp_neutral) .+
        s.dhdv[s.surfp[j], s.surfp[j]] .* s.dhp_coup)

        # @views s.F .= s.F .+ s.dhdv[s.surfp[j], s.surfp[j]] .* s.dhp_coup
    end
    s.F .= s.F .+ s.dhp_neutral
end

#---------------------simulation main time loop for n=2:tsteps-----------------
function simulate!(s::Simulation)
    for n in 2:tsteps
        do_we_stop = v_z_condition_check(s)
        if do_we_stop == -1
            break
        end
        # Record minimial distance of z approach and orientational angle
        @views min_z_no = minimum(s.x[1:2, 3])
        if min_z_no <= s.trajzmin[1]
            s.trajzmin .= min_z_no
            s.trajtheta .= acos(s.Δ_no[3]/norm(s.Δ_no))
        end
        t_prog = round(n/tsteps * 100.0; digits=3)
        str_t_prog = string(t_prog)*"%"
        traj_prog = round(min_z_no/Å; digits=3)
        str_trajzmin_prog = string(traj_prog)
        println(str_t_prog*", "*str_trajzmin_prog)
        # get state corresponding to current surface
        @inbounds for j in 1:Ne
            s.ϕ[:, j] = view(s.Γ, :, s.surfp[j])
        end

        #Calculate nonadiabatic coupling matrix DM between adiabatic orbitals
        get_dm!(s::Simulation)
        # Get occupied (particle), unoccupied (hole) states corresponding to current surface specified by surfp
        # surfp is Ne x 1 array of particle states
        # surfh is (Ms+1-Ne) x 1 array of hole states
        get_occupied_unoccpied_states!(s)
        # Get b_{kl} where k is current many-electron adiabatic state and l is possible new
        # 	many-electron adiabatic state
        # k corresponds to occupation of states specified by surfp
        # l corresponds to occupation of states specified by a one electron-hole-pair excitation out of surfp
        get_blk2_pbmaxest!(s)
        # Calculate upper bound of maximum hopping probability, Pbmaxest

        #beginn hopping
        #Generate random number for surface hopping
        hoprand = rand()
        #! Only attempt hop if hoprand < Pbmaxest to avoid expensive calculation
        #! of real hopping elements blk
        if hoprand < s.Pbmaxest[1]
            hopping!(s, hoprand, n)
        end
        # Propagate electronic Hamiltonian
        uu = MVector{Ms+1, ComplexF64}(zeros(ComplexF64, Ms+1))
        uuu = zeros(ComplexF64, Ms+1, Ne)
        propagate_hamiltonian!(s, uu, uuu)
        #storage
        if logopt == 1
            storage_sim1!(s, n)
        end
        #propagate nuclear motion
        propagate_nuclear!(s)
        #get adiabatic eigenvalues, eigenvectors
        get_eigen_sim!(s)
        # Get dH / d Ep matrices where Ep is parameter in Newns-Anderson Hamiltonian
        # (i.e. Ep = V or Ep = E_a = (EI-EN) )
        get_dhdea_dhdv_loop!(s)

        # Calculate forces
        s.F .= 0
        update_forces!(s)
        s.F[399:530, :] .= 0
        s.v[399:530, :] .= 0
        s.F .= -s.F

        s.v .= s.vtemp .+ 0.5 ./m_spread .* s.F * dt
        if logopt == 1
            @views s.storage_xno[1:3, 1] = s.x[1, :]
            @views s.storage_xno[4:6, 1] = s.x[2, :]
            @views s.storage_vno[1:3, 1] = s.v[1, :]
            @views s.storage_vno[4:6, 1] = s.v[2, :]#check seems okay
        end
        s.nf .= s.nf .+ one(Int64)
    end
    s.storage_e[5, :] = sum(s.storage_e[1:4, :], dims=1) .- sum(s.storage_e[1:4, 1])

end
#-------------------------------------------------------------------------------

#check if the projectile has reached simulation boundaries
function v_z_condition_check(s::Simulation) #should be ok
    @views v_cond = (s.v[1, 3]*m_N + s.v[2, 3]*m_O)/(m_N + m_O)
    @views z_cond = (s.x[1, 3] + s.x[2, 3])/2.0
    if z_cond > z_end && v_cond > 0.0
        println("Boundary reached:End Simulation")
        return -1
    end
    return 1
end

@inbounds @inline function get_dm!(s::Simulation) #should be ok
    # Calculate nonadiabatic coupling elements DM between adiabatic orbitals s.dm
    # phipsi = <phi_k|psi> is overlap between "current" adiabatic state |phi_k> and electronic state |psi>
    temp_1 = transpose(s.ϕ) * s.ψ   #preallocation
    s.phipsi[1] = det(temp_1)
    # Calculate akk as defined in Tully JCP 1990
    s.akk[1] = abs(s.phipsi[1])^2 #corrected error

    # Calculate nonadiabatic coupling elements DM between adiabatic orbitals

    temp_2 = s.dhp_ion .- s.dhp_neutral #preallocation
    temp_2 .= temp_2 .* s.v
    temp_2_sum = sum(temp_2)
    s.dm .= temp_2_sum .* s.dhdea

    temp_3 = s.dhp_coup .* s.v
    temp_3_sum = sum(temp_3)
    s.dm .= s.dm .+ temp_3_sum .* s.dhdv

    temp_rep1 = repeat(s.λ, 1, Ms+1)
    temp_rep2 = temp_rep1' .- temp_rep1
    s.dm .= s.dm ./ temp_rep2
    for j in 1:Ms+1
        s.dm[j,j] = zero(eltype(s.dm))
    end #should be ok
end

@inline @inbounds function get_occupied_unoccpied_states!(s::Simulation) #should be ok
    s.occnum .= 0
    for j in 1:Ne
        s.occnum[s.surfp[j]] = 1
    end
    k = 1
    for j in 1:Ms+1
        if s.occnum[j] == 0
            s.surfh[k] = j
            k = k+1
        end
    end #should be ok
end

@inline @inbounds function get_blk2_pbmaxest!(s::Simulation) #should be ok
    rtemp = abs(s.phipsi[1] * s.phipsi[1])
    if rtemp > one(eltype(rtemp))
        rtemp = 1.0
    end
    for jp in 1:Ne
        for jh in 1:Ms+1-Ne
            s.blk2[jp, jh] = 2.0 * abs(s.phipsi[1] * s.dm[s.surfp[jp], s.surfh[jh]])
        end
    end

    s.Pbmaxest[1] = sqrt(sum(s.blk2.^2)) * sqrt(1.0 - rtemp)/s.akk[1] * dt * Float64(thop) #should be ok
end

#---------------------start hopping! hoprand < pbmaxest -----------------------
function hopping!(s::Simulation, hoprand::Float64, n::Int64)
    hopping!_pbmaxest!(s, hoprand, n)
end

function hopping!_pbmaxest!(s::Simulation, hoprand::Float64, n::Int64)
    get_blk_akl!(s)
    get_pb!(s)

    s.storage_phop[n] = s.Pb[Ne*(Ms+1-Ne)]
    #attempt hop conditioned on random number
    if hoprand < s.Pb[Ne*(Ms+1-Ne)]
        print("Hopping! ")
        hopping!_pbmaxest!_attempt!(s, hoprand, n)
    end
end

@inline function get_blk_akl!(s::Simulation) #seems ok
    ctemp1 = Matrix{ComplexF64}(undef, Ne, Ne)
    for jp in 1:Ne
        for jh in 1:Ms+1-Ne
            s.surfpnew .= s.surfp
            s.surfpnew[jp] = s.surfh[jh]
            for j in 1:Ne
                @views s.ϕ[:, j] = s.Γ[:, s.surfpnew[j]]
            end
            #Ctemp = <phi_l|psi> is overlap between "new" adiabatic state
            #|phi_l> and electronic state |psi>
            mul!(ctemp1, transpose(s.ϕ) , s.ψ)
            ctemp = det(ctemp1)
            s.akl[1] = s.phipsi[1] * conj(ctemp)
            @views s.blk[jp, jh] = 2.0 * real(s.akl[1] * s.dm[s.surfp[jp], s.surfh[jh]])
        end
    end
end

@inline function get_pb!(s::Simulation) #seems ok
    rtemp = 0.0
    for jh in 1:Ms+1-Ne
        for jp in 1:Ne
            s.Pb[(jh-1) * Ne + jp] = rtemp
            @views temp = s.blk[jp,jh]
            if temp <= 0.0
                temp = 0.0
            end
            s.Pb[(jh-1) * Ne + jp] = s.Pb[(jh-1)*Ne + jp] + temp
            rtemp = s.Pb[(jh - 1) * Ne + jp]
        end
    end
    s.Pb .= s.Pb ./ s.akk[1] * dt * Float64(thop)
end
#-----------------start hopping! hoprand < pbmaxest && hoprand < Pbmax ---------
function hopping!_pbmaxest!_attempt!(s::Simulation, hoprand::Float64, n::Int64)
    hopstate = which_state(s, hoprand) #ok
    jh= which_orbital_jh(hopstate) #ok
    jp= which_orbital_jp(hopstate) #ok
    rhop = which_direction(s, jp, jh) #needs durther testing, not sure about rescaling
    e_el_old, e_el_new, k_rhop, k_temp = check_energy1(s, jp, jh, rhop)
    bool_e = check_energy2(e_el_old, e_el_new, k_rhop, k_temp)
    s.attnum .= s.attnum .+ one(Int64)

    if bool_e == true
        update_vel_populations!(s, jp, jh, rhop, k_temp, k_rhop, e_el_old, e_el_new, n) #needs more storage
    end
end

@inline function which_state(s::Simulation, hoprand::Float64)::Int64 #ok
    hopstate = one(Int64)
    while s.Pb[hopstate] <= hoprand
        hopstate = hopstate + one(Int64)
    end
    return hopstate
end

@inline function which_orbital_jh(hopstate::Int64)::Int64 #ok
    jh = Int(floor((hopstate - 1) / Ne)) + one(Int64)
    return jh
end
@inline function which_orbital_jp(hopstate::Int64)::Int64 #ok
    jp = ((hopstate - one(Int64)) % Ne) + one(Int64)
    return jp
end

@inline function which_direction(s::Simulation, jp::Int64, jh::Int64)::float_array #needs checking
    @views temp = s.dhdea[s.surfp[jp], s.surfh[jh]] .* (s.dhp_ion .- s.dhp_neutral)
    @views temp .= temp .+ s.dhdv[s.surfp[jp], s.surfh[jh]] .* s.dhp_coup
    temp2 = sqrt(sum(temp.^2)/sum(s.v .^2))
    temp .= temp ./ temp2
    temp3 = sum(temp .* s.v)
    temp .= temp .* copysign(1.0, temp3)
    return temp
end

@inline function check_energy1(s::Simulation, jp::Int64, jh::Int64, rhop::float_array) # seems ok
    s.surfpnew .= s.surfp
    s.surfpnew[jp] = s.surfh[jh]
    e_el_old = sum(s.λ[s.surfp])
    e_el_new = sum(s.λ[s.surfpnew])
    k_rhop = 0.5 * sum(1.0 ./ m_spread .* rhop .^2)
    k_temp = sum(s.v .* rhop)

    return e_el_old, e_el_new, k_rhop, k_temp
end

@inline function check_energy2(e_el_old::Float64, e_el_new::Float64, k_rhop::Float64, k_temp::Float64)::Bool
    bool_e = k_temp^2 - 4.0*k_rhop * (e_el_new - e_el_old) > 0.0
    return bool_e
end

@inline function update_vel_populations!(s::Simulation, jp::Int64, jh::Int64, rhop::float_array,
    k_temp::Float64, k_rhop::Float64, e_el_old::Float64, e_el_new::Float64, n::Int64)
    rtemp = s.vscale[1]
    s.vscale[1] = (k_temp - sqrt(k_temp^2 - 4.0 * k_rhop * (e_el_new - e_el_old)))/2.0/k_rhop
    @views k_no = 0.50 * sum(view(m_spread, 1:2, :).* view(s.v, 1:2, :).^2)
    @views k_au = 0.50 * sum(view(m_spread, 3:N+2, :) .* view(s.v, 3:N+2, :).^2)
    s.v .= s.v .- s.vscale[1] .* (1.0 ./m_spread .* rhop)
    s.storage_deltaKNO[s.exnum[1]+1] = k_no - 0.5*sum(view(m_spread, 1:2, :) .* view(s.v, 1:2, :) .^2)
    s.storage_deltaKAu[s.exnum[1]+1] = k_au - 0.5*sum(view(m_spread, 3:N+2, :) .* view(s.v, 3:N+2, :) .^2)
    s.storage_state[2*s.exnum[1] + 1] = s.surfp[jp]
    s.storage_state[2*s.exnum[1] + 2] = s.surfh[jh]
    s.surfp .= s.surfpnew
    s.storage_hoptimes[s.exnum[1]+1] = n
    s.exnum .= s.exnum .+ one(Int64)
    #endif
end
#-----------------END hopping! hoprand < pbmaxest && hoprand < Pbmax ---------

@inline function propagate_hamiltonian!(s::Simulation, uu::MVector{Ms+1, ComplexF64},
     uuu::Matrix{ComplexF64}) # seems ok
    s.ψ .= transpose(s.Γ) * s.ψ

    uu .= exp.(-1im * s.λ * dt/hbar)
    uuu .= repeat(uu, 1, Ne)
    s.ψ .= uuu .* s.ψ
    s.ψ .= s.Γ * s.ψ
end

#storage_op_e
@inline function storage_sim1!(s::Simulation, n::Int64) #seems ok
    get_temperature!(s, n)
    s.storage_K[n] = 0.5 * sum(s.v .^2 .* m_spread)
    s.storage_P[n] = s.hp[1] + sum(s.λ[s.surfp])
    @views s.storage_aop[n] = sum(abs.(s.ψ[1, :]).^2) #rhoa
    store_temp = sum(abs.(transpose(s.Γ) * s.ψ).^2, dims=2)
    s.storage_op[:, n] .= store_temp[:, 1] #psip

    rtemp = norm(s.Δ_no)
    s.storage_e[1, n] = 0.5 * (mass_arr[1] + mass_arr[2]) *
     sum(((mass_arr[1].*s.v[1, :] + mass_arr[2].*s.v[2, :])/(mass_arr[1] + mass_arr[2])).^2)
    @views T_vib = 0.50 * μ * sum(((s.v[1, :] .- s.v[2, :]).*s.Δ_no/rtemp).^2)
    U_vib = F_n *(1.0 - exp(-δ_n * (rtemp -  r_0_NO)))^2
    E_vib = U_vib + T_vib
    s.storage_e[2, n] = E_vib
    @views T_tot = 0.5 * mass_arr[1] * sum(s.v[1, :].^2) + 0.5*mass_arr[2]*sum(s.v[2, :].^2)
    s.storage_e[3, n] = T_tot - T_vib - s.storage_e[1, n]
    s.storage_e[4, n] = sum(s.λ[s.surfp]) - sum(s.λ[s.surfpinit])

end
@inline function propagate_nuclear!(s::Simulation)
    s.vdot .= 1.0./m_spread .* s.F
    s.x .= s.x .+ s.v * dt .+ 0.50*s.vdot * dt^2
    s.Δ_no .=  s.x[1, :] - s.x[2, :]
    s.vtemp .= s.v .+ 0.50*s.vdot*dt #ARRAY ALLOCATION !!!! solve later
end

@inline function get_eigen_sim!(s::Simulation)
    s.hp .= get_E_all(s)    #update energy surface
    get_F_all(s)            #update partial derivatives dhp_i
    s.H .= h0 .+ vm .* s.hp[3] ./sqrt_de
    s.H[1, 1] = s.H[1, 1] + s.hp[2] - s.hp[1]
    Γ_hold = deepcopy(s.Γ) #ARRAY ALLOCATION !!!! solve later
    s.λ .= eigvals(s.H)
    s.Γ .= eigvecs(s.H)
    temp = zeros(Float64, Ms+1)
    @inbounds for j in 1:Ms+1
        @views temp .= abs.(transpose(Γ_hold) * s.Γ[:, j])
        Γmaxloc = argmax(temp)
        temp2 = dot(Γ_hold[:, Γmaxloc], s.Γ[:, j])
        s.Γ[:, j] = s.Γ[:, j]/copysign(1.0, temp2)
    end
end

function get_dhdea_dhdv_loop!(s::Simulation)
    s.blk .= 0.0
    @inbounds for j in 1:Ms+1
        s.dhdea[:, j] = s.Γ[1, :] .* s.Γ[1, j]
    end
    temp = (vm * s.Γ) #ARRAY ALLOCATION !!!! solve later
    mul!(s.dhdv, transpose(s.Γ), temp)
    s.dhdv .= s.dhdv ./sqrt_de
end

function get_temperature!(s::Simulation, n::Int)
    speed = [norm(s.v[i, :]) for i in 3:398]
    s.storage_temp[n] = 2*m_au*sum(speed.^2)/(3*(N*0.75-3)*kb)
end


text:
function model_uptake(
    conv::Function = expconv;
    t::AbstractVector,
    params::NamedTuple,
    ca::AbstractVector,
)
    @extract (fp, ps, vp) params
    Tp = vp / (fp + ps)
    E = ps / (fp + ps)
    ca_conv_1 = zeros(length(t))
    for idx = 2:length(ca)
        ca_conv_1[idx] = ca_conv_1[idx-1] + (t[idx] - t[idx-1]) * ca[idx]
    end
    ct = fp .* ((1 - E) .* expconv(ca, 1 / Tp, t) .+ E .* ca_conv_1)
    return ct
end

function fit_uptake_nls(
    conv::Function = expconv;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    fp, ps, vp = (fill(NaN, volume_size...) for _ = 1:3)
    model(x, p) = model_uptake(; t = x, ca, params = (fp = p[1], ps = p[2], vp = p[3]))
    lls_est = fit_uptake_lls(; t, ca, ct, mask).est
    init_fp, init_ps, init_vp = select(lls_est, (:fp, :ps, :vp))
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        initialvalue = [init_fp[idx], init_ps[idx], init_vp[idx]]
        (fp[idx], ps[idx], vp[idx]) = curve_fit(model, t, ct[idx, :], initialvalue).param
    end
    return (; est = (; fp, ps, vp))
end

function fit_uptake_lls(;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    x1, fp, x3, ps, vp = (fill(NaN, volume_size...) for _ = 1:5)

    M = zeros(num_timepoints, 3)
    M[:, 2] = cumul_integrate(t, ca, TrapezoidalFast())
    M[:, 3] = cumul_integrate(t, M[:, 2], TrapezoidalFast())
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        M[:, 1] = -cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        (x1[idx], fp[idx], x3[idx]) = M \ ct[idx, :]
    end
    denum = @. x1 * fp - x3
    @. ps = fp * x3 / denum
    @. vp = fp^2 / denum
    @. ps[denum==0] = 0
    @. vp[denum==0] = 0
    return (; est = (; fp, ps, vp))
end


text:
module G0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    -3.02912483102257e-6*l1.^5 - 1.78242802467006e-5*l1.^4.*l2 + 3.29870469288301e-5*l1.^4.*l3 + 9.37009405046133e-8*l1.^4 - 2.7515585012037e-5*l1.^3.*l2.^2 + 0.000126439665488924*l1.^3.*l2.*l3 + 1.41766156856612e-7*l1.^3.*l2 - 0.00012928389889518*l1.^3.*l3.^2 - 5.17114922403405e-7*l1.^3.*l3 + 3.82562374888515e-10*l1.^3 - 4.68312335895801e-7*l1.^2.*l2.^3 + 8.41756449126506e-5*l1.^2.*l2.^2.*l3 - 3.5556688146475e-7*l1.^2.*l2.^2 - 0.00027422615773196*l1.^2.*l2.*l3.^2 + 2.84603410436225e-7*l1.^2.*l2.*l3 + 8.24471484365873e-10*l1.^2.*l2 + 0.000220912928627533*l1.^2.*l3.^3 + 6.34804124641222e-7*l1.^2.*l3.^2 - 1.97385538081975e-9*l1.^2.*l3 + 4.72058579046329e-13*l1.^2 + 1.03644793616222e-5*l1.*l2.^4 - 4.03814334300281e-5*l1.*l2.^3.*l3 - 2.54281526051452e-8*l1.*l2.^3 - 2.4037236653465e-5*l1.*l2.^2.*l3.^2 + 7.88483312239122e-7*l1.*l2.^2.*l3 - 7.16189502159426e-10*l1.*l2.^2 + 0.00019939210057147*l1.*l2.*l3.^3 - 1.07291993224609e-6*l1.*l2.*l3.^2 - 2.21211033699355e-10*l1.*l2.*l3 + 1.17459514671635e-12*l1.*l2 - 0.00016055210411278*l1.*l3.^4 - 6.65739999913559e-8*l1.*l3.^3 + 2.0885116396423e-9*l1.*l3.^2 - 2.1200173246156e-12*l1.*l3 - 2.90488121299558e-5*l1 + 2.92654312504004e-6*l2.^5 - 2.49701624716648e-5*l2.^4.*l3 + 6.91116586537353e-8*l2.^4 + 7.00070469556942e-5*l2.^3.*l3.^2 - 2.50342189017214e-7*l2.^3.*l3 - 7.11777798043519e-11*l2.^3 - 6.1726059406604e-5*l2.^2.*l3.^3 - 2.02753578297076e-8*l2.^2.*l3.^2 + 9.30013738403587e-10*l2.^2.*l3 - 4.12072853349556e-13*l2.^2 - 1.9256742517163e-5*l2.*l3.^4 + 3.72132609316059e-7*l2.*l3.^3 - 8.17394230662383e-10*l2.*l3.^2 - 3.53791896171088e-13*l2.*l3 + 1.77973142556746e-5*l2 + 3.60656416727372e-5*l3.^5 - 7.63807770899228e-8*l3.^4 - 4.2573148239698e-10*l3.^3 + 1.23925501369454e-12*l3.^2 + 1.14158920149061e-5*l3 - 2.4957191258213e-8
end
end


text:
using QContinuousMeasurement
using Test
using Random

@testset "kind=0" begin
    params = ModelParameters(Nj=2, kind=0, kcoll=0)
    model = LocalDephasingModel(params)
    state = blockdiag_css(params.Nj)

    model2 = CollectiveDephasingModel(params)
    state2 = fixedj_css(params.Nj)

    seed = 1

    Random.seed!(seed)
    res = simulate_trajectory(model, state)

    Random.seed!(seed)
    res2 = simulate_trajectory(model2, state2)

    for k in keys(res)
        @test getfield(res, k) ≈ getfield(res2, k)
    end
end


text:
"""
    taxonomy(tx::NCBITaxon)

    Builds a taxonomy with `tx` at the root by running
    breadth-first search.

"""
function taxonomy(tx::AbstractTaxon, txnmytype::Type{T}) where {T <: AbstractTaxonomy}
    root = TaxonNode(tx, TaxonomicRank(rank(tx.val),0)) # no good

    txnmy = txnmytype(root)
    _taxonomy!(root, txnmy)

    return txnmy
end

"""
    _taxonomy!(root::TaxonNode{T}, tree::T)

    all natural, organic, free-range depth-first search.
"""
function _taxonomy!(root::TaxonNode{T}, tax::FixedRankTaxonomy) where {T <: NCBITaxonWrapper}
    currentnode = root
    for child in NCBITaxonomy.children(taxon(currentnode))
        childrank = rank(child)
        if _isfixedrank(childrank, ranks(tax))
            newnode = _addnode!(NCBITaxonWrapper(child), childrank, currentnode)
            _taxonomy!(newnode, tax)
        end
    end
end

function _isfixedrank(rank, ranks)
    retval =  rank ∈ [r.name for r in ranks] # slow
    return retval
end

function _taxonomy!(root::TaxonNode{T}, tax::Type{FlexibleRankTaxonomy}) where {T <: NCBITaxonWrapper}
    currentnode = root
    for child in NCBITaxonomy.children(taxon(currentnode))
        newnode = _addnode!(NCBITaxonWrapper(child), currentnode)
        _taxonomy!(newnode, tax)
    end
end


text:
abstract type AbstractKernel end

function covariance!(b, k, θ, X, Y)
    distance!(b, k, θ, X, Y)
    b .= (x -> radial(k, θ, x)).(b)
    return b
end

Zygote.@adjoint function covariance!(b, k, θ, X, Y)
    tmp, back = pullback(b, k, θ, X, Y) do b, k, θ, X, Y
        d = distance!(b, k, θ, X, Y)
        (x -> radial(k, θ, x)).(d)
    end
    b .= tmp
    return b, c -> back(c)
end

function jitter!(b, ϵ; rescale=true)
    for i in diagind(b)
        b[i] += ϵ
    end
    if rescale
        b ./= b[1, 1]
    end
    return b
end

Zygote.@adjoint function jitter!(b, ϵ; rescale=true)
    tmp, back = pullback(b, ϵ) do b, ϵ
        c = b + UniformScaling(ϵ)
        if rescale
            c ./ c[1]
        else
            c
        end
    end
    b .= tmp
    return b, c -> back(c)
end

sqeuclidean!(b, X, Y, ℓ) = sqeuclidean!(b, X ./ ℓ, Y ./ ℓ)
function sqeuclidean!(b, X, Y)
    mul!(b, X', Y, -2, 0)
    r1 = sum!(abs2, similar(X, 1, size(X, 2)), X)
    r2 = sum!(abs2, similar(Y, 1, size(Y, 2)), Y)
    b .+= r1' .+ r2
    return b
end

Zygote.@adjoint function sqeuclidean!(b, X, Y)
    tmp, back = pullback(X, Y) do A, B
        pairwise(SqEuclidean(), A, B; dims=2)
    end
    b .= tmp
    return b, c::AbstractMatrix -> (nothing, back(c)...)
end

struct Gaussian <: AbstractKernel end
distance!(b, ::Gaussian, θ, X, Y) = sqeuclidean!(b, X, Y, θ[1])
radial(::Gaussian, _, τ) = exp(-τ / 2)

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ParallelCoverallsSetup
using Test

@testset "ParallelCoverallsSetup.jl" begin
    @test_nowarn i_am_fully_tested()
    @test_nowarn i_am_partially_tested()
    @test_throws AssertionError i_will_throw()
end


text:
using Coverage

coverage = process_folder()

print("Overall coverage: ")
println(get_summary(coverage))

covered_lines, total_lines = get_summary(coverage)

for c in coverage
  print(c.filename * ": ")
  println(get_summary(c))
end


text:
using DataArrays, DataFrames, ReTerms
using Base.Test

include(joinpath(dirname(@__FILE__),"data.jl"))

include("paramlowertriangular.jl")
include("scalarReTerm.jl")
include("vectorReTerm.jl")
include("pls.jl")
include("simulation.jl")


text:
@testset "Metadata" begin
    #HalfMBB - 2x2

    # 7   8   9
    #  x--x--x      x--x--x
    #  |  |  |      |C3|C4|
    #  x--x--x      x--x--x
    # 4| 5|  |6     |C1|C2|
    #  x--x--x      x--x--x
    # 1   2   3

    #Node coords:
    #Ferrite.Node{2,Float64}([0.0, 0.0])
    #Ferrite.Node{2,Float64}([1.0, 0.0])
    #Ferrite.Node{2,Float64}([2.0, 0.0])
    #Ferrite.Node{2,Float64}([0.0, 1.0])
    #Ferrite.Node{2,Float64}([1.0, 1.0])
    #Ferrite.Node{2,Float64}([2.0, 1.0])
    #Ferrite.Node{2,Float64}([0.0, 2.0])
    #Ferrite.Node{2,Float64}([1.0, 2.0])
    #Ferrite.Node{2,Float64}([2.0, 2.0])

    #Cells
    #Ferrite.Cell{2,4,4}((1, 2, 5, 4))
    #Ferrite.Cell{2,4,4}((2, 3, 6, 5))
    #Ferrite.Cell{2,4,4}((4, 5, 8, 7))
    #Ferrite.Cell{2,4,4}((5, 6, 9, 8))

    # 7   8   9
    #  x--x--x      x--x--x
    #  |  |  |      |C3|C4|
    #  x--x--x      x--x--x
    # 4| 5|  |6     |C1|C2|
    #  x--x--x      x--x--x
    # 1   2   3

    #Node dofs
    ##1 #2 #3  #4 #5 #6  #7  #8  #9
    #1  3   9  7  5  11  15  13  17
    #2  4  10  8  6  12  16  14  18

    #Cell dofs
    # #1  #2  #3  #4
    # 1   3   7   5
    # 2   4   8   6
    # 3   9   5  11
    # 4  10   6  12
    # 5  11  13  17
    # 6  12  14  18
    # 7   5  15  13
    # 8   6  16  14

    problem = HalfMBB(Val{:Linear}, (2, 2), (1.0, 1.0), 1.0, 0.3, 1.0)

    coords = [
        (0.0, 0.0),
        (1.0, 0.0),
        (2.0, 0.0),
        (0.0, 1.0),
        (1.0, 1.0),
        (2.0, 1.0),
        (0.0, 2.0),
        (1.0, 2.0),
        (2.0, 2.0),
    ]

    for (i, node) in enumerate(problem.ch.dh.grid.nodes)
        @test node.x.data == coords[i]
    end

    cells = [
        Ferrite.Cell{2,4,4}((1, 2, 5, 4)),
        Ferrite.Cell{2,4,4}((2, 3, 6, 5)),
        Ferrite.Cell{2,4,4}((4, 5, 8, 7)),
        Ferrite.Cell{2,4,4}((5, 6, 9, 8)),
    ]
    @test problem.ch.dh.grid.cells == cells

    node_dofs = [
        1 3 9 7 5 11 15 13 17
        2 4 10 8 6 12 16 14 18
    ]
    @test problem.metadata.node_dofs == node_dofs

    cell_dofs = [
        1 3 7 5
        2 4 8 6
        3 9 5 11
        4 10 6 12
        5 11 13 17
        6 12 14 18
        7 5 15 13
        8 6 16 14
    ]
    @test problem.metadata.cell_dofs == cell_dofs

    dof_cells = problem.metadata.dof_cells
    cell_dofs = problem.metadata.cell_dofs
    for i in 1:Ferrite.ndofs(problem.ch.dh)
        d_cells = dof_cells[i]
        for c in d_cells
            (cellid, localdof) = c
            @test i == cell_dofs[localdof, cellid]
        end
    end

    node_first_cells = [
        (1, 1), (1, 2), (2, 2), (1, 4), (1, 3), (2, 3), (3, 4), (3, 3), (4, 3)
    ]
    # First node is the first cell's first node.
    @test problem.metadata.node_cells[1] == [(1, 1)]
    # Second node is the first cell's second node,
    # and the second cell's first node.
    @test problem.metadata.node_cells[2] == [(1, 2), (2, 1)]
    # Third node is the second cell's second node.
    @test problem.metadata.node_cells[3] == [(2, 2)]
    # Fourth node is the first cell's fourth node,
    # and the third cell's first node.
    @test problem.metadata.node_cells[4] == [(1, 4), (3, 1)]
    # Fifth node is the first cell's third node, the second cell's fourth node,
    # the third cell's second node and the fourth cell's first node.
    @test problem.metadata.node_cells[5] == [(1, 3), (2, 4), (3, 2), (4, 1)]
    # Sixth node is the second cell's third node,
    # and the fourth cell's second node.
    @test problem.metadata.node_cells[6] == [(2, 3), (4, 2)]
    # Seventh node is the third cell's fourth node.
    @test problem.metadata.node_cells[7] == [(3, 4)]
    # Eigth node is the third cell's third node,
    # and the fourth cell's fourth node.
    @test problem.metadata.node_cells[8] == [(3, 3), (4, 4)]
    # Ninth node is the fourth cell's third node.
    @test problem.metadata.node_cells[9] == [(4, 3)]
end


text:
"
A multivariate distribution made up of independent component distributions.
"
immutable MvIndependent <: ContinuousMultivariateDistribution
    components::Array{ContinuousUnivariateDistribution, 1}
end

length(d::MvIndependent) = length(d.components)
_logpdf{T<:Real}(d::MvIndependent, x::AbstractVector{T}) = sum(Float64[logpdf(d.components[i], x[i]) for i in 1:length(d)])
function _rand!{T<:Real}(d::MvIndependent, x::AbstractVector{T})
    x = T[rand(c) for c in d.components]
end


text:
eachcoordinate(vsr::VerySimpleRaster) = product(vsr.xs, vsr.ys)
coordinates(vsr::VerySimpleRaster) = collect(flatten(eachcoordinate(vsr)))


text:
using LoweredCodeUtils
using Test

# @testset "Ambiguity" begin
#     @test isempty(detect_ambiguities(LoweredCodeUtils, LoweredCodeUtils.JuliaInterpreter, Base, Core))
# end

include("signatures.jl")
include("codeedges.jl")


text:
module MUMIMO

using LinearAlgebra
using LLLplus
using Plots
using Printf
using Statistics

export
    mimoUplink,
    scodes,
    mimo_slice,
    siso_demod

include("mimoUplink.jl")
include("scodes.jl")

end # LLLplus module

max_stars_repo_name: christianpeel/MUMIMO

max_issues_repo_name: christianpeel/MUMIMO

max_forks_repo_name: christianpeel/MUMIMO


text:
mutable struct BaseEnumeration{T} <: AbstractSolver{T}
    xstar::Array{T, 2}
    vstar::T
    status::Status
    colindexs::Array{Int, 2}
    removedLine::BitVector
    n::Int
    function BaseEnumeration(lp::AbstractLP{T}) where T
        n = nvar(lp)
        m = ncons(lp)
        xstar = Array{T, 2}(undef, n + 1, 0) #a made up variable
        vstar = typemax(T)
        status = Unknown
        colindexs = Array{Int, 2}(undef, m + 1, 0) #made up constraint
        removedLine = falses(m)
        return new{T}(xstar, vstar, status, colindexs, removedLine, n)
    end
end
const BE{T} = BaseEnumeration{T}
function (be::BE{T})(lp::AbstractLP{T}; verbose::Bool = false) where T
    A = getA(lp)
    b = getb(lp)
    c = getc(lp)
    r1 = rank(A)
    r2 = rank([A b])
    if r1 != r2
        verbose && println("Infeasible since rank(A) != rank([A b])")
        #this is Infeasible
        be.status = lp.status = Infeasible
    end
    if r1 < size(A, 1)
        #there are redundent constraint
        removedLine = reduceproblem(A, b)
        @assert sum(removedLine) == r1 "Impossible to reduce the problem to a feasible form"
        be.removedLine[:] = removedLine
        verbose && println("A not full rank, $(sum(!, be.removedLine)) rows removed from problem.")
        A = A[removedLine, :]
        b = b[removedLine]
    end
    m, n = size(A)
    A = [A zeros(T, m); c' one(T)]
    b = [b ; abignumber(T)]
    c = [c; zero(T)]
    verbose && @show A
    verbose && @show b
    verbose && @show c
    m, n = size(A)
    C =  combinations(1:n, m)
    verbose && println("will loop over the $(length(C)) possible basis.")
    for (iter, comb) in enumerate(C)
        B = @view A[:, comb]
        if rank(B) == m
            xB = B\b
            if all(xB .>= - tol(T))
                cB = c[comb]
                vtemp = dot(cB, xB)
                if vtemp < be.vstar - tol(T)
                    #found a better solution
                    be.vstar = vtemp
                    be.xstar = zeros(T, be.n + 1, 1)
                    be.xstar[comb, end] = xB
                    be.colindexs = reshape(comb, m, 1)
                elseif be.vstar - tol(T) <= vtemp <= be.vstar + tol(T)
                    #found a similar solution
                    be.xstar = [be.xstar zeros(T, be.n + 1)]
                    be.xstar[comb, end] = xB
                    be.colindexs = [be.colindexs comb]
                end
                verbose && println("basis $iter feasible, solution $vtemp")
            else
                verbose && println("basis $iter full rank but not feasible xB = $xB")
            end
        else
            verbose && println("basis $iter not full rank")
        end
    end
    if be.vstar == typemax(T)
        be.status = Infeasible
        lp.status = Infeasible
    elseif be.vstar == -abignumber(T)
        be.status = Unbounded
        lp.status = Unbounded
    else
        be.status = Optimal
        lp.status = Optimal
        lp.xstar = be.xstar[1:end - 1, 1]
        lp.vstar = be.vstar
    end
    return be.status
end

function reduceproblem(A::Matrix{T}, b::Vector{T})::BitVector where T
    M = [A b]
    m, n = size(M)
    reducedprob = trues(m)
    for i in 1:m
        j = findfirst(!iszero, M[i, :])
        @assert j != n "unexpected error in reduceproblem, line $i, M = $M"
        isnothing(j) ? (reducedprob[i] = false) : pivot!(M, i, j)
    end
    return reducedprob
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
"""
    local_clustering_coefficient(g, v)
    local_clustering_coefficient(g, vs)

Return the [local clustering coefficient](https://en.wikipedia.org/wiki/Clustering_coefficient)
for node `v` in graph `g`. If a list of vertices `vs` is specified, return a vector
of coefficients for each node in the list.
"""
function local_clustering_coefficient(g::AbstractGraph, v::Integer)
    ntriang, alltriang = local_clustering(g, v)

    return alltriang == 0 ? 0. : ntriang * 1.0 / alltriang
end
local_clustering_coefficient(g::AbstractGraph, vs = vertices(g)) =
    [local_clustering_coefficient(g, v) for v in vs]

@doc_str """
    local_clustering(g, v)
    local_clustering(g, vs)

Return a tuple `(a, b)`, where `a` is the number of triangles in the neighborhood
of `v` and `b` is the maximum number of possible triangles. If a list of vertices
`vs` is specified, return two vectors representing the number of triangles and
the maximum number of possible triangles, respectively, for each node in the list.

This function is related to the local clustering coefficient `r` by ``r=\frac{a}{b}``.
"""
function local_clustering(g::AbstractGraph, v::Integer)
    k = degree(g, v)
    k <= 1 && return (0, 0)
    neighs = neighbors(g, v)
    c = 0
    for i in neighs, j in neighs
        i == j && continue
        if has_edge(g, i, j)
            c += 1
        end
    end
    return is_directed(g) ? (c , k*(k-1)) : (div(c,2) , div(k*(k-1),2))
end
function local_clustering(g::AbstractGraph, vs = vertices(g))
    ntriang = zeros(Int, length(vs))
    nalltriang = zeros(Int, length(vs))
    i = 0
    for (i, v) in enumerate(vs)
        ntriang[i], nalltriang[i] = local_clustering(g, v)
    end
    return ntriang, nalltriang
end

"""
    triangles(g[, v])
    triangles(g, vs)

Return the number of triangles in the neighborhood of node `v` in graph `g`.
If a list of vertices `vs` is specified, return a vector of number of triangles
for each node in the list. If no vertices are specified, return the number
of triangles for each node in the graph.
"""
triangles(g::AbstractGraph, v::Integer) = local_clustering(g, v)[1]
triangles(g::AbstractGraph, vs = vertices(g)) = local_clustering(g, vs)[1]

"""
    global_clustering_coefficient(g)

Return the [global clustering coefficient](https://en.wikipedia.org/wiki/Clustering_coefficient)
of graph `g`.
"""
function global_clustering_coefficient(g::AbstractGraph)
    c = 0
    ntriangles = 0
    for v in vertices(g)
        neighs = neighbors(g, v)
        for i in neighs, j in neighs
            i == j && continue
            if has_edge(g, i, j)
                c += 1
            end
        end
        k = degree(g, v)
        ntriangles += k*(k-1)
    end
    ntriangles == 0 && return 1.
    return c / ntriangles
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
struct Rewriter <: AbstractRewriter
    rewriters::Dict{AbstractTheory,AbstractRewriter}
    Rewriter() = new(Dict{AbstractTheory,AbstractRewriter}())
end

function Base.push!(rw::Rewriter, (p, b)::Pair{<:AbstractTerm})
    th = theory(p)
    haskey(rw.rewriters, th) || (rw.rewriters[th] = rewriter(th))
    push!(rw.rewriters[th], p => b)
    rw
end

function compile(rw::Rewriter)
    fn_name = gensym(:rewrite)

    rewriter_exprs = Expr[]

    rewriters = Dict()
    for (th, rewriter) ∈ rw.rewriters
        fn, expr = compile(rewriter)
        rewriters[th] = fn
        push!(rewriter_exprs, expr)
    end

    tree = gensym(:rewriter_TREE)

    fn_name, quote
        $(rewriter_exprs...)
        $tree = Dict($((:($th => $fn) for (th, fn) ∈ rewriters)...))
        function $fn_name(t)
            while true
                th = $theory(t)
                t = map($fn_name, t)
                $haskey($tree, th) || return t

                t′ = $tree[th](t)
                t′ === nothing && return t
                t = t′
            end
        end
    end
end


text:
using StatsBase
using ClusterValidityIndices
using Logging
using Plots

# Plotting options
theme(:dark)    # Plotting style
dpi = 300       # Plotting dots-per-inch
# gr()            # GR backend (default for Plots.jl)
unicodeplots()

# Include the library definitions
# include(projectdir("julia/lib_sim.jl"))

# Set the logging level to Info and standardize the random seed
LogLevel(Logging.Info)
# Random.seed!(0)
# data_dir(args...) = projectdir("work/data/meta_icvi", args...)
# results_dir(args...) = projectdir("work/results/meta_icvi", args...)
data_dir(args...) = joinpath("../data/training", args...)
results_dir(args...) = joinpath("../data/results", args...)

# Setup the data
data_paths = [
    data_dir("correct_partition.csv"),
    data_dir("over_partition.csv"),
    data_dir("under_partition.csv")
]
n_samples = 2000
n_data = 3

# Setup ICVIs
cvis = [
    # cSIL(),
    # PS()
    PS(),
    GD43()
]
n_cvis = length(cvis)

# Setup the correlation
n = 100                         # Window size
n_corr = n_samples - n + 1      # Number of correlations
corrs = zeros(n_corr, n_data)   # Preallocation of correlations
iter_corr = n:n_samples         # Iterator for the correlation generation

# Iterate over all partition types
for dx = 1:n_data
    data_path = data_paths[dx]

    # Load the training data
    data, labels = get_cvi_data(data_path)
    labels = relabel_cvi_data(labels)
    data_name = splitext(basename(data_path))[1]    # Get the data_name

    # Instantiate the icvis
    local_cvis = deepcopy(cvis)
    criterion_values = zeros(n_samples, n_cvis)

    # Run the ICVIs iteratively
    for i = 1:n_samples
        sample = data[:, i]
        label = labels[i]
        for j = 1:n_cvis
            criterion_values[i, j] = get_icvi!(local_cvis[j], sample, label)
        end
    end

    # Calculate the windowed correlations
    for i = iter_corr
        corrs[i-n+1, dx] = corspearman(criterion_values[i-n+1:i, 1], criterion_values[i-n+1:i, 2])
    end
end

p = plot(dpi=dpi, legend=:bottomleft)
for i = 1:n_data
    plot!(p, 1:n_corr, corrs[:, i], label=basename(data_paths[i]), size=(800,500))
end
# plot!(p, 1:n_samples, labels, label="Labels")
# plot!(p, 1:n_samples, data.train_y, label="True Labels")
# plot!(p, 1:n_corr, criterion_values_p, label="Porcelain")
# title!("Spearman, Data: " * basename(data_path))
title!("Partition ICVI Correlations")
xlabel!("Sample Index")
ylabel!("Spearman")
# xlims!(1, n_corr)
# ylims!(0, Inf)

# Save and show the plot
# png(p, results_dir("2_correlation_partition"))
display(p)
dest_file = results_dir("2_correlation_partition.txt")
savefig(p, dest_file)


text:
"""
    function identify_relevant_dps(ahpq::AHPQdata, query::AbstractArray)

Computes the inner products between the query and the input pre-clustering vectors.
Returns the datapoints assigned to the closest `config[:b]` centers and the corresponding distances.
"""

function identify_relevant_dps(ahpq::AHPQdata, query::AbstractArray)
    n_centers = size(ahpq.clusterdata.centers)[2]
    distances = zeros(n_centers)
    for i in 1:n_centers
        @inbounds distances[i] = dot(query, @inbounds @view ahpq.clusterdata.centers[:,i])
    end
    cluster_indexes = partialsortperm(distances, 1:ahpq.config[:b], rev=true)
    dp_to_inspect = ahpq.clusterdata.assignments[cluster_indexes]
    return dp_to_inspect, distances[cluster_indexes]
end

struct cluster_iterator
    codebook_i::Int
    n_centers::Int
    n_dims_center::Int
end

function Base.iterate(f::cluster_iterator,stepid=0)
    # state the ending conditions first
    stepid += 1
    if f.n_centers < stepid
        return nothing
    end
    return ((stepid,(f.codebook_i-1)*f.n_centers*f.n_dims_center + (stepid-1)*f.n_dims_center + 1),stepid)
end

"""
    function create_complete_lookup_table(qd::QuantizerData, query::AbstractArray)

Helper function of the main `MIPS` function.\n
Computes the inner products between all codebook centers and the query. Distances are stored in a complete lookup table and is returned to speed up searching.
"""
function create_complete_lookup_table(qd::QuantizerData, query::AbstractArray)
    lookup_table = zeros(qd.n_codebooks, qd.n_centers)
    for i in 1:qd.n_codebooks
        clust_iterator = cluster_iterator(i, qd.n_centers, qd.n_dims_center)
        @inbounds for (j, c_i) in clust_iterator
            @fastmath lookup_table[i,j] = dot(
                (@view(query[((i-1)*qd.n_dims_center+1):i*qd.n_dims_center])),
                (@view(qd.C[c_i:(c_i+qd.n_dims_center-1)])))
        end
    end
    return lookup_table
end
"""
    function compute_distances(qd::QuantizerData, dp_ids::AbstractArray, LUT::AbstractMatrix, distances::Int, norms::AbstractArray)

Computes the approximate inner products between the datapoints and the query, based on preclustering distance and the PQ distance.
Returns a vector of distances for all datapoints in `dp_ids`, using the lookup table `LUT` and preclustering distances `distances`.
"""
function compute_distances(qd::QuantizerData, dp_ids::AbstractArray, LUT::AbstractMatrix, distances::AbstractArray, norms::AbstractArray)
    distances_new = Float32[]
    dp_i = 0
    for (distance_i, bin) in enumerate(dp_ids)
        for dp in bin
            dp_i += 1
            push!(distances_new, distances[distance_i])
            j = UInt8(0)

            @inbounds for k in @inbounds @view qd.I.assignments[:,dp]
                j+= UInt8(1)
                distances_new[dp_i] += norms[dp]*LUT[j,k]
            end
        end
    end
    distances_new
end

function compute_distances(qd::QuantizerData, dp_ids::AbstractArray, LUT::AbstractMatrix, distances::Int, norms::AbstractArray)
    n_dp = length(dp_ids)
    distances = zeros(n_dp)
    for (i,ii) in enumerate(dp_ids)
        j = UInt8(0)
        @inbounds for k in @inbounds @view qd.I.assignments[:,ii]
            j+= UInt8(1)
            @inbounds distances[i] += @inbounds LUT[j,k]
        end
    end
    return distances
end

"""
    function MIPS(ahpq::AHPQdata, query::AbstractArray, k::Int)

Searches the most `k` max inner products for the input query/queries for the
given indexes and configurations `ahpq`. Takes a single query-vector or a set of queries (matrix).

# Example
```julia
traindata = rand(d, n)
queries = rand(d, m)
ahpq = builder(traindata, 0.2)

k_neares_neighbours = MIPS(ahpq, queries, 10)
k_neares_neighbours = MIPS(ahpq, queries[:,1], 10)
"""
function MIPS(ahpq::AHPQdata, query::AbstractArray, k::Int)

     # Step 1
    if ahpq.config.a > 0
        dp_ids, distances = identify_relevant_dps(ahpq, query)
    else
        dp_ids = 1:ahpq.qd.n_dp
        distances = 0
    end

    # Step 2
    LUT = create_complete_lookup_table(ahpq.qd, query)
    distances = compute_distances(ahpq.qd, dp_ids, LUT, distances, ahpq.norms)

    # Step 3 -- Using reordering settings
    if ahpq.config.reorder == 0
        return @inbounds collect(Iterators.flatten(dp_ids))[partialsortperm(distances, 1:k, rev=true)]
    else
        dp_ids_searcher = partialsortperm(distances, 1:ahpq.config[:reorder], rev=true)
        dp_real_ids = collect(Iterators.flatten(dp_ids))
        exact_distances = dropdims(query'ahpq.residuals[:,dp_real_ids[dp_ids_searcher]],dims=1)
        distances[dp_ids_searcher] .+= exact_distances
        return @inbounds @view dp_real_ids[partialsortperm(distances, 1:k, rev=true)]
    end
end

function MIPS(ahpq::AHPQdata, queries::AbstractMatrix, n_neighbors::Int)
    n_queries = size(queries)[2]
    ranking = Array{Int64,2}(undef, n_neighbors, n_queries)
    Threads.@threads for j in 1:n_queries
        ranking[:,j] = MIPS(ahpq, @inbounds(@view(queries[:,j])), n_neighbors)
    end
    return ranking
end


text:
Base.Enums.basetype(::Type{CE}) where {CE<:Cenum} = Base.Enums.basetype(bitstype(CE))

Base.Enums.namemap(ce::Cenum) = Base.Enums.namemap(typeof(ce))
function Base.Enums.namemap(::Type{CE}) where {CE<:Cenum}
	result = ()
	for (n, v) in values(bitstype(CE))
		result = (result..., (n, v))
	end
	return result
end

(::Type{CE})(val::Integer = 0) where {CE<:Cenum} = Core.Intrinsics.bitcast(bitstype(CE), convert(Base.Enums.basetype(CE), val))::CE
(::Type{T})(ce::CE) where {T<:Integer, CE<:Cenum} = convert(T, Core.Intrinsics.bitcast(Base.Enums.basetype(CE), ce))::T

Base.sizeof(::Type{CE}) where {CE<:Cenum} = invoke(sizeof, Tuple{Any}, bitstype(CE))
Base.zero(::Type{CE}) where {CE<:Cenum} = CE()

Base.convert(::Type{CE}, val::Integer) where {CE<:Cenum} = CE(val)
Base.convert(::Type{T}, ce::CE) where {CE<:Cenum, T<:Integer} = T(ce)

Base.promote_rule(::Type{T}, ::Type{CE}) where {T<:Integer, CE<:Cenum} = promote_type(T, Base.Enums.basetype(CE))
Base.:(~)(ce::CE) where {CE<:Cenum} = ~convert(Base.Enums.basetype(CE), ce)
for op in (:+, :-, :&, :|, :xor, :(==))
	@eval begin
		Base.$(op)(x::CE, y::CE) where {CE<:Cenum} = ($op)(Core.Intrinsics.bitcast(Base.Enums.basetype(CE), x), Core.Intrinsics.bitcast(Base.Enums.basetype(CE), y))
		Base.$(op)(x::CE, y::T) where {CE<:Cenum, T<:Integer} = ($op)(promote_type(T, CE)(x), promote_type(T, CE)(y))
		Base.$(op)(x::T, y::CE) where {T<:Integer, CE<:Cenum} = ($op)(promote_type(T, CE)(x), promote_type(T, CE)(y))
	end
end

function Base.typemin(::Type{CE}) where {CE<:Cenum}
	x = nothing
	for (n, v) in Base.Enums.namemap(CE)
		x = x === nothing || v < x ? v : x
	end
	return x === nothing ? CE() : CE(x)
end

function Base.typemax(::Type{CE}) where {CE<:Cenum}
	x = nothing
	for (n, v) in Base.Enums.namemap(CE)
		x = x === nothing || v >= x ? v : x
	end
	return x === nothing ? CE() : CE(x)
end

function Base.instances(::Type{CE}) where {CE<:Cenum}
	result = ()
	for (n, v) in Base.Enums.namemap(CE)
		result = (result..., CE(v))
	end
	return result
end

function Base.Symbol(ce::Cenum)
	for (n, v) in Base.Enums.namemap(ce)
		typeof(v)(ce) == v && return n::Symbol
	end
	error("Value of $(typeof(ce)) is not one of $(join(map(first, Base.Enums.namemap(ce)), ", "))")
end

function Base.show(io::IO, x::Cenum)
    sym = Symbol(x)
    if !get(io, :compact, false)
        from = get(io, :module, Main)
        def = typeof(x).name.module
        if from === nothing || !Base.isvisible(sym, def, from)
            show(io, def)
            print(io, ".")
        end
    end
    print(io, sym)
end

function Base.show(io::IO, ::MIME"text/plain", x::Cenum)
    print(io, x, "::")
    show(IOContext(io, :compact => true), typeof(x))
    print(io, " = ")
    show(io, Integer(x))
end

function Base.show(io::IO, m::MIME"text/plain", t::Type{<:Cenum})
    if isconcretetype(t)
        print(io, "Cenum ")
        Base.show_datatype(io, t)
        print(io, ":")
        for x in instances(t)
            print(io, "\n", Symbol(x), " = ")
            show(io, Integer(x))
        end
    else
        invoke(show, Tuple{IO, MIME"text/plain", Type}, io, m, t)
    end
end


text:
using Consecutive
using Test

# Test some normal iterations.
@testset "Testing Normal Iteration" begin
    y = (1,2,3,4,5)
    I = consecutive(y)
    @test length(I) == length(y) - 1
    @test Base.IteratorSize(I) == Base.IteratorSize(y)
    @test Base.IteratorEltype(I) == Base.IteratorEltype(y)

    for (index, item) in enumerate(I)
        @test item == (y[index], y[index + 1])
    end
    @test collect(consecutive(y)) == [(1,2), (2,3), (3,4), (4,5)]

    # Test "restrict"
    @test Consecutive.restrict(Base.HasShape{4}()) == Base.HasLength()
    @test Consecutive.restrict(Base.HasLength()) == Base.HasLength()
    @test Consecutive.restrict(Base.SizeUnknown()) == Base.SizeUnknown()
end

@testset "Testing Vectors" begin
    y = [[1 2];[3 4]]
    I = consecutive(y)

    @test length(I) == length(y) - 1
    @test eltype(I) == Tuple{Int,Int}
    @test Base.IteratorSize(I) == Base.HasLength()

    expected = ((1,3), (3,2), (2,4))
    for (a, b) in zip(I, expected)
        @test a == b
    end

    # Compare against a base function
    A = rand(10000)
    mydiff(x) = [b-a for (a,b) in consecutive(x)]
    @test diff(A) == mydiff(A)
end

@testset "Testing IO Iteration" begin

    parts = [
        "The quick",
        "brown fox jumped over",
        "the lazy",
        "cat"
    ]
    str = join(parts, "\n")

    iobuf = IOBuffer(str)

    for (index, (lna, lnb)) in enumerate(consecutive(eachline(iobuf)))
        @test lna == parts[index]
        @test lnb == parts[index + 1]
    end
end


text:
@symbol_func function omega_pe(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_n = n_rho(cur_reactor, cur_rho)

  cur_omega = 5.64e11

  cur_omega *= sqrt(cur_n)

  cur_omega
end


text:
## some unit tests for internal functions ##

using Base.Test

## translate()
@test MCMC.translate(12)                  == 12
@test MCMC.translate(:(b = a+6))          == :(b = a+6)
@test MCMC.translate(:(copy!(b,z)))       == :(copy!(b,z))
@test MCMC.translate(:(z ~ Normal(1,2)))  == :(__acc += logpdf(Normal(1,2),z))
@test MCMC.translate([:(b = Array(Float64,10,10)), :(z ~ Normal(1,2))]) ==
	[:(b = Array(Float64,10,10)), :(__acc += logpdf(Normal(1,2),z))]

## vec2var()
@test MCMC.vec2var(x=3.)                  == [:(x=__beta[1])]
@test MCMC.vec2var(z=[1., 2, 3])          == [:(z=__beta[1:3])]
@test MCMC.vec2var(A=[12. 45 0; 18 1 2])  == [:(A=reshape(__beta[1:6],2,3))]
@test MCMC.vec2var(A=[12. 45 0; 18 1 2], x=3., z=[1., 2, 3]) ==
	[:(A=reshape(__beta[1:6],2,3)), :(x=__beta[7]), :(z=__beta[8:10])]

## var2vec()
dp(v) = MCMC.ReverseDiffSource.dprefix(v)
@test MCMC.var2vec(x=3.)                  == :([$(dp(:x))])
@test MCMC.var2vec(z=[1., 2, 3])          == :([$(dp(:z))])
@test MCMC.var2vec(A=[12. 45 0; 18 1 2])  == :([vec($(dp(:A)))])
@test MCMC.var2vec(z=[1., 2, 3], A=[12. 45 0; 18 1 2], x=3.) == :([$(dp(:z)), vec($(dp(:A))), $(dp(:x))])

## modelVars()
@test MCMC.modelVars(x=3.)                    == (1, [:x => (1,())]      , [3.])
@test MCMC.modelVars(x=3)                     == (1, [:x => (1,())]      , [3.])
@test MCMC.modelVars(z=[1., 2, 3])            == (3, [:z => (1,(3,))]    , [1.,2,3])
@test MCMC.modelVars(z=[true, false, false])  == (3, [:z => (1,(3,))]    , [1.,0,0])
@test MCMC.modelVars(A=[12. 45 0; 18 1 2])    == (6, [:A => (1,(2,3))]   , [12.,18,45,1,0,2])
@test MCMC.modelVars(z=[1., 2, 3], x=3., A=[12. 45 0; 18 1 2]) ==
	(10, [:A => (5,(2,3)), :z => (1,(3,)), :x => (4,())]   , [1.,2,3,3,12,18,45,1,0,2])


text:
include("../roblib.jl")

function mark_perimeter(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))

    side = Up
    for i in 0:1:3
        putmarkers!(r, side)
        side = right(side)
    end

    movements!(r, (Up, Right), reverse(num_steps))
end


text:
module ASCIIPlots
    using Printf, Statistics

    export ASCIIPlot
    export imagesc, lineplot, scatterplot

    struct ASCIIPlot
        s::String
    end

    Base.show(io::IO, p::ASCIIPlot) = print(io, p.s)

    include("imagesc.jl")
    include("lineplot.jl")
    include("scatterplot.jl")
end