text:
daysofweek = ["Mo", "Tu", "We", "Th", "Fr", "Sa","Su"]
specialdays = ["PH", "SH"]

"""Parse a conditional tag of the following structure:
<restriction-value> @ <condition>[;<restriction-value> @ <condition>]

Rturns a dictionary of complonents.
"""
function parseconditionaltag(tag::AbstractString)::Dict{Symbol, Dict}
    # regex for values: find value-parts of string
    regexvalue = r"\;?\s?(\-?[[:alnum:]\_\.]*)\s?@"
    values = collect(m.captures[1] for m = eachmatch(regexvalue, tag))

    # regex for condition: find codition-parts of string
    regexcondi1 = r"@\s?\(([[:alnum:]\s\;\:\,\.\-\<\>]+)\)\;?"
    regexcondi2 = r"@\s?([[:alnum:]\:\,\.\-\<\>]+\s?[[:alnum:]\:\,\.\-\<\>]+)\;?"
    condis = [collect(m.captures[1] for m = eachmatch(regexcondi1, tag));
              collect(m.captures[1] for m = eachmatch(regexcondi2, tag))]
    values = [string(i) for i = values]
    condis = [string(i) for i = condis]

    # check if result is plausible
    # todo: allow turning this off, as it is costly
    if length(values) == length(condis)
        norg = length(tag)
        norgspe = sum(length.(findall.(isequal.(['@', ';', '(', ')', ' ']), tag)))
        nres = sum([length.(i) for i in [values; condis]])
        nresspe = sum([sum(length.(findall.(isequal.([';', ' ']), i))) for i = [values; condis]])
        mengthmatch = norg - (nres + norgspe - nresspe)
        if mengthmatch != 0
            @error "The tokenizer of the conditional tags resulted in an implausible character count" condition=tag tokens=[values, condis]
        end
    else
        @error "The number of value and condition tokens don't match" values=values conditions=condis
    end

    # gether parsed data into a dictionary
    out = Dict{Symbol, Dict}()
    count = 1
    for (v, c) = zip(values, condis)
        out[Symbol("condition" * string(count))] =  Dict{Symbol, Union{Dict, AbstractString, Integer}}(
            :value => v,
            :rules => parsecondition(c)
        )

        count += 1
    end

    return out
end

""" Parse condition part of a conditional tag.  Returns a rules as dictionary.
"""
function parsecondition(condition::AbstractString)::Dict{Symbol, Dict}
    condis = strip.(split(condition, ';'))

    out = Dict{Symbol, Dict}()
    count = 1
    for token = condis
        # todo: capture AND !!!!!!!
        # haven't seen it in the tags of intrest, but it is part of the syntax
        days = ruledays(token)
        hours = rulehours(token)
        daytimeword = Dict{Symbol, Union{AbstractArray, Dict}}(
            :words => rulewords(token),
            :hoursofweek => hoursofweek(hours, days[1], days[2])
        )
        out[Symbol("rule" * string(count))] = daytimeword

        count += 1
    end

    return out
end

""" Parse time (hour of day) component of the condition token as array of int.
"""
function rulehours(token::AbstractString)::AbstractArray{UInt8}
    token = replace(token, "24:00" => "00:00")
    # one or more time intervals
    regextime = r"(\d{2}:\d{2})-(\d{2}:\d{2})"
    timeinterval = collect(string.(m.captures) for m = eachmatch(regextime, token))
    if length(timeinterval) > 0
        starttimes, endtimes = [[Dates.hour(Dates.Time.(i)) for i = j] for j = zip(timeinterval...)]
        endtimes = [i == 0 ? 24 : i for i = endtimes]
        if length(starttimes) > 1 &&
                issorted([i for j = zip(starttimes, endtimes) for i = j])
            hoursofday =  UInt8[i for (s, e) = zip(starttimes, endtimes) for i in collect(s:(e-1))]
        elseif length(starttimes) == 1
            if starttimes[1] > endtimes[1]
                hoursofday = [collect(UInt8, 0:(endtimes[1]-1)); collect(UInt8, starttimes[1]:23)]
            else
                hoursofday = collect(UInt8, starttimes[1]:(endtimes[1]-1))
            end
        else
            throw(ErrorException("The time part of the condition can't be interpreted: $token"))
        end

    else
        hoursofday = UInt8[]
    end

    return hoursofday
end

""" Parse the day component of the condition token. Return dict containng a array
of `daysofweek` as numbers between 1  (Mo) and 7 (Su) and `specialdays` array of strings.
"""
function ruledays(token::AbstractString)::AbstractArray
    regexdayint = r"(Mo|Tu|We|Th|Fr|Sa|Su)-(Mo|Tu|We|Th|Fr|Sa|Su)"
    dayinterval = match(regexdayint, token)
    if dayinterval != nothing
        firstday = findfirst(dayinterval.captures[1] .== daysofweek)
        lastday = findfirst(dayinterval.captures[2] .== daysofweek)
        if firstday <= lastday
            days = collect(UInt8, firstday:lastday)
            sdays = String[]
        else
            days = [collect(UInt8, 1:lastday); collect(UInt8, firstday:7)]
            sdays = String[]
        end
    else
        regexdaycol = r"(Mo|Tu|We|Th|Fr|Sa|Su|PH)"
        dayscolect = collect(string(m.captures[1]) for m = eachmatch(regexdaycol, token))
        days = UInt8[]
        sdays = String[]
        if length(dayscolect) > 0
            for dc = dayscolect
                d = findfirst(dc .== daysofweek)
                if d != nothing
                    days = UInt8[days; d]
                end
                sd = specialdays[dc .== specialdays]
                if sd != nothing
                    sdays = [sdays; sd]
                end
            end
        end
    end

    return[days, sdays]
end

""" Find special word in condition token. Return array of Strings
"""
function rulewords(
        token::AbstractString
    )::AbstractArray{AbstractString, 1}
    regexword = r"(off|signal|wet)"
    words = collect(string(m.captures[1]) for m = eachmatch(regexword, token))

    return words
end

""" Use day of week and hour of day and return the hour of week
"""
function hoursofweek(
        hours::Array{<:Integer,1}=UInt8[],
        days::Array{<:Integer,1}=UInt8[],
        specialdays::Array{<:AbstractString,1}=String[]
    )::Array{Integer,1}
    nd = length(days)
    nh = length(hours)
    nsd = length(specialdays)

    # if days and hours are empty return empty array
    # if days are empty but specialdays are not, return empty array
    if nd == 0 && nh == 0 || nd == 0 && nsd != 0
        return Array{UInt8,1}[]
    # if only days are empty generate all days
    elseif nd == 0 && nsd == 0
        days = collect(UInt8, 1:7)
        nd = length(days)
    # if hours are empty generate all hours
    elseif nh == 0
        hours = collect(UInt8, 0:23)
        nh = length(hours)
    end
    daysashours = (days .- 1) .* 24

    # hours per day
    out = reshape((daysashours' .+ hours), nd * nh)

    return out
end


text:
#=
import Pkg; Pkg.activate(".")
=#
using Distributions
import Random
using RCall
import LinearAlgebra

include("MCMC.jl")

prettymat(X) = begin show(stdout, "text/plain", X); println(); end

### MAIN ###
Random.seed!(1)

# easy = true  # bivariate normal
easy = false  # 250-variate normal

if easy
  K = 3
  S = MCMC.eye(K)
  S[1,2] = S[2,1] = -0.8
  S[1,3] = S[3,1] = 0.5
else
  K = 50
  S = rand(InverseWishart(K, MCMC.eye(K)))
end

mvn = MvNormal(randn(K) * 10, S)
logprob(x::Vector{Float64}) = logpdf(mvn, x)

init = randn(K)

lp(s::Vector{Float64}) = 0.0
out, propcov = MCMC.mcmc(init, logprob, lp, thin=10, batchsize=100,
                         temper_power=3,
                         propcov_init=MCMC.eye(length(init)) * 1e-3,
                         nburn=10000, nmcmc=1000, max_temper=1000.0,
                         verbose=1)

acc_rate = size(unique(out, dims=1), 1) / size(out, 1)

ll = [logprob(out[i, :]) for i in 1:size(out, 1)]
R"plot($ll, type='l')"

G = 8
println("cov(out):")
prettymat(round.(cov(out[:, 1:G]), digits=3))
println("cov(mvn):")
prettymat(round.(cov(mvn)[1:G, 1:G], digits=3))

[mean(out, dims=1)' mean(mvn)]

R"rcommon::plotPosts($(out)[, 1:5])";
println("Acceptance rate: $(acc_rate)")

ext: jl


text:
"""
    IntersectionProperty{N<:Real} <: Property

Type that represents a safety property characterized with a set of bad states.
A safety property is satisfied by a given set of states if the intersection with
the set of bad states is empty.

### Fields

- ``bad`` -- convex set representing the bad states
- ``witness`` -- witness point (empty vector if not set)
"""
mutable struct IntersectionProperty{N<:Real} <: Property
    bad::LazySet
    witness::Vector{N}

    IntersectionProperty{N}(bad::LazySet) where {N<:Real} = new(bad, N[])
end

# type-less convenience constructor
IntersectionProperty(bad::LazySet{N}) where {N} = IntersectionProperty{N}(bad)

"""
    inout_map_property(prop::IntersectionProperty,
                       partition::AbstractVector{<:AbstractVector{Int}},
                       blocks::AbstractVector{Int},
                       n::Int
                      )::IntersectionProperty

Map an `IntersectionProperty` to the dimensions of analyzed blocks.

### Input

- `prop`      -- property
- `partition` -- block partition; elements are start and end indices of a block
- `blocks`    -- list of all output block indices in the partition
- `n`         -- total number of input dimensions

### Output

A new property of reduced dimension.

### Notes

If the dimension is not reduced, we keep the original set.
Otherwise, the dimension reduction is achieved with a `LinearMap`.
"""
function inout_map_property(prop::IntersectionProperty{N},
                            partition::AbstractVector{<:AbstractVector{Int}},
                            blocks::AbstractVector{Int},
                            n::Int
                           )::IntersectionProperty{N} where {N<:Real}
    @assert dim(prop.bad) == n "the property has dimension $(dim(prop.bad)) but should have dimension $n"
    proj = projection_map(partition, blocks)
    if length(proj) == n
        # no change in the dimension, copy the old property (keep the set)
        return IntersectionProperty(prop.bad)
    else
        M = sparse(proj, proj, ones(N, length(proj)), n, n)
        return IntersectionProperty(M * prop.bad)
    end
end

"""
    check_property(set::LazySet, prop::IntersectionProperty)::Bool

Checks whether a convex set satisfies a safety property.

### Input

- ``set``  -- convex set
- ``prop`` -- safety property with bad states

### Output

`true` iff the safety property is satisfied for the given set of states.
This is the case iff the set of states does not intersect with the bad states.
"""
@inline function check_property(set::LazySet, prop::IntersectionProperty)::Bool
    nonempty_intersection, witness = is_intersection_empty(set, prop.bad, true)
    if nonempty_intersection
        # store violation witness
        prop.witness = witness
    end
    return !nonempty_intersection
end


text:
Base.@deprecate lla2ecef ECEF
Base.@deprecate ecef2lla LLA
Base.@deprecate ecef2enu ENU
Base.@deprecate lla2enu ENU
Base.@deprecate centerBounds(bounds::Bounds) center(bounds)


text:
# Copyright (c) <2021>, <Sven Leyffer & Jongeun Kim>
# All rights reserved.

# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

using Pkg

function install_packages()
    pkgs = ["CSV", "DataFrames", "DataStructures", "Glob", "JuMP", "MathOptInterface", "SCIP"]
    Pkg.add(pkgs)
    Pkg.status()
end

install_packages()

max_stars_repo_name: jongeunkim/STreCH

max_issues_repo_name: jongeunkim/STreCH

max_forks_repo_name: jongeunkim/STreCH


text:
## ExoplanetsSysSim/examples/setup_abc_christiansen.jl
## (c) 2016 Eric B. Ford & Danley C. Hsu

module EvalSysSimModel
  export setup, get_param_vector, get_ss_obs #, evaluate_model
  export gen_data, calc_summary_stats, calc_distance, is_valid
  using ExoplanetsSysSim
  include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "christiansen_func.jl"))
  include(joinpath(pwd(), "param_file.jl"))

  sim_param_closure = SimParam()
  summary_stat_ref_closure =  CatalogSummaryStatistics()

    function is_valid(param_vector::Vector{Float64})
      global sim_param_closure
      update_sim_param_from_vector!(param_vector,sim_param_closure)
      const sigma_log_r = get_real(sim_param_closure,"sigma_log_planet_radius")
      const sigma_log_P = get_real(sim_param_closure,"sigma_log_planet_period")
      const rate_tab::Array{Float64,2} = get_any(sim_param_closure, "obs_par", Array{Float64,2})
      const lambda = sum_kbn(rate_tab)
    # println("Lambda = ", lambda)
      if sigma_log_r <= 0. || sigma_log_P<=0. || lambda > 10. || any(x -> x < 0., rate_tab)
         return false
      end
      return true
    end

#=
    function norm_christiansen(param_vector::Vector{Float64})
      log_z_sum = logsumexp(param_vector)
      param_vector = exp(param_vector-z_sum)
    end
=#

    function gen_data(param_vector::Vector{Float64})
      global sim_param_closure
      update_sim_param_from_vector!(param_vector,sim_param_closure)
      cat_phys = generate_kepler_physical_catalog(sim_param_closure)
      cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys, sim_param_closure)
      cat_obs = ExoplanetsSysSim.observe_kepler_targets_single_obs(cat_phys_cut, sim_param_closure)
      return cat_obs
    end

    function calc_summary_stats(cat::KeplerObsCatalog )
      global sim_param_closure
      #sum_stat = calc_summary_stats_obs_binned_rates(cat[1], sim_param_closure)
      sum_stat = calc_summary_stats_obs_binned_rates(cat, sim_param_closure)
      return sum_stat
    end

    function calc_distance(sum_stat_obs::CatalogSummaryStatistics,sum_stat_sim::CatalogSummaryStatistics, n::Integer = 0)
      global sim_param_closure
      dist1 = calc_distance_vector_binned(sum_stat_obs,sum_stat_sim, 1, sim_param_closure)
      num_available = length(dist1)
      num_to_use = n>0 ? min(n,num_available) : num_available
      return calc_scalar_distance(dist1[1:num_to_use])
    end

  function setup()
    global sim_param_closure = setup_sim_param_christiansen()

    sim_param_closure = set_test_param(sim_param_closure)
    add_param_fixed(sim_param_closure,"stellar_catalog","q1_q12_christiansen.jld")

    ### Use simulated planet candidate catalog data
    add_param_fixed(sim_param_closure,"num_kepler_targets",150969)  # For "observed" catalog
    cat_obs = simulated_read_kepler_observations(sim_param_closure)
    #inv_det_prob(cat_obs, sim_param_closure)
    #println("------------------------")
    global summary_stat_ref_closure = calc_summary_stats_obs_binned_rates(cat_obs,sim_param_closure, trueobs_cat = true)
    ###

    ### Use real planet candidate catalog data
    #add_param_fixed(sim_param_closure,"koi_catalog","q1_q16_koi_cand.csv")
    #cat_obs = setup_actual_planet_candidate_catalog(setup_star_table_christiansen(sim_param_closure), sim_param_closure)
    #global summary_stat_ref_closure = calc_summary_stats_obs_binned_rates(cat_obs,sim_param_closure, trueobs_cat = true)
    ###

    num_targ = get_int(sim_param_closure,"num_kepler_targets")
    add_param_fixed(sim_param_closure,"num_targets_sim_pass_one",num_targ)  # Set universal simulated catalog size

  end

  get_param_vector() = make_vector_of_sim_param(sim_param_closure)
  get_ss_obs() = summary_stat_ref_closure

  function set_simparam_ss(sim_param::ExoplanetsSysSim.SimParam, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics)
    global sim_param_closure = sim_param
    global summary_stat_ref_closure = ss_true
  end

