text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a1e77490-ba51-11eb-3cd8-1d82e8b5f7e4
begin
	using Pkg
	Pkg.activate(joinpath(@__DIR__, ".."))
end

# ╔═╡ 9bb18432-873d-4979-9e2c-a71e1007b2d7
using Agents, InteractiveDynamics, CairoMakie, PlutoUI

# ╔═╡ f3d42d7d-3e92-4f78-9284-689e2ac50017
using Setfield

# ╔═╡ 6ca2d3d4-d9ac-4e9d-9b1b-e8a66f00c7e6
using Random

# ╔═╡ d4ed1cd4-89a1-49f3-8c18-19215b4c2e8e
using DataFrames

# ╔═╡ 89d0deb2-13c1-45ca-b2ec-ee9704ae6e15
using Base.Threads

# ╔═╡ 3097bf32-e7ba-4b36-8674-3ddd37e1f395
mutable struct BPH <: AbstractAgent
	id::Int
	pos::Dims{2}
	energy
	age
	nb_reproduce
end

# ╔═╡ d073c77d-390c-437b-8ede-630f07e5d9d9
# module MyModel
# using Ag
# @agent BPH GridAgent{2} begin
# 	energy
# 	age
# 	nb_reproduce
# end
# end

# ╔═╡ 4abaea9e-d7ff-4d00-9433-0dd6e3d6da8d
# BPH = MyModel.BPH

# ╔═╡ d948d629-73bc-48b2-9275-51c0bd82a29a
function calculate_positions(pos, nb_flower, field_size)
	positions = if pos === :corner
		Iterators.product(1:5, 1:5)
	elseif pos === :border
		_, field_height = field_size
		Iterators.product(1:2, 1:field_height)
	elseif pos === :random
		field_width, field_height = field_size
		flower_start = (field_width ÷ 2 - nb_flower ÷ 2)
		Iterators.product(1:flower_start, 1:field_height)
	end
	collect(positions)
end

# ╔═╡ 60c02b02-92d3-4c58-b9f1-6e71dfaa28fc
function init_model(nb_bph_init, nb_flower, init_position, field_size;
		energy_miss = 0.025,
		age_init = 168,
		age_reproduce = 504,
		age_old = 600,
		age_die = 720,
		pr_move = 0.99,
		pr_reproduce = 0.07,
		pr_egg_death = 0.0025,
		pr_old_death = 0.04,
		offspring_max = 12,
		offspring_min = 5,
		energy_max = 1.0,
		energy_transfer = 0.1,
		energy_consume = 0.025,
		energy_move = 0.2,
		energy_reproduce = 0.8,
		seed = 1)
	field_width, field_height = field_size
	properties = (
		age_init = 168,
		age_reproduce = 504,
		age_old = 600,
		age_die = 720,
		pr_move = 0.99,
		pr_reproduce = 0.7,
		pr_egg_death = 0.0025,
		pr_old_death = 0.04,
		offspring_max = 12,
		offspring_min = 5,
		energy_miss = 0.025,
		energy_max = 1.0,
		energy_transfer = 0.1,
		energy_consume = 0.025,
		energy_move = 0.2,
		energy_reproduce = 0.8,
		flower_start = field_width ÷ 2 - nb_flower ÷ 2,
		flower_end = field_width ÷ 2 + nb_flower ÷ 2,
		food = fill(0.4, field_size)
	)
	field = GridSpace(field_size; periodic=false)
	rng = MersenneTwister(seed)
	scheduler = Schedulers.randomly
	model = ABM(BPH, field; properties, rng, scheduler)

	init_positions = calculate_positions(init_position, nb_flower, field_size)
	for _ in 1:nb_bph_init
		id = nextid(model)
		pos = rand(model.rng, init_positions)
		energy = rand(model.rng, 0.4:0.01:0.6)
		bph = BPH(id, pos, energy, rand(model.rng, 0:299), rand(model.rng, 1:21))
		add_agent_pos!(bph, model)
	end

	return model
end

# ╔═╡ cca251da-24fe-4b97-87db-337cc2d9490d
function bph_reflex!(bph, model)
	# older
	bph.age = bph.age + 1

	# move around
	directions = Iterators.product([-1, -2, 1, 2], [-1, -2, 1, 2])
	if bph.age > model.age_init && rand(model.rng) < model.pr_move
		walk!(bph, rand(model.rng, directions |> collect), model)
	end

	# loss energy normally
	if bph.age ≥ model.age_init
		bph.energy = bph.energy - model.energy_consume
	end

	# loss energy in flower
	x, y = bph.pos
	if model.flower_start ≤ x < model.flower_end
		bph.energy -= model.energy_miss
	elseif model.food[x, y] > 0 && bph.age ≥ model.age_init
		transfer = min(model.energy_transfer, model.food[x, y])
		model.food[x, y] -= transfer
		bph.energy = min(bph.energy + transfer, model.energy_max)
	end

	# reproduce
	if (bph.age ≥ model.age_reproduce &&
		bph.energy ≥ model.energy_reproduce &&
		bph.nb_reproduce < 21 &&
		rand(model.rng) ≤ model.pr_reproduce)
		#
		bph.nb_reproduce = bph.nb_reproduce + 1
		nb_offspring = rand(model.rng, model.offspring_min:model.offspring_max)
		for _ in 1:nb_offspring
			id = nextid(model)
			bph = BPH(id, bph.pos, 0.4, 0, 0)
			add_agent_pos!(bph, model)
		end
	end

	# die
	if bph.energy ≤ 0 ||
			(bph.age ≥ model.age_die) ||
			(model.age_die ≥ bph.age ≥ model.age_old && rand(model.rng) ≤ model.pr_old_death)
			(bph.age < model.age_init && rand(model.rng) ≤ model.pr_egg_death)
		kill_agent!(bph, model)
	end
end

# ╔═╡ af4042f1-7679-4cc5-ad20-b23937ffa40b
function rice_grow!(model)
	alive = @. model.food > 0
	@. model.food = min(model.food + model.food * 0.008 * alive , 1.0)
end

# ╔═╡ e8c37ce6-9376-4215-ae47-a3e3df26276a
bph_color(bph) = if bph.age < 168
	RGBf0(0, 0, 255)
else
	RGBf0(255, 0, 0)
end

# ╔═╡ 7e042ab7-6d6e-4df0-a9ef-b521fdf9f1de
rice_color(model) = map(positions(model)) do pos
	x, y = pos
	if model.flower_start ≤ x ≤ model.flower_end
		# RGBf0(255, 255, 0)
		-1
	else
		# RGBAf0(0, 255, 0, model.food[x, y])
		max(model.food[x, y], 0)
	end
end

# ╔═╡ 9ff14f09-e64b-4898-935a-c425ef1d3eaf
plotkwargs = (ac = bph_color, heatarray = rice_color)

# ╔═╡ 456e997d-c424-4647-9ba9-cd98a4c43ed2
heatkwargs = (nan_color = RGBf0(1.0, 1.0, 0),
	colormap = [RGBAf0(0, 1.0, 0, i) for i in 0:0.01:1]
)

# ╔═╡ 6f4dc104-c0f1-4c9f-af54-5de98a6e62e9
md"""
Generate Video $(@bind should_video PlutoUI.CheckBox())

Collect Data $(@bind should_data PlutoUI.CheckBox())

Flower $(@bind flower PlutoUI.Slider(12:20; show_value=true))
"""

# ╔═╡ 2da972ea-1e4b-4096-81b2-1e86b4eb4a10
@bind seed PlutoUI.NumberField(1:2000)