end  # module EvalSysSimModel

include(joinpath(Pkg.dir("ABC"),"src/composite.jl"))

module SysSimABC
  export setup_abc, run_abc
  import ABC
  import Distributions
  using CompositeDistributions
  import ExoplanetsSysSim
  import EvalSysSimModel
  include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "christiansen_func.jl"))
  include(joinpath(pwd(), "param_file.jl"))

  function setup_abc(num_dist::Integer = 0)
    EvalSysSimModel.setup()
    #println("# setup_abc: ",EvalSysSimModel.sim_param_closure)
    theta_true = EvalSysSimModel.get_param_vector()
    param_prior = CompositeDist( Distributions.ContinuousDistribution[Distributions.Uniform(0., 1.) for x in 1:length(theta_true)] )
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats,calc_distance_ltd, param_prior, make_proposal_dist=ABC.make_proposal_dist_gaussian_diag_covar, is_valid=EvalSysSimModel.is_valid,
                                     num_part=40, num_max_attempt=50, num_max_times=200, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel, adaptive_quantiles = false, epsilon_reduction_factor=0.9, tau_factor=2.0);
  end

  function run_abc_largegen(pop::ABC.abc_population_type, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics, epshist_targ::Float64, npart::Integer = 1000, num_dist::Integer = 0)
    sim_param_closure = setup_sim_param_christiansen()
    sim_param_closure = set_test_param(sim_param_closure)
    #add_param_fixed(sim_param_closure,"num_kepler_targets",150000)  # For "observed" catalog
    EvalSysSimModel.set_simparam_ss(sim_param_closure, ss_true)

    theta_true = EvalSysSimModel.get_param_vector()
    param_prior = CompositeDist( Distributions.ContinuousDistribution[Distributions.Uniform(0., 1.) for x in 1:length(theta_true)] )
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats,calc_distance_ltd, param_prior, is_valid=EvalSysSimModel.is_valid,
                                     num_part=npart, num_max_attempt=50, num_max_times=1, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel);

    println("# run_abc_largegen: ",EvalSysSimModel.sim_param_closure)
    sampler_plot = abc_plan.make_proposal_dist(pop, abc_plan.tau_factor)
    theta_plot = Array{Float64}(size(pop.theta, 1), npart)
    for i in 1:npart
      theta_plot[:,i], dist_plot, attempts_plot = ABC.generate_theta(abc_plan, sampler_plot, ss_true, epshist_targ)
    end
    return theta_plot
  end

  function run_abc(abc_plan::ABC.abc_pmc_plan_type)
    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    #println("True catalog SS: ", ss_true)
    pop_out = ABC.run_abc(abc_plan,ss_true;verbose=true)
  end

  function run_abc(abc_plan::ABC.abc_pmc_plan_type, pop::ABC.abc_population_type)
    dist_threshold = maximum(pop.dist)
    EvalSysSimModel.add_param_fixed(EvalSysSimModel.sim_param_closure,"minimum ABC dist skip pass 2",dist_threshold)
    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    pop_out = ABC.run_abc(abc_plan,ss_true,pop;verbose=true)
  end

end # module SysSimABC


text:
"""
`z` is required to be 0 or 1.
With probability `probFlip`, change z to 1 or 0.
"""
function flip_bit(z::Bool, probFlip::Float64)::Bool
  return probFlip > rand() ? !z : z
end

"""
`Z` is required to be AbstractArray{Bool}.
Flip a random selection of `num_bits` bits.
The resulting `Z` should be different from the original `Z` by
`num_bits` bits.
"""
function flip_bits(Z::AbstractArray{Bool}, num_bits::Integer)
  idx = CartesianIndices(Z)
  idx_to_flip = Distributions.sample(idx, num_bits, replace=false)
  Z_new = deepcopy(Z)
  Z_new[idx_to_flip] .= .!Z_new[idx_to_flip]
  return Z_new
end

"""
default similarity function used in computing log probability of Z_repFAM.
= sum(abs.(z1 .- z2))
"""
function similarity_default(z1::Vector{Bool}, z2::Vector{Bool})::Float64
  exp(-sum(abs.(z1 .- z2)))
end

"""
log penalty term as related to column k only.
"""
function log_penalty_repFAM(k::Integer, Z::Matrix{Bool}, similarity::Function)
  K = size(Z, 2)

  log_penalty = 0.0

  for q in 1:K
    if q != k
      log_penalty += MCMC.log1m(similarity(Z[:, k], Z[:, q]))
    end
  end

  return log_penalty
end

"""
log penalty term in repFAM.
"""
function log_penalty_repFAM(Z::Matrix{Bool}, similarity::Function)::Float64
  K = size(Z, 2)

  log_penalty = 0.0

  for k1 in 1:(K - 1)
    for k2 in (k1 + 1):K
      log_penalty += MCMC.log1m(similarity(Z[:, k1], Z[:, k2]))
    end
  end

  return log_penalty
end

"""
log probability of Z ~ repFam_K(v, C), WITHOUT NORMALIZING CONSTANT
where v ~ Beta(a/K, 1),
and similarity(z_{k1}, z_{k2}) computes the similarity of binary vectors z_{k1}
and z_{k2}.

We require similarity(⋅, ⋅) ∈ (0, 1). And when z_{k1} == z_{k2} exactly,
similarity = 1. Similarly, when distance between z_{k1} and z_{k2} approaches
∞, similarity = 0.
"""
function logprob_Z_repFAM(Z::Matrix{Bool}, v::Vector{Float64},
                          similarity::Function)::Float64
  J, K = size(Z)

  # IBP component
  lp = sum(Z * log.(v) + (1 .- Z) * MCMC.log1m.( v))

  # Repulsive component
  log_penalty = log_penalty_repFAM(Z, similarity)

  return lp + log_penalty
end

function logfullcond_Z_repfam(Z::Matrix{Bool}, s::State, c::Constants,
                              d::Data, t::Tuners, sb_ibp::Bool)
  v = sb_ibp ? cumprod(s.v) : s.v

  # Log prior
  lp = logprob_Z_repFAM(Z, v, c.similarity_Z)

  # log likelihood
  ll = 0.0

  if lp > -Inf
    for i in 1:d.I
      for j in 1:d.J
        for n in 1:d.N[i]
          # Cell type
          k = s.lam[i][n]
          if k > 0  # if not noisy cell type
            z = Z[j, k]
            ll += log(dmixture(z, i, n, j, s, c, d))
          end
        end
      end
    end
  end

  return lp + ll
end

function update_Z_repFAM!(cand_Z::Matrix{Bool}, s::State, c::Constants,
                          d::Data, t::Tuners, sb_ibp::Bool; debug::Bool=false)
  # Current Z
  curr_Z = s.Z

  # Log full conditional for candidate Z and current Z
  log_fc_cand = logfullcond_Z_repfam(cand_Z, s, c, d, t, sb_ibp)
  log_fc_curr = logfullcond_Z_repfam(curr_Z, s, c, d, t, sb_ibp)

  # Whether or not to accept proposed Z
  accept = log_fc_cand - log_fc_curr > log(rand())

  if accept
    if debug
      println("Joint-update of Z resulted in a move.")
    end
    s.Z .= cand_Z
  end

  # Update tuning param
  # MCMC.update_tuning_param_default(t.Z, accept)
end

function update_Z_repFAM!(j::Integer, k::Integer,
                          s::State, c::Constants, d::Data, t::Tuners,
                          sb_ibp::Bool; debug::Bool=false)
  # Gibbs:
  Z0 = deepcopy(s.Z)
  Z0[j, k] = false

  Z1 = deepcopy(s.Z)
  Z1[j, k] = true

  log_prob_0 = logfullcond_Z_repfam(Z0, s, c, d, t, sb_ibp)
  log_prob_1 = logfullcond_Z_repfam(Z1, s, c, d, t, sb_ibp)

  p1 = 1 / (1 + exp(log_prob_0 - log_prob_1))

  s.Z[j, k] = p1 > rand()

  # Or Metropolis?
  # cand_Z = deepcopy(s.Z)
  # cand_Z[j, k] = !s.Z[j, k]
  # update_Z_repFAM!(cand_Z, s, c, d, t, sb_ibp, debug=debug)
end

function update_Z_repFAM!(s::State, c::Constants, d::Data, t::Tuners,
                          sb_ibp::Bool; debug::Bool=false)
  if c.probFlip_Z == 1  # Update each Z[j, k] sequentially, instead of jointly.
    for j in 1:d.J
      for k in 1:c.K
        update_Z_repFAM!(j, k, s, c, d, t, sb_ibp, debug=debug)
      end
    end
    # Update a bunch at a time
    cand_Z = flip_bits(s.Z, 3)
    update_Z_repFAM!(cand_Z, s, c, d, t, sb_ibp, debug=debug)
  else  # Flip only some bits to obtain proposed Z.
    num_bits_to_flip = round(Int, c.probFlip_Z * d.J * c.K)
    num_bits_to_flip = clamp(num_bits_to_flip, 1, d.J * c.K)

    # cand_Z = flip_bit.(s.Z, MCMC.logit(t.Z.value, a=0.0, b=1.0))
    # cand_Z = Matrix{Bool}(flip_bit.(s.Z, c.probFlip_Z))
    cand_Z = flip_bits(s.Z, num_bits_to_flip)
    update_Z_repFAM!(cand_Z, s, c, d, t, sb_ibp, debug=debug)
  end
end

# function update_Z_repFAM!(s::State, c::Constants, d::Data, t::Tuners,
#                           sb_ibp::Bool; config::Any=nothing)
#
# end

#= TODO:
Let R_t = number of rows to update at MCMC iteration t.
Require that
- P(R_1 = J) = 0.95
- P(R_∞ = j) = .05
- R_t ~ Bin(J - 1; exp(-(t - 1) / 1000)) + 1

1. Row-by-row updates for R_t rows, updating multiple elements at a time.
2. Row-by-row updates, updating one element at at a time.
=#

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
@testset "DataSplits" begin
    @testset "calc_mask" begin
        pools = fill(:U, 5)
        pools[1] = :Lin
        pools[2] = :Lout
        init_mask = trues(5)
        init_mask[5] = false
        @test calc_mask(FullSplitStrat(), init_mask, pools) == init_mask
        @test calc_mask(UnlabeledSplitStrat(), init_mask, pools) == [false, false, true, true, false]
        @test calc_mask(UnlabeledAndLabeledInlierSplitStrat(), init_mask, pools) == [true, false, true, true, false]
        @test calc_mask(LabeledSplitStrat(), init_mask, pools) == [true, true, false, false, false]
        @test calc_mask(LabeledInlierSplitStrat(), init_mask, pools) == [true, false, false, false, false]
        @test calc_mask(LabeledOutlierSplitStrat(), init_mask, pools) == [false, true, false, false, false]
    end

    @testset "get_local_idx" begin
        pools = fill(:U, 6)
        train = BitArray([true, true, false, false, true, true])
        ds = DataSplits(train, .~train, FullSplitStrat())

        @test OneClassActiveLearning.get_local_idx(1, ds, pools, Val(:train)) == 1
        @test OneClassActiveLearning.get_local_idx(5, ds, pools, Val(:train)) == 3
        @test OneClassActiveLearning.get_local_idx(6, ds, pools, Val(:train)) == 4
        @test_throws ErrorException OneClassActiveLearning.get_local_idx(3, ds, pools, Val(:train))

        @test OneClassActiveLearning.get_local_idx(3, ds, pools, Val(:test)) == 1
        @test_throws ErrorException OneClassActiveLearning.get_local_idx(5, ds, pools, Val(:test))

        pools[3:6] .= :Lin
        ds = DataSplits(train, .~train, LabeledInlierSplitStrat())

        @test_throws ErrorException OneClassActiveLearning.get_local_idx(1, ds, pools, Val(:train))
        @test OneClassActiveLearning.get_local_idx(5, ds, pools, Val(:train)) == 1
    end

    @testset "filter_query_id" begin
        train = vcat(trues(5), falses(5))
        query_ids = [1,3,5]
        query_labels = [:Lin, :Lout, :U]
        @test OneClassActiveLearning.filter_query_id(query_ids, DataSplits(train, .~train, FullSplitStrat()), query_labels, Val(:train)) == [1,3,5]
        @test OneClassActiveLearning.filter_query_id(query_ids, DataSplits(train, .~train, LabeledInlierSplitStrat()), query_labels, Val(:train)) == [1]
        @test OneClassActiveLearning.filter_query_id(query_ids, DataSplits(train, .~train, LabeledOutlierSplitStrat()), query_labels, Val(:train)) == [3]
        @test OneClassActiveLearning.filter_query_id(query_ids, DataSplits(train, .~train, UnlabeledAndLabeledInlierSplitStrat()), query_labels, Val(:train)) == [1,5]
    end

    @testset "no holdout" begin
        data = reshape(collect(1:10), (2, 5))
        pools = fill(:U, 5)
        pools[1] = :Lin
        pools[2] = :Lout
        train = trues(5)
        ds = DataSplits(train, UnlabeledSplitStrat())
        train_data, train_pools, train_idx = get_train(ds, data, pools)
        @test size(train_data, 2) == length(train_pools) == length(train_idx)
        @test train_data == data[:, pools .== :U]
        @test all(train_pools .== :U)
        @test train_idx == [3, 4, 5]
    end

    @testset "with holdout" begin
        data = zeros(2, 5)
        pools = fill(:U, 5)
        train = BitArray([true, true, true, false, false])
        test = BitArray([false, false, false, true, true])
        @assert all(train .⊻ test)
        ds = DataSplits(train, test, FullSplitStrat())
        _, _, train_idx = get_train(ds, data, pools)
        _, _, test_idx = get_test(ds, data, pools)
        _, _, query_idx = get_query(ds, data, pools)
        @test train_idx == [1, 2, 3]
        @test test_idx == [4, 5]
        @test query_idx == [1, 2, 3]
    end
end

@testset "splits and init pools" begin
    data = zeros(2, 100)
    labels = fill(:inlier, 100)
    labels[1:10] .= :outlier
    @test_throws ArgumentError OneClassActiveLearning.get_splits_and_init_pools(data, labels, "foo", "Pu")
    @test_throws ArgumentError OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "foo")
    @testset "full split strategy" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pu")
        @test ds.train == ds.test == trues(size(data, 2))
        @test ip == fill(:U, size(data, 2))
        for p in [0.1, 0.2]
            ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pp"; p=p)
            @test ds.train == ds.test == trues(size(data, 2))
            @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == round(Int, size(data, 2) * p)
        end
        for n in [10, 20, 30]
            ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pn"; n=n)
            @test ds.train == ds.test == trues(size(data, 2))
            @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == n
        end
        for n in [11, 22, 33]
            ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pnin"; n=n)
            @test ds.train == ds.test == trues(size(data, 2))
            @test length(findall(x -> x == :Lin, ip)) == n
        end
        for x in [5, 10, 20]
            ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pa"; x=x)
            @test ds.train == ds.test == trues(size(data, 2))
            @test length(findall(x -> x == :Lin, ip)) == size(data, 1) + x
            ds, ip = OneClassActiveLearning.get_splits_and_init_pools(zeros(150, 100), labels, "Sf", "Pa"; x=x)
            @test ds.train == ds.test == trues(size(data, 2))
            @test length(findall(x -> x == :Lin, ip)) == length(findall(x -> x == :inlier, labels))
        end
    end
    @testset "holdout" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pu")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pu"; holdout_p=0.1)
        @test (sum(ds.train) == 90) && (sum(ds.test) == 10) && all(ds.train .⊻ ds.test)
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pp")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip[ds.train])) == sum(ds.train) * 0.1
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pn")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip[ds.train])) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pnin")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> x == :Lin, ip)) == 20
        @test length(findall(x -> x == :Lin, ip[ds.train])) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pa")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> x == :Lin, ip[ds.train])) == size(data, 1) + 10
    end
    @testset "inff" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pu")
        @test ds.train == ds.test == trues(size(data, 2))
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pp")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == sum(ds.train) * 0.1
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pn")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pnin")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pa")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == size(data, 1) + 10
    end
    @testset "lff" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sl", "Pu")
        @test ds.train == ds.test == trues(size(data, 2))
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sl", "Pp")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == sum(ds.train) * 0.1
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sl", "Pn")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sl", "Pnin")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sl", "Pa")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == size(data, 1) + 10
    end
end


text:
# * Relativistic orbital

"""
    κ2ℓ(κ::Integer) -> Integer

Calculate the `ℓ` quantum number corresponding to the `κ` quantum number.

Note: `κ` and `ℓ` values are always integers.
"""
function κ2ℓ(κ::Integer)
    κ == zero(κ) && throw(ArgumentError("κ can not be zero"))
    (κ < 0) ? -(κ + 1) : κ
end

"""
    κ2j(κ::Integer) -> HalfInteger

Calculate the `j` quantum number corresponding to the `κ` quantum number.

Note: `κ` is always an integer, but `j` will be a half-integer value.
"""
function κ2j(kappa::Integer)
    kappa == zero(kappa) && throw(ArgumentError("κ can not be zero"))
    half(2*abs(kappa) - 1)
end

"""
    ℓj2κ(ℓ::Integer, j::Real) -> Integer

Converts a valid `(ℓ, j)` pair to the corresponding `κ` value.

Note: there is a one-to-one correspondence between valid `(ℓ,j)` pairs and `κ` values
such that for `j = ℓ ± 1/2`, `κ = ∓(j + 1/2)`.
"""
function ℓj2κ(ℓ::Integer, j::Real)
    assert_orbital_ℓj(ℓ, j)
    (j < ℓ) ? ℓ : -(ℓ + 1)
end

function assert_orbital_ℓj(ℓ::Integer, j::Real)
    j = HalfInteger(j)
    s = half(1)
    (ℓ == j + s) || (ℓ == j - s) ||
        throw(ArgumentError("Invalid (ℓ, j) = $(ℓ), $(j) pair, expected j = ℓ ± 1/2."))
    return
end

"""
    struct RelativisticOrbital{N <: AtomicLevels.MQ} <: AbstractOrbital

Label for an atomic orbital with a principal quantum number `n::N` and well-defined total
angular momentum ``j``. The angular component of the orbital is labelled by the ``(\\ell, j)``
pair, conventionally written as ``\\ell_j`` (e.g. ``p_{3/2}``).

The ``\\ell`` and ``j`` can not be arbitrary, but must satisfy ``j = \\ell \\pm 1/2``.
Internally, the ``\\kappa`` quantum number, which is a unique integer corresponding to every
physical ``(\\ell, j)`` pair, is used to label each allowed pair.
When ``j = \\ell \\pm 1/2``, the corresponding ``\\kappa = \\mp(j + 1/2)``.

When printing and parsing `RelativisticOrbital`s, the notation `nℓ` and `nℓ-` is used (e.g.
`2p` and `2p-`), corresponding to the orbitals with ``j = \\ell + 1/2`` and
``j = \\ell - 1/2``, respectively.

The type parameter `N` has to be such that it can represent a proper principal quantum number
(i.e. a subtype of [`AtomicLevels.MQ`](@ref)).

# Properties

The following properties are part of the public API:

* `.n :: N` -- principal quantum number ``n``
* `.κ :: Int` -- ``\\kappa`` quantum number
* `.ℓ :: Int` -- the orbital angular momentum label ``\\ell``
* `.j :: HalfInteger` -- total angular momentum ``j``

```jldoctest
julia> orb = ro"5g-"
5g-

julia> orb.n
5

julia> orb.j
7/2

julia> orb.ℓ
4
```

# Constructors

    RelativisticOrbital(n::Integer, κ::Integer)
    RelativisticOrbital(n::Symbol, κ::Integer)
    RelativisticOrbital(n, ℓ::Integer, j::Real)

Construct an orbital label with the quantum numbers `n` and `κ`.
If the principal quantum number `n` is an integer, it has to positive and the orbital angular
momentum must satisfy `0 <= ℓ < n`.
Instead of `κ`, valid `ℓ` and `j` values can also be specified instead.

```jldoctest
julia> RelativisticOrbital(1, 0, 1//2)
1s

julia> RelativisticOrbital(2, -1)
2s

julia> RelativisticOrbital(:K, 2, 3//2)
Kd-
```
"""
struct RelativisticOrbital{N<:MQ} <: AbstractOrbital
    n::N
    κ::Int
    function RelativisticOrbital(n::Integer, κ::Integer)
        n ≥ 1 || throw(ArgumentError("Invalid principal quantum number $(n)"))
        κ == zero(κ) && throw(ArgumentError("κ can not be zero"))
        ℓ = κ2ℓ(κ)
        0 ≤ ℓ && ℓ < n || throw(ArgumentError("Angular quantum number has to be ∈ [0,$(n-1)] when n = $(n)"))
        new{Int}(n, κ)
    end
    function RelativisticOrbital(n::Symbol, κ::Integer)
        κ == zero(κ) && throw(ArgumentError("κ can not be zero"))
        new{Symbol}(n, κ)
    end
end
RelativisticOrbital(n::MQ, ℓ::Integer, j::Real) = RelativisticOrbital(n, ℓj2κ(ℓ, j))

Base.:(==)(a::RelativisticOrbital, b::RelativisticOrbital) =
    a.n == b.n && a.κ == b.κ

Base.hash(o::RelativisticOrbital, h::UInt) = hash(o.n, hash(o.κ, h))

"""
    mqtype(::RelativisticOrbital{MQ}) = MQ

Returns the main quantum number type of a [`RelativisticOrbital`](@ref).
"""
mqtype(::RelativisticOrbital{MQ}) where MQ = MQ

Base.propertynames(::RelativisticOrbital) = (fieldnames(RelativisticOrbital)..., :j, :ℓ)
function Base.getproperty(o::RelativisticOrbital, s::Symbol)
    s === :j ? κ2j(o.κ) :
    s === :ℓ ? κ2ℓ(o.κ) :
    getfield(o, s)
end

function Base.show(io::IO, orb::RelativisticOrbital)
    write(io, "$(orb.n)$(spectroscopic_label(κ2ℓ(orb.κ)))")
    orb.κ > 0 && write(io, "-")
end

function flip_j(orb::RelativisticOrbital)
    orb.κ == -1 && return RelativisticOrbital(orb.n, -1) # nothing to flip for s-orbitals
    RelativisticOrbital(orb.n, orb.κ < 0 ? abs(orb.κ) - 1 : -(orb.κ + 1))
end

degeneracy(orb::RelativisticOrbital{N}) where N = 2*abs(orb.κ) # 2j + 1 = 2|κ|

function Base.isless(a::RelativisticOrbital, b::RelativisticOrbital)
    nisless(a.n, b.n) && return true
    aℓ, bℓ = κ2ℓ(a.κ), κ2ℓ(b.κ)
    a.n == b.n && aℓ < bℓ && return true
    a.n == b.n && aℓ == bℓ && abs(a.κ) < abs(b.κ) && return true
    false
end

parity(orb::RelativisticOrbital) = p"odd"^κ2ℓ(orb.κ)
symmetry(orb::RelativisticOrbital) = orb.κ

isbound(::RelativisticOrbital{Int}) = true
isbound(::RelativisticOrbital{Symbol}) = false

"""
    angular_momenta(orbital)

Returns the angular momentum quantum numbers of `orbital`.

# Examples
```jldoctest
julia> angular_momenta(ro"2p-")
(1/2,)

julia> angular_momenta(ro"3d")
(5/2,)
```
"""
angular_momenta(orbital::RelativisticOrbital) = (orbital.j,)
angular_momentum_labels(::RelativisticOrbital) = ("j",)

# ** Saving/loading