# ╔═╡ 48fb77fb-ab35-4a4b-948b-ce3599f7ec60
if should_video
	vid = "bph.mp4"
	let model = init_model(200, flower, :corner, (100, 100); seed = seed)
		abm_video(vid, model, bph_reflex!, rice_grow!; frames=1200 - 1, framerate=10, heatkwargs, plotkwargs...)
		LocalResource(vid)
	end
end

# ╔═╡ 149636d6-95b3-4ead-a607-f37c4de5ed52
1234

# ╔═╡ f0ddf1cc-12e5-4444-a5f4-d4085b8056b3
# @time max_times = map(1:100) do seed
# 	bph(x) = true
# 	model = init_model(200, 16, :corner, (100, 100); seed = seed)
# 	adf, _ = run!(model, bph_reflex!, rice_grow!, 2880; adata = [(bph, count)])
# 	seed => maximum(adf.step)
# end

# ╔═╡ 9e235da3-3f0f-42de-af15-10d940b3304c
@bind rerun_btn PlutoUI.Button("Re-run")

# ╔═╡ 3492f6a0-a867-4403-b32f-4dff1f706142
begin
	rerun_btn
	model = init_model(200, flower, :corner, (100, 100); seed)
end

# ╔═╡ 063e3c58-fb8f-46d7-8f37-677af66598bf
md"## Run & Collect data"

# ╔═╡ 09cf52d0-d6f5-4dad-ad43-68aea0f44cdd
# adata = [x -> x.pos, x -> x.age]

# ╔═╡ 736ae5e6-3c3f-4e24-9b93-6fc578d378db
adata = let bph_A1(bph) = bph.age ≥ model.age_init && bph.pos[1] < model.flower_start
	bph_E1(bph) = bph.age < model.age_init && bph.pos[1] < model.flower_start
	bph_A2(bph) = bph.age ≥ model.age_init && bph.pos[1] ≥ model.flower_end
	bph_E2(bph) = bph.age < model.age_init && bph.pos[1] ≥ model.flower_end
	[(bph_E1, count), (bph_E2, count), (bph_A1, count), (bph_A2, count)]
end

# ╔═╡ ff33a882-4092-42e5-a272-043cc9c2dbf0
# max_step = maximum(adfa.step)

# ╔═╡ bc177477-8bb5-4821-9076-517b07c0c0f7
# GC.gc()

# ╔═╡ 6835bd93-4824-4768-a575-912fdc4a94ee
let model = init_model(200, flower, :random, (100, 100))
	GC.gc()
	PlutoUI.with_terminal() do
		@time run!(model, bph_reflex!, rice_grow!, 2880; adata=adata);
	end
end

# ╔═╡ 8a7027e6-7959-43e9-b0b6-245ffbeeec74
adf, mdf = run!(model, bph_reflex!, rice_grow!, 2880; adata=adata);

# ╔═╡ 85e8ab56-3345-4e31-a6f0-cff58c42ac6c
# adf2 = rename(adf, [:step, :id, :pos, :age]);

# ╔═╡ 1471bab5-baab-48f3-8f13-975b2aabcd7b
adf

# ╔═╡ 4a9553e5-9229-47ec-8248-e72be094fa34

# ╔═╡ 088a6bab-0be7-4f0e-96bf-7497b6dd7209
let f = Figure()
	Axis(f[1, 1])
	step = adf.step[adf.count_bph_A1 .> 0]
	bph1 = adf.count_bph_A1[adf.count_bph_A1 .> 0]
	lines!(step, bph1)
	lines!(adf.step ./ 24, adf.count_bph_A1)
	lines!(adf.step ./ 24, adf.count_bph_A2)
	f
end

# ╔═╡ 96336b96-44d5-4356-ab9d-ab8cbd9ea1bb
let f = Figure()
	Axis(f[1, 1])
	@threads for seed = 1:5
		model = init_model(20, 18, :corner, (100, 100); seed = seed);
		adf, mdf = run!(model, bph_reflex!, rice_grow!, 2880; adata=adata);
		bph1 = adf.count_bph_A1[1:1700]
		steps = adf.step[1:1700]
		lines!(steps, bph1)
		# lines!(adf.step ./ 24, adf.count_bph_A1)
		# lines!(adf.step ./ 24, adf.count_bph_A2)
	end
	f
end

# ╔═╡ 67366c81-42fd-407c-bf21-1fe04cc4dda0
# adfa = let adult(age) = age .≥ model.age_init
# 	egg(age) = age .< model.age_init
# 	field1(pos) = getindex.(pos, 1) .< model.flower_start
# 	field2(pos) = getindex.(pos, 1) .≥ model.flower_end
# 	combine(groupby(adf2, :step),
# 		[:pos, :age] => ((p, a) -> count(egg(a) .& field1(p))) => :e1,
# 		[:pos, :age] => ((p, a) -> count(adult(a) .& field1(p))) => :a1,
# 		[:pos, :age] => ((p, a) -> count(egg(a) .& field2(p))) => :e2,
# 		[:pos, :age] => ((p, a) -> count(adult(a) .& field2(p))) => :a2,
# 	)
# end;

# ╔═╡ 5a8a235a-64ff-41a3-981b-684a1424fa07
# let f = Figure()
# 	Axis(f[1, 1])
# 	# lines!(adfa.step / 24, adfa.e1)
# 	# lines!(adfa.step / 24, adfa.e2)
# 	lines!(adfa.step, adfa.a1)
# 	lines!(adfa.step, adfa.a2)
# 	f
# end

# ╔═╡ d6502588-0196-4953-86f0-5244b519f9ee

# ╔═╡ Cell order:
# ╠═a1e77490-ba51-11eb-3cd8-1d82e8b5f7e4
# ╠═9bb18432-873d-4979-9e2c-a71e1007b2d7
# ╠═f3d42d7d-3e92-4f78-9284-689e2ac50017
# ╠═3097bf32-e7ba-4b36-8674-3ddd37e1f395
# ╠═d073c77d-390c-437b-8ede-630f07e5d9d9
# ╠═4abaea9e-d7ff-4d00-9433-0dd6e3d6da8d
# ╠═d948d629-73bc-48b2-9275-51c0bd82a29a
# ╠═60c02b02-92d3-4c58-b9f1-6e71dfaa28fc
# ╠═6ca2d3d4-d9ac-4e9d-9b1b-e8a66f00c7e6
# ╠═3492f6a0-a867-4403-b32f-4dff1f706142
# ╠═cca251da-24fe-4b97-87db-337cc2d9490d
# ╠═af4042f1-7679-4cc5-ad20-b23937ffa40b
# ╠═e8c37ce6-9376-4215-ae47-a3e3df26276a
# ╠═7e042ab7-6d6e-4df0-a9ef-b521fdf9f1de
# ╠═9ff14f09-e64b-4898-935a-c425ef1d3eaf
# ╠═456e997d-c424-4647-9ba9-cd98a4c43ed2
# ╠═48fb77fb-ab35-4a4b-948b-ce3599f7ec60
# ╠═6f4dc104-c0f1-4c9f-af54-5de98a6e62e9
# ╠═2da972ea-1e4b-4096-81b2-1e86b4eb4a10
# ╠═149636d6-95b3-4ead-a607-f37c4de5ed52
# ╠═f0ddf1cc-12e5-4444-a5f4-d4085b8056b3
# ╟─9e235da3-3f0f-42de-af15-10d940b3304c
# ╟─063e3c58-fb8f-46d7-8f37-677af66598bf
# ╠═09cf52d0-d6f5-4dad-ad43-68aea0f44cdd
# ╠═736ae5e6-3c3f-4e24-9b93-6fc578d378db
# ╠═ff33a882-4092-42e5-a272-043cc9c2dbf0
# ╠═bc177477-8bb5-4821-9076-517b07c0c0f7
# ╠═6835bd93-4824-4768-a575-912fdc4a94ee
# ╠═8a7027e6-7959-43e9-b0b6-245ffbeeec74
# ╠═d4ed1cd4-89a1-49f3-8c18-19215b4c2e8e
# ╠═85e8ab56-3345-4e31-a6f0-cff58c42ac6c
# ╠═1471bab5-baab-48f3-8f13-975b2aabcd7b
# ╠═4a9553e5-9229-47ec-8248-e72be094fa34
# ╠═088a6bab-0be7-4f0e-96bf-7497b6dd7209
# ╠═89d0deb2-13c1-45ca-b2ec-ee9704ae6e15
# ╠═96336b96-44d5-4356-ab9d-ab8cbd9ea1bb
# ╠═67366c81-42fd-407c-bf21-1fe04cc4dda0
# ╠═5a8a235a-64ff-41a3-981b-684a1424fa07
# ╠═d6502588-0196-4953-86f0-5244b519f9ee

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
function iterative_probabilistic_improvement(tuning_run::Run,
                                             channel::RemoteChannel;
                                             threshold::AbstractFloat = 2.)
    function iterate(tuning_run::Run)
        return probabilistic_improvement(tuning_run,
                                         threshold = threshold)
    end

    technique(tuning_run,
              channel,
              iterate,
              name = "Iterative Probabilistic Improvement")