Base.write(io::IO, o::RelativisticOrbital{Int}) = write(io, 'i', o.n, o.κ)
Base.write(io::IO, o::RelativisticOrbital{Symbol}) = write(io, 's', sizeof(o.n), o.n, o.κ)

function Base.read(io::IO, ::Type{RelativisticOrbital})
    kind = read(io, Char)
    n = if kind == 'i'
        read(io, Int)
    elseif kind == 's'
        b = Vector{UInt8}(undef, read(io, Int))
        readbytes!(io, b)
        Symbol(b)
    else
        error("Unknown RelativisticOrbital type $(kind)")
    end
    κ = read(io, Int)
    RelativisticOrbital(n, κ)
end

# * Orbital construction from strings

function Base.parse(::Type{<:RelativisticOrbital}, orb_str)
    m = match(r"^([0-9]+|.)([a-z]|\[[0-9]+\])([-]{0,1})$", orb_str)
    isnothing(m) && throw(ArgumentError("Invalid orbital string: $(orb_str)"))
    n = parse_orbital_n(m)
    ℓ = parse_orbital_ℓ(m)
    j = ℓ + (m[3] == "-" ? -1 : 1)*1//2
    RelativisticOrbital(n, ℓ, j)
end

"""
    @ro_str -> RelativisticOrbital

A string macro to construct an [`RelativisticOrbital`](@ref) from the canonical string
representation.

```jldoctest
julia> ro"1s"
1s

julia> ro"2p-"
2p-

julia> ro"Kf-"
Kf-
```
"""
macro ro_str(orb_str)
    parse(RelativisticOrbital, orb_str)
end

"""
    @ros_str -> Vector{RelativisticOrbital}

Can be used to easily construct a list of [`RelativisticOrbital`](@ref)s.

# Examples

```jldoctest
julia> ros"2[s-p] 3[p] k[0-d]"
10-element Array{RelativisticOrbital,1}:
 2s
 2p-
 2p
 3p-
 3p
 ks
 kp-
 kp
 kd-
 kd
```
"""
macro ros_str(orbs_str)
    orbitals_from_string(RelativisticOrbital, orbs_str)
end

"""
    str2κ(s) -> Int

A function to convert the canonical string representation of a ``\\ell_j`` angular label
(i.e. `ℓ-` or `ℓ`) into the corresponding ``\\kappa`` quantum number.

```jldoctest
julia> str2κ.(["s", "p-", "p"])
3-element Array{Int64,1}:
 -1
  1
 -2
```
"""
function str2κ(κ_str)
    m = match(r"^([a-z]|\[[0-9]+\])([-]{0,1})$", κ_str)
    m === nothing && throw(ArgumentError("Invalid κ string: $(κ_str)"))
    ℓ = parse_orbital_ℓ(m, 1)
    j = ℓ + half(m[2] == "-" ? -1 : 1)
    ℓj2κ(ℓ, j)
end

"""
    @κ_str -> Int

A string macro to convert the canonical string representation of a ``\\ell_j`` angular label
(i.e. `ℓ-` or `ℓ`) into the corresponding ``\\kappa`` quantum number.

```jldoctest
julia> κ"s", κ"p-", κ"p"
(-1, 1, -2)
```
"""
macro κ_str(κ_str)
    str2κ(κ_str)
end

"""
    nonrelorbital(o)

Return the non-relativistic orbital corresponding to `o`.

# Examples

```jldoctest
julia> nonrelorbital(o"2p")
2p

julia> nonrelorbital(ro"2p-")
2p
```
"""
nonrelorbital(o::Orbital) = o
nonrelorbital(o::RelativisticOrbital) = Orbital(o.n, o.ℓ)

export RelativisticOrbital, @ro_str, @ros_str, @κ_str, nonrelorbital, str2κ, κ2ℓ, κ2j, ℓj2κ


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

mutable struct StorageProfile <: SwaggerModel
    imageReference::Any # spec type: Union{ Nothing, ImageReference } # spec name: imageReference
    osDisk::Any # spec type: Union{ Nothing, OSDisk } # spec name: osDisk
    dataDisks::Any # spec type: Union{ Nothing, Vector{DataDisk} } # spec name: dataDisks

    function StorageProfile(;imageReference=nothing, osDisk=nothing, dataDisks=nothing)
        o = new()
        validate_property(StorageProfile, Symbol("imageReference"), imageReference)
        setfield!(o, Symbol("imageReference"), imageReference)
        validate_property(StorageProfile, Symbol("osDisk"), osDisk)
        setfield!(o, Symbol("osDisk"), osDisk)
        validate_property(StorageProfile, Symbol("dataDisks"), dataDisks)
        setfield!(o, Symbol("dataDisks"), dataDisks)
        o
    end
end # type StorageProfile

const _property_map_StorageProfile = Dict{Symbol,Symbol}(Symbol("imageReference")=>Symbol("imageReference"), Symbol("osDisk")=>Symbol("osDisk"), Symbol("dataDisks")=>Symbol("dataDisks"))
const _property_types_StorageProfile = Dict{Symbol,String}(Symbol("imageReference")=>"ImageReference", Symbol("osDisk")=>"OSDisk", Symbol("dataDisks")=>"Vector{DataDisk}")
Base.propertynames(::Type{ StorageProfile }) = collect(keys(_property_map_StorageProfile))
Swagger.property_type(::Type{ StorageProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_StorageProfile[name]))}
Swagger.field_name(::Type{ StorageProfile }, property_name::Symbol) =  _property_map_StorageProfile[property_name]

function check_required(o::StorageProfile)
    true
end

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


text:
# ContVar structure simulation with horizontal transfer
export cont_var_simulation, fitness
using Statistics
using DataStructures
try    # "using Printf" fails in julia v6, but avoids a deprecation error in julia v7
  using Printf
catch
end
#using DataStructures
import DataStructures.counter
import DataStructures.Accumulator
#=
Recommended command line to run:
>  julia run.jl configs/example1
=#

@doc """ function cont_var_simulation()
  Wright-Fisher model simulation (as opposed to Moran model)
  Parameters:
    N     MetaPopulation size
    m     number of subpopulations   # for now, subpopulation size = N/m
    ngens number of generations after burn in
    num_attributes   number of quantitative attributes of a variant
    variant_table Keeps track fitnesses and variant parent and innovation ancestor
"""
function cont_var_simulation( simrecord::cont_var_result_type; bins_to_vect::Bool=:false, bins_cutoff::Float64=1.0/simrecord.N )
  fit_diff_counter = counter(Int64)
  variant_table = Dict{Int64,variant_type}()
  #int_burn_in = Int(round(simrecord.burn_in*simrecord.N+50.0))  # reduce for testing
  id = Int[1]
  n = Int(floor(simrecord.N/simrecord.num_subpops))    # size of subpopulations
  #println("N: ",simrecord.N,"  mutation_stddev: ",simrecord.mutation_stddev,"  num_attributes: ",simrecord.num_attributes,"  int_burn_in: ",simrecord.int_burn_in)
  cumm_fitness_means = zeros(Float64,simrecord.num_subpops)
  cumm_fitness_medians = zeros(Float64,simrecord.num_subpops)
  cumm_fitness_coef_vars = zeros(Float64,simrecord.num_subpops)
  cumm_attr_means = [ zeros(Float64,simrecord.num_attributes) for i in 1:simrecord.num_subpops]
  cumm_attr_medians = [ zeros(Float64,simrecord.num_attributes) for i in 1:simrecord.num_subpops]
  cumm_attr_coef_vars = [ zeros(Float64,simrecord.num_attributes) for i in 1:simrecord.num_subpops]
  count_gens = 0
  subpops = PopList()
  for j = 1:simrecord.num_subpops
    Base.push!( subpops, Population() )
    for i = 1:n
      Base.push!( subpops[j], new_innovation( id, simrecord.ideal, simrecord.num_attributes, variant_table, simrecord.neutral ) )
    end
  end
  previous_variant_id = 1
  current_variant_id = id[1]
  previous_subpops = deepcopy(subpops)
  for g = 1:simrecord.ngens+simrecord.int_burn_in
    #println("before g: ",g,"  pop: ",subpops[1],"  pop attr: ",[ variant_table[subpops[1][i]].attributes[1] for i = 1:n ])
    after_burn_in = g > simrecord.int_burn_in
    previous_previous_variant_id = previous_variant_id
    previous_variant_id = current_variant_id
    current_variant_id = id[1]
    for j = 1:simrecord.num_subpops
      for i = 1:n
        subpops[j][i] = mutate_attributes( previous_subpops[j][i], id, variant_table, simrecord, after_burn_in, fit_diff_counter, bins_to_vect=bins_to_vect, bins_cutoff=bins_cutoff )
      end
      if simrecord.neutral
        # Equivalent to proportional selection with equal fitnesses
        r = rand(1:n,n)
        subpops[j] = subpops[j][ r ]
      else
        subpops[j] = propsel( subpops[j], n, variant_table )
      end
    end
    previous_subpops = deepcopy(subpops)
    if after_burn_in
      cumm_fitness_means += [ mean( [variant_table[v].fitness for v in s]) for s in subpops]
      cumm_fitness_medians += [ median( [variant_table[v].fitness for v in s]) for s in subpops]
      cumm_fitness_coef_vars += [ coef_var( [variant_table[v].fitness for v in s]) for s in subpops]
      # cumm_attr_means[s][i] is the mean of attribute i for subpop s, where the mean is over elements of s
      cumm_attr_means += [ [ mean( [ variant_table[v].attributes[i] for v in s]) for i =1:simrecord.num_attributes ] for s in subpops]
      # cumm_attr_means[s][i] is the median of attribute i for subpop s, where the median is over elements of s
      cumm_attr_medians += [ [ median( [ variant_table[v].attributes[i] for v in s]) for i =1:simrecord.num_attributes ] for s in subpops]
      # cumm_attr_coef_vars[s][i] is the coefficient of variation of attribute i for subpop s, where the mean is over elements of s
      cumm_attr_coef_vars += [ [ coef_var( [ variant_table[v].attributes[i] for v in s]) for i =1:simrecord.num_attributes ] for s in subpops]
      #println("fitness_mean: ", [ mean( [variant_table[v].fitness for v in s]) for s in subpops])
      #println("attr_coef_var: ",  [ [ coef_var( [ variant_table[v].attributes[i] for v in s]) for i =1:simrecord.num_attributes ] for s in subpops])
      #println("attr mean: ", [ [ mean( [ variant_table[v].attributes[i] for v in s]) for i =1:simrecord.num_attributes ] for s in subpops])
      count_gens += 1
    end
    clean_up_variant_table(previous_previous_variant_id,previous_variant_id,variant_table)
  end  # for g
  @assert count_gens == simrecord.ngens
  cumm_fitness_means /= simrecord.ngens
  cumm_fitness_medians /= simrecord.ngens
  cumm_fitness_coef_vars /= simrecord.ngens
  cumm_attr_means /= simrecord.ngens
  cumm_attr_medians /= simrecord.ngens
  cumm_attr_coef_vars /= simrecord.ngens
  # The next 3 means are over subpops
  simrecord.fitness_mean = mean(cumm_fitness_means)
  simrecord.fitness_median = mean(cumm_fitness_medians)
  simrecord.fitness_coef_var = mean(cumm_fitness_coef_vars)
  # The next 3 means are over attributes and subpops
  simrecord.attribute_mean = mean(mean(cumm_attr_means))
  simrecord.attribute_median = mean(mean(cumm_attr_medians))
  simrecord.attribute_coef_var = mean(mean(cumm_attr_coef_vars))
  if !bins_to_vect
    (simrecord.neg_count, simrecord.neg_neutral, simrecord.pos_neutral, simrecord.pos_count ) = summarize_bins( fit_diff_counter )
    return simrecord
  else
    (bins_vector,lower_bounds_vect) = bins_to_vector( fit_diff_counter, bins_cutoff )
    return (simrecord,bins_vector,lower_bounds_vect)
  end
end

function fitness( attributes::Vector{Float64}, ideal::Vector{Float64}, neutral::Bool )
  fit_slope = 1.0
  if length(attributes) != length(ideal)
    error("length(attributes) must equal length(ideal) in fitness")
  end
  if neutral
    return 1.0
  end
  dis = 0.0
  for k = 1:length(attributes)
    dis += abs( attributes[k] - ideal[k] )
  end
  if fit_slope == 0.0  # use the older linear method of computing fitness
    result = 1.0-dis/length(attributes)
    if result < 0.0
      #println("negative fitness")
      #println("fitness: attributes: ",attributes,"  ideal: ",ideal," fit: ",result)
      result = 0.0
    end
    @assert result >= 0.0
  else  # new method of computing fitness added on 11/14/17
    result = 1.0/(fit_slope*dis+1.0)
  end
  return result
end

@doc """ function new_innovation()
  Create a new variant_table record whose attributes are all set to ideal (specificied as an input parameter)
"""
function new_innovation( id::Vector{Int64}, ideal::Float64, num_attributes::Int64, variant_table::Dict{Int64,variant_type}, neutral::Bool )
  i = id[1]
  variant_table[i] = variant_type( 0.0, fill( ideal, num_attributes ) )
  variant_table[i].fitness = fitness( variant_table[i].attributes, fill( ideal, num_attributes), neutral )
  id[1] += 1
  i
end

@doc """ function mutate_attributes()
  Mutate the attributes corresponding the variant_table[v].
  The mutation is actually done by the function mutate().
"""
function mutate_attributes( v::Int64, id::Vector{Int64}, variant_table::Dict{Int64,variant_type},
    simrecord::cont_var_result_type, after_burn_in::Bool, fit_diff_counter::Accumulator{Int64,Int64};
    bins_to_vect::Bool=:false, bins_cutoff::Float64=1.0/simrecord.N )
  i = id[1]
  vt = variant_table[v]
  new_attributes = mutate( vt.attributes, simrecord.mutation_stddev )
  new_fit = fitness( new_attributes, fill( simrecord.ideal, simrecord.num_attributes), simrecord.neutral )
  #println("new attr: ",new_attributes,"  new_fit: ",new_fit)
  if after_burn_in
    if !bins_to_vect
      increment_bins( fit_diff_counter, new_fit-vt.fitness, 1.0/simrecord.N )
    else
      increment_bins( fit_diff_counter, new_fit-vt.fitness, bins_cutoff )
    end
  end
  variant_table[i] = deepcopy(vt)
  variant_table[i].fitness = new_fit
  variant_table[i].attributes = new_attributes
  id[1] += 1
  return i
end

@doc """ function mutate()
  Mutate the attributes array.
"""
function mutate( attributes::Vector{Float64}, mutation_stddev::Float64 )
  new_attributes = deepcopy(attributes)
      for i = 1:length(new_attributes)
        if new_attributes[i] <= 0.0
          println("neg attribute: ",new_attributes[i])
          new_attributes[i] = 1.0e-6
        end
        multiplier = (1.0+mutation_stddev*randn())
        #println("multiplier: ",multiplier)
        while multiplier <= 1.0e-6
          #println("neg multiplier")
          multiplier = (1.0+mutation_stddev*randn())
        end
        new_attributes[i] *= multiplier
        if new_attributes[i] < 0.0
          println("negative attribute with i=",i,": attribute: ",new_attributes[i])
        end
      end
    return new_attributes
end

function print_subpop( subpop::Vector{Int64}, variant_table::Dict{Int64,variant_type} )
  "[ $(join([ @sprintf(" %5d:%5.4f",vt,variant_table[vt].fitness) for vt in subpop ]))]"
end

function print_pop( stream::IO, subpops::PopList, variant_table::Dict{Int64,variant_type} )
  for sp in subpops
    print(stream,print_subpop(sp,variant_table))
  end
  println(stream)
end

# compute and save statistics about subpopulations and populations

function means( subpops::PopList, variant_table::Dict{Int64,variant_type} )
  fit(v) = variant_table[v].fitness
  means = [ mean(map(fit,s)) for s in subpops ]
  vars = [ var(map(fit,s)) for s in subpops ]
  return means, vars
end

function attr_means( subpops::PopList, variant_table::Dict{Int64,variant_type}, num_attributes::Int64 )
  ave_means = zeros(Float64,length(subpops))
  i = 1
  for s in subpops
    att_means = [ mean([variant_table[v].attributes[j] for v in s]) for j =1:num_attributes]
    ave_means[i] = mean(att_means)
    i += 1
  end
  #println("ave_means: ",ave_means)
  return ave_means
end

function attr_vars( subpops::PopList, variant_table::Dict{Int64,variant_type}, num_attributes::Int64 )
  ave_vars = zeros(Float64,length(subpops))
  i = 1
  for s in subpops
    att_vars = [ var([variant_table[v].attributes[j] for v in s]) for j =1:num_attributes]
    ave_vars[i] = mean(att_vars)
    i += 1
  end
  return ave_vars
end

function coef_var( lst )
  return std(lst)/mean(lst)
end

function clean_up_variant_table( previous_variant_id::Int64, previous_previous_variant_id::Int64,
    variant_table::Dict{Int64,variant_type} )
  #println("clean up:  ppv: ",previous_previous_variant_id,"  pv: ",previous_variant_id)
  for v = previous_variant_id:previous_previous_variant_id-1
    delete!(variant_table,v)
  end
end


text:
#Written By Johnathan Bizzano
exppow(x) = MathConstants.e ^ x
softmax(s, x) = x / s

"multinomial logistic regression and is often used as the last activation function of a neural network to
normalize the output of a network to a probability distribution over predicted output classes - Wikipedia"
function Softmax()
    return LayerGenerator(
        function(pb, sett, input_shape)
                LayerDesign(Sig(:SoftMaxLayer), input_shape, input_shape, (0), nothing,
                    function (ins, outs, cons)
                            return function ()
                                        broadcast!(exppow, ins, ins)
                                        broadcast!(softmax, outs, sum(ins), ins)
                                   end
                    end)
        end)
end

export SoftmaxLayer
SoftMaxLayer() = SoftMax()


text:
module OnionSockets

using Sockets
using DiffieHellman

import Serialization
import SecureIO
import Multiplexers

import Multiplexers: Line
import SecureIO: SecureSerializer

serialize(io::Union{TCPSocket,IOBuffer},msg) = Serialization.serialize(io,msg)
deserialize(io::Union{TCPSocket,IOBuffer}) = Serialization.deserialize(io)

serialize(io::Line,msg) = Multiplexers.serialize(io,msg)
deserialize(io::Line) = Multiplexers.deserialize(io)

serialize(io::SecureSerializer,msg) = SecureIO.serialize(io,msg)
deserialize(io::SecureSerializer) = SecureIO.deserialize(io)

BallotIOs = Union{TCPSocket,IOBuffer,Line,SecureSerializer}

Multiplexers.serialize(io::BallotIOs,msg) = serialize(io,msg)
Multiplexers.deserialize(io::BallotIOs) = deserialize(io)

SecureIO.serialize(io::BallotIOs,msg) = serialize(io,msg)
SecureIO.deserialize(io::BallotIOs) = deserialize(io)

import Sockets.connect
import Sockets.accept
import Sockets.TCPServer

function connect(socket::BallotIOs,id,dh::DH)
    send = x-> serialize(socket,x)
    get = () -> deserialize(socket)

    key,ballotid = diffie(send,get,dh)
    sroutersocket = SecureSerializer(socket,key)
    return sroutersocket
end

connect(port,id,dh::DH) = connect(connect(port),id,dh)

function accept(socket::BallotIOs,members,dh::DH)
    send = x -> serialize(socket,x)
    get = () -> deserialize(socket)

    key,unknownid = hellman(send,get,dh)
    #@assert unknownid==nothing

    securesocket = SecureSerializer(socket,key)
    return securesocket
end

accept(socket::BallotIOs,dh::DH) = accept(socket,nothing,dh)
accept(server::TCPServer,members,dh::DH) = accept(accept(server),members,dh)

export connect, accept, serialize, deserialize

end # module


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this = parse(Int64, this)
# num_workers, this = 10,1
println("Worker " * string(this) * " of " * string(num_workers))

include("../src/init.jl")
include("../src/fct.jl")
include("../src/fitting-fct.jl")
include("../src/worker_fct.jl")
data_type = "cal"
dataset = "v07.01"
datasets = [
    [53,64],
    [65,79],
    [80,92],
    [93,113]
]

channels = []
for ch in 0:1:36
    if channel_to_bege[ch] && ch != 6
        push!(channels, ch)
    end
end

channels = get_share_for_worker(channels, num_workers, this)
# channels = [0]
@info("")

for ch in channels
    ch_str = lpad(ch, 2, "0");
    plots_path = "../plots/data/"
    base_path  = "../../waveforms/calib/v07.01/"
    base_path_AE  = "../../waveforms/data/calib_AE/"

    calib_filepath = "../dicts/calib.json"
    calib = isfile(calib_filepath) ? JSON.parsefile(calib_filepath) : Dict()
    det_lib = JSON.parsefile("../dicts/det_lib.json")

    ch_str = lpad(ch, 2, "0");
    ds = det_lib[ctn[ch]]["ds"]

    high_cut = 4
    set = "data";
    if !haskey(calib, ctn[ch])
        calib[ctn[ch]] = Dict()
        calib[ctn[ch]]["data"] = Dict()
        calib[ctn[ch]]["sim"] = Dict()
    end

    filepath = joinpath(base_path_AE, ch_str * "-" * ctn[ch] * "-20MWA.h5")
    !isdir(dirname(filepath)) ? mkpath(dirname(filepath)) : ""
    if !isfile(filepath)

        files = glob(joinpath(base_path, ch_str * "-" * ctn[ch] * "/*.h5"));

        numofele = 5
        BackDelta5 = div(numofele,2)
        ForwardDelta5 = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1
        numofele = 201
        BackDelta201 = div(numofele,2)
        ForwardDelta201 = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1

        all_data = Table(energy = [], run = [], channel = [], AoEvetoed = [], datasetID = [], AoEclassifier = [], A = [], E = [])
        @showprogress 1 "Read data for Ch" * ch_str * "..." for file in files
            tmp = HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end
            tmp = tmp |> @filter(datasets[ds][1] <= _.run <= datasets[ds][2]) |> Table
            if channel_to_bege[ch] && ch != 6
                if 0 in unique(tmp.datasetID)
                    tmp = tmp |> @filter(_.datasetID == 0) |> Table
                    A = []
                    E = []
                    for wf in tmp.waveform
                        push!(A, maximum(movingaverage(diff(wf.value),5,BackDelta5,ForwardDelta5,3)))
                        push!(E, maximum(movingaverage(diff(wf.value),201,BackDelta201,ForwardDelta201,20)))
                    end
                    if size(tmp,1) > 1
                        append!(all_data, Table(energy = tmp.energy, run = tmp.run, channel = tmp.channel, AoEvetoed = tmp.AoEvetoed, datasetID = tmp.datasetID, AoEclassifier = tmp.AoEclassifier, A = A, E = E))
                    end
                end
            else
                A = []
                E = []
                for wf in tmp.waveform
                    push!(A, maximum(movingaverage(diff(wf.value),5,BackDelta5,ForwardDelta5,3)))
                    push!(E, maximum(movingaverage(diff(wf.value),201,BackDelta201,ForwardDelta201,20)))
                end
                if size(tmp,1) > 1
                    append!(all_data, Table(energy = tmp.energy, run = tmp.run, channel = tmp.channel, AoEvetoed = tmp.AoEvetoed, datasetID = tmp.datasetID, AoEclassifier = tmp.AoEclassifier, A = A, E = E))
                end
            end
        end

        HDF5.h5open(filepath, "w") do h5f
            LegendHDF5IO.writedata(h5f, "data", Table(
                A = float.(all_data.A),
                E = float.(all_data.E),
                run = Array{Int64, 1}(all_data.run),
                AoEvetoed = Array{Int64, 1}(all_data.AoEvetoed),
                datasetID = Array{Int64, 1}(all_data.datasetID),
                AoEclassifier = float.(all_data.AoEclassifier),
                energy = float.(all_data.energy)
            ))
        end
    end
end


text:
dosomething() = sleep(abs(randn()))

function runNsecondsworthofjobs(N)
    times = Vector{Float64}()
    totaltime = 0
    runcount = 0
    while totaltime < N
        t = @elapsed(dosomething())
        push!(times, t)
        totaltime += t
        runcount += 1
    end
    println("Ran job $runcount times, for total time of $totaltime seconds.")
    println("Average time per run was $(sum(times)/length(times)) seconds.")
    println("Individual times of the jobs in seconds were:")
    for t in times
        println("    $t")
    end
end

runNsecondsworthofjobs(5)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
dummyvalue(::Type{T}) where T<:AbstractFloat = T(NaN)
dummyvalue(::Type{T}) where T = typemax(T)

isdummy(val::T) where T = isequal(val, dummyvalue(T))

"""
    qnthresh(N)

Return the minimum number of points needed to specify the quasi-Newton quadratic model
in `N` dimensions.
"""
qnthresh(N) = ((N+1)*(N+2))÷2

"""
   xvert, fvert, qcoef = qfit(xm=>fm, x0=>f0, xp=>fp)

Given three points `xm < x0 < xp ` and three corresponding
values `fm`, `f0`, and `fp`, fit a quadratic. Returns the position `xvert` of the vertex,
the quadratic's value `fvert` at `xvert`, and the coefficient `qcoef` of the quadratic term.
`xvert` is a minimum if `qcoef > 0`.

Note that if the three points lie on a line, `qcoef == 0` and both `xvert` and `fvert` will
be infinite.
"""
function qfit(xfm, xf0, xfp)
    cm, c0, cp = lagrangecoefs(xfm, xf0, xfp)
    xm, fm = xfm
    x0, f0 = xf0
    xp, fp = xfp
    qvalue(x) = cm*(x-x0)*(x-xp) + c0*(x-xm)*(x-xp) + cp*(x-xm)*(x-x0)
    qcoef = cm+c0+cp
    if fm == f0 == fp
        return x0, f0, zero(qcoef) # when it's flat, use the middle point as the "vertex"
    end
    xvert = (cm*(x0+xp) + c0*(xm+xp) + cp*(xm+x0))/(2*qcoef)
    return xvert, qvalue(xvert), qcoef
end