end

ext: jl


text:
using Documenter
using CarrierCapture

makedocs(
    sitename = "CarrierCapture",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [CarrierCapture],
    pages = [
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Library" => Any[
            "Public" => "lib/public.md",
            "Brooglie" => "lib/brooglie.md",
            "Plotter" => "lib/plotter.md"
        ],
    ]
)

deploydocs(
    target = "build",
    repo   = "github.com/WMD-group/CarrierCapture.jl.git",
)


text:
using TextAnalysis
function loadYearsAndTitles(fileName::String)::Vector{Tuple{Int,String}}
    yearTitles = Vector{Tuple{Int,String}}()
    file = open(fileName)
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line, ". ")
        @assert length(fields) == 2
        ymd = split(fields[1], "-")
        year = parse(Int, ymd[1])
        title = fields[2]
        push!(yearTitles, (year, title))
    end
    yearTitles
end

function unhyphenAndStem(stemmer::Stemmer, toks::Vector{String})::Vector{String}
    result = Vector{String}()
    for tok in toks
        if tok == ""
            continue
        end
        try
            tok = replace(tok, "/" => "-or-")
            while startswith(tok, "-")
                tok = tok[2:end]
            end
            while endswith(tok, "-")
                tok = tok[1:end - 1]
            end
            subToks = split(tok, "-")
            stemmedSubToks = stem(stemmer, subToks)
            for subTok in stemmedSubToks
                push!(result, subTok)
            end
        catch
            push!(result, tok)
        end
    end
    result
end

YearTokens = Tuple{Int,String,Vector{String}}
function tokenizeYearTitles(yearTitles::Vector{Tuple{Int,String}})::Vector{YearTokens}
    numCPUs = length(Sys.cpu_info())
    jobs = Channel{Tuple{Int,Int}}(numCPUs)
    Segments = Vector{Tuple{Int,Int}}
    SegmentedYearTokens = Tuple{Segments,Vector{Vector{YearTokens}}}
    jobOutputs = Channel{SegmentedYearTokens}(numCPUs)

    function runJob(iCPU::Int)
        stemmer = Stemmer("english")
        mySegments = Segments()
        myYearTokens = Vector{Vector{YearTokens}}()
        for job in jobs
            push!(mySegments, job)
            jobYearTokens = Vector{YearTokens}()
            for yearTitle in yearTitles[job[1]:job[2]]
                year = yearTitle[1]
                title = replace(lowercase(yearTitle[2]), "—" => "-")
                toksOfTitle = tokenize(title)
                toksOfTitle = unhyphenAndStem(stemmer, toksOfTitle)
                push!(jobYearTokens, (year, title, toksOfTitle))
            end
            push!(myYearTokens, jobYearTokens)
        end
        put!(jobOutputs, (mySegments, myYearTokens))
    end

    function makeJobs(batchSize::Int)
        n = length(yearTitles)
        for i in 1:batchSize:n
            put!(jobs, (i, min(i + batchSize - 1, n)))
        end
    end
    bind(jobs, @async makeJobs(1000))
    for iCPU in 1:numCPUs
        Threads.@spawn runJob(iCPU)
    end

    yearTitleTokens = Vector{YearTokens}(undef, length(yearTitles))
    for iCPU in 1:numCPUs
        segmentedYearTokens = take!(jobOutputs)
        segments = segmentedYearTokens[1]
        for (segID, segment) in enumerate(segments)
            yearTitleTokens[segment[1]:segment[2]] = segmentedYearTokens[2][segID]
        end
    end
    yearTitleTokens
end

function saveYearTitleTokens(fileName::String, yearTitleTokens::Vector{YearTokens})
    file = open(fileName, "w")
    for yearTokens in yearTitleTokens
        year = yearTokens[1]
        title = yearTokens[2]
        tokens = yearTokens[3]
        print(file, "$year. $title.")
        for token in tokens
            print(file, " $token")
        end
        print(file, "\n")
    end
    close(file)
end

for i in 0:106
    println("processing title block $i ...")
    yearTitles = loadYearsAndTitles("SCC2021/MAGTitles/MAGEnTitles.$i.txt")
    yearTitleTokens = tokenizeYearTitles(yearTitles)
    saveYearTitleTokens("TokenizedMAGTitles/MAGTitles.$i.txt", yearTitleTokens)
end

max_stars_repo_path: src/00 tokenize_mag_titles.jl

max_issues_repo_path: src/00 tokenize_mag_titles.jl

max_forks_repo_path: src/00 tokenize_mag_titles.jl


text:
abstract type CumulantSchedule end
function update! end
function get_cumulant(cs::CumulantSchedule, goal) end
function get_cumulant_eval_values end

#cumulants for TMaze
include("environments/tmaze_cumulants.jl")

# TMaze environments
include("environments/tabular_tmaze.jl")
include("environments/1d-tmaze.jl")

# Cont 2d Worlds environment.
include("environments/cont_2d_worlds_spec.jl")

# Mountain Car
include("environments/mountain_car.jl")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function initializemdp(version)

    app = initializeapp(version);
    config = initializeconfig(app);

    return app, config;

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using DataDeps