@inline function lagrangecoefs(xfm, xf0, xfp)
    xm, fm = xfm.first, xfm.second
    x0, f0 = xf0.first, xf0.second
    xp, fp = xfp.first, xfp.second
    @assert(xp > x0 && x0 > xm && isfinite(xm) && isfinite(xp))
    cm = fm/((xm-x0)*(xm-xp))  # coefficients of Lagrange polynomial
    c0 = f0/((x0-xm)*(x0-xp))
    cp = fp/((xp-xm)*(xp-x0))
    cm, c0, cp
end

"""
    Δf = qdelta(box)

Return the difference `fmin - fbox`, where `fbox` is the value of `f` at the evaluation
point of `box` and `fmin` is the minimum of a one-dimensional quadratic fit of `f` (using
the data in `box.parent`) over the bounds of `box`.

Note that `Δf` might be -Inf, if `box` is unbounded and the quadratic estimate
is not convex.
"""
function qdelta(box::Box)
    xv, fv = box.parent.xvalues, box.parent.fvalues
    xm, x0, xp = xv[1], xv[2], xv[3]
    fm, f0, fp = fv[1], fv[2], fv[3]
    fbox = box.parent.fvalues[box.parent_cindex]
    cm, c0, cp = lagrangecoefs(xm=>fm, x0=>f0, xp=>fp)
    qvalue(x) = cm*(x-x0)*(x-xp) + c0*(x-xm)*(x-xp) + cp*(x-xm)*(x-x0)
    bb = boxbounds(box)
    qcoef = cm+c0+cp
    xvert = (cm*(x0+xp) + c0*(xm+xp) + cp*(xm+x0))/(2*qcoef)
    if qcoef > 0 && bb[1] <= xvert <= bb[2]
        # Convex and the vertex is inside the box
        return qvalue(xvert) - fbox
    end
    # Otherwise, the minimum is achieved at one of the edges
    # This needs careful evaluation in the case of infinite boxes to avoid Inf - Inf == NaN.
    if isinf(bb[1]) || isinf(bb[2])
        lcoef = -cm*(x0+xp) - c0*(xm+xp) - cp*(xm+x0)
        ccoef = cm*x0*xp + c0*xm*xp + cp*xm*x0
        if qcoef == 0  # the function is linear
            let lcoef = lcoef, ccoef = ccoef
                lvalue(x) = lcoef*x + ccoef
                return min(lvalue(bb[1]), lvalue(bb[2])) - fbox
            end
        else
            qvalue_inf(x) = isinf(x) ? abs(x)*sign(qcoef) : qvalue(x)
            return min(qvalue_inf(bb[1]), qvalue_inf(bb[2])) - fbox
        end
    end
    return min(qvalue(bb[1]), qvalue(bb[2])) - fbox
end

function qdelta(box::Box{T}, splitdim::Integer) where T
    p = find_parent_with_splitdim(box, splitdim)
    return p.parent.splitdim == splitdim ? qdelta(p) : zero(T)
end

function is_diag_convex(box)
    # Test whether we're likely to be in a convex patch. This only checks the
    # diagonals, because that's quick.
    isdiagconvex = true
    for i = 1:ndims(box)
        p = find_parent_with_splitdim(box, i)
        if isroot(p)
            isdiagconvex = false
        else
            xs, fs = p.parent.xvalues, p.parent.fvalues
            xvert, fvert, qcoef = qfit(xs[1]=>fs[1], xs[2]=>fs[2], xs[3]=>fs[3])
            isdiagconvex &= qcoef > 0
        end
        isdiagconvex || break
    end
    return isdiagconvex
end

## Minimum Edge List utilities
Base.empty!(mel::MELink) = (mel.next = mel; return mel)

function dropnext!(prev, next)
    if next != next.next
        # Drop the next item from the list
        next = next.next
        prev.next = next
    else
        # Drop the last item from the list
        prev.next = prev
        next = prev
    end
    return next
end

"""
    prev, next = trim!(mel::MELink, w, label=>fvalue)

Remove entries from `mel` that are worse than `(w, label=>fvalue)`.
Returns linked-list positions on either side of the putative new value.
"""
function trim!(mel::MELink, w, lf::Pair)
    l, f = lf.first, lf.second
    prev, next = mel, mel.next
    while prev != next && w > next.w
        if f <= next.f
            next = dropnext!(prev, next)
        else
            prev = next
            next = next.next
        end
    end
    if w == next.w && f < next.f
        next = dropnext!(prev, next)
    end
    return prev, next
end

function Base.insert!(mel::MELink, w, lf::Pair)
    l, f = lf
    prev, next = trim!(mel, w, lf)
    # Perform the insertion
    if prev == next
        # we're at the end of the list
        prev.next = typeof(mel)(l, w, f)
    else
        if f < next.f
            prev.next = typeof(mel)(l, w, f, next)
        end
    end
    return mel
end

function Base.iterate(mel::MELink)
    mel == mel.next && return nothing
    return (mel.next, mel.next)
end
function Base.iterate(mel::MELink, state::MELink)
    state == state.next && return nothing
    return (state.next, state.next)
end

function popfirst!(mel::MELink)
    item = mel.next
    mel.next = item.next == item ? mel : item.next
    return item
end

Base.length(mel::MELink) = count(x->true, mel)

function Base.show(io::IO, mel::MELink)
    print(io, "List(")
    next = mel.next
    while mel != next
        print(io, '(', next.w, ", ", next.l, "=>", next.f, "), ")
        mel = next
        next = next.next
    end
    print(io, ')')
end

### Box utilities
function Base.show(io::IO, box::Box)
    x = fill(NaN, ndims(box))
    position!(x, box)
    val = isroot(box) ? "Root" : value(box)
    print(io, "Box$val@", x)
end

function value(box::Box)
    isroot(box) && error("root box does not have a unique value")
    box.parent.fvalues[box.parent_cindex]
end
value_safe(box::Box{T}) where T = isroot(box) ? typemax(T) : value(box)

Base.isless(box1::Box, box2::Box) = isless(value_safe(box1), value_safe(box2))

function pick_other(xvalues, fvalues, idx)
    j = 1
    if j == idx j += 1 end
    xf1 = xvalues[j] => fvalues[j]
    j += 1
    if j == idx j += 1 end
    xf2 = xvalues[j] => fvalues[j]
    return xf1, xf2
end

function treeprint(io::IO, f::Function, root::Box)
    show(io, root)
    y = f(root)
    y != nothing && print(io, y)
    if !isleaf(root)
        print(io, '(')
        treeprint(io, f, root.children[1])
        print(io, ", ")
        treeprint(io, f, root.children[2])
        print(io, ", ")
        treeprint(io, f, root.children[3])
        print(io, ')')
    end
end
treeprint(io::IO, root::Box) = treeprint(io, x->nothing, root)

function add_children!(parent::Box{T}, splitdim, xvalues, fvalues, u::Real, v::Real) where T
    isleaf(parent) || error("cannot add children to non-leaf node")
    (length(xvalues) == 3 && xvalues[1] < xvalues[2] < xvalues[3]) || throw(ArgumentError("xvalues must be monotonic, got $xvalues"))
    parent.splitdim = splitdim
    p = find_parent_with_splitdim(parent, splitdim)
    if isroot(p)
        parent.minmax = (u, v)
    else
        parent.minmax = boxbounds(p)
    end
    for i = 1:3
        @assert(parent.minmax[1] <= xvalues[i] <= parent.minmax[2])
    end
    minsep = eps(T)*(xvalues[3]-xvalues[1])
    @assert(xvalues[2]-xvalues[1] > minsep)
    @assert(xvalues[3]-xvalues[2] > minsep)
    parent.xvalues = xvalues
    parent.fvalues = fvalues
    for i = 1:3
        Box(parent, i)  # creates the children of parent
    end
    parent
end

function cycle_free(box)
    p = parent(box)
    while !isroot(p)
        p == box && return false
        p = p.parent
    end
    return true
end

function isparent(parent, child)
    parent == child && return true
    while !isroot(child)
        child = child.parent
        parent == child && return true
    end
    return false
end

"""
    boxp = find_parent_with_splitdim(box, splitdim::Integer)

Return the first node at or above `box` who's parent box was split
along dimension `splitdim`. If `box` has not yet been split along
`splitdim`, returns the root box.
"""
function find_parent_with_splitdim(box::Box, splitdim::Integer)
    while !isroot(box)
        p = parent(box)
        if p.splitdim == splitdim
            return box
        end
        box = p
    end
    return box
end

"""
    box = greedy_smallest_child_leaf(root)

Walk the tree recursively, choosing the child with smallest function value at each stage.
`box` will be a leaf node.
"""
function greedy_smallest_child_leaf(box::Box)
    # Not guaranteed to be the smallest function value, it's the smallest that can be
    # reached stepwise
    while !isleaf(box)
        idx = argmin(box.fvalues)
        box = box.children[idx]
    end
    box
end

"""
    box = find_smallest_child_leaf(root)

Return the node below `root` with smallest function value.
"""
function find_smallest_child_leaf(box::Box)
    vmin, boxmin = value(box), box
    for leaf in leaves(box)
        vleaf = value(leaf)
        if vleaf <= vmin
            vmin, boxmin = vleaf, leaf
        end
    end
    return boxmin
end

function unique_smallest_leaves(boxes)
    uboxes = Set{eltype(boxes)}()
    for box in boxes
        push!(uboxes, find_smallest_child_leaf(box))
    end
    return uboxes
end

"""
    box = find_leaf_at(root, x)

Return the leaf-node `box` that contains `x`.
"""
function find_leaf_at(root::Box, x)
    isleaf(root) && return root
    while !isleaf(root)
        i = root.splitdim
        found = false
        for box in (root.children[1], root.children[2], root.children[3])
            bb = boxbounds(box)
            if bb[1] <= x[i] <= bb[2]
                root = box
                found = true
                break
            end
        end
        found || error("$(x[i]) not within $(root.minmax)")
    end
    root
end

"""
    box, success = find_leaf_at_edge(root, x, splitdim, dir)

Return the node `box` that contains `x` with an edge at `x[splitdim]`.
If `dir > 0`, a box to the right of the edge will be returned; if `dir < 0`, a box to the
left will be returned. `box` will be a leaf-node unless `success` is false;
the usual explanation for this is that the chosen edge is at the edge of the allowed
domain, and hence there isn't a node in that direction.

This is a useful utility for finding the neighbor of a given box. Example:

    # Let's find the neighbors of `box` along its parent's splitdim
    x = position(box, x0)
    i = box.parent.splitdim
    bb = boxbounds(box)
    # Right neighbor
    x[i] = bb[2]
    rnbr = find_leaf_at_edge(root, x, i, +1)
    # Left neighbor
    x[i] = bb[1]
    lnbr = find_leaf_at_edge(root, x, i, -1)
"""
function find_leaf_at_edge(root::Box, x, splitdim::Integer, dir::Signed)
    isleaf(root) && return root, false
    while !isleaf(root)
        i = root.splitdim
        found = false
        for box in (root.children[1], root.children[2], root.children[3])
            bb = boxbounds(box)
            if within(x[i], bb, dir)
                root = box
                found = true
                break
            end
        end
        found || return (root, false)
    end
    return (root, true)
end

"""
    x = position(box)
    x = position(box, x0)

Return the n-dimensional position vector `x` at which this box was evaluated
when it was a leaf. Some entries of `x` might be `NaN`, if `box` is sufficiently
near the root and not all dimensions have been split.
The variant supplying `x0` fills in those dimensions with the corresponding values
from `x0`.
"""
Base.position(box::Box) = position!(fill(NaN, ndims(box)), box)

function Base.position(box::Box, x0::AbstractVector)
    x = similar(x0)
    flag = falses(length(x0))
    position!(x, flag, box, x0)
end
function position!(x, flag, box::Box, x0::AbstractVector)
    copyto!(x, x0)
    position!(x, flag, box)
end
function position!(x, box::Box)
    flag = falses(length(x))
    position!(x, flag, box)
    return x
end
function position!(x, flag, box::Box)
    fill!(flag, false)
    nfilled = 0
    while !isroot(box) && nfilled < length(x)
        i = box.parent.splitdim
        if !flag[i]
            x[i] = box.parent.xvalues[box.parent_cindex]
            flag[i] = true
            nfilled += 1
        end
        box = box.parent
    end
    x
end
function default_position!(x, flag, xdefault)
    length(x) == length(flag) == length(xdefault) || throw(DimensionMismatch("all three inputs must have the same length"))
    for i = 1:length(x)
        if !flag[i]
            x[i] = xdefault[i]
        end
    end
    x
end

"""
    left, right = boxbounds(box)

Compute the bounds of `box` along the `splitdim` of `box`'s parent.
This throws an error for the root box.
"""
function boxbounds(box::Box)
    isroot(box) && error("cannot compute bounds on root Box")
    p = parent(box)
    if box.parent_cindex == 1
        return (p.minmax[1], (p.xvalues[1]+p.xvalues[2])/2)
    elseif box.parent_cindex == 2
        return ((p.xvalues[1]+p.xvalues[2])/2, (p.xvalues[2]+p.xvalues[3])/2)
    elseif box.parent_cindex == 3
        return ((p.xvalues[2]+p.xvalues[3])/2, p.minmax[2])
    end
    error("invalid parent_cindex $(box.parent_cindex)")
end

"""
    left, right = boxbounds(box, lower::Real, upper::Real)

Compute the bounds of `box` along the `splitdim` of `box`'s parent.
For the root box, returns `(lower, upper)`.
"""
function boxbounds(box::Box, lower::Real, upper::Real)
    isroot(box) && return (lower, upper)
    return boxbounds(box)
end

"""
    bb = boxbounds(box, lower::AbstractVector, upper::AbstractVector)

Compute the bounds of `box` along all dimensions.
"""
function boxbounds(box::Box{T}, lower::AbstractVector, upper::AbstractVector) where T
    length(lower) == length(upper) == ndims(box) || throw(DimensionMismatch("lower and upper must match dimensions of box"))
    bb = [(T(lower[i]), T(upper[i])) for i = 1:ndims(box)]
    boxbounds!(bb, box)
end

"""
    bb = boxbounds(box, splitdim::Integer, lower::AbstractVector, upper::AbstractVector)

Compute the bounds of `box` along dimension `splitdim`.
"""
function boxbounds(box::Box{T}, splitdim::Integer, lower::AbstractVector, upper::AbstractVector) where T
    p = find_parent_with_splitdim(box, splitdim)
    boxbounds(p, lower[splitdim], upper[splitdim])
end

function boxbounds!(bb, box::Box)
    flag = falses(ndims(box))
    boxbounds!(bb, flag, box)
    return bb
end
function boxbounds!(bb, flag, box::Box)
    fill!(flag, false)
    if isleaf(box)
        bb[box.parent.splitdim] = boxbounds(box)
        flag[box.parent.splitdim] = true
    else
        bb[box.splitdim] = box.minmax
        flag[box.splitdim] = true
    end
    nfilled = 1
    while !isroot(box) && nfilled < ndims(box)
        i = box.parent.splitdim
        if !flag[i]
            bb[i] = boxbounds(box)
            flag[i] = true
            nfilled += 1
        end
        box = box.parent
    end
    bb
end
function boxbounds!(bb, flag, box::Box, lower, upper)
    for i = 1:ndims(box)
        bb[i] = (lower[i], upper[i])
    end
    boxbounds!(bb, flag, box)
end

"""
    scale = boxscale(box, splits)

Return a vector containing a robust measure of the "scale" of `box` along
each coordinate axis. Specifically, it is typically related to the gap between
evaluation points of its parent boxes, falling back on the initial user-supplied `splits`
if it hasn't been split along a particular axis.

Note that a box that extends to infinity still has a finite `scale`. Moreover, a box
where one evaluation point is at the edge has a scale bigger than zero.
"""
function boxscale(box::Box{T,N}, splits) where {T,N}
    bxscale(s1, s2, s3) = s1 == s2 ? s3 - s2 :
                          s2 == s3 ? s2 - s1 :
                          min(s2-s1, s3-s2)
    bxscale(s) = bxscale(s[1], s[2], s[3])
    scale = Vector{T}(undef, N)
    for i = 1:N
        p = find_parent_with_splitdim(box, i)
        splitsi = splits[i]
        if isroot(p)
            scale[i] = bxscale(splitsi)
        else
            scale[i] = bxscale(p.parent.xvalues)
        end
    end
    return scale
end

function width(box::Box, splitdim::Integer, xdefault::Real, lower::Real, upper::Real)
    p = find_parent_with_splitdim(box, splitdim)
    bb = boxbounds(p, lower, upper)
    x = isroot(p) ? xdefault : p.parent.xvalues[p.parent_cindex]
    max(x-bb[1], bb[2]-x)
end
width(box::Box, splitdim::Integer, xdefault, lower, upper) =
    width(box, splitdim, xdefault[splitdim], lower[splitdim], upper[splitdim])

function isinside(x, lower, upper)
    ret = true
    for i = 1:length(x)
        ret &= lower[i] <= x[i] <= upper[i]
    end
    ret
end
function isinside(x, bb::Vector{Tuple{T,T}}) where T
    ret = true
    for i = 1:length(x)
        bbi = bb[i]
        ret &= bbi[1] <= x[i] <= bbi[2]
    end
    ret
end

"""
    within(x, (left, right), dir)

Return `true` if `x` lies between `left` and `right`. If `x` is on the edge,
`dir` must point towards the interior (positive if `x==left`, negative if `x==right`).
"""
function within(x::Real, bb::Tuple{Real,Real}, dir)
    if !(bb[1] <= x <= bb[2])
        return false
    end
    ((x == bb[1]) & (dir < 0)) && return false
    ((x == bb[2]) & (dir > 0)) && return false
    return true
end

function epswidth(bb::Tuple{T,T}) where T<:Real
    w1 = isfinite(bb[1]) ? eps(bb[1]) : T(0)
    w2 = isfinite(bb[2]) ? eps(bb[2]) : T(0)
    return 10*min(w1, w2)
end

function count_splits(box::Box)
    nsplits = Vector{Int}(undef, ndims(box))
    count_splits!(nsplits, box)
end

function count_splits!(nsplits, box::Box)
    fill!(nsplits, 0)
    if box.splitdim == 0
        box = box.parent
        box.splitdim == 0 && return nsplits  # an unsplit tree
    end
    nsplits[box.splitdim] += 1
    while !isroot(box)
        box = box.parent
        nsplits[box.splitdim] += 1
    end
    return nsplits
end

function Base.extrema(root::Box)
    isleaf(root) && error("tree is empty")
    minv, maxv = extrema(root.fvalues)
    for bx in root
        isleaf(bx) && continue
        mn, mx = extrema(bx.fvalues)
        minv = min(minv, mn)
        maxv = max(maxv, mx)
    end
    minv, maxv
end

## Utilities for experimenting with topology of the tree
"""
    splitprint([io::IO], box)

Print a representation of all boxes below `box`, using parentheses prefaced by a number `n`
to denote a split along dimension `n`, and `l` to represent a leaf.

See also [`parse`](@ref) for the inverse: converting a string representation to a tree of boxes.
"""
function splitprint(io::IO, box::Box)
    if isleaf(box)
        print(io, 'l')
    else
        print(io, box.splitdim, '(')
        splitprint(io, box.children[1])
        print(io, ", ")
        splitprint(io, box.children[2])
        print(io, ", ")
        splitprint(io, box.children[3])
        print(io, ')')
    end
end
splitprint(box::Box) = splitprint(stdout, box)

"""
    splitprint_colored([io::IO], box, innerbox)

Like [`splitprint`](@ref), except that `innerbox` is highlighted in red, and the chain
of parents of `innerbox` are highlighted in cyan.
"""
function splitprint_colored(io::IO, box::Box, thisbox::Box, allparents=get_allparents(thisbox))
    if isleaf(box)
        box == thisbox ? printstyled(io, 'l', color=:light_red) : print(io, 'l')
    else
        if box == thisbox
            printstyled(io, box.splitdim, color=:light_red)
        elseif box ∈ allparents
            printstyled(io, box.splitdim, color=:cyan)
        else
            print(io, box.splitdim)
        end
        print(io, '(')
        splitprint_colored(io, box.children[1], thisbox, allparents)
        print(io, ", ")
        splitprint_colored(io, box.children[2], thisbox, allparents)
        print(io, ", ")
        splitprint_colored(io, box.children[3], thisbox, allparents)
        print(io, ')')
    end
end
splitprint_colored(box::Box, thisbox::Box) = splitprint_colored(stdout, box, thisbox)

function get_allparents(box)
    allparents = Set{typeof(box)}()
    p = box
    while !isroot(p)
        p = parent(p)
        push!(allparents, p)
    end
    allparents
end

"""
    root = parse(Box{T,N}, string)

Parse a `string`, in the output format of [`splitprint`](@ref), and generate
a tree of boxes with that structure.
"""
function Base.parse(::Type{B}, str::AbstractString) where B<:Box
    b = B()
    splitbox!(b, str)
end

# splitbox! uses integer-valued positions and sets all function values to 0
function splitbox!(box::Box{T,N}, dim) where {T,N}
    x = position(box, zeros(N))
    xd = x[dim]
    add_children!(box, dim, [xd,xd+1,xd+2], zeros(3), -Inf, Inf)
    box
end

function splitbox!(box::Box, str::AbstractString)
    str == "l" && return
    m = match(r"([0-9]*)\((.*)\)", str)
    dim = parse(Int, m.captures[1])
    dimstr = m.captures[2]
    splitbox!(box, dim)
    commapos = [0,0]
    commaidx = 0
    open = 0
    i = firstindex(dimstr)
    while i <= ncodeunits(dimstr)
        c, i = iterate(dimstr, i)
        if c == '('
            open += 1
        elseif c == ')'
            open -= 1
        elseif c == ',' && open == 0
            commapos[commaidx+=1] = prevind(dimstr, i)
        end
    end
    splitbox!(box.children[1], strip(dimstr[1:prevind(dimstr, commapos[1])]))
    splitbox!(box.children[2], strip(dimstr[nextind(dimstr, commapos[1]):prevind(dimstr, commapos[2])]))
    splitbox!(box.children[3], strip(dimstr[nextind(dimstr, commapos[2]):end]))
    return box
end

## Tree traversal

"""
    root = get_root(box)

Return the root node for `box`.
"""
function get_root(box::Box)
    while !isroot(box)
        box = parent(box)
    end
    box
end

abstract type DepthFirstIterator end
Base.IteratorSize(::Type{<:DepthFirstIterator}) = Base.SizeUnknown()

struct DepthFirstLeafIterator{B<:Box} <: DepthFirstIterator
    root::B
end

function leaves(root::Box)
    DepthFirstLeafIterator(root)
end

function Base.iterate(iter::DepthFirstLeafIterator)
    state = find_next_leaf(iter, iter.root)
    state === nothing && return nothing
    return state, state
end
Base.iterate(root::Box) = root, root

function Base.iterate(iter::DepthFirstLeafIterator, state::Box)
    @assert(isleaf(state))
    state = find_next_leaf(iter, state)
    state === nothing && return nothing
    return (state, state)
end
function find_next_leaf(iter::DepthFirstLeafIterator, state::Box)
    ret = iterate(iter.root, state)
    ret === nothing && return nothing
    while !isleaf(ret[1])
        ret = iterate(iter.root, ret[2])
        ret === nothing && return nothing
    end
    return ret[2]
end

function Base.iterate(root::Box, state::Box)
    item = state  # old item
    if isleaf(item)
        box, i = up(item, root)
        if i <= length(box.children)
            item = box.children[i]
            return (item, item)
        end
        @assert(box == root)
        return nothing
    end
    item = item.children[1]
    return (item, item)
end

function up(box, root)
    local i
    box == root && return (box, length(box.children)+1)
    while true
        box, i = box.parent, box.parent_cindex+1
        box == root && return (box, i)
        i <= length(box.children) && break
    end
    return (box, i)
end

function Base.length(iter::DepthFirstLeafIterator)
    ret = iterate(iter)
    len = 0
    while ret !== nothing
        item, state = ret
        ret = iterate(iter, state)
        len += 1
    end
    return len
end

## Utilities for working with both mutable and immutable vectors
replacecoordinate!(x, i::Integer, val) = (x[i] = val; x)

replacecoordinate!(x::SVector{N,T}, i::Integer, val) where {N,T} =
    SVector{N,T}(_rpc(Tuple(x), i-1, T(val)))
@inline _rpc(t, i, val) = (ifelse(i == 0, val, t[1]), _rpc(Base.tail(t), i-1, val)...)
_rps(::Tuple{}, i, val) = ()

ipcopy!(dest, src) = copyto!(dest, src)
ipcopy!(dest::SVector, src) = src

## Other utilities
lohi(x, y) = x <= y ? (x, y) : (y, x)
function lohi(x, y, z)
    @assert(x <= y)
    z <= x && return z, x, y
    z <= y && return x, z, y
    return x, y, z
end

function order_pairs(xf1, xf2)
    x1, f1 = xf1
    x2, f2 = xf2
    return x1 <= x2 ? (xf1, xf2) : (xf2, xf1)
end

function order_pairs(xf1, xf2, xf3)
    xf1, xf2 = order_pairs(xf1, xf2)
    xf2, xf3 = order_pairs(xf2, xf3)
    xf1, xf2 = order_pairs(xf1, xf2)
    return xf1, xf2, xf3
end

function biggest_interval(a, b, c, d)
    ab, bc, cd = b-a, c-b, d-c
    if ab <= bc && ab <= cd
        return (a, b)
    elseif bc <= ab && bc <= cd
        return (b, c)
    end
    return (c, d)
end

function ensure_distinct(x::T, xref, bb::Tuple{Real,Real}; minfrac = 0.1) where T
    Δx = min(xref - bb[1], bb[2] - xref)
    if !isfinite(Δx)
        x != xref && return x
        return xref + 1
    end
    Δxmin = T(minfrac*Δx)
    if abs(x - xref) < Δxmin
        s = x == xref ? (bb[2] - xref > xref - bb[1] ? 1 : -1) : sign(x-xref)
        x = T(xref + Δxmin*s)
    end
    return x
end

function ensure_distinct(x::T, x1, x2, bb::Tuple{Real,Real}; minfrac = 0.1) where T
    x1, x2 = lohi(x1, x2)
    Δxmin = minfrac*min(x2-x1, bb[2] == x2 ? T(Inf) : bb[2]-x2, bb[1] == x1 ? T(Inf) : x1-bb[1])
    @assert(Δxmin > 0)
    if abs(x - x1) < Δxmin
        s = x == x1 ? 1 : sign(x-x1)
        x = T(max(bb[1], x1 + Δxmin*s))
    elseif abs(x - x2) < Δxmin
        s = x == x2 ? -1 : sign(x-x2)
        x = T(min(bb[2], x2 + Δxmin*s))
    end
    return x