register(DataDep(
  "immigrant-salaries",
  "http://www.stat.ufl.edu/~winner/data/immwork.txt",
  "http://www.stat.ufl.edu/~winner/data/immwork.dat",
  "d4e517a5725a613bf30b224b53a5b3b4509cf7a126f7813a12b3d2769de6e470",
  post_fetch_method=(path -> begin
    # IDs have spaces in them which is the delimiter
    lines = readlines(path)
    map!(x -> lstrip(x[30:end]), lines, lines)
    write(path, join(lines, '\n'))

    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:4,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module CUFFT
using CUDArt
using Compat

import Base: convert
import CUDArt: destroy
export Plan, compatibility, plan, tie, RCpair, RCfree
# "Public" but not exported: version()

include("libcufft.jl")
const lib = LibCUFFT

type Plan{From,To,N}
    p::Cint
end
convert(::Type{lib.cufftHandle}, p::Plan) = p.p

function destroy(p::Plan)
#     println("Destroying ", p)
#     Base.show_backtrace(STDOUT, backtrace())
    lib.cufftDestroy(p)
end

plan_dict = @compat Dict(
    (Float32,Complex64) => lib.CUFFT_R2C,
    (Complex64,Float32) => lib.CUFFT_C2R,
    (Complex64,Complex64) => lib.CUFFT_C2C,
    (Float64,Complex128) => lib.CUFFT_D2Z,
    (Complex128,Float64) => lib.CUFFT_Z2D,
    (Complex128,Complex128) => lib.CUFFT_Z2Z
)

# For in-place R2C and C2R transforms
function RCpair{T<:AbstractFloat}(realtype::Type{T}, realsize) # TODO?: add dims
    csize = [realsize...]
    csize[1] = div(realsize[1],2) + 1
    C = CudaPitchedArray(Complex{T}, csize...)
    R = reinterpret(T, C, tuple(realsize...))
    R, C
end

function RCpair{T<:AbstractFloat}(A::Array{T}) # TODO?: add dims
    R, C = RCpair(eltype(A), size(A))
    copy!(R, A)
    R, C
end

RCfree{T<:AbstractFloat}(R::CudaPitchedArray{T}, C::CudaPitchedArray{Complex{T}}) = free(C)

function plan_size(dest, src)
    nd = ndims(dest)
    ndims(src) == nd || throw(DimensionMismatch())
    sz = Array(Cint, nd)
    for i = 2:nd
        if size(dest,i) != size(src,i)
            throw(DimensionMismatch())
        end
        sz[i] = size(dest,i)
    end
    local plantype
    if eltype(dest) == eltype(src)
        size(dest,1) == size(src,1) || throw(DimensionMismatch())
        sz[1] = size(dest,1)
    else
        szbig = max(size(dest,1), size(src,1))
        szsmall = min(size(dest,1), size(src,1))
        szsmall == div(szbig,2)+1 || throw(DimensionMismatch())
        sz[1] = szbig
    end
    reverse(sz)
end

# Returns a function dofft!(dest, src, forward).
# Unlike FFTW's plan, this does not destroy the inputs
# TODO?: add dims
function plan(dest::AbstractCudaArray, src::AbstractCudaArray; compat::Symbol = :native, stream=null_stream)
    p = Cint[0]
    sz = plan_size(dest, src)
    inembed = reverse(Cint[size(src)...])
    onembed = reverse(Cint[size(dest)...])
    inembed[end] = pitchel(src)
    onembed[end] = pitchel(dest)
    plantype = plan_dict[(eltype(src),eltype(dest))]
    lib.cufftPlanMany(p, ndims(dest), sz, inembed, 1, 1, onembed, 1, 1, plantype, 1)
    pl = Plan{eltype(src),eltype(dest),ndims(dest)}(p[1])
    compatibility(pl, compat)
    tie(pl, stream)
    (dest,src,forward) -> exec!(pl, src, dest, forward)
end

direction(forward::Bool) = forward ? lib.CUFFT_FORWARD : lib.CUFFT_INVERSE

exec!(p::Plan{Complex64,Complex64}, input, output, forward::Bool = true) =
    lib.cufftExecC2C(p, input, output, direction(forward))
exec!(p::Plan{Float32,Complex64}, input, output, forward::Bool = true) =
    lib.cufftExecR2C(p, input, output)
exec!(p::Plan{Complex64,Float32}, input, output, forward::Bool = true) =
    lib.cufftExecC2R(p, input, output)
exec!(p::Plan{Complex128,Complex128}, input, output, forward::Bool = true) =
    lib.cufftExecZ2Z(p, input, output, direction(forward))
exec!(p::Plan{Float64,Complex128}, input, output, forward::Bool = true) =
    lib.cufftExecD2Z(p, input, output)
exec!(p::Plan{Complex128,Float64}, input, output, forward::Bool = true) =
    lib.cufftExecZ2D(p, input, output)

version() = (ret = Cint[0]; lib.cufftGetVersion(ret); ret[1])

modedict = @compat Dict(
    :native     => lib.CUFFT_COMPATIBILITY_NATIVE,
    :padding    => lib.CUFFT_COMPATIBILITY_FFTW_PADDING,
    :asymmetric => lib.CUFFT_COMPATIBILITY_FFTW_ASYMMETRIC,
    :all        => lib.CUFFT_COMPATIBILITY_FFTW_ALL
)

compatibility(p::Plan, mode::Symbol) = lib.cufftSetCompatibilityMode(p, modedict[mode])

tie(p::Plan, s::Stream) = lib.cufftSetStream(p, s)

end


text:
module SIAMFANLEquations
using LinearAlgebra
using LinearAlgebra.BLAS
using SparseArrays
using SuiteSparse
using BandedMatrices
using Printf

export nsolsc
export ptcsolsc
export ptcsol
export ptcsoli
export nsol
export nsoli
export aasol
export nofact
export secant
export armijosc
export kl_gmres
export kl_bicgstab
export Katv
export kstore
export Orthogonalize!
export EvalF!
export solhistinit

include("Tools/armijo.jl")
include("Tools/PrintError.jl")
include("Tools/FunctionJacobianEvals.jl")
include("Tools/ManageStats.jl")
include("Tools/IterationInit.jl")
include("Tools/ErrorTest.jl")
include("Tools/NewtonKrylov_Tools.jl")
include("Tools/PTCTools.jl")
include("Tools/PTCToolsi.jl")
include("Tools/AA_Tools.jl")
include("Solvers/Chapter1/nsolsc.jl")
include("Solvers/Chapter1/ptcsolsc.jl")
include("Solvers/Chapter1/secant.jl")
include("Solvers/ptcsol.jl")
include("Solvers/ptcsoli.jl")
include("Solvers/nsol.jl")
include("Solvers/nsoli.jl")
include("Solvers/aasol.jl")
include("Solvers/aa_qr_update.jl")
include("Solvers/LinearSolvers/kl_gmres.jl")
include("Solvers/LinearSolvers/kl_bicgstab.jl")
include("Solvers/LinearSolvers/Orthogonalize!.jl")
include("PlotsTables/printhist.jl")

module TestProblems
using LinearAlgebra
using SparseArrays
using SuiteSparse
using BandedMatrices
using AbstractFFTs
using FFTW
using Printf

export
    #Functions
    fcos,
    fpatan,
    spitchfork,
    linatan,
    sptestp,
    sptest,
    ftanx,
    ftanxp,
    heqinit,
    heqf!,
    heqJ!,
    HeqFix!,
    simple!,
    jsimple!,
    JVsimple,
    heqbos!,
    setc!,
    chandprint,
    bvpinit,
    Fbvp!,
    Jbvp!,
    FBeam!,
    FBeamtd!,
    BeamJ!,
    BeamtdJ!,
    beaminit,
    ptctest,
    pdeF!,
    pdeJ!,
    Jvec2d,
    pdeinit,
    pdegminit,
    fishinit,
    fish2d,
    sintv,
    isintv,
    Pfish2d,
    Pvec2d,
    Lap2d,
    Dx2d,
    Dy2d,
    solexact,
    l2dexact,
    dxexact,
    dyexact,
    hardleft!,
    hardleftFix!

include("TestProblems/Scalars/fcos.jl")
include("TestProblems/Scalars/fpatan.jl")
include("TestProblems/Scalars/spitchfork.jl")
include("TestProblems/Scalars/linatan.jl")
include("TestProblems/Scalars/ftanx.jl")
include("TestProblems/Systems/simple!.jl")
include("TestProblems/Systems/Fbvp!.jl")
include("TestProblems/Systems/FBeam!.jl")
include("TestProblems/Systems/Hequation.jl")
include("TestProblems/Systems/EllipticPDE.jl")
include("TestProblems/Systems/PDE_Tools.jl")
end

module Examples
using SIAMFANLEquations
using SIAMFANLEquations.TestProblems
using LinearAlgebra
using BandedMatrices

export ptciBeam
export ptcBeam
export ivpBeam
export BVP_solve
export nsolheq
export NsolPDE
export NsoliPDE
export PDE_aa

include("Examples/ptciBeam.jl")
include("Examples/ptcBeam.jl")
include("Examples/ivpBeam.jl")
include("Examples/BVP_solve.jl")
include("Examples/NsolPDE.jl")
include("Examples/NsoliPDE.jl")
include("Examples/PDE_aa.jl")
include("Examples/Internal/nsolheq.jl")
end

end # module


text:
#---------------------------------------------------------
# # [SPECTrecon ML-EM](@id 4-mlem)
#---------------------------------------------------------

# This page illustrates ML-EM reconstruction with the Julia package
# [`SPECTrecon`](https://github.com/JeffFessler/SPECTrecon.jl).

# ### Setup

# Packages needed here.

using SPECTrecon
using MIRTjim: jim, prompt
using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# Maximum-likelihood expectation-maximization (ML-EM)
# is a classic algorithm for performing SPECT image reconstruction.

# ### Simulation data

nx,ny,nz = 64,64,50
T = Float32
xtrue = zeros(T, nx,ny,nz)
xtrue[(1nx÷4):(2nx÷3), 1ny÷5:(3ny÷5), 2nz÷6:(3nz÷6)] .= 1
xtrue[(2nx÷5):(3nx÷5), 1ny÷5:(2ny÷5), 4nz÷6:(5nz÷6)] .= 2

average(x) = sum(x) / length(x)
function mid3(x::AbstractArray{T,3}) where {T}
    (nx,ny,nz) = size(x)
    xy = x[:,:,ceil(Int, nz÷2)]
    xz = x[:,ceil(Int,end/2),:]
    zy = x[ceil(Int, nx÷2),:,:]'
    return [xy xz; zy fill(average(xy), nz, nz)]
end
jim(mid3(xtrue), "Middle slices of xtrue")

# ### PSF

# Create a synthetic depth-dependent PSF for a single view
px = 11
psf1 = psf_gauss( ; ny, px)
jim(psf1, "PSF for each of $ny planes")

# In general the PSF can vary from view to view
# due to non-circular detector orbits.
# For simplicity, here we illustrate the case
# where the PSF is the same for every view.

nview = 60
psfs = repeat(psf1, 1, 1, 1, nview)
size(psfs)

# ### SPECT system model using `LinearMapAA`

dy = 8 # transaxial pixel size in mm
mumap = zeros(T, size(xtrue)) # zero μ-map just for illustration here
plan = SPECTplan(mumap, psfs, dy; T)

using LinearMapsAA: LinearMapAA
using LinearAlgebra: mul!
forw! = (y,x) -> project!(y, x, plan)
back! = (x,y) -> backproject!(x, y, plan)
idim = (nx,ny,nz)
odim = (nx,nz,nview)
A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)

# ### Basic Expectation-Maximization (EM) algorithm

# Noisy data
using Distributions: Poisson

if !@isdefined(ynoisy) # generate (scaled) Poisson data
    ytrue = A * xtrue
    target_mean = 20 # aim for mean of 20 counts per ray
    scale = target_mean / average(ytrue)
    scatter_fraction = 0.1 # 10% uniform scatter for illustration
    scatter_mean = scatter_fraction * average(ytrue) # uniform for simplicity
    ynoisy = rand.(Poisson.(scale * (ytrue .+ scatter_mean))) / scale
end
jim(ynoisy, "$nview noisy projection views")

# ### ML-EM algorithm - basic version

# This basic ML-EM version uses the linear map, but it is still allocating.

function mlem(x0, ynoisy, background, A; niter::Int = 20)
    all(>(0), background) || throw("need background > 0")
    x = copy(x0)
    asum = A' * ones(eltype(ynoisy), size(ynoisy))
	time0 = time()
    for iter = 1:niter
        @show iter, extrema(x), time() - time0
        ybar = A * x .+ background # forward model
        x .*= (A' * (ynoisy ./ ybar)) ./ asum # multiplicative update
    end
    return x
end

# This preferable ML-EM version modifies the input `x`,
# so no memory allocation is needed within the loop!

function mlem!(x, ynoisy, background, A; niter::Int = 20)
    all(>(0), background) || throw("need background > 0")
    asum = A' * ones(eltype(ynoisy), size(ynoisy)) # this allocates
    ybar = similar(ynoisy)
    yratio = similar(ynoisy)
    back = similar(x)
	time0 = time()
    for iter = 1:niter
        @show iter, extrema(x), time() - time0
        mul!(ybar, A, x)
        @. yratio = ynoisy / (ybar + background) # coalesce broadcast!
        mul!(back, A', yratio) # back = A' * (ynoisy / ybar)
        @. x *= back / asum # multiplicative update
    end
    return x
end

# Apply both versions of ML-EM to this simulated data

x0 = ones(T, nx, ny, nz) # initial uniform image

niter = 30
if !@isdefined(xhat1)
    xhat1 = mlem(x0, ynoisy, scatter_mean, A; niter)
end

if !@isdefined(xhat2)
    xhat2 = copy(x0)
    mlem!(xhat2, ynoisy, scatter_mean, A; niter)
end

@assert xhat1 ≈ xhat2

jim(mid3(xhat2), "ML-EM at $niter iterations")


text:
const MESSAGES_REGEX = r"^/channels/\d+/messages/\d+$"
const ENDS_MAJOR_ID_REGEX = r"(?:channels|guilds|webhooks)/\d+$"
const ENDS_ID_REGEX = r"/\d+$"
const EXCEPT_TRAILING_ID_REGEX = r"(.*?)/\d+$"

mutable struct Bucket <: Threads.AbstractLock
    remaining::Union{Int, Nothing}
    reset::Union{DateTime, Nothing}  # UTC.
    sem::Base.Semaphore

    Bucket() = new(nothing, nothing, Base.Semaphore(1))
    Bucket(remaining::Int, reset::DateTime) = new(remaining, reset, Base.Semaphore(1))
end

mutable struct Limiter
    reset::Union{DateTime, Nothing}  # API-wide limit.
    buckets::Dict{AbstractString, Bucket}
    lock::Threads.AbstractLock

    Limiter() = new(nothing, Dict(), Threads.SpinLock())
end

function Bucket(l::Limiter, method::Symbol, endpoint::AbstractString)
    locked(l.lock) do
        endpoint = parse_endpoint(endpoint, method)
        if !haskey(l.buckets, endpoint)
            l.buckets[endpoint] = Bucket()
        end
        return l.buckets[endpoint]
    end
end

Base.lock(b::Bucket) = Base.acquire(b.sem)
Base.unlock(b::Bucket) = Base.release(b.sem)

function reset!(b::Bucket)
    b.remaining = nothing
    b.reset = nothing
end

# Note: These DateTime operations only work if your system clock is accurate.

function Base.wait(l::Limiter, b::Bucket)
    n = now(UTC)
    if l.reset !== nothing && l.reset > n
        sleep(l.reset - n)
    end

    if b.reset !== nothing && n < b.reset
        sleep(b.reset - n)
        reset!(b)
    end
end

function update!(l::Limiter, b::Bucket, r::HTTP.Messages.Response)
    if r.status == 429
        d = JSON.parse(String(copy(r.body)))
        if get(d, "global", false)
            l.reset = now(UTC) + Millisecond(get(d, "retry_after", 0))
            return
        end
    end

    headers = Dict(r.headers)

    haskey(headers, "X-RateLimit-Remaining") || return
    haskey(headers, "X-RateLimit-Reset") || return

    b.remaining = parse(Int, headers["X-RateLimit-Remaining"])
    b.reset = unix2datetime(parse(Int, headers["X-RateLimit-Reset"]))
end

function islimited(l::Limiter, b::Bucket)
    n = now(UTC)

    if l.reset !== nothing
        if n < l.reset
            return true
        else
            l.reset = nothing
        end
    end

    (b.remaining === nothing || b.reset === nothing) && return false
    if n > b.reset
        reset!(b)
        return false
    end

    return b.remaining == 0
end

function parse_endpoint(endpoint::AbstractString, method::Symbol)
    return if method === :DELETE && match(MESSAGES_REGEX, endpoint) !== nothing
        # Special case 1: Deleting messages has its own rate limit.
        first(match(EXCEPT_TRAILING_ID_REGEX, endpoint).captures) * " $method"
    elseif startswith(endpoint, "/invites/")
        # Special case 2: Invites are identified by a non-numeric code.
        "/invites"
    elseif match(ENDS_MAJOR_ID_REGEX, endpoint) !== nothing
        endpoint
    elseif match(ENDS_ID_REGEX, endpoint) !== nothing
        first(match(EXCEPT_TRAILING_ID_REGEX, endpoint).captures)
    else
        endpoint
    end
end


text:
using OpenQuantumBase, Test

@test σx*PauliVec[1][1] == PauliVec[1][1]
@test σx*PauliVec[1][2] == -PauliVec[1][2]
@test σy*PauliVec[2][1] == PauliVec[2][1]
@test σy*PauliVec[2][2] == -PauliVec[2][2]
@test σz*PauliVec[3][1] == PauliVec[3][1]
@test σz*PauliVec[3][2] == -PauliVec[3][2]

@test bloch_to_state(π/2, 0.0) ≈ PauliVec[1][1]
@test bloch_to_state(0, 0) ≈ PauliVec[3][1]
@test bloch_to_state(π/2, π/2) ≈ PauliVec[2][1]
@test_throws ArgumentError bloch_to_state(2π, 0)
@test_throws ArgumentError bloch_to_state(π, 3π)

@test creation_operator(3) ≈ [0 0 0; 1 0 0; 0 sqrt(2) 0]
@test annihilation_operator(3) ≈ [0 1 0; 0 0 sqrt(2); 0 0 0]
@test number_operator(3) ≈ [0 0 0; 0 1 0; 0 0 2]


text:
# this files downloads files for reproducing from Zenodo dataset
# https://zenodo.org/record/3823283
using HTTP, DelimitedFiles

zenodo = expanduser("~/SeisNoise_paper/DATA/zenodo.txt")
files = readdlm(zenodo,String)[:]
zenodourl = "https://zenodo.org/record/3823283/files/"

# make output directory
if !isdir(expanduser("~/SeisNoise_paper/DATA/YEARSAC"))
    mkpath(expanduser("~/SeisNoise_paper/DATA/YEARSAC"))
end

# download files from zenodo
for (ii,f) in enumerate(files)
    dlpath = joinpath(zenodourl,f * "?download=1")
    outpath = joinpath(expanduser("~/SeisNoise_paper/DATA/YEARSAC"), f)
    println("Downloading file $f, $ii of $(length(files))")
    HTTP.download(dlpath,outpath)
end


text:
module CSpPlots

using CoordinateSplittingPTrees
using Makie, GeometryTypes, Colors, PerceptualColourMaps

export plotboxes

get_finite(x, default) = isfinite(x) ? x : default
outerbounds(b1, b2) = (min(b1[1], b2[1]), max(b1[2], b2[2]))
outerbounds(b1, b2::Real) = (min(b1[1], b2), max(b1[2], b2))
outerbounds_finite(b1, b2) = (min(b1[1], get_finite(b2[1], b1[1])), max(b1[2], get_finite(b2[2], b1[2])))
widenbounds(b) = (Δx = 0.1*(b[2]-b[1]); return (b[1]-Δx, b[2]+Δx))

function finitebounds(root::Box{p,T}) where {p,T}
    @assert(ndims(root) == 2)
    xbounds, ybounds = boxbounds(root)
    if all(isfinite, xbounds) && all(isfinite, ybounds)
        return widenbounds(xbounds), widenbounds(ybounds)
    end
    # Get bounds needed to contain all finite positions inside `root`
    xbounds, ybounds = (Inf, -Inf), (Inf, -Inf)
    lower, upper = root.world.lower, root.world.upper
    bb = Vector{Tuple{T,T}}(undef, 2)
    x = [NaN, NaN]
    lfilled, rfilled = [false, false], [false, false]
    for box in leaves(root)
        for i = 1:2
            bb[i] = (lower[i], upper[i])
        end
        CoordinateSplittingPTrees.boxbounds!(bb, lfilled, rfilled, box)
        xbounds = outerbounds_finite(xbounds, bb[1])
        ybounds = outerbounds_finite(ybounds, bb[2])
        x = CoordinateSplittingPTrees.position!(x, lfilled, box)
        xbounds = outerbounds(xbounds, x[1])
        ybounds = outerbounds(ybounds, x[2])
    end
    widenbounds(xbounds), widenbounds(ybounds)
end

function plotboxes(f::Function, root::Box, xbounds::Tuple{Any,Any}, ybounds::Tuple{Any,Any};
                   clim=f.(extrema(root)),
                   cs::AbstractVector{<:Colorant}=cmap("RAINBOW3"))
    @assert(ndims(root) == 2)
    scene = Scene()
    rects, cols, points = Rectangle{Float32}[], eltype(cs)[], Point{2,Float32}[]
    for bx in leaves(root)
        addrect!(f, rects, cols, points, bx, clim, cs, xbounds, ybounds)
    end
    poly!(scene, rects, color=cols, linecolor=:black)
    scatter!(scene, points, color=:black, markersize=0.03)
    center!(scene)
    return scene
end
plotboxes(root::Box; kwargs...) =
    plotboxes(root, finitebounds(root)...; kwargs...)
plotboxes(root::Box, xbounds::Tuple{Any,Any}, ybounds::Tuple{Any,Any}; kwargs...) =
    plotboxes(value, root, xbounds, ybounds; kwargs...)
plotboxes(f::Function, root::Box; kwargs...) =
    plotboxes(f, root, finitebounds(root)...; kwargs...)

function addrect!(f, rects, cols, points, box::Box, clim, cs, xbounds, ybounds)
    @assert(isleaf(box))
    @assert(ndims(box) == 2)
    x = position(box)
    bbx, bby = boxbounds(box)
    bbx, bby = clamp.(bbx, xbounds...), clamp.(bby, ybounds...)
    rect = Rectangle{Float32}(bbx[1], bby[1], bbx[2]-bbx[1], bby[2]-bby[1])
    fval = f(box)
    if isfinite(fval)
        push!(rects, rect)
        cnorm = (clamp(fval, clim...) - clim[1])/(clim[2] - clim[1])
        col = cs[round(Int, (length(cs)-1)*cnorm) + 1]
        push!(cols, col)
    end
    push!(points, Point{2,Float32}(x[1], x[2]))
    return rects, cols, points
end

end


text:
using Interpolations: LinearInterpolation, Flat, lbounds, ubounds
using SpecialFunctions: gamma
using HDF5

"""
    line_absorption(linelist, λs, temp, nₑ, n_densities, partition_fns, ξ
                   ; α_cntm=nothing, cutoff_threshold=1e-3, window_size=20.0*1e-8)

Calculate the opacity coefficient, α, in units of cm^-1 from all lines in `linelist`, at wavelengths
`λs` [cm^-1].

other arguments:
- `temp` the temerature in K
- `n_densities`, a Dict mapping species to absolute number density [cm^-3].
- `partition_fns`, a Dict containing the partition function of each species
- `ξ` is the microturbulent velocity in cm/s (n.b. NOT km/s)

The three keyword arguments specify how the size of the window over which each line is calculated
is chosen.
- If `α_cntm` is passed as a callable returning the continuum opacity as a function of
  wavelength, the line window will extend to 4 doppler widths or the wavelength at which the Lorentz
  wings of the line are at `cutoff_threshold * α_cntm[line.wl]`, whichever is greater.
   `cutoff_threshold` defaults to 1e-3.
- if `α_cntm` is not passed, defaults to `window_size`, which is 2e-7 (in cm, i.e. 20 Å) unless
  otherwise specified
"""
function line_absorption!(α, linelist, λs, temp, nₑ, n_densities::Dict, partition_fns::Dict, ξ
                         ; α_cntm=nothing, cutoff_threshold=1e-3, window_size=20.0*1e-8)
    if length(linelist) == 0
        return zeros(length(λs))
    end

    #lb and ub are the indices to the upper and lower wavelengths in the "window", i.e. the shortest
    #and longest wavelengths which feel the effect of each line
    lb = 1
    ub = 1
    β = @. 1/(kboltz_eV * temp)
    for line in linelist
        m = get_mass(line.species)

        #doppler-broadening parameter
        σ = doppler_width.(line.wl, temp, m, ξ)

        #get all damping params from linelist.  There may be better sources for this.
        Γ = line.gamma_rad
        if !ismolecule(line.species)
            Γ = Γ .+ (nₑ .* scaled_stark.(line.gamma_stark, temp) +
                      n_densities[species"H_I"] .* [scaled_vdW(line.vdW, m, T) for T in temp])
        end
        #calculate the lorentz broadenign parameter in in wavelength. Doing this involves an
        #implicit aproximation that λ(ν) is linear over the line window.
        γ = @. Γ * line.wl^2 / c_cgs

        E_upper = line.E_lower + c_cgs * hplanck_eV / line.wl
        levels_factor = @. (exp(-β*line.E_lower) - exp(-β*E_upper)) ./ partition_fns[line.species].(temp)

        #total wl-integrated absorption coefficient
        amplitude = @. 10.0^line.log_gf*n_densities[line.species]*sigma_line(line.wl)*levels_factor

        if !isnothing(α_cntm)
            ρ_crit = [cntm(line.wl) * cutoff_threshold for cntm in α_cntm] ./ amplitude
            Δλ_D = maximum(inverse_gaussian_density.(ρ_crit, σ))
            Δλ_L = maximum(inverse_lorentz_density.(ρ_crit, γ))
            window_size = max(Δλ_D, Δλ_L)
        end
        lb, ub = move_bounds(λs, lb, ub, line.wl, window_size)
        if lb >= ub
            continue
        end

        invσ = 1.0./σ
        @inbounds view(α, :, lb:ub) .+= line_profile.(line.wl, invσ, γ, amplitude, view(λs, lb:ub)')
    end
end

"""
    inverse_gaussian_density(ρ, σ)

Calculate the inverse of a (0-centered) Gaussian PDF with standard deviation `σ`, i.e. the value of
`x` for which `ρ = exp(-0.5 x^2/σ^2}) / √[2π]`, which is given by `σ √[-2 log (√[2π]σρ)]`.  Returns
0 when ρ is larger than any value taken on by the PDF.

See also: [`inverse_lorentz_density`](@ref).
"""
function inverse_gaussian_density(ρ, σ)
    if ρ > 1/(sqrt(2π) * σ)
        0.0
    else
        σ * sqrt(-2log(sqrt(2π) * σ * ρ))
    end
end

"""
    inverse_lorentz_density(ρ, γ)

Calculate the inverse of a (0-centered) Lorentz PDF with width `γ`, i.e. the value of `x` for which
`ρ = 1 / (π γ (1 + x^2/γ^2))`, which is given by `√[γ/(πρ) - γ^2]`. Returns 0 when ρ is larger than
any value taken on by the PDF.

See also: [`inverse_gaussian_density`](@ref).
"""
function inverse_lorentz_density(ρ, γ)
    if ρ > 1/(π*γ)
        0.0
    else
        sqrt(γ/(π * ρ) - γ^2)
    end
end

"""
    load_hydrogen_stark_profiles()

Load the (doppler-convolved) Stark-broadening profiles from disk.
"""
function setup_hydrogen_stark_profiles(fname=joinpath(_data_dir,
                                                      "Stehle-Hutchson-hydrogen-profiles.h5"))
    hline_stark_profiles = h5open(fname, "r") do fid
        map(keys(fid)) do transition
            temps = read(fid[transition], "temps")
            nes = read(fid[transition], "electron_number_densities")
            delta_nu_over_F0 = read(fid[transition], "delta_nu_over_F0")

            siz = length(delta_nu_over_F0)
            P = read(fid[transition], "profile")

            (temps=temps,
             electron_number_densities=nes,
             #flat BCs to allow wls very close to the line center
             profile=LinearInterpolation((temps, nes, [floatmin() ; log.(delta_nu_over_F0[2:end])]),
                                         log.(P); extrapolation_bc=Flat()),
             lower = HDF5.read_attribute(fid[transition], "lower"),
             upper = HDF5.read_attribute(fid[transition], "upper"),
             Kalpha = HDF5.read_attribute(fid[transition], "Kalpha"),
             log_gf = HDF5.read_attribute(fid[transition], "log_gf"),
             λ0 = LinearInterpolation((temps, nes), read(fid[transition], "lambda0") * 1e-8)
            )
        end
    end
end

#used in hydrogen_line_absorption
_zero2epsilon(x) = x + (x == 0) * floatmin()

"""
    hydrogen_line_absorption(λs, T, nₑ, nH_I, UH_I, hline_stark_profiles, ξ;
                                  stark_window_size=3e-7, self_window_size=1e-6)

Calculate contribution to the the absorption coefficient, α, from hydrogen lines in units of cm^-1,
at wavelengths `λs`.

Arguments:
- `T`: temperature [K]
- `nₑ`: electron number density [cm^-3]
- `nH_I`: neutral hydrogen number density [cm^-3]
- `hline_stark_profiles`: (returned by setup_hline_stark_profiles`)
- `ξ`: microturbulent velocity [cm/s]. This is only applied to Hα-Hγ.  Other hydrogen lines profiles
   are dominated by stark broadening, and the stark broadened profiles are pre-convolved with a
   doppler profile.
- `stark_window_size`: the max distance from each line center [cm] at which to calculate the line
   absorption for lines besides Hα-Hγ (those dominated by stark broadening).
- `self_window_size`: the max distance from each line center [cm] at which to calculate the line
   absorption for Hα-Hγ (those dominated by self-broadening).
"""
function hydrogen_line_absorption(λs, T, nₑ, nH_I, UH_I, hline_stark_profiles, ξ;
                                  stark_window_size=3e-7, self_window_size=1e-6)
    νs = c_cgs ./ λs
    dνdλ = c_cgs ./ λs.^2
    #This is the Holtzmark field, by which the frequency-unit-detunings are divided for the
    #interpolated stark profiles
    F0 = 1.25e-9 * nₑ^(2/3)
    α_hlines = zeros(length(λs))
    for line in hline_stark_profiles
        if !all(lbounds(line.λ0.itp)[1:2] .< (T, nₑ) .< ubounds(line.λ0.itp)[1:2])
            continue #transitions to high levels are omitted for high nₑ and T
        end
        λ₀ = line.λ0(T, nₑ)

        Elo = RydbergH_eV * (1 - 1/line.lower^2)
        Eup = RydbergH_eV * (1 - 1/line.upper^2)
        β = 1/(kboltz_eV * T)
        levels_factor = (exp(-β*Elo) - exp(-β*Eup)) / UH_I(T)
        amplitude = 10.0^line.log_gf * nH_I * sigma_line(λ₀) * levels_factor

        #compute profile with either self or stark broadening
        Hmass = get_mass(Formula("H"))
        if line.lower == 2 && line.upper in [3, 4, 5]
            #ABO params and line center
            λ₀, σ, α = if line.upper == 3
                6.56460998e-5, 1180.0, 0.677
            elseif line.upper == 4
                4.8626810200000004e-5, 2320.0, 0.455
            elseif line.upper == 5
                4.34168232e-5, 4208.0, 0.380
            end

            #use Barklem+ 2000 p-d approximation for resonant broadening of first 3 balmer lines
            #λ₀ may not be the most appropriate choice here?
            Δλ_D = doppler_width(λ₀, T, Hmass, ξ)

            Γ = scaled_vdW((σ*bohr_radius_cgs^2, α), Hmass, T) * nH_I
            Δλ_L = Γ * λ₀^2 / c_cgs

            lb, ub = move_bounds(λs, 0, 0, λ₀, self_window_size)
            if lb == ub
                continue
            end
            @inbounds view(α_hlines, lb:ub) .+= line_profile.(λ₀, 1.0/Δλ_D, Δλ_L, amplitude,
                                                             view(λs, lb:ub))
        else #use Stehle+ 1999 Stark-broadened profiles
            lb, ub = move_bounds(λs, 0, 0, λ₀, stark_window_size)
            if lb == ub
                continue
            end

            ν₀ = c_cgs / (λ₀)
            scaled_Δν = _zero2epsilon.(abs.(view(νs,lb:ub) .- ν₀) ./ F0)
            dIdν = exp.(line.profile.(T, nₑ, log.(scaled_Δν)))
            @inbounds view(α_hlines,lb:ub) .+= dIdν .* view(dνdλ, lb:ub) .* amplitude
        end
    end
    α_hlines
end

"the width of the doppler-broadening profile"
doppler_width(λ₀, T, m, ξ) = λ₀ * sqrt(2kboltz_cgs*T / m + ξ^2) / c_cgs

"the stark broadening gamma scaled acording to its temperature dependence"
scaled_stark(γstark, T; T₀=10_000) = γstark * (T/T₀)^(1/6)

"""
the vdW broadening gamma scaled acording to its temperature dependence, using either simple scaling
or ABO
"""
scaled_vdW(vdW::Real, m, T, T₀=10_000) = vdW * (T/T₀)^0.3
function scaled_vdW(vdW::Tuple{F, F}, m, T) where F <: Real
    v₀ = 1e6 #σ is given at 10_000 m/s = 10^6 cm/s
    σ = vdW[1]
    α = vdW[2]

    invμ = 1/(1.008*amu_cgs) + 1/m #inverse reduced mass
    vbar = sqrt(8 * kboltz_cgs * T / π * invμ) #relative velocity
    #n.b. "gamma" is the gamma function, not a broadening parameter
    2 * (4/π)^(α/2) * gamma((4-α)/2) * v₀ * σ * (vbar/v₀)^(1-α)
end

function move_bounds(λs::AbstractRange, lb, ub, λ₀, window_size)
    len = length(λs)
    lb = clamp(Int(cld(λ₀ - window_size - λs[1], step(λs)) + 1), 1, len)
    ub = clamp(Int(fld(λ₀ + window_size - λs[1], step(λs)) + 1), 1, len)
    lb,ub
end
#walk lb and ub to be window_size away from λ₀. assumes λs is sorted
function move_bounds(λs, lb, ub, λ₀, window_size)
    while lb+1 < length(λs) && λs[lb] < λ₀ - window_size
        lb += 1
    end
    while lb > 1 && λs[lb-1] > λ₀ - window_size
        lb -= 1
    end
    while ub < length(λs) && λs[ub+1] < λ₀ + window_size
        ub += 1
    end
    while ub > 1 && λs[ub] > λ₀ + window_size
        ub -= 1
    end
    lb, ub
end

"""
    sigma_line(wl)

The cross-section (divided by gf) at wavelength `wl` in Ångstroms of a transition for which the
product of the degeneracy and oscillator strength is `10^log_gf`.
"""
function sigma_line(λ::Real)
    #work in cgs
    e  = electron_charge_cgs
    mₑ = electron_mass_cgs
    c  = c_cgs

    #the factor of |dλ/dν| = λ²/c is because we are working in wavelength rather than frequency
    (π*e^2/mₑ/c) * (λ^2/c)
end

"""
    line_profile(λ₀, invΔλ_D, Δλ_L, line_amplitude, λ)

A voigt profile centered on λ₀ with Doppler width 1/`invΔλ_D` and Lorentz width `Δλ_L`
evaluated at `λ` (cm).  Note that this returns values in units of cm^-1.
"""
function line_profile(λ₀::Real, invΔλ_D::Real, Δλ_L::Real, line_amplitude::Real, λ::Real)
    _line_profile(λ₀, invΔλ_D, Δλ_L*invΔλ_D/(4π), line_amplitude*invΔλ_D/sqrt(π), λ)
end
function _line_profile(λ₀::Real, invΔλ_D::Real, Δλ_L_invΔλ_D_div_4π::Real,
                       amplitude_invΔλ_D_div_sqrt_π::Real, λ::Real)
    voigt(Δλ_L_invΔλ_D_div_4π, abs(λ-λ₀) * invΔλ_D) * amplitude_invΔλ_D_div_sqrt_π
end

@inline function harris_series(v) # assume v < 5
    v2 = v*v
    H₀ = exp(-(v2))
    H₁ = if (v < 1.3)
        -1.12470432 + (-0.15516677 + (3.288675912 +(-2.34357915 + 0.42139162*v)*v)*v)*v
    elseif v < 2.4
        -4.48480194 + (9.39456063 + (-6.61487486 + (1.98919585  - 0.22041650*v)*v)*v)*v
    else #v < 5
        ((0.554153432 + (0.278711796 + (-0.1883256872 + (0.042991293 - 0.003278278*v)*v)*v)*v) /
         (v2 - 3/2))
    end
    H₂ = (1 - 2v2) * H₀
    H₀, H₁, H₂
end

"""
    voigt(α, v)

The [voigt function](https://en.wikipedia.org/wiki/Voigt_profile#Voigt_functions), ``H``.
Approximation from Hunger 1965.
"""
function voigt(α, v)
    v2 = v*v
    if α <= 0.2 && (v >= 5)
        invv2 = (1/v2)
        (α/sqrt(π) * invv2) * (1 + 1.5invv2 + 3.75*invv2^2)
    elseif α <= 0.2 #v < 5
        H₀, H₁, H₂ = harris_series(v)
        H₀ + (H₁ + H₂*α)*α
    elseif (α <= 1.4) && (α + v < 3.2)
        #modified harris series: M_i is H'_i in the source text
        H₀, H₁, H₂ = harris_series(v)
        M₀ = H₀
        M₁ = H₁ + 2/sqrt(π) * M₀
        M₂ = H₂ - M₀ + 2/sqrt(π) * M₁
        M₃ = 2/(3sqrt(π))*(1-H₂) - (2/3)*v2*M₁ + (2/sqrt(π))*M₂
        M₄ = 2/3 * v2*v2 * M₀ - 2/(3sqrt(π)) * M₁ + 2/sqrt(π) * M₃
        ψ =  0.979895023 + (-0.962846325 + (0.532770573 - 0.122727278*α)*α)*α
        ψ * (M₀ + (M₁ + (M₂ + (M₃ + M₄*α)*α)*α)*α)
    else #α > 1.4 or (α > 0.2 and α + v > 3.2)
        r2 = (v2)/(α*α)
        α_invu = 1/sqrt(2) / ((r2 + 1) * α)
        α2_invu2 = α_invu * α_invu
        sqrt(2/π) * α_invu * (1 + (3*r2 - 1 + ((r2-2)*15*r2+2)*α2_invu2)*α2_invu2)
    end
end


text:
using CompileBot

# using runtests:
snoop_bot(BotConfig("TestPackage2", os = ["windows", "linux"], tmin =0.0))