end

"""
    t, exitdim = pathlength_box_exit(x0, dx, bb)

Given a ray `x0 + t*dx`, compute the value of `t` at which the ray exits the box
delimited by `bb` (a vector of `(lo, hi)` tuples). Also return the coordinate dimension
along which the exit occurs.
"""
function pathlength_box_exit(x0, dx, bb)
    t = oftype((bb[1][1] - x0[1])/dx[1], Inf)
    exitdim = 0
    for i = 1:length(x0)
        xi, dxi, bbi = x0[i], dx[i], bb[i]
        ti = (ifelse(dxi >= 0, bbi[2], bbi[1]) - xi)/dxi
        if ti < t
            t = ti
            exitdim = i
        end
    end
    t, exitdim
end

"""
    t, intersectdim = pathlength_hyperplane_intersect(x0, dx, xtarget, tmax)

Compute the maximum pathlength `t` (up to a value of `tmax`) at which the ray `x0 + t*dx`
intersects one of the hyperplanes specified by `x[i] = xtarget[i]` for any dimension `i`.
"""
function pathlength_hyperplane_intersect(x0, dx, xtarget, tmax)
    t = zero(typeof((xtarget[1] - x0[1])/dx[1]))
    intersectdim = 0
    for i = 1:length(x0)
        xi, dxi, xti = x0[i], dx[i], xtarget[i]
        ti = (xti - xi)/dxi
        if ti <= tmax && ti > t
            t = ti
            intersectdim = i
        end
    end
    t, intersectdim
end

# function different_basins(boxes::AbstractVector{B}, x0, lower, upper) where B<:Box
#     basinboxes = B[]
#     for box in boxes
#         ubasin = true
#         for bbox in basinboxes
#             if !is_different_basin(box, bbox, x0, lower, upper)
#                 ubasin = false
#                 break
#             end
#         end
#         if ubasin
#             push!(basinboxes, box)
#         end
#     end
#     return basinboxes
# end

# function is_different_basin(box1, box2, x0, lower, upper)
#     # This convexity test has its limits: we're comparing the maximum along the secant
#     # within the box to a function value calculated at a point that isn't along the secant.
#     # False positives happen, but this is better than false negatives.
#     root = get_root(box1)
#     v1, v2 = value(box1), value(box2)
#     x1, x2 = position(box1, x0), position(box2, x0)
#     dx = x2 - x1
#     bb = boxbounds(box1, lower, upper)
#     flag = Vector{Bool}(undef, length(lower))
#     leaf = box1
#     t, exitdim = pathlength_box_exit(x1, dx, bb)
#     # For consistency (e.g., commutivity with box1 and box2), we have to check
#     # exit condition even at first box, even though it's likely a false positive.
#     vmax = v1 + t*(v2-v1)
#     qdtot = oftype(vmax, 0)
#     # for i = 1:ndims(leaf)   # commented out to avoid false negatives
#     #     qdtot += qdelta(leaf, i)
#     # end
#     if value(leaf)+qdtot > vmax
#         return true
#     end
#     while t < 1
#         x2[:] .= x1 .+ t.*dx
#         x2[exitdim] = bb[exitdim][dx[exitdim] > 0 ? 2 : 1]  # avoid roundoff error in the critical coordinate
#         leaf_old = leaf
#         leaf, success = find_leaf_at_edge(root, x2, exitdim, dx[exitdim] > 0 ? +1 : -1)
#         success || break
#         boxbounds!(bb, flag, leaf, lower, upper)
#         tnext, exitdim = pathlength_box_exit(x1, dx, bb)
#         tnext = max(tnext, t)
#         while tnext == t  # must have hit a corner
#             tnext += oftype(t, 1e-4)
#             x2[:] .= x1 .+ tnext.*dx
#             bxtmp = find_leaf_at(root, x2)
#             boxbounds!(bb, flag, bxtmp, lower, upper)
#             tnext, exitdim = pathlength_box_exit(x1, dx, bb)
#         end
#         vmax = v1 + t*(v2-v1)
#         if tnext < 1
#             vmax = max(vmax, v1 + tnext*(v2-v1))
#         end
#         qdtot = oftype(vmax, 0)
#         # for i = 1:ndims(leaf)
#         #     qdtot += qdelta(leaf, i)
#         # end
#         if value(leaf)+qdtot > vmax
#             return true
#         end
#         t = tnext
#     end
#     return false
# end

"""
    a, b, c = pick3(a, b, (lower::Real, upper::Real))

Returns an ordered triple `a, b, c`, with two agreeing with the input `a` and `b`,
and the third point bisecting the largest interval between `a`, `b`, and the edges
`lower`, `upper`.
"""
function pick3(a, b, bb)
    a, b = lohi(a, b)
    imin, imax = biggest_interval(bb[1], a, b, bb[2])
    if isinf(imin)
        return a-2*(b-a), a, b
    elseif isinf(imax)
        return a, b, b+2*(b-a)
    end
    return a, b, c = lohi(a, b, (imin+imax)/2)
end

function issame(x1, x2, scale, rtol=sqrt(eps(eltype(x1))))
    same = true
    for i = 1:length(x1)
        same &= abs(x1[i] - x2[i]) < rtol*scale[i]
    end
    return same
end


text:
abstract type SpecialAngle{T<:Real} <: Real end
struct Piby2Radians{T<:Real} <: SpecialAngle{T} end
abstract type ScaledPi{T<:Real} <: SpecialAngle{T} end
struct ZeroRadians{T<:Real} <: ScaledPi{T} end
struct TwoNPiRadians{T} <: ScaledPi{T}
	n :: Int

	function TwoNPiRadians{T}(n::Int) where {T}
		new{T}(n)
	end
end
struct TwoNPlusOnePiRadians{T} <: ScaledPi{T}
	n :: Int

	function TwoNPlusOnePiRadians{T}(n::Int) where {T}
		new{T}(n)
	end
end

for DT in [:ZeroRadians, :Piby2Radians]
	@eval $DT() = $DT{Float64}()
	@eval $DT{T}(::$DT) where {T} = $DT{T}()
end

for DT in [:TwoNPiRadians, :TwoNPlusOnePiRadians]
	@eval $DT(n::Int) = $DT{Float64}(n)
	@eval $DT{T}(x::$DT) where {T} = $DT{T}(x.n)
end

const Zero = ZeroRadians()
const Pi = TwoNPlusOnePiRadians(0)
const TwoPi = TwoNPiRadians(1)
const ThreePi = TwoNPlusOnePiRadians(1)
const FourPi = TwoNPiRadians(2)
const Piby2 = Piby2Radians()

Base.:(-)(a::SpecialAngle) = -float(a)
Base.:(-)(a::ZeroRadians) = a
Base.:(-)(a::T) where {T<:TwoNPiRadians} = T(-a.n)
Base.:(-)(a::T) where {T<:TwoNPlusOnePiRadians} = T(-a.n - 1)

function Base.:(+)(::ZeroRadians{T1}, ::ZeroRadians{T2}) where {T1,T2}
	ZeroRadians{promote_type(T1,T2)}()
end
for DT in [:Piby2Radians]
	@eval function Base.:(+)(::ZeroRadians{T1}, ::$DT{T2}) where {T1,T2}
		$DT{promote_type(T1,T2)}()
	end
	@eval function Base.:(+)(::$DT{T2}, ::ZeroRadians{T1}) where {T1,T2}
		$DT{promote_type(T1,T2)}()
	end
end

for DT in [:TwoNPiRadians, :TwoNPlusOnePiRadians]
	@eval function Base.:(+)(::ZeroRadians{T1}, x::$DT{T2}) where {T1,T2}
		$DT{promote_type(T1,T2)}(x.n)
	end
end

function Base.:(+)(::Piby2Radians{T1}, ::Piby2Radians{T2}) where {T1,T2}
	TwoNPlusOnePiRadians{promote_type(T1,T2)}(0)
end
function Base.:(+)(a::TwoNPiRadians{T1}, b::TwoNPiRadians{T2}) where {T1,T2}
	TwoNPiRadians{promote_type(T1,T2)}(a.n + b.n)
end

function Base.:(+)(a::TwoNPiRadians{T1}, b::TwoNPlusOnePiRadians{T2}) where {T1,T2}
	TwoNPlusOnePiRadians{promote_type(T1,T2)}(a.n + b.n)
end
function Base.:(+)(a::TwoNPlusOnePiRadians{T2}, b::TwoNPiRadians{T1}) where {T1,T2}
	TwoNPlusOnePiRadians{promote_type(T1,T2)}(a.n + b.n)
end
function Base.:(+)(a::TwoNPlusOnePiRadians{T2}, b::TwoNPlusOnePiRadians{T1}) where {T1,T2}
	TwoNPiRadians{promote_type(T1,T2)}(a.n + b.n + 1)
end

# Returns exp(i α π/2) = cos(α π/2) + i*sin(α π/2)
# if α is a half-integer, we set α = m/2 for an odd integer m
# exp(i m π/4) = cos(m π/4) + i*sin(m π/4)
function cis_special(::Type{T}, α::HalfInt, p::Piby2Radians) where {T}
	if isinteger(α)
		αmod4 = mod(floor(Int,α),4)
		if αmod4 == 0
			res = Complex(one(T),zero(T))
		elseif αmod4 == 1
			res = Complex(zero(T),one(T))
		elseif αmod4 == 2
			res = Complex(-one(T),zero(T))
		elseif αmod4 == 3
			res = Complex(zero(T),-one(T))
		end
	else
		invsqrt2 = 1/√(T(2))
		m = mod(floor(Int,twice(α)),8)
		if m == 1
			res = Complex(invsqrt2, invsqrt2)
		elseif m == 3
			res = Complex(-invsqrt2, invsqrt2)
		elseif m == 5
			res = Complex(-invsqrt2, -invsqrt2)
		elseif m == 7
			res = Complex(invsqrt2, -invsqrt2)
		end
	end
	return res
end

cis_special(::Type{T}, α::HalfInt, ::ZeroRadians) where {T} = one(Complex{T})

# Returns exp(i α 2πn) = cos(2πn α) + i*sin(2πn α) = (-1)^(n * 2α)
function cis_special(::Type{T}, α::HalfInt, x::TwoNPiRadians) where {T}
	isodd(x.n * twice(α)) ? Complex(-one(T)) : Complex(one(T))
end

# Returns exp(i α (2n+1)π) = cos((2n+1)π α) + i*sin((2n+1)π α)
# if α is an integer, we obtain exp(i α (2n+1)π) = cos((2n+1)π α) = (-1)^((2n+1)α)
# if α is a half-integer = m/2 for odd m, we obtain
# exp(i α (2n+1)π) = i*sin((2n+1) * m * π/2)
function cis_special(::Type{T}, α::HalfInt, x::TwoNPlusOnePiRadians) where {T}
	r, i = zero(T), zero(T)
	if isinteger(α)
		if isodd(floor(Int,α))
			r = -one(T)
		else
			r = one(T)
		end
	else
		m = twice(α)
		p = (2x.n + 1)*m
		if mod(p, 4) == 1
			i = one(T)
		else
			i = -one(T)
		end
	end
	Complex{T}(r, i)
end

cis_special(α::Real, p::SpecialAngle{T}) where {T} = cis_special(T, HalfInt(α), p)

Base.cos(::ZeroRadians{T}) where {T} = one(T)
Base.cos(::TwoNPiRadians{T}) where {T} = one(T)
Base.cos(::TwoNPlusOnePiRadians{T}) where {T} = -one(T)
Base.cos(::Piby2Radians{T}) where {T} = zero(T)
Base.sin(::Piby2Radians{T}) where {T} = one(T)
Base.sin(::ScaledPi{T}) where {T} = zero(T)

for T in [:Float16, :Float32, :Float64, :BigFloat]
	@eval Base.$T(::Piby2Radians) = $T(π)/2
	@eval Base.$T(::ZeroRadians) = zero($T)
	@eval Base.$T(x::TwoNPiRadians) = $T(π)*2x.n
	@eval Base.$T(x::TwoNPlusOnePiRadians) = $T(π)*(2x.n + 1)
end
Base.AbstractFloat(p::SpecialAngle{T}) where {T} = T(p)

Base.one(::SpecialAngle{T}) where {T} = one(T)
Base.zero(::SpecialAngle{T}) where {T} = zero(T)

Base.promote_rule(::Type{<:SpecialAngle{T}}, R::Type{<:Real}) where {T} = promote_type(T,R)

Base.show(io::IO, x::ZeroRadians) = print(io, "0")
Base.show(io::IO, x::Piby2Radians) = print(io, "π/2")
function Base.show(io::IO, x::TwoNPiRadians)
	if iszero(x.n)
		print(io, 0)
	else
		print(io, 2x.n,"π")
	end
end
function Base.show(io::IO, x::TwoNPlusOnePiRadians)
	if iszero(x.n)
		print(io, "π")
	elseif x.n == -1
		print(io, "-π")
	else
		print(io, 2x.n + 1,"π")
	end
end