text:
#using Distributions
#include("constants.jl")
#include("orbit.jl")
#include("planet.jl")

#if !@isdefined PlanetarySystemAbstract
  @compat abstract type PlanetarySystemAbstract end

  struct PlanetarySystem{StarT<:StarAbstract} <: PlanetarySystemAbstract
    star::StarT
    planet::Vector{Planet}
    orbit::Vector{Orbit}

      # TODO DETAIL: Setup inner constructor to enforce equal number of planets & orbits
      function PlanetarySystem{StarT}(s::StarT, p::AbstractVector{Planet}, o::AbstractVector{Orbit}) where {StarT<:StarAbstract}
        @assert(length(p)==length(o)) # else error(string("Number of planets must match number of orbits: Np= ",length(p)," No= ",length(o)))
        new(s,p,o)
      end
  end

  const PlanetarySystemSingleStar = PlanetarySystem{SingleStar}

#end

function PlanetarySystem(s::StarT) where {StarT<:StarAbstract}
     PlanetarySystem(s,Vector{Planet}(undef,0),Vector{Orbit}(undef,0))  # Constructor for a Planetary System with no planets
end

function PlanetarySystem(s::StarT, p::Planet, o::Orbit) where {StarT<:StarAbstract}
   PlanetarySystem(s,[p],[o])  # Constructor for a single Planet System
end

function PlanetarySystem(s::StarT, p::AbstractVector{Planet}, o::AbstractVector{Orbit}) where {StarT<:StarAbstract}
   PlanetarySystem{StarT}(s,p,o)  # Constructor for a single Planet System
end

function PlanetarySystem(ps::PlanetarySystem{StarT}, keep::AbstractVector{Int64})  where {StarT<:StarAbstract} # Why doesn't this work?
   PlanetarySystem{StarT}(ps.star,ps.planet[keep],ps.orbit[keep])
end

function star( ps::PlanetarySystem{StarT} )::StarT where {StarT<:StarAbstract}
  return ps.star
end

function planets( ps::PlanetarySystem{StarT} )::Vector{Planet} where {StarT<:StarAbstract}
  return ps.planet
end

function orbits( ps::PlanetarySystem{StarT} )::Vector{Orbit} where {StarT<:StarAbstract}
  return ps.orbit
end

#function PlanetarySystemSingleStar(ps::PlanetarySystemSingleStar, keep::Vector{Int64})
#function PlanetarySystem{StarT<:StarAbstract}(ps::PlanetarySystem{StarT}, keep::Vector{Int64})
#   PlanetarySystem(ps.star,ps.planet[keep],ps.orbit[keep])
#end

flux(ps::PlanetarySystem{StarT}) where {StarT<:StarAbstract} = flux(star(ps))
#flux(ps::PlanetarySystem{Star}) = flux(ps.star)
#flux(ps::PlanetarySystem{BinaryStar}) = flux(ps.star)
#flux(ps::PlanetarySystem{MultipleStar}) = flux(ps.star)

function num_planets(s::PlanetarySystem{StarT}) where {StarT<:StarAbstract}
  @assert( length(planets(s)) == length(orbits(s)) )    # TODO OPT: Deactivate inner assert's like this for speed once tested
  return length(planets(s))
end

function calc_hill_sphere(a::Float64, mu::Float64)
    return a*(mu/3)^(1//3)
end

function calc_mutual_hill_radii(ps::PlanetarySystem{StarT}, pl1::Int64, pl2::Int64) where StarT <: StarAbstract
    mu = (ps.planet[pl1].mass + ps.planet[pl2].mass)/ps.star.mass
    a = 0.5*(ps.orbit[pl1].a + ps.orbit[pl2].a)
    return calc_hill_sphere(a, mu)
end

function test_stability_circular(P::AbstractVector{Float64}, mass::AbstractVector{Float64}, star_mass::Float64, sim_param::SimParam)
    @assert length(P) == length(mass)
    min_num_mutual_hill_radii = get_real(sim_param, "num_mutual_hill_radii")
    found_instability = false
    order = sortperm(P)
    a2 = semimajor_axis(P[order[1]], mass[order[1]]+star_mass)
    for pl in 1:(length(P)-1)
        a1 = a2   # semimajor_axis(P[order[pl]], mass[order[pl]]+star_mass)
        a2 = semimajor_axis(P[order[pl+1]], mass[order[pl+1]]+star_mass)
        a = 0.5*(a1+a2)
        mu = (mass[order[pl]] + mass[order[pl+1]])/star_mass
        mutual_hill_radius = calc_hill_sphere(a, mu)
        if a2-a1  < min_num_mutual_hill_radii*mutual_hill_radius
            found_instability = true
            break
        end
    end # loop over neighboring planet pairs within cluster
    return !found_instability
end

function test_stability(P::AbstractVector{Float64}, mass::AbstractVector{Float64}, star_mass::Float64, sim_param::SimParam; ecc::AbstractVector{Float64}=zeros(length(P)))
    @assert length(P) == length(mass) == length(ecc)
    min_num_mutual_hill_radii = get_real(sim_param, "num_mutual_hill_radii")
    found_instability = false
    order = sortperm(P)
    a2 = semimajor_axis(P[order[1]], mass[order[1]]+star_mass)
    for pl in 1:(length(P)-1)
        a1 = a2   # semimajor_axis(P[order[pl]], mass[order[pl]]+star_mass)
        a2 = semimajor_axis(P[order[pl+1]], mass[order[pl+1]]+star_mass)
        a = 0.5*(a1+a2)
        mu = (mass[order[pl]] + mass[order[pl+1]])/star_mass
        mutual_hill_radius = calc_hill_sphere(a, mu)
        e1 = ecc[order[pl]]
        e2 = ecc[order[pl+1]]
        if a2*(1-e2)-a1*(1+e1) < min_num_mutual_hill_radii*mutual_hill_radius
            found_instability = true
            break
        end
    end # loop over neighboring planet pairs within cluster
    return !found_instability
end

function is_period_ratio_near_resonance(period_ratio::Float64, sim_param::SimParam)
    resonance_width = get_real(sim_param, "resonance_width")
    resonance_width_factor = 1+resonance_width
    period_ratios_to_check = get_any(sim_param, "period_ratios_mmr", Array{Float64,1})
    result = false
    for period_ratio_mmr in period_ratios_to_check
        if period_ratio_mmr <= period_ratio <= period_ratio_mmr*resonance_width_factor
            result = true
            break
        end
    end
    return result
end

function calc_if_near_resonance(P::AbstractVector{Float64}, sim_param::SimParam)
    @assert issorted(P)   # TODO: OPT: Could remove once know it is safe
    result = falses(length(P))
    if length(P) >= 2
        for i in 1:(length(P)-1)
            if is_period_ratio_near_resonance(P[i+1]/P[i], sim_param)
                result[i] = true
                result[i+1] = true
            end # near mmr
        end # planets
    end # at least two planets
    return result
end

# Code to generate simple planetary systems.  For more sophisticated algorithms, see clustered model in He et al. 2019

function generate_planet_mass_from_radius_powerlaw(r::Float64, sim_param::SimParam)
  mr_power_index::Float64 = get_real(sim_param,"mr_power_index")
  mr_const::Float64 = get_real(sim_param,"mr_const")
  mr_max_mass::Float64 = get_real(sim_param,"mr_max_mass")
  m = mr_const*earth_mass*(r/earth_radius)^mr_power_index
  if m > mr_max_mass
     m = mr_max_mass
  end
  return m
end

function generate_planet_mass_from_radius_powerlaw(r::Float64, s::Star, o::Orbit, sim_param::SimParam)  # TODO USER SCI: This is importnat if you are using a stability criteria.  In that case, this should be replacde w/ better M-R relationship.  See Matthias's clustered  model for example.
  generate_planet_mass_from_radius_powerlaw(r,sim_param)
end

function generate_num_planets_poisson(lambda::Real, max_planets::Integer; min_planets::Integer = 0)
  ##### Note: this function produces odd behaviour if lambda < min_planets due to bugs in Distributions.Truncated(); for example, Distributions.Truncated(Distributions.Poisson(lambda),min_planets,max_planets) returns only values >=2 if min_planets=1 and lambda<~0.95
  if lambda < min_planets*1e-3
      return min_planets
  end
  bug_fixed = false # TODO OPT: true case should work, but Danley found bug in Distributions package.  Revert once fixed for speed.
  local n
  if bug_fixed
     d = Distributions.Truncated(Distributions.Poisson(lambda),min_planets,max_planets)
     n = rand(d)
  else
     if min_planets == 0
        min_planets = -1
     end
     d = Distributions.Truncated(Distributions.Poisson(lambda),min_planets,max_planets)
     n = rand(d)
     #=
     n = -1
     while !(min_planets<=n<=max_planets)
        n = rand(Distributions.Poisson(lambda))
     end
     =#
  end
  return n
end

function draw_truncated_poisson(lambda::Real; min::Integer=0, max::Integer=20, n::Integer=1)
    pmf = [(lambda^k*exp(-lambda))/factorial(k) for k in min:max]
    pmf ./= sum(pmf)
    cmf = zeros(Float64,length(pmf))
    for i in 1:length(pmf)
        for j in i:-1:1
            cmf[i] += pmf[j]
        end
    end
    result = Array{Int64}(undef,n)
    for i in 1:n
        u = rand()
        result[i] = findfirst(x-> x>=u,cmf) + min -1
    end
    return result
end

function generate_num_planets_poisson(s::Star, sim_param::SimParam)
  lambda::Float64 = exp(get_real(sim_param,"log_eta_pl"))
  max_tranets_in_sys::Int64 = get_int(sim_param,"max_tranets_in_sys")
  generate_num_planets_poisson(lambda,max_tranets_in_sys)
end

function generate_period_and_sizes_log_normal(s::Star, sim_param::SimParam; num_pl::Integer = 1)  # TODO USER SCI:  User should make sure planetary properties are being drawn appropriately for their scientific purposes
    mu_log_r::Float64 = get_real(sim_param,"mean_log_planet_radius")
    sigma_log_r::Float64 = get_real(sim_param,"sigma_log_planet_radius")
    mu_log_P::Float64 = get_real(sim_param,"mean_log_planet_period")
    sigma_log_P::Float64 = get_real(sim_param,"sigma_log_planet_period")
    min_period::Float64 = get_real(sim_param,"min_period")
    max_period::Float64 = get_real(sim_param,"max_period")
    min_radius::Float64 = get_real(sim_param,"min_radius")
    max_radius::Float64 = get_real(sim_param,"max_radius")
    max_draws::Int64 = 100

    if   sigma_log_r <= 0. || sigma_log_P<=0.
     println("# mu_log_r= ", mu_log_r, " sigma_log_r= ", sigma_log_r, " mu_log_P= ", mu_log_P, " sigma_log_P= ", sigma_log_P)
    end
    rdist = LogNormal(mu_log_r,sigma_log_r)
    Pdist = LogNormal(mu_log_P,sigma_log_P)
    #Rlist = rand(rdist,num_pl)
    #Plist = rand(Pdist,num_pl)
    #idx_keep = find(i->(min_radius<=Rlist[i]<=max_radius) && (min_period<=Plist[i]<=max_period), 1:num_pl )
    #return Plist[idx_keep], Rlist[idx_keep]  # replaced because want to return exactly num_pl.  Could use Truncated to restore above code.
    Rlist = zeros(num_pl)
    Plist = zeros(num_pl)
    for i in 1:num_pl
      j = 0
      while ! (min_radius<Rlist[i]<=max_radius) && j<max_draws
            Rlist[i] = rand(rdist)
            j+=1
      end
      if j>=max_draws
         println("# Struggled to draw size for: ",mu_log_r, " ", sigma_log_r)
      end
      j = 0
      while ! (min_period<Plist[i]<=max_period) && j<max_draws
            Plist[i] = rand(Pdist)
            j+=1
      end
      if j>=max_draws
         println("# Struggled to draw period for: ",mu_log_P, " ", sigma_log_P)
      end
    end
    return Plist, Rlist
end

function draw_power_law(n::Real, x0::Real, x1::Real, num_pl::Integer)
    if n != -1
        return ((x1^(n+1) - x0^(n+1)).*rand(num_pl) .+ x0^(n+1)).^(1/(n+1))
    else #if n == -1
        return x0*(x1/x0).^rand(num_pl)
    end
end

function draw_power_law!(out::Array, n::Real, x0::Real, x1::Real, num_pl::Integer)
    @assert length(out)==num_pl
    if n != -1
        return out .= ((x1^(n+1) - x0^(n+1)).*rand(num_pl) .+ x0^(n+1)).^(1/(n+1))
    else #if n == -1
        return out .= x0*(x1/x0).^rand(num_pl)
    end
end

function draw_broken_power_law(n1::Real, n2::Real, x0::Real, x1::Real, xb::Real, num_pl::Integer)
    #x0 and x1 are the lower and upper truncation limits, and xb is the break point, i.e. x0 <= xb <= x1 (all must be positive)
    #n1 and n2 are the power law indices between x0 and xb, and xb and x1, respectively (can be positive or negative)
    @assert(x0 <= xb <= x1)
    @assert(num_pl >= 1)

    u_draws = rand(num_pl) #'num_pl' draws from the uniform distribution between 0 and 1
    x_draws = zeros(num_pl)

    if (n1 != -1) & (n2 != -1)
        C1 = 1.0/(((xb^(n1+1) - x0^(n1+1))/(n1+1)) + ((xb^(n1-n2)*(x1^(n2+1) - xb^(n2+1)))/(n2+1))) #normalization constant
        ub = (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1) #break point in u, between 0 and 1
        for (i,u) in enumerate(u_draws)
            if u <= ub
                x_draws[i] = (((n1+1)*u)/C1 + x0^(n1+1))^(1/(n1+1))
            else #if u > ub
                x_draws[i] = (((n2+1)/(C1*xb^(n1-n2)))*(u - (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1)) + xb^(n2+1))^(1/(n2+1))
            end
        end
    elseif (n1 == -1) & (n2 != -1)
        C1 = 1.0/(log(xb/x0) + ((xb^(-1-n2))*(x1^(n2+1)) - 1)/(n2+1)) #normalization constant
        ub = C1*log(xb/x0) #break point in u, between 0 and 1
        for (i,u) in enumerate(u_draws)
            if u <= ub
                x_draws[i] = x0*exp(u/C1)
            else #if u > ub
                x_draws[i] = (((n2+1)/(C1*xb^(-1-n2)))*(u - C1*log(xb/x0)) + xb^(n2+1))^(1/(n2+1))
            end
        end
    elseif (n1 != -1) & (n2 == -1)
        C1 = 1.0/(((xb^(n1+1) - x0^(n1+1))/(n1+1)) + (xb^(n1+1))*log(x1/xb)) #normalization constant
        ub = (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1) #break point in u, between 0 and 1
        for (i,u) in enumerate(u_draws)
            if u <= ub
                x_draws[i] = (((n1+1)*u)/C1 + x0^(n1+1))^(1/(n1+1))
            else #if u > ub
                x_draws[i] = xb*exp((1/(C1*xb^(n1+1)))*(u - (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1)))
            end
        end
    else #if n1 == -1 and n2 == -1 (i.e. it is a single power-law with index of -1)
        for (i,u) in enumerate(u_draws)
            x_draws[i] = x0*exp(u*log(x1/x0))
        end
    end

    return x_draws
end

function generate_periods_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    power_law_P::Float64 = get_real(sim_param,"power_law_P")
    min_period::Float64 = get_real(sim_param,"min_period")
    max_period::Float64 = get_real(sim_param,"max_period")
    Plist = draw_power_law(power_law_P,min_period,max_period, num_pl)
    return Plist
end

function generate_sizes_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    power_law_r::Float64 = get_real(sim_param,"power_law_r")
    min_radius::Float64 = get_real(sim_param,"min_radius")
    max_radius::Float64 = get_real(sim_param,"max_radius")
    Rlist = draw_power_law(power_law_r,min_radius,max_radius, num_pl)
    return Rlist
end

function generate_sizes_broken_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    power_law_r1::Float64 = get_real(sim_param,"power_law_r1")
    power_law_r2::Float64 = get_real(sim_param,"power_law_r2")
    min_radius::Float64 = get_real(sim_param,"min_radius")
    max_radius::Float64 = get_real(sim_param,"max_radius")
    break_radius::Float64 = get_real(sim_param,"break_radius")
    Rlist = draw_broken_power_law(power_law_r1,power_law_r2,min_radius,max_radius,break_radius, num_pl)
    return Rlist
end

function generate_period_and_sizes_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    return (generate_periods_power_law(s, sim_param, num_pl=num_pl), generate_sizes_power_law(s, sim_param, num_pl=num_pl))
end

using SpecialFunctions

cdf_lognormal(x::Float64; μ::Float64=0., σ::Float64=1.) = 0.5*erfc(-(log(x) - μ)/(σ*sqrt(2.)))

function invert_cdf_lognormal(y::Float64; μ::Float64=0., σ::Float64=1.)
    @assert 0<=y<=1
    return exp(μ - sqrt(2.)*σ*erfcinv(2*y))
end

function cdf_power_law(x::Float64; x0::Float64, x1::Float64, α::Float64)
    @assert x0 <= x <= x1
    if α != -1
        return (x^(α+1) - x0^(α+1))/(x1^(α+1) - x0^(α+1))
    else # if α == -1
        return log(x/x0)/log(x1/x0)
    end
end

function invert_cdf_power_law(y::Float64; x0::Float64, x1::Float64, α::Float64)
    @assert x0 < x1
    if α != -1
        return (x0^(α+1) + (x1^(α+1) - x0^(α+1))*y)^(1/(α+1))
    else # if α == -1
        return x0*(x1/x0)^y
    end
end

"""
    draw_segmented_uniform(segments)

Draw a random uniform variable given segments in the interval (0,1).

# Arguments:
- `segments::Array{Tuple{Float64,Float64},1}`: array of segments (a,b) between 0 and 1.
NOTE: all segments must be between [0,1] and be non-overlapping!

# Returns:
A value (float) contained in one of the given segments, drawn uniformly.
"""
function draw_segmented_uniform(segments::Array{Tuple{Float64,Float64},1})
    # Note: "segments" must be non-overlapping
    n_segs = length(segments) # number of segments
    @assert n_segs > 0
    for (i,seg) in enumerate(segments)
        @assert 0 <= seg[1] < seg[2] <= 1 # check that all segments are valid
        if i != n_segs
            @assert segments[i][2] < segments[i+1][1] # check to make sure there are no overlapping segments
        end
    end

    seg_lengths = [seg[2]-seg[1] for seg in segments]
    sum_lengths = sum(seg_lengths)

    u_draw = rand()
    for (i,seg) in enumerate(segments)
        # Map the union of the segments uniformly to the domain [0,1]:
        sum_l = sum(seg_lengths[1:i-1])
        if sum_l <= u_draw*sum_lengths < sum_l + seg_lengths[i]
            return rand(Uniform(seg[1], seg[2]))
        end
    end
end

"""
    compute_unstable_regions_periods_given_planets(P, mass, insert_pl_mass, star_mass, sim_param; ecc=zeros(length(P)), insert_pl_ecc=0., use_mutualHill=true, verbose=false)

Compute the unstable intervals in orbital period given the existing planets and the mass and eccentricity of the planet we are trying to insert.

# Arguments:
- `P::AbstractVector{Float64}`: periods of the existing planets.
- `mass::AbstractVector{Float64}`: masses of the existing planets.
- `insert_pl_mass::Float64`: mass of the planet we are trying to insert.
- `star_mass::Float64`: mass of the star.
- `sim_param::SimParam`: a SimParam object containing various simulation parameters.
- `ecc::AbstractVector{Float64}=zeros(length(P))`: eccentricities of the existing planets (defaults to circular).
- `insert_pl_ecc::Float64=0.`: eccentricity of the planet we are trying to insert (defaults to circular).
- `use_mutualHill::Bool=true`: whether to use the minimum separation in terms of the mutual Hill radii (default) or the Hill radii (if false).
- `verbose::Bool=false`: whether to print various messages.
NOTE: the existing planets do not have to be sorted in period.

# Returns:
`P_segments_unstable::Array{Tuple{Float64,Float64},1}`: an array of segments. Inserting the planet somewhere inside one of these segments would deem the system unstable.

NOTE: these segments may be overlapping!
"""
function compute_unstable_regions_periods_given_planets(P::AbstractVector{Float64}, mass::AbstractVector{Float64}, insert_pl_mass::Float64, star_mass::Float64, sim_param::SimParam; ecc::AbstractVector{Float64}=zeros(length(P)), insert_pl_ecc::Float64=0., use_mutualHill::Bool=true, verbose::Bool=false)
    @assert length(P) == length(mass) == length(ecc)
    min_num_mutual_hill_radii = get_real(sim_param, "num_mutual_hill_radii")
    order = sortperm(P)

    P_segments_unstable = Array{Tuple{Float64,Float64},1}(undef,length(P))
    for pl in 1:length(P)
        a = semimajor_axis(P[order[pl]], mass[order[pl]]+star_mass)
        if use_mutualHill
            mu = (mass[order[pl]] + insert_pl_mass)/star_mass
            term = (min_num_mutual_hill_radii/2.)*(mu/3.)^(1//3)
            a_lower = a*(1. - ecc[order[pl]] - term)/(1. + insert_pl_ecc + term)
            a_upper = (0 < 1-insert_pl_ecc-term) ? a*(1. + ecc[order[pl]] + term)/(1. - insert_pl_ecc - term) : Inf # 0<(1-insert_pl_ecc-term) must be true or else inequality flips sign and there are no physical solutions for a_upper
        else # use Hill radii
            mu = mass[order[pl]]/star_mass
            hill_radius = calc_hill_sphere(a, mu)
            a_lower, a_upper = a*(1. - ecc[order[pl]]) - hill_radius*min_num_mutual_hill_radii, a*(1. + ecc[order[pl]]) + hill_radius*min_num_mutual_hill_radii
        end
        a_lower = max(a_lower, 0.) # set a_lower=0 if negative
        @assert 0 <= a_lower < a_upper
        P_lower, P_upper = period_given_semimajor_axis(a_lower, insert_pl_mass+star_mass), period_given_semimajor_axis(a_upper, insert_pl_mass+star_mass)
        if verbose
            println("P_blocked: ", (P_lower, P_upper))
        end
        P_segments_unstable[pl] = (P_lower, P_upper)
    end
    return P_segments_unstable
end

"""
    compute_unstable_mutualHill_regions_periodscales_given_clusters(pl_per_cl, P_cl, mass_cl, insert_cl_ρ, insert_cl_mass, star_mass, sim_param; ecc_cl=zeros(length(P_cl)), insert_cl_ecc=zeros(length(insert_cl_mass)), verbose=false)

Compute the unstable intervals in period scale given the existing clusters and the cluster we are trying to insert, requiring that no clusters overlap.

# Arguments:
- `pl_per_cl::Vector{Int64}`: number of planets in each existing cluster.
- `P_cl::AbstractVector{Float64}`: periods of the planets in the existing clusters.
- `mass_cl::AbstractVector{Float64}`: masses of the planets in the existing clusters.
- `insert_cl_ρ::AbstractVector{Float64}`: unscaled periods of the planets in the cluster we are trying to insert.
- `insert_cl_mass::AbstractVector{Float64}`: masses of the planets in the cluster we are trying to insert.
- `star_mass::Float64`: mass of the star.
- `sim_param::SimParam`: a SimParam object containing various simulation parameters.
- `ecc_cl::AbstractVector{Float64}=zeros(length(P_cl))`: eccentricities of the planets in the existing clusters (defaults to circular).
- `insert_cl_ecc::AbstractVector{Float64}=zeros(length(insert_cl_mass))`: eccentricities of the planets in the cluster we are trying to insert (defaults to circular).
- `verbose::Bool=false`: whether to print various messages.
NOTE: the existing clusters must be grouped (and in the order implied by `pl_per_cl`), but the planets do not have to be sorted in period.

# Returns:
`Pc_segments_unstable_nonoverlapping::Array{Tuple{Float64,Float64},1}`: an array of (non-overlapping) segments. Inserting the cluster somewhere inside one of these segments would either deem the system unstable, or cause clusters to overlap.
"""
function compute_unstable_mutualHill_regions_periodscales_given_clusters(pl_per_cl::Vector{Int64}, P_cl::AbstractVector{Float64}, mass_cl::AbstractVector{Float64}, insert_cl_ρ::AbstractVector{Float64}, insert_cl_mass::AbstractVector{Float64}, star_mass::Float64, sim_param::SimParam; ecc_cl::AbstractVector{Float64}=zeros(length(P_cl)), insert_cl_ecc::AbstractVector{Float64}=zeros(length(insert_cl_mass)), verbose::Bool=false)
    n_cl = length(pl_per_cl)
    @assert all(pl_per_cl .> 0)
    @assert sum(pl_per_cl) == length(P_cl) == length(mass_cl) == length(ecc_cl)
    @assert length(insert_cl_ρ) == length(insert_cl_mass) == length(insert_cl_ecc)
    min_num_mutual_hill_radii = get_real(sim_param, "num_mutual_hill_radii")

    # Get the inner-most and outer-most planets in the cluster we are trying to insert:
    order_insert_cl = sortperm(insert_cl_ρ) # ρ denotes unscaled periods in this function
    insert_pl_in_ρ, insert_pl_in_mass, insert_pl_in_ecc = insert_cl_ρ[order_insert_cl[1]], insert_cl_mass[order_insert_cl[1]], insert_cl_ecc[order_insert_cl[1]]
    insert_pl_out_ρ, insert_pl_out_mass, insert_pl_out_ecc = insert_cl_ρ[order_insert_cl[end]], insert_cl_mass[order_insert_cl[end]], insert_cl_ecc[order_insert_cl[end]]
    max_ratio = insert_pl_out_ρ/insert_pl_in_ρ

    Pc_segments_unstable = Array{Tuple{Float64,Float64},1}(undef,n_cl)
    pl_start = 1
    pl_stop = 0
    for c in 1:n_cl
        pl_stop += pl_per_cl[c]
        order = (pl_start - 1) .+ sortperm(P_cl[pl_start:pl_stop])

        # Compare the inner-most planet in this cluster to the outer-most planet in the cluster we are trying to fit in to calculate the lower bound for the unstable region:
        a_in = semimajor_axis(P_cl[order[1]], mass_cl[order[1]]+star_mass)
        mu = (mass_cl[order[1]] + insert_pl_out_mass)/star_mass
        term = (min_num_mutual_hill_radii/2.)*(mu/3.)^(1//3)
        a_lower = a_in*(1. - ecc_cl[order[1]] - term)/(1. + insert_pl_out_ecc + term)
        a_lower = max(a_lower, 0.) # set a_lower=0 if negative
        P_lower = period_given_semimajor_axis(a_lower, insert_pl_out_mass+star_mass)

        # Compare the outer-most planet in this cluster to the inner-most planet in the cluster we are trying to fit in to calculate the upper bound for the unstable region:
        a_out = semimajor_axis(P_cl[order[end]], mass_cl[order[end]]+star_mass)
        mu = (mass_cl[order[end]] + insert_pl_in_mass)/star_mass
        term = (min_num_mutual_hill_radii/2.)*(mu/3.)^(1//3)
        a_upper = (0 < 1-insert_pl_in_ecc-term) ? a_out*(1. + ecc_cl[order[end]] + term)/(1. - insert_pl_in_ecc - term) : Inf # 0<(1-insert_pl_ecc-term) must be true or else inequality flips sign and there are no physical solutions for a_upper
        P_upper = period_given_semimajor_axis(a_upper, insert_pl_in_mass+star_mass)

        # Now, convert P_lower and P_upper to bounds for the period scale (which depend on the unscaled periods of the cluster we are trying to insert):
        Pc_lower, Pc_upper = P_lower/insert_pl_out_ρ, P_upper/insert_pl_in_ρ # Pc_lower should decrease if insert_pl_out_ρ > 1, and likewise Pc_upper should increase if insert_pl_in_ρ < 1
        if verbose
            println("Pc_blocked: ", (Pc_lower, Pc_upper))
        end
        Pc_segments_unstable[c] = (Pc_lower, Pc_upper)

        pl_start += pl_per_cl[c]
    end
    Pc_segments_unstable = sort(Pc_segments_unstable)

    Pc_segments_unstable_nonoverlapping = Tuple{Float64,Float64}[]
    x_start, x_stop = Pc_segments_unstable[1]
    for i in 1:length(Pc_segments_unstable)-1
        if x_stop > Pc_segments_unstable[i+1][1] # segments overlap
            x_stop = Pc_segments_unstable[i+1][2]
        else # segments do not overlap
            if Pc_segments_unstable[i+1][1]/x_stop > max_ratio
                push!(Pc_segments_unstable_nonoverlapping, (x_start, x_stop))
                x_start, x_stop = Pc_segments_unstable[i+1]
            else # segments do not overlap, but there is not enough room between the two segments to insert the cluster
                x_stop = Pc_segments_unstable[i+1][2]
            end
        end
    end
    push!(Pc_segments_unstable_nonoverlapping, (x_start, x_stop))
    if verbose
        println("Pc_blocked with no room: ", Pc_segments_unstable_nonoverlapping)
    end
    return Pc_segments_unstable_nonoverlapping
end

"""
    compute_regions_non_overlapping(segments_blocked)

Compute the non-overlapping segments given an array of (potentially overlapping) segments.

# Arguments:
- `segments_blocked::Array{Tuple{Float64,Float64},1}`: array of potentially overlapping segments.
NOTE: the input segments must be sorted in increasing order!

# Returns:
`segments_blocked_nonoverlapping::Array{Tuple{Float64,Float64},1}`: array of non-overlapping segments.
"""
function compute_regions_non_overlapping(segments_blocked::Array{Tuple{Float64,Float64},1})
    # Make a list of non-overlapping segments:
    segments_blocked_nonoverlapping = Tuple{Float64,Float64}[]
    x_start, x_stop = segments_blocked[1]
    for i in 1:length(segments_blocked)-1
        if x_stop > segments_blocked[i+1][1] # segments overlap
            x_stop = segments_blocked[i+1][2]
        else # segments do not overlap
            push!(segments_blocked_nonoverlapping, (x_start, x_stop))
            x_start, x_stop = segments_blocked[i+1]
        end
    end
    push!(segments_blocked_nonoverlapping, (x_start, x_stop))
    return segments_blocked_nonoverlapping
end

"""
    compute_allowed_regions_cdf_lognormal(segments_blocked; μ=0., σ=1., x_min=0., x_max=Inf, ϵ=1e-12)

Compute the allowed regions in the cumulative distribution function (CDF) of a lognormal distribution, given blocked segments in the domain of the variable.

# Arguments:
- `segments_blocked::Array{Tuple{Float64,Float64},1}`: array of blocked segments.
- `μ::Float64=0.`: mean parameter for the lognormal distribution.
- `σ::Float64=1.`: standard deviation parameter for the lognormal distribution.
- `x_min::Float64=0.`: minimum value of the variable.
- `x_max::Float64=Inf`: maximum value of the variable.
- `ϵ::Float64=1e-12`: minimum segment length in CDF to be included as an allowed segment, to reduce the rate of error due to machine precision.
NOTE: the array of blocked segments do not have to be non-overlapping; this function will first compute the non-overlapping segments.

# Returns:
`cdf_segments_allowed::Array{Tuple{Float64,Float64},1}`: array of (non-overlapping) segments in the CDF that are not blocked.

NOTE: returns an empty array (length zero) if the entire domain between `x_min` and `x_max` is blocked!
"""
function compute_allowed_regions_cdf_lognormal(segments_blocked::Array{Tuple{Float64,Float64},1}; μ::Float64=0., σ::Float64=1., x_min::Float64=0., x_max::Float64=Inf, ϵ::Float64=1e-12)
    @assert 0 <= x_min < x_max <= Inf

    # Make a list of non-overlapping segments in period deemed unstable:
    segments_blocked_nonoverlapping = compute_regions_non_overlapping(segments_blocked)

    # Make a list of non-overlapping, allowed segments in the cdf:
    cdf_segments_allowed = Tuple{Float64,Float64}[]
    n_seg_max = length(segments_blocked_nonoverlapping)+1 # maximum number of allowed segments if all the unstable segments fit between x_min and x_max
    for i in 1:n_seg_max
        if i==1 && x_min < segments_blocked_nonoverlapping[i][1]
            cdf_start = cdf_lognormal(x_min; μ=μ, σ=σ)
            cdf_stop = cdf_lognormal(segments_blocked_nonoverlapping[i][1]; μ=μ, σ=σ)
            if cdf_stop - cdf_start > ϵ
                push!(cdf_segments_allowed, (cdf_start, cdf_stop))
            end
        elseif i==n_seg_max && x_max > segments_blocked_nonoverlapping[i-1][2]
            cdf_start = cdf_lognormal(segments_blocked_nonoverlapping[i-1][2]; μ=μ, σ=σ)
            cdf_stop = cdf_lognormal(x_max; μ=μ, σ=σ)
            if cdf_stop - cdf_start > ϵ
                push!(cdf_segments_allowed, (cdf_start, cdf_stop))
            end
        elseif 1 < i < n_seg_max
            cdf_start = cdf_lognormal(segments_blocked_nonoverlapping[i-1][2]; μ=μ, σ=σ)
            cdf_stop = cdf_lognormal(segments_blocked_nonoverlapping[i][1]; μ=μ, σ=σ)
            if cdf_stop - cdf_start > ϵ
                push!(cdf_segments_allowed, (cdf_start, cdf_stop))
            end
        end
    end

    return cdf_segments_allowed
end

"""
    compute_allowed_regions_cdf_power_law(segments_blocked; x0, x1, α, ϵ=1e-12)

Compute the allowed regions in the cumulative distribution function (CDF) of a power-law distribution, given blocked segments in the domain of the variable.

# Arguments:
- `segments_blocked::Array{Tuple{Float64,Float64},1}`: array of blocked segments.
- `x0::Float64`: minimum value of the variable.
- `x1::Float64`: maximum value of the variable.
- `α::Float64`: power-law index.
- `ϵ::Float64=1e-12`: minimum segment length in CDF to be included as an allowed segment, to reduce the rate of error due to machine precision.
NOTE: the array of blocked segments do not have to be non-overlapping; this function will first compute the non-overlapping segments.

# Returns:
`cdf_segments_allowed::Array{Tuple{Float64,Float64},1}`: array of (non-overlapping) segments in the CDF that are not blocked.

NOTE: returns an empty array (length zero) if the entire domain between `x0` and `x1` is blocked!
"""
function compute_allowed_regions_cdf_power_law(segments_blocked::Array{Tuple{Float64,Float64},1}; x0::Float64, x1::Float64, α::Float64, ϵ::Float64=1e-12)
    @assert x0 < x1

    # Make a list of non-overlapping segments in period scale deemed unstable:
    segments_blocked_nonoverlapping = compute_regions_non_overlapping(segments_blocked)

    # Make a list of non-overlapping, allowed segments in the cdf:
    cdf_segments_allowed = Tuple{Float64,Float64}[]
    n_seg_max = length(segments_blocked_nonoverlapping)+1 # maximum number of allowed segments if all the unstable segments fit between x0 and x1
    for i in 1:n_seg_max
        if i==1 && x0 < segments_blocked_nonoverlapping[i][1]
            cdf_start = 0.
            cdf_stop = cdf_power_law(segments_blocked_nonoverlapping[i][1]; x0=x0, x1=x1, α=α)
            if cdf_stop - cdf_start > ϵ
                push!(cdf_segments_allowed, (cdf_start, cdf_stop))
            end
        elseif i==n_seg_max && x1 > segments_blocked_nonoverlapping[i-1][2]
            cdf_start = cdf_power_law(segments_blocked_nonoverlapping[i-1][2]; x0=x0, x1=x1, α=α)
            cdf_stop = 1.
            if cdf_stop - cdf_start > ϵ
                push!(cdf_segments_allowed, (cdf_start, cdf_stop))
            end
        elseif 1 < i < n_seg_max
            cdf_start = cdf_power_law(segments_blocked_nonoverlapping[i-1][2]; x0=x0, x1=x1, α=α)
            cdf_stop = cdf_power_law(segments_blocked_nonoverlapping[i][1]; x0=x0, x1=x1, α=α)
            if cdf_stop - cdf_start > ϵ
                push!(cdf_segments_allowed, (cdf_start, cdf_stop))
            end
        end
    end

    return cdf_segments_allowed
end

"""
    draw_lognormal_allowed_regions(segments_blocked; μ=0., σ=1., x_min=0., x_max=Inf, ϵ=1e-12, verbose=false)

Draw a random value that is not inside a blocked region, from a lognormal distribution.

# Arguments:
- `segments_blocked::Array{Tuple{Float64,Float64},1}`: array of blocked segments.
- `μ::Float64=0.`: mean parameter for the lognormal distribution.
- `σ::Float64=1.`: standard deviation parameter for the lognormal distribution.
- `x_min::Float64=0.`: minimum value of the variable.
- `x_max::Float64=Inf`: maximum value of the variable.
- `ϵ::Float64=1e-12`: minimum segment length in CDF to be included as an allowed segment, to reduce the rate of error due to machine precision.
- `verbose::Bool=false`: whether to print various messages.
NOTE: the array of blocked segments do not have to be non-overlapping.

# Returns:
A value (float) drawn from the lognormal distribution, that is not inside a blocked region.

NOTE: returns NaN if the entire region between `x_min` and `x_max` is blocked!
"""
function draw_lognormal_allowed_regions(segments_blocked::Array{Tuple{Float64,Float64},1}; μ::Float64=0., σ::Float64=1., x_min::Float64=0., x_max::Float64=Inf, ϵ::Float64=1e-12, verbose::Bool=false)
    cdf_segments_allowed = compute_allowed_regions_cdf_lognormal(segments_blocked; μ=μ, σ=σ, x_min=x_min, x_max=x_max, ϵ=ϵ)
    if length(cdf_segments_allowed) == 0
        if verbose
            println("No allowed regions left to draw from; returning NaN.")
        end
        return NaN
    else
        y = draw_segmented_uniform(cdf_segments_allowed)
        return invert_cdf_lognormal(y; μ=μ, σ=σ)
    end
end

"""
    draw_power_law_allowed_regions(segments_blocked; x0, x1, α, ϵ=1e-12, verbose=false)

Draw a random value that is not inside a blocked region, from a power-law distribution.

# Arguments:
- `segments_blocked::Array{Tuple{Float64,Float64},1}`: array of blocked segments.
- `x0::Float64`: minimum value of the variable.
- `x1::Float64`: maximum value of the variable.
- `α::Float64`: power-law index.
- `ϵ::Float64=1e-12`: minimum segment length in CDF to be included as an allowed segment, to reduce the rate of error due to machine precision.
- `verbose::Bool=false`: whether to print various messages.
NOTE: the array of blocked segments do not have to be non-overlapping.

# Returns:
A value (float) drawn from the power-law distribution, that is not inside a blocked region.

NOTE: returns NaN if the entire region between `x0` and `x1` is blocked!
"""
function draw_power_law_allowed_regions(segments_blocked::Array{Tuple{Float64,Float64},1}; x0::Float64, x1::Float64, α::Float64, ϵ::Float64=1e-12, verbose::Bool=false)
    cdf_segments_allowed = compute_allowed_regions_cdf_power_law(segments_blocked; x0=x0, x1=x1, α=α, ϵ=ϵ)
    if length(cdf_segments_allowed) == 0
        if verbose
            println("No allowed regions left to draw from; returning NaN.")
        end
        return NaN
    else
        y = draw_segmented_uniform(cdf_segments_allowed)
        return invert_cdf_power_law(y; x0=x0, x1=x1, α=α)
    end
end

"""
    draw_period_lognormal_allowed_regions(P, mass, insert_pl_mass, star_mass, sim_param; μ=0., σ=1., x_min=0., x_max=Inf, ecc=zeros(length(P)), insert_pl_ecc=0., use_mutualHill=true, ϵ=1e-12, verbose=false)

Draw a period from a lognormal distribution, that would allow the system to be stable if the planet is inserted there.

# Arguments:
- `P::AbstractVector{Float64}`: periods of the existing planets.
- `mass::AbstractVector{Float64}`: masses of the existing planets.
- `insert_pl_mass::Float64`: mass of the planet we are trying to insert.
- `star_mass::Float64`: mass of the star.
- `sim_param::SimParam`: a SimParam object containing various simulation parameters.
- `μ::Float64=0.`: mean parameter for the lognormal distribution.
- `σ::Float64=1.`: standard deviation parameter for the lognormal distribution.
- `x_min::Float64=0.`: minimum value of the variable.
- `x_max::Float64=Inf`: maximum value of the variable.
- `ecc::AbstractVector{Float64}=zeros(length(P))`: eccentricities of the existing planets (defaults to circular).
- `insert_pl_ecc::Float64=0.`: eccentricity of the planet we are trying to insert (defaults to circular).
- `use_mutualHill::Bool=true`: whether to use the minimum separation in terms of the mutual Hill radii (default) or the Hill radii (if false).
- `ϵ::Float64=1e-12`: minimum segment length in CDF to be included as an allowed segment, to reduce the rate of error due to machine precision.
- `verbose::Bool=false`: whether to print various messages.
NOTE: the existing planets do not have to be sorted in period.

# Returns:
`P_draw::Float64`: period drawn from the lognormal distribution, that would allow the system to be stable if the planet is inserted there.

NOTE: returns NaN if any period between `x_min` and `x_max` would make the system unstable!
"""
function draw_period_lognormal_allowed_regions(P::AbstractVector{Float64}, mass::AbstractVector{Float64}, insert_pl_mass::Float64, star_mass::Float64, sim_param::SimParam; μ::Float64=0., σ::Float64=1., x_min::Float64=0., x_max::Float64=Inf, ecc::AbstractVector{Float64}=zeros(length(P)), insert_pl_ecc::Float64=0., use_mutualHill::Bool=true, ϵ::Float64=1e-12, verbose::Bool=false)
    @assert length(P) == length(mass) == length(ecc)
    @assert 0 <= insert_pl_mass < star_mass
    @assert all(P .> 0)
    @assert all(mass .>= 0)
    @assert all(1 .> ecc .>= 0)
    @assert 0 <= insert_pl_ecc < 1

    if length(P) > 0
        P_segments_unstable = compute_unstable_regions_periods_given_planets(P, mass, insert_pl_mass, star_mass, sim_param; ecc=ecc, insert_pl_ecc=insert_pl_ecc, use_mutualHill=use_mutualHill, verbose=verbose)
        P_draw = draw_lognormal_allowed_regions(P_segments_unstable; μ=μ, σ=σ, x_min=x_min, x_max=x_max, ϵ=ϵ, verbose=verbose)
    else
        P_draw = invert_cdf_lognormal(rand(Uniform(cdf_lognormal(x_min; μ=μ, σ=σ), cdf_lognormal(x_max; μ=μ, σ=σ))); μ=μ, σ=σ)
    end
    return P_draw
end

"""
    draw_period_lognormal_allowed_regions_Hill(P, mass, insert_pl_mass, star_mass, sim_param; μ=0., σ=1., x_min=0., x_max=Inf, ecc=zeros(length(P)), insert_pl_ecc=0., ϵ=1e-12, verbose=false)

Draw a period that is not inside a blocked region using Hill radii, from a lognormal distribution, by calling `draw_period_lognormal_allowed_regions`.

NOTE: returns NaN if any period between `x_min` and `x_max` would make the system unstable!
"""
draw_period_lognormal_allowed_regions_Hill(P::AbstractVector{Float64}, mass::AbstractVector{Float64}, insert_pl_mass::Float64, star_mass::Float64, sim_param::SimParam; μ::Float64=0., σ::Float64=1., x_min::Float64=0., x_max::Float64=Inf, ecc::AbstractVector{Float64}=zeros(length(P)), insert_pl_ecc::Float64=0., ϵ::Float64=1e-12, verbose::Bool=false) = draw_period_lognormal_allowed_regions(P, mass, insert_pl_mass, star_mass, sim_param; μ=μ, σ=σ, x_min=x_min, x_max=x_max, ecc=ecc, insert_pl_ecc=insert_pl_ecc, use_mutualHill=false, ϵ=ϵ, verbose=verbose)

"""
    draw_period_lognormal_allowed_regions_mutualHill(P, mass, insert_pl_mass, star_mass, sim_param; μ=0., σ=1., x_min=0., x_max:=Inf, ecc=zeros(length(P)), insert_pl_ecc=0., ϵ=1e-12, verbose=false)

Draw a period that is not inside a blocked region using mutual Hill radii, from a lognormal distribution, by calling `draw_period_lognormal_allowed_regions`.

NOTE: returns NaN if any period between `x_min` and `x_max` would make the system unstable!
"""
draw_period_lognormal_allowed_regions_mutualHill(P::AbstractVector{Float64}, mass::AbstractVector{Float64}, insert_pl_mass::Float64, star_mass::Float64, sim_param::SimParam; μ::Float64=0., σ::Float64=1., x_min::Float64=0., x_max::Float64=Inf, ecc::AbstractVector{Float64}=zeros(length(P)), insert_pl_ecc::Float64=0., ϵ::Float64=1e-12, verbose::Bool=false) = draw_period_lognormal_allowed_regions(P, mass, insert_pl_mass, star_mass, sim_param; μ=μ, σ=σ, x_min=x_min, x_max=x_max, ecc=ecc, insert_pl_ecc=insert_pl_ecc, use_mutualHill=true, ϵ=ϵ, verbose=verbose)

"""
    draw_periodscale_power_law_allowed_regions_mutualHill(pl_per_cl, P_cl:, mass_cl:, insert_cl_ρ, insert_cl_mass, star_mass, sim_param; x0, x1, α, ecc_cl=zeros(length(P_cl)), insert_cl_ecc=zeros(length(insert_cl_mass)), ϵ=1e-12, verbose=false)

Draw a period scale from a power-law distribution, that would allow the system to be stable if the cluster is inserted there (and cause no overlapping clusters).

# Arguments:
- `pl_per_cl::Vector{Int64}`: number of planets in each existing cluster.
- `P_cl::AbstractVector{Float64}`: periods of the planets in the existing clusters.
- `mass_cl::AbstractVector{Float64}`: masses of the planets in the existing clusters.
- `insert_cl_ρ::AbstractVector{Float64}`: unscaled periods of the planets in the cluster we are trying to insert.
- `insert_cl_mass::AbstractVector{Float64}`: masses of the planets in the cluster we are trying to insert.
- `star_mass::Float64`: mass of the star.
- `sim_param::SimParam`: a SimParam object containing various simulation parameters.
- `x0::Float64`: minimum value of the variable.
- `x1::Float64`: maximum value of the variable.
- `α::Float64`: power-law index.
- `ecc_cl::AbstractVector{Float64}=zeros(length(P_cl))`: eccentricities of the planets in the existing clusters (defaults to circular).
- `insert_cl_ecc::AbstractVector{Float64}=zeros(length(insert_cl_mass))`: eccentricities of the planets in the cluster we are trying to insert (defaults to circular).
- `ϵ::Float64=1e-12`: minimum segment length in CDF to be included as an allowed segment, to reduce the rate of error due to machine precision.
- `verbose::Bool=false`: whether to print various messages.
NOTE: the existing clusters must be grouped (and in the order implied by `pl_per_cl`), but the planets do not have to be sorted in period.

# Returns:
`Pc_draw::Float64`: period scale drawn from the power-law distribution, that would allow the system to be stable if the cluster is inserted there.

NOTE: returns NaN if any period scale between `x0` and `x1` would make the system unstable or cause clusters to overlap!
"""
function draw_periodscale_power_law_allowed_regions_mutualHill(pl_per_cl::Vector{Int64}, P_cl::AbstractVector{Float64}, mass_cl::AbstractVector{Float64}, insert_cl_ρ::AbstractVector{Float64}, insert_cl_mass::AbstractVector{Float64}, star_mass::Float64, sim_param::SimParam; x0::Float64, x1::Float64, α::Float64, ecc_cl::AbstractVector{Float64}=zeros(length(P_cl)), insert_cl_ecc::AbstractVector{Float64}=zeros(length(insert_cl_mass)), ϵ::Float64=1e-12, verbose::Bool=false)
    @assert x0 < x1
    @assert sum(pl_per_cl) == length(P_cl) == length(mass_cl) == length(ecc_cl)
    @assert all(P_cl .> 0)
    @assert all(mass_cl .> 0)
    @assert all(0 .<= ecc_cl .< 1)
    @assert length(insert_cl_ρ) == length(insert_cl_mass) == length(insert_cl_ecc)
    @assert all(0 .< insert_cl_ρ)
    @assert all(0 .< insert_cl_mass .< star_mass)
    @assert all(0 .<= insert_cl_ecc .< 1)

    if length(P_cl) > 0
        Pc_segments_unstable = compute_unstable_mutualHill_regions_periodscales_given_clusters(pl_per_cl, P_cl, mass_cl, insert_cl_ρ, insert_cl_mass, star_mass, sim_param; ecc_cl=ecc_cl, insert_cl_ecc=insert_cl_ecc, verbose=verbose)
        Pc_draw = draw_power_law_allowed_regions(Pc_segments_unstable; x0=x0, x1=x1, α=α, ϵ=ϵ, verbose=verbose)
    else
        Pc_draw = invert_cdf_power_law(rand(); x0=x0, x1=x1, α=α)
    end
    return Pc_draw
end

function TruncatedUpper(d::Distributions.UnivariateDistribution, u::Float64)
    zero(u) < u || error("lower bound should be less than upper bound.")
    lcdf = zero(u)
    ucdf = isinf(u) ? one(u) : cdf(d, u)
    tp = ucdf - lcdf
    Distributions.Truncated{typeof(d),Distributions.value_support(typeof(d))}(d, zero(u), u, lcdf, ucdf, tp, log(tp))
end

function generate_e_omega_rayleigh_direct(sigma_hk::Float64; max_e::Float64 = 1.0)
  @assert(0<max_e<=1.0)
  ecc::Float64 = rand( TruncatedUpper(Rayleigh(sigma_hk),max_e) )
  w::Float64 = 2pi*rand()
  return ecc, w
end

function generate_e_omega_rayleigh_two_gaussians(sigma_hk::Float64; max_e::Float64 = 1.0)
  @assert(0<max_e<=1.0)
  h = k = 1.0
  while h*h+k*k >= max_e*max_e
    h = sigma_hk*randn()
    k = sigma_hk*randn()
  end
  ecc::Float64 = sqrt(h*h+k*k)
  #w::Float64 = atan2(k,h)
  w::Float64 = atan(k,h)
  return ecc, w
end

function generate_e_omega_rayleigh(sigma_hk::Float64; max_e::Float64 = 1.0)
   if max_e > sigma_hk
      return generate_e_omega_rayleigh_two_gaussians(sigma_hk,max_e=max_e)
   else
      return generate_e_omega_rayleigh_direct(sigma_hk,max_e=max_e)
   end
end

function generate_e_omega_rayleigh(sim_param::SimParam; max_e::Float64 = 1.0)
  sigma_hk::Float64 = get_real(sim_param,"sigma_hk")
  generate_e_omega_rayleigh(sigma_hk, max_e=max_e)
end

struct ParamsTriangle
    id_xy::Tuple{Int64,Int64} # indices of the two transformed parameters in a vector of parameters
    # A,B,C are vertices of the triangle for parameters (x,y), in no particular order
    A::Tuple{Float64,Float64}
    B::Tuple{Float64,Float64}
    C::Tuple{Float64,Float64}
end

"""
    map_square_to_triangle(r1, r2, pt)

Map a point (r1, r2) in the unit square to a point P=(x,y) in the triangle defined by vertices A,B,C. If r1,r2 are uniformly drawn in [0,1], then the point P is also uniformly drawn in the triangle; see http://www.cs.princeton.edu/~funk/tog02.pdf (Section 4.2) for a reference.

# Arguments:
- `r1r2::Tuple{Float64,Float64}`: a point in the unit square representing two transformed parameters.
- `pt::ParamsTriangle`: an object containing the indices of the two transformed parameters and vertices of the triangle defining the bounds for the two parameters.

# Returns:
- `P::Tuple{Float64,Float64}`: a point for the two parameters in the triangle corresponding to the point given by `r1r2` in the unit square.
"""
function map_square_to_triangle(r1r2::Tuple{Float64,Float64}, pt::ParamsTriangle)
    @assert all(0 .<= r1r2 .<= 1)
    r1, r2 = r1r2
    P = (1. - sqrt(r1)) .* pt.A .+ (sqrt(r1)*(1. - r2)) .* pt.B .+ (sqrt(r1)*r2) .* pt.C
    return P
end

function generate_planetary_system_hardcoded_example(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # in this version we specify fixed functions that are known at compile time, allowing for additional optimizations (~0.6 second faster per Kepler catalog out of ~3.6 sec on my laptop w/ 1 core)
  generate_planet_mass_from_radius = generate_planet_mass_from_radius_powerlaw
  generate_num_planets = generate_num_planets_poisson
  generate_period_and_sizes = generate_period_and_sizes_log_normal

  generate_e_omega =  generate_e_omega_rayleigh

  # generate_star = get_function(sim_param,"generate_star")
  # star::StarAbstract = generate_star(sim_param)
  num_pl::Int64 = generate_num_planets(star, sim_param)

  if( num_pl==0 )
    return PlanetarySystem(star)
  else
    (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
    idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred

    min_a_in_rstar = 2.0
    min_P_orbit = day_in_year*sqrt((min_a_in_rstar*star.radius*rsol_in_au)^3 / star.mass) # minimum semi-major axis of two stellar radii
    idx = idx[Plist[idx] .> min_P_orbit]
    if( length(idx)==0 )
        return PlanetarySystem(star)
    end

    pl = Array{Planet}(length(idx))
    orbit = Array{Orbit}(length(idx))
    a = map(i->semimajor_axis(Plist[i],star.mass),idx)
    max_e = ones(length(idx))
    max_e_factor = 0.999 # A factor just less than 1 to prevent numerical issues with near-crossing orbits
    if length(a)>=2
       for i in 1:length(a)
         if i==1
           max_e[i] = max_e_factor*(1-a[i]/a[i+1])/(1+a[i]/a[i+1])
         elseif i==length(a)
           max_e[i] = max_e_factor*(1-a[i-1]/a[i])/(1+a[i-1]/a[i])
         else
           max_e[i] = max_e_factor*min( (1-a[i]/a[i+1])/(1+a[i]/a[i+1]), (1-a[i-1]/a[i])/(1+a[i-1]/a[i]) )
         end
       end
    end
    for i in 1:length(idx)
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      Rpl = Rlist[idx[i]]
      (ecc::Float64,  omega::Float64) = generate_e_omega(sim_param, max_e=max_e[i])
      incl::Float64 = acos(rand())
      orbit[i] = Orbit(P,ecc,incl,omega,2pi*rand(),2pi*rand())
      mass::Float64 = generate_planet_mass_from_radius(Rpl, star, orbit[i], sim_param)
      pl[i] = Planet( Rpl,  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

function generate_planetary_system_uncorrelated_incl(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # load functions to use for drawing parameters
   generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
   generate_num_planets = get_function(sim_param,"generate_num_planets")
   generate_period_and_sizes = get_function(sim_param,"generate_period_and_sizes")
   generate_e_omega = get_function(sim_param,"generate_e_omega")

  #  generate_star = get_function(sim_param,"generate_star")
  # star::StarAbstract = generate_star(sim_param)
  num_pl::Int64 = generate_num_planets(star, sim_param)::Int64
  sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

  if( num_pl==0 )
    return PlanetarySystem(star)::PlanetarySystem
  else
     (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
    idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred

    min_a_in_rstar = 2.0
    min_P_orbit = day_in_year*sqrt((min_a_in_rstar*star.radius*rsol_in_au)^3 / star.mass) # minimum semi-major axis of two stellar radii
    idx = idx[Plist[idx] .> min_P_orbit]
    if( length(idx)==0 )
        return PlanetarySystem(star)
    end

    pl = Array{Planet}(undef,length(idx))
    orbit = Array{Orbit}(undef,length(idx))
    a = map(i->semimajor_axis(Plist[i],star.mass),idx)
    max_e = ones(length(idx))
    max_e_factor = 0.999 # A factor just less than 1 to prevent numerical issues with near-crossing orbits
    if length(a)>=2
       for i in 1:length(a)
          if i==1
            max_e[i] = max_e_factor*(1-a[i]/a[i+1])/(1+a[i]/a[i+1])
          elseif i==length(a)
            max_e[i] = max_e_factor*(1-a[i-1]/a[i])/(1+a[i-1]/a[i])
          else
            max_e[i] = max_e_factor*min( (1-a[i]/a[i+1])/(1+a[i]/a[i+1]), (1-a[i-1]/a[i])/(1+a[i-1]/a[i]) )
          end
        end
    end
    for i in 1:length(idx)
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      Rpl = Rlist[idx[i]]
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      (ecc::Float64,  omega::Float64) = generate_e_omega(sim_param, max_e=max_e[i])
      incl::Float64 = acos(rand())
      orbit[i] = Orbit(P,ecc,incl,omega,2pi*rand(),2pi*rand())
      # set!(orbit[idx[i]],P,ecc,incl,omega,2pi*rand(),2pi*rand())
      mass::Float64 = generate_planet_mass_from_radius(Rpl, star, orbit[i], sim_param)
      pl[i] = Planet( Rpl,  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

# This version generates more systems roughly near a common plane, but until incorporate CORBITS data, ABC can't match input param
function generate_planetary_system_simple(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # load functions to use for drawing parameters
    generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
    generate_num_planets = get_function(sim_param,"generate_num_planets")
  #  generate_num_planets = generate_num_planets_christiansen
  #  generate_period_and_sizes = generate_period_and_sizes_christiansen
    generate_period_and_sizes = get_function(sim_param,"generate_period_and_sizes")
    generate_e_omega = get_function(sim_param,"generate_e_omega")
    sigma_incl = deg2rad(get_real(sim_param,"sigma_incl"))

  #   generate_star = get_function(sim_param,"generate_star")
  #  star::StarAbstract = generate_star(sim_param)
   num_pl = generate_num_planets(star, sim_param)::Int64
   sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

  if( num_pl==0 )
    return PlanetarySystem(star)
  else
    (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
    idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred
    incl_sys = acos(rand())

    min_a_in_rstar = 2.0
    min_P_orbit = day_in_year*sqrt((min_a_in_rstar*star.radius*rsol_in_au)^3 / star.mass) # minimum semi-major axis of two stellar radii
    idx = idx[Plist[idx] .> min_P_orbit]
    if( length(idx)==0 )
        return PlanetarySystem(star)
    end

    pl = Array{Planet}(undef,length(idx))
    orbit = Array{Orbit}(undef,length(idx))
    for i in 1:length(idx)
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      Rpl = Rlist[idx[i]]
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      (ecc,  omega) = generate_e_omega(sim_param)::Tuple{Float64,Float64}
      incl_mut = sigma_incl*sqrt(randn()^2+randn()^2) # rand(Distributions.Rayleigh(sigma_incl)) # sigma_incl*randn()
      asc_node = 2pi*rand()
      mean_anom = 2pi*rand()
      #incl = incl_sys + sigma_incl*randn()
      incl =  incl_mut!=zero(incl_mut) ? acos( cos(incl_sys)*cos(incl_mut) + sin(incl_sys)*sin(incl_mut)*cos(asc_node) ) : incl_sys
      orbit[i] = Orbit(P,ecc,incl,omega,asc_node,mean_anom)
      mass = generate_planet_mass_from_radius(Rpl, star, orbit[i], sim_param)::Float64
      pl[i] = Planet( Rpl,  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

function test_planetary_system_constructors(sim_param::SimParam)
  generate_star = get_function(sim_param,"generate_star")
  star = generate_star(sim_param)
  empty_sys = PlanetarySystem(star)
  earth = Planet(earth_radius,earth_mass)
  earth_orbit = Orbit(365.2425,0.0167,0.5*pi,0.0,0.0,0.0)
  solar_sys = PlanetarySystem(star, earth,earth_orbit)
  m = generate_planet_mass_from_radius_powerlaw(0.02,star,earth_orbit,sim_param)/earth_mass
  generate_planetary_system_simple(star,sim_param,verbose=true)
end


text:
@testset "GAS recursion" begin
    @testset "Normal Distribution 0 scaling" begin
        # Normal distribution
        # ω = [1.0; 0.1]
        # A = Diagonal([0.0; 0.5])
        # B = Diagonal([0.9; 0.5])
        # dist = Normal()
        # scaling = 0.0
        # This was generated on the R ScoreDrivenModels package where Normal is parametrized with sigma^2
        obs = [9.38058; 9.76652; 11.4661; 10.086; 10.1325; 11.6142; 10.5992; 8.6258; 9.18878; 9.52337; 11.204; 10.3993; 10.398; 10.1044; 9.49838; 11.6369; 10.6589; 7.83496; 11.1071; 9.41395; 8.86579; 9.75966; 8.99453; 9.23643; 9.36574; 8.34479; 11.0988; 10.1779; 8.85946; 11.3691; 10.5032; 9.6958; 10.8524; 10.8788; 10.8396; 10.7027; 10.5508; 9.94048; 9.72156; 9.6593; 9.37897; 9.80628; 8.86069; 12.3324; 12.0924; 8.35426; 9.47027; 9.49339; 10.7748; 9.91682; 10.2349; 9.97429; 9.96225; 11.1917; 9.75298; 11.481; 8.10715; 10.8092; 10.1411; 10.2142; 10.3529; 9.54254; 9.69554; 9.08409; 8.92653; 10.3253; 10.4355;10.0497; 10.8289; 12.0053; 9.23809; 7.24977; 11.9831; 8.95158; 9.17356; 11.1064; 9.68704; 8.80066; 10.2009; 9.86381; 10.0053; 10.343; 9.66946; 10.5743; 9.7966; 10.2974; 10.9768; 10.4428; 9.68766; 11.0573; 11.0428; 10.5897; 10.2385; 9.41362; 11.2815; 9.31883; 12.2613; 12.6291; 9.61586; 8.77585]
        parameters = [10.0 1.2214; 10.0 1.02894; 10.0 0.884717; 10.0 1.48609; 10.0 1.05055; 10.0 0.88589; 10.0 1.69009; 10.0 1.17998; 10.0 1.39493; 10.0 1.14381; 10.0 0.96738; 10.0 1.23123; 10.0 0.986465; 10.0 0.889888;10.0 0.81443; 10.0 0.839127; 10.0 1.75167; 10.0 1.21197; 10.0 2.49184; 10.0 1.53647; 10.0 1.1282; 10.0 1.21577; 10.0 0.960372; 10.0 1.09741; 10.0 1.02973; 10.0 0.963018; 10.0 1.72012; 10.0 1.34537; 10.0 1.00422; 10.0 1.19237; 10.0 1.39234; 10.0 1.06286; 10.0 0.906873; 10.0 1.00144; 10.0 1.04446; 10.0 1.04133; 10.0 0.988868; 10.0 0.92414; 10.0 0.828211; 10.0 0.801845; 10.0 0.799131; 10.0 0.868089; 10.0 0.810647; 10.0 1.15644; 10.0 3.00046; 10.0 2.14722; 10.0 1.72881; 10.0 1.17856; 10.0 0.986679; 10.0 0.995392; 10.0 0.860216; 10.0 0.811198; 10.0 0.775368; 10.0 0.758244; 10.0 1.19709; 10.0 0.953791; 10.0 1.49371; 10.0 1.9161; 10.0 1.29769; 10.0 0.984255; 10.0 0.863918; 10.0 0.829355; 10.0 0.834877; 10.0 0.808577; 10.0 1.00314; 10.0 1.14884; 10.0 0.944035; 10.0 0.879349; 10.0 0.807684; 10.0 0.956813; 10.0 2.40762; 10.0 1.4185; 10.0 3.88788; 10.0 2.18542; 10.0 1.44289; 10.0 1.16377; 10.0 1.20778; 10.0 0.965281; 10.0 1.22737; 10.0 0.961417; 10.0 0.84802; 10.0 0.792615; 10.0 0.79525; 10.0 0.794373; 10.0 0.85104; 10.0 0.803728; 10.0 0.793162; 10.0 1.03552; 10.0 0.918325; 10.0 0.847009; 10.0 1.10177; 10.0 1.1563; 10.0 0.997798; 10.0 0.872097; 10.0 0.887044; 10.0 1.28777; 10.0 1.0688;10.0 2.94284; 10.0 2.65624; 10.0 1.4224; 10.0 1.33584]

        gas = Model(1, 1, Normal, 0.0)
        gas.ω = [1.0; 0.1]
        gas.A[1] = convert(Matrix{Float64}, Diagonal([0.0; 0.5]))
        gas.B[1] = convert(Matrix{Float64}, Diagonal([0.9; 0.5]))
        param = score_driven_recursion(gas, obs)

        @test all(parameters - param .< 1e-3)
    end
end

ext: jl


text:
f(x)=cos(x)
g(x)=-sin(x)
h(x)=-cos(x)

x=2.0

x=optimize_newton(x,f,g,h)

@test x≈π rtol=1e-4


text:
export simplify

add_verbose_scope(:Simplify)

@doc Markdown.doc"""
    simplify(K::AnticNumberField; canonical::Bool = false) -> AnticNumberField, NfToNfMor

Tries to find an isomorphic field $L$ given by a "simpler" defining polynomial.
By default, "simple" is defined to be of smaller index, testing is done only
using a LLL-basis of the maximal order.

If `canonical` is set to `true`, then a canonical defining polynomial is found,
where canonical is using the definition of PARI's `polredabs`, which is described in
http://beta.lmfdb.org/knowledge/show/nf.polredabs.

Both versions require a LLL reduced basis for the maximal order.
"""
function simplify(K::AnticNumberField; canonical::Bool = false, cached::Bool = true, save_LLL_basis::Bool = true)
  Qx, x = PolynomialRing(FlintQQ, "x")

  if degree(K) == 1
    L = number_field(x - 1, cached = cached, check = false)[1]
    return L, hom(L, K, gen(K), check = false)
  end
  if canonical
    if !is_defining_polynomial_nice(K)
      K1, mK1 = simplify(K, cached = false, save_LLL_basis = false)
      K2, mK2 = simplify(K1, cached = cached, save_LLL_basis = save_LLL_basis, canonical = true)
      return K2, mK2*mK1
    end
    a, f1 = polredabs(K)
    f = Qx(f1)
    L = NumberField(f, cached = cached, check = false)[1]
    m = hom(L, K, a, check = false)
    return L, m
  end
  n = degree(K)
  OK = maximal_order(K)
  if isdefined(OK, :lllO)
    @vprint :Simplify 1 "LLL basis was already there\n"
    ZK = OK.lllO::typeof(OK)
  else
    b = _simplify(OK)
    if b != gen(K)
      @vprint :Simplify 1 "The basis of the maximal order contains a better primitive element\n"
      f1 = Qx(minpoly(representation_matrix(OK(b))))
      L1 = NumberField(f1, cached = cached, check = false)[1]
      #Before calling again the simplify on L1, we need to define the maximal order of L1
      mp = hom(L1, K, b, check = false)
      _assure_has_inverse_data(mp)
      B = basis(OK, K)
      BOL1 = Vector{nf_elem}(undef, degree(L1))
      for i = 1:degree(L1)
        BOL1[i] = mp\(B[i])
      end
      OL1 = NfOrd(BOL1, false)
      OL1.is_maximal = 1
      set_attribute!(L1, :maximal_order => OL1)
      @vprint :Simplify 3 "Trying to simplify $(L1.pol)\n"
      L2, mL2 = simplify(L1, cached = cached, save_LLL_basis = save_LLL_basis)
      h = mL2 * mp
      return L2, mL2*mp
    end
    prec = 100 + 25*div(n, 3) + Int(round(log(abs(discriminant(OK)))))
    @vtime :Simplify 3 ZK = lll(OK, prec = prec)
    OK.lllO = ZK
  end
  @vtime :Simplify 3 a = _simplify(ZK)
  if a == gen(K)
    f = K.pol
  else
    @vtime :Simplify 3 f = Qx(minpoly(representation_matrix(OK(a))))
  end
  L = NumberField(f, cached = cached, check = false)[1]
  m = hom(L, K, a, check = false)
  if save_LLL_basis
    _assure_has_inverse_data(m)
    B = basis(ZK, K)
    BOL = Vector{nf_elem}(undef, degree(L))
    for i = 1:degree(L)
      BOL[i] = m\(B[i])
    end
    OL = NfOrd(BOL, false)
    if isdefined(ZK, :disc)
      OL.disc = ZK.disc
      if is_defining_polynomial_nice(L)
        OL.index = root(divexact(numerator(discriminant(L.pol)), OL.disc), 2)
      end
    end
    OL.is_maximal = 1
    set_attribute!(L, :maximal_order => OL)
  end
  if cached
    embed(m)
    embed(inv(m))
  end
  return L, m
end

function _simplify(O::NfAbsOrd)
  K = nf(O)

  B = basis(O, K, copy = false)
  nrep = min(3, degree(K))
  Bnew = elem_type(K)[]
  for i = 1:length(B)
    push!(Bnew, B[i])
    for j = 1:nrep
      push!(Bnew, B[i]+B[j])
      push!(Bnew, B[i]-B[j])
    end
  end
  #First, we search for elements that are primitive using block systems in the simple case.
  B1 = _sieve_primitive_elements(Bnew)

  #Now, we select the one of smallest T2 norm
  a = primitive_element(K)
  d = denominator(a, O)
  if !isone(d)
    a *= d
  end
  I = t2(a)
  for i = 1:length(B1)
    t2n = t2(B1[i])
    if t2n < I
      a = B1[i]
      I = t2n
    end
  end
  return a
end

function primitive_element(K::AnticNumberField)
  return gen(K)
end

function _sieve_primitive_elements(B::Vector{NfAbsNSElem})
  K = parent(B[1])
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  pols = [Zx(to_univariate(Globals.Qx, x)) for x in K.pol]
  p, d = _find_prime(pols)
  F = FlintFiniteField(p, d, "w", cached = false)[1]
  Fp = GF(p, cached = false)
  Fpt = PolynomialRing(Fp, ngens(K))[1]
  Ft = PolynomialRing(F, "t", cached = false)[1]
  rt = Vector{Vector{fq_nmod}}(undef, ngens(K))
  for i = 1:length(pols)
    rt[i] = roots(pols[i], F)
  end
  rt_all = Vector{Vector{fq_nmod}}(undef, degree(K))
  ind = 1
  it = cartesian_product_iterator([1:degrees(K)[i] for i in 1:ngens(K)], inplace = true)
  for i in it
    rt_all[ind] = fq_nmod[rt[j][i[j]] for j = 1:length(rt)]
    ind += 1
  end
  indices = Int[]
  for i = 1:length(B)
    if length(vars(data(B[i]))) != ngens(K)
      continue
    end
    if isone(denominator(B[i]))
      continue
    end
    if _is_primitive_via_block(B[i], rt_all, Fpt)
      push!(indices, i)
    end
  end
  return B[indices]
end

function _is_primitive_via_block(el::NfAbsNSElem, rt::Vector{Vector{fq_nmod}}, Rt::MPolyRing)
  K = parent(el)
  fR = map_coefficients(base_ring(Rt), data(el), parent = Rt)
  s = Set{fq_nmod}()
  for x in rt
    val = evaluate(fR, x)
    if val in s
      return false
    end
    push!(s, val)
    if length(s) > div(degree(K), 2)
      return true
    end
  end
  error("Something went wrong")
end

function _block(el::NfAbsNSElem, rt::Vector{Vector{fq_nmod}}, R::GaloisField)
  fR = map_coefficients(R, data(el))
  s = fq_nmod[evaluate(fR, x) for x in rt]
  b = Vector{Int}[]
  a = BitSet()
  i = 0
  n = length(rt)
  while i < n
    i += 1
    if i in a
      continue
    end
    z = s[i]
    push!(b, findall(x->s[x] == z, 1:n))
    for j in b[end]
      push!(a, j)
    end
  end
  return b
end

function AbstractAlgebra.map_coefficients(F::GaloisField, f::fmpq_mpoly; parent = PolynomialRing(F, nvars(parent(f)), cached = false)[1])
  dF = denominator(f)
  d = F(dF)
  if iszero(d)
    error("Denominator divisible by p!")
  end
  m = inv(d)
  ctx = MPolyBuildCtx(parent)
  for x in zip(coefficients(f), exponent_vectors(f))
    el = numerator(x[1]*dF)
    push_term!(ctx, F(el)*m, x[2])
  end
  return finish(ctx)
end

function _sieve_primitive_elements(B::Vector{nf_elem})
  K = parent(B[1])
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  f = Zx(K.pol*denominator(K.pol))
  a = gen(K)*denominator(K.pol)

  p, d = _find_prime(fmpz_poly[f])

  F = FlintFiniteField(p, d, "w", cached = false)[1]
  Ft = PolynomialRing(F, "t", cached = false)[1]
  ap = zero(Ft)
  fit!(ap, degree(K)+1)
  rt = roots(f, F)

  n = degree(K)
  indices = Int[]
  for i = 1:length(B)
    if isone(denominator(B[i]))
      continue
    end
    b = _block(B[i], rt, ap)
    if length(b) == n
      push!(indices, i)
    end
  end
  return B[indices]

end

 #a block is a partition of 1:n
 #given by the subfield of parent(a) defined by a
 #the embeddings used are in R
 #K = parent(a)
 # then K has embeddings into the finite field (parent of R[1])
 # given by the roots (in R) of the minpoly of K
 #integers in 1:n are in the same block iff a(R[i]) == a(R[j])
 #the length of such a block (system) is the degree of Q(a):Q, the length
 # of a block is the degree K:Q(a)
 # a is primitive iff the block system has length n
function _block(a::nf_elem, R::Vector{fq_nmod}, ap::fq_nmod_poly)
  # TODO:
  # Maybe this _tmp business has to be moved out of this function too
  _R = GF(Int(characteristic(base_ring(ap))), cached = false)
  _Ry, _ = PolynomialRing(_R, "y", cached = false)
  _tmp = _Ry()
  nf_elem_to_gfp_poly!(_tmp, a, false) # ignore denominator
  set_length!(ap, length(_tmp))
  for i in 0:(length(_tmp) - 1)
    setcoeff!(ap, i, base_ring(ap)(_get_coeff_raw(_tmp, i)))
  end
#  ap = Ft(Zx(a*denominator(a)))
  s = fq_nmod[evaluate(ap, x) for x = R]
  b = Vector{Int}[]
  a = BitSet()
  i = 0
  n = length(R)
  while i < n
    i += 1
    if i in a
      continue
    end
    z = s[i]
    push!(b, findall(x->s[x] == z, 1:n))
    for j in b[end]
      push!(a, j)
    end
  end
  return b
end

#given 2 block systems b1, b2 for elements a1, a2, this computes the
#system for Q(a1, a2), the compositum of Q(a1) and Q(a2) as subfields of K
function _meet(b1::Vector{Vector{Int}}, b2::Vector{Vector{Int}})
  b = Vector{Int}[]
  for i=b1
    for j = i
      for h = b2
        if j in h
          s = intersect(i, h)
          if !(s in b)
            push!(b, s)
          end
        end
      end
    end
  end
  return b
end

function _find_prime(v::Vector{fmpz_poly})
  p = 2^10
  total_deg = prod(degree(x) for x in v)
  n_attempts = min(total_deg, 10)
  candidates = Vector{Tuple{Int, Int}}(undef, n_attempts)
  i = 1
  polsR = Vector{gfp_poly}(undef, length(v))
  while i < n_attempts+1
    p = next_prime(p)
    R = GF(p, cached=false)
    Rt = PolynomialRing(R, "t", cached = false)[1]
    found_bad = false
    for j = 1:length(v)
      fR = map_coefficients(R, v[j], parent = Rt)
      if degree(fR) != degree(v[j]) || !is_squarefree(fR)
        found_bad = true
        break
      end
      polsR[j] = fR
    end
    if found_bad
      continue
    end
    d = 1
    for j = 1:length(polsR)
      fR = polsR[j]
      FS = factor_shape(fR)
      d1 = lcm(Int[x for (x, v) in FS])
      d = lcm(d1, d)
    end
    if d <= total_deg^2
      candidates[i] = (p, d)
      i += 1
    end
  end
  res = candidates[1]
  for j = 2:n_attempts
    if candidates[j][2] < res[2]
      res = candidates[j]
    end
  end
  return res[1], res[2]
end

function polredabs(K::AnticNumberField)
  #intended to implement
  # http://beta.lmfdb.org/knowledge/show/nf.polredabs
  #as in pari
  #TODO: figure out the separation of T2-norms....
  ZK = lll(maximal_order(K))
  I = index(ZK)^2
  D = discriminant(ZK)
  B = basis(ZK, copy = false)
  Zx = FlintZZ["x"][1]
  f = Zx(K.pol)
  p, d = _find_prime(fmpz_poly[f])

  F = FlintFiniteField(p, d, "w", cached = false)[1]
  Ft = PolynomialRing(F, "t", cached = false)[1]
  ap = zero(Ft)
  fit!(ap, degree(K)+1)
  rt = roots(f, F)

  n = degree(K)

  b = _block(B[1].elem_in_nf, rt, ap)
  i = 2
  while length(b) < degree(K)
    bb = _block(B[i].elem_in_nf, rt, ap)
    b = _meet(b, bb)
    i += 1
  end
  i -= 1
#  println("need to use at least the first $i basis elements...")
  pr = 100
  old = precision(BigFloat)
  local E
  while true
    try
      setprecision(BigFloat, pr)
      E = enum_ctx_from_ideal(ideal(ZK, 1), zero_matrix(FlintZZ, 1, 1), prec = pr, TU = BigFloat, TC = BigFloat)

      if E.C[end] + 0.0001 == E.C[end]  # very very crude...
        pr *= 2
        continue
      end
      break
    catch e
      if isa(e, InexactError) || isa(e, LowPrecisionLLL) || isa(e, LowPrecisionCholesky)
        pr *= 2
        continue
      end
      rethrow(e)
    end
  end

  l = zeros(FlintZZ, n)
  l[i] = 1

  scale = 1.0
  enum_ctx_start(E, matrix(FlintZZ, 1, n, l), eps = 1.01)

  a = gen(K)
  all_a = nf_elem[a]
  la = length(a)*BigFloat(E.t_den^2)
  Ec = BigFloat(E.c//E.d)
  eps = BigFloat(E.d)^(1//2)

  found_pe = false
  first = true
  while !found_pe
    while first || enum_ctx_next(E)
      first = false
#      @show E.x
      M = E.x*E.t
      q = elem_from_mat_row(K, M, 1, E.t_den)
      bb = _block(q, rt, ap)
      if length(bb) < n
        continue
      end
      found_pe = true
#  @show    llq = length(q)
#  @show sum(E.C[i,i]*(BigFloat(E.x[1,i]) + E.tail[i])^2 for i=1:E.limit)/BigInt(E.t_den^2)
      lq = Ec - (E.l[1] - E.C[1, 1]*(BigFloat(E.x[1,1]) + E.tail[1])^2)
#      @show lq/E.t_den^2

      if lq < la + eps
        if lq > la - eps
          push!(all_a, q)
#          @show "new one", q, minpoly(q), bb
        else
          a = q
          all_a = nf_elem[a]
          if lq/la < 0.8
#            @show "re-init"
            enum_ctx_start(E, E.x, eps = 1.01)  #update upperbound
            first = true
            Ec = BigFloat(E.c//E.d)
          end
          la = lq
  #        @show Float64(la/E.t_den^2)
        end
      end
    end
    scale *= 2
    enum_ctx_start(E, matrix(FlintZZ, 1, n, l), eps = scale)
    first = true
    Ec = BigFloat(E.c//E.d)
  end

  setprecision(BigFloat, old)
  #try to find the T2 shortest element
  #the precision management here needs a revision once we figure out
  #how....
  #examples that require this are Gunters:
  #=
  die drei Polynome

[ 10834375376002294480896, x^18 - x^16 - 6*x^14 - 4*x^12 - 4*x^10 + 2*x^8 +
6*x^6 - 4*x^4 + 3*x^2 - 1 ],
[ 10834375376002294480896, x^18 - 3*x^16 + 4*x^14 - 6*x^12 - 2*x^10 + 4*x^8 +
4*x^6 + 6*x^4 + x^2 - 1 ],
[ 10834375376002294480896, x^18 + x^16 - x^14 - 8*x^12 - 3*x^8 + 27*x^6 -
25*x^4 + 8*x^2 - 1 ],

werden alle als 'canonical' ausgegeben, obwohl sie isomorphe
K"orper definieren ??
=#
  sort!(all_a, lt = (a,b) -> length(a) < length(b))
  i = length(all_a)
  la1 = length(all_a[1])
  while i >= 1 && la1 <= length(all_a[i]) - 1e-10
    i -= 1
  end
  all_a = all_a[1:i]

  all_f = Tuple{nf_elem, fmpq_poly}[(x, minpoly(x)) for x=all_a]
  all_d = fmpq[abs(discriminant(x[2])) for x= all_f]
  m = minimum(all_d)
  L1 = Tuple{nf_elem, fmpq_poly}[]
  for i = 1:length(all_f)
    if all_d[i] == m
      push!(L1, all_f[i])
    end
  end
  L2 = Tuple{nf_elem, fmpq_poly}[minQ(x) for x=L1]
  if length(L2) == 1
    return L2[1]
  end
  L3 = sort(L2, lt = il)

  return L3[1]
end

function Q1Q2(f::PolyElem)
  q1 = parent(f)()
  q2 = parent(f)()
  g = gen(parent(f))
  for j=0:degree(f)
    if isodd(j)
      q2 += coeff(f, j)*g^div(j, 2)
    else
      q1 += coeff(f, j)*g^div(j, 2)
    end
  end
  return q1, q2
end

function minQ(A::Tuple{nf_elem, fmpq_poly})
  a = A[1]
  f = A[2]
  q1, q2 = Q1Q2(f)
  if leading_coefficient(q1)>0 && leading_coefficient(q2) > 0
    return (-A[1], f(-gen(parent(f)))*(-1)^degree(f))
  else
    return (A[1], f)
  end
end

function int_cmp(a, b)
  if a==b
    return 0
  end
  if abs(a) == abs(b)
    if a>b
      return 1
    else
      return -1
    end
  end
  return cmp(abs(a), abs(b))
end

function il(F, G)
  f = F[2]
  g = G[2]
  i = degree(f)
  while i>0 && int_cmp(coeff(f, i), coeff(g, i))==0
    i -= 1
  end
  return int_cmp(coeff(f, i), coeff(g, i))<0
end


text:
using Documenter, Pkg, PlanktonIndividuals, Plots
import PlutoSliderServer
Pkg.precompile()

examples = ["vertical_2D_example.jl", "horizontal_2D_example.jl", "surface_mixing_3D_example.jl",
            "0D_experiment.jl", "global_ocean_2D_example.jl", "global_ocean_3D_example.jl"]

model_description = [
    # "Table Of Content" => "model_description.md",
    "Phytoplankton Physiology" => "phyto_equations.md",
    "Biogeochemistry" => "bgc_equations.md",
]

makedocs(;
         modules = [PlanktonIndividuals],
         format = Documenter.HTML(collapselevel = 1, mathengine = MathJax3()),
         pages = [
             "Home" => "index.md",
             "Model Description" => model_description,
             "Model Configuration" => "model_setup.md",
             "Model Simulation" => "model_run.md",
             "GPU Support" => "gpu_support.md",
             "Examples"  => "examples.md",
             "Benchmarks"  => "benchmarks.md",
             "Library" => "library.md",
             "Function index" => "function_index.md"
             ],
         repo="https://github.com/JuliaOcean/PlanktonIndividuals.jl/blob/{commit}{path}#L{line}",
         sitename = "PlanktonIndividuals.jl",
         authors="ZhenWu <zhenwu@mit.edu>",
)

for i in examples
    fil_in=joinpath(@__DIR__,"..", "examples",i)
    fil_out=joinpath(@__DIR__,"build", "examples",i[1:end-2]*"html")
    PlutoSliderServer.export_notebook(fil_in)
    mv(fil_in[1:end-2]*"html",fil_out)
    cp(fil_in,fil_out[1:end-4]*"jl")
end

deploydocs(repo="github.com/JuliaOcean/PlanktonIndividuals.jl.git")


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Statistics
using ArraysOfArrays, ElasticArrays, StatsBase
using Distributions
using StableRNGs

@testset "autocorr" begin
    stblrng() = StableRNG(789990641)

    # Primitive MCMC of flat distribution between [0, 0, 0] and [10, 5, 8]:
    rng = stblrng()
    v = nestedview(ElasticArray{Float64, 2}(undef, 3, 0))
    push!(v, [0, 0, 0])
    for i in 1:(10^4 - 1)
        push!(v, clamp.(last(v) .+ randn(rng, 3), [0, 0, 0], [10, 5, 8]))
    end

    v1 = flatview(v)[1, :]

    @testset "BAT._ac_next_pow_two" begin
        @test @inferred(BAT._ac_next_pow_two(0)) == 1
        @test @inferred(BAT._ac_next_pow_two(1)) == 1
        @test BAT._ac_next_pow_two(7) == 8
        @test BAT._ac_next_pow_two(8) == 8
    end

    @testset "bat_integrated_autocorr_len" begin
        dist = product_distribution([Normal() for i in 1:3])
        data = @inferred(nestedview(rand(dist, 10^4)))
        @test isapprox(@inferred(bat_integrated_autocorr_len(data)).result, ones(length(dist)), rtol=0.2)
    end

    @testset "BAT.fft_autocov" begin
        @test @inferred(BAT.fft_autocov(v)) isa ArrayOfSimilarArrays{Float64,1,1}
        @test flatview(BAT.fft_autocov(v)[1:20]) ≈ StatsBase.autocov(flatview(v)', 0:19)'

        @test @inferred(BAT.fft_autocov(v1)) isa Vector{Float64}
        @test BAT.fft_autocov(v1)[1:20] ≈ StatsBase.autocov(v1, 0:19)
    end

    @testset "BAT.fft_autocor" begin
        @test @inferred(BAT.fft_autocor(v)) isa ArrayOfSimilarArrays{Float64,1,1}
        @test flatview(BAT.fft_autocor(v)[1:20]) ≈ StatsBase.autocor(flatview(v)', 0:19)'

        @test @inferred(BAT.fft_autocor(v1)) isa Vector{Float64}
        @test BAT.fft_autocor(v1)[1:20] ≈ StatsBase.autocor(v1, 0:19)
    end
end


text:
"""
    ZeroOneNormalization(low::T, high::T)

This is a functor that normalizes a vector to be in the range [0,1]. Initial upper and lower bounds for each element are needed.

See also: [`PosNegNormalization`](@ref), [`GaussianNormalization`](@ref)

# Examples
```jldoctest
julia> low = [0.0, -1.0];

julia> high = [3.0, 0.5];

julia> nrm = ZeroOneNormalization(low, high);

julia> x = [1.0, 0.0];

julia> feats = nrm(x)
2-element Array{Float64,1}:
 0.3333333333333333
 0.6666666666666666

julia> y = zero(x);  # create buffer to prevent allocations

julia> feats = nrm(y, x);  # no allocation return

julia> feats = nrm(y, x, fit=true);  # update upper and lower bounds if x is outside the range

```
"""
struct ZeroOneNormalization{T} <: Any where {T}
    low::T
    high::T

    function ZeroOneNormalization(ranges)
        low = ranges[:, 1]
        high = ranges[:, 2]
        return new{typeof(low)}(low, high)
    end

    function ZeroOneNormalization(low::T, high::T) where {T}
        return new{T}(deepcopy(low), deepcopy(high))
    end
end

"""
    PosNegNormalization(low::T, high::T)

This is a functor that normalizes a vector to be in the range [-1,1]. Initial upper and lower bounds for each element are needed.

See also: [`ZeroOneNormalization`](@ref), [`GaussianNormalization`](@ref)

# Examples
```jldoctest
julia> low = [0.0, -1.0];

julia> high = [3.0, 0.5];

julia> nrm = PosNegNormalization(low, high);

julia> x = [1.0, 0.0];

julia> feats = nrm(x)
2-element Array{Float64,1}:
 -0.33333333333333337
  0.33333333333333326

julia> y = zero(x);  # create buffer to prevent allocations

julia> feats = nrm(y, x);  # no allocation return

julia> feats = nrm(y, x, fit=true);  # update upper and lower bounds if x is outside the range

```
"""
struct PosNegNormalization{T} <: Any where {T}
    low::T
    high::T

    function PosNegNormalization(ranges)
        low = ranges[:, 1]
        high = ranges[:, 2]
        return new{typeof(low)}(low, high)
    end

    function PosNegNormalization(low::T, high::T) where {T}
        return new{T}(deepcopy(low), deepcopy(high))
    end
end

"""
    GaussianNormalization(num_features::Int[, weight])

This is a functor that normalizes each element of a vector to be centered around a mean and have a variance of one.
OnlineStats.KahanVariance is used for computing and tracking the mean and variance. Likewise any OnlineStats.Weight can
be used to weigh each sample.

See also: [`ZeroOneNormalization`](@ref), [`PosNegNormalization`](@ref)

# Examples
```jldoctest
julia> nrm = GaussianNormalization(2, 1e-4);

julia> x = [1.0, 0.0];

julia> feats = nrm(x)
2-element Array{Float64,1}:
 1.0
 0.0

julia> y = zero(x);  # create buffer to prevent allocations

julia> feats = nrm(y, x);  # no allocation return

julia> feats = nrm(y, x, fit=true)  # update mean and variance using x
2-element Array{Float64,1}:
 0.9999
 0.0

```
"""
struct GaussianNormalization{T} <: Any where {T}
    v::T

    function GaussianNormalization(num_features::Int)
        v = num_features*KahanVariance()
        return new{typeof(v)}(v)
    end

    function GaussianNormalization(num_features::Int, weight::T) where {T<:Real}
        v = num_features*KahanVariance(weight=ExponentialWeight(weight))
        return new{typeof(v)}(v)
    end

    function GaussianNormalization(num_features::Int, weight::OnlineStats.Weight)
        v = num_features*KahanVariance(weight=weight)
        return new{typeof(v)}(v)
    end
end

function (f::ZeroOneNormalization)(x; fit=false)
    if fit
        @. f.low = min(x, f.low)
        @. f.high = max(x, f.high)
    end
    y = zero(x)
    @. y = (x - f.low) / (f.high - f.low)
    return y
end

function (f::ZeroOneNormalization)(buff::T, x::T; fit=false) where {T}
    if fit
        @. f.low = min(x, f.low)
        @. f.high = max(x, f.high)
    end
    @. buff = (x - f.low) / (f.high - f.low)
    return buff
end

function (f::PosNegNormalization)(x; fit=false)
    if fit
        @. f.low = min(x, f.low)
        @. f.high = max(x, f.high)
    end
    y = zero(x)
    T2 = eltype(x)
    a = convert(T2, 2.0)
    b = convert(T2, 1.0)
    @. y = a * (x - f.low) / (f.high - f.low) - b
    return y
end

function (f::PosNegNormalization)(buff::T, x::T; fit=false) where {T}
    if fit
        @. f.low = min(x, f.low)
        @. f.high = max(x, f.high)
    end
    T2 = eltype(x)
    a = convert(T2, 2.0)
    b = convert(T2, 1.0)
    @. buff = a * (x - f.low) / (f.high - f.low) - b
    return buff
end

function (f::GaussianNormalization)(x; fit=false)
    if fit
        fit!(f.v, x)
    end
    y = zero(x)
    T2 = eltype(x)
    @. y = (x - convert(T2, mean(f.v.stats))) / convert(T2, std(f.v.stats))
    return y
end

function (f::GaussianNormalization)(buff::T, x::T; fit=false) where {T}
    if fit
        fit!(f.v, x)
    end
    T2 = eltype(x)
    @. buff = (x - convert(T2, mean(f.v.stats))) / convert(T2, std(f.v.stats))
    return buff
end


text:
"""
Calculates the exact boundingbox of a Scene/Plot, without considering any transformation
"""
function raw_boundingbox(x::Atomic)
    bb = data_limits(x)
end
rootparent(x) = rootparent(parent(x))
rootparent(x::Scene) = x

function raw_boundingbox(x::Annotations)
    bb = raw_boundingbox(x.plots)
    inv(modelmatrix(rootparent(x))) * bb
end
function raw_boundingbox(x::Combined)
    raw_boundingbox(x.plots)
end
function boundingbox(x)
    raw_boundingbox(x)
end

function combined_modelmatrix(x)
    m = Mat4f0(I)
    while true
        m = modelmatrix(x) * m
        if parent(x) !== nothing && parent(x) isa Combined
            x = parent(x)
        else
            break
        end
    end
    m
end

function modelmatrix(x)
    t = transformation(x)
    transformationmatrix(t.translation[], t.scale[], t.rotation[])
end

function boundingbox(x::Atomic)
    bb = raw_boundingbox(x)
    combined_modelmatrix(x) * bb
end

boundingbox(scene::Scene) = raw_boundingbox(scene)
function raw_boundingbox(scene::Scene)
    if scene[Axis] !== nothing
        return raw_boundingbox(scene[Axis])
    elseif scene.limits[] !== automatic
        return scene_limits(scene)
    elseif cameracontrols(scene) == EmptyCamera()
        # Empty camera means this is a parent scene that itself doesn't display anything
        return raw_boundingbox(scene.children)
    else
        plots = plots_from_camera(scene)
        children = filter(scene.children) do child
            child.camera == scene.camera
        end
        return raw_boundingbox([plots; children])
    end
end

function raw_boundingbox(plots::Vector)
    isempty(plots) && return FRect3D()
    plot_idx = iterate(plots)
    bb = FRect3D()
    while plot_idx !== nothing
        plot, idx = plot_idx
        plot_idx = iterate(plots, idx)
        # isvisible(plot) || continue
        bb2 = boundingbox(plot)
        isfinite(bb) || (bb = bb2)
        isfinite(bb2) || continue
        bb = union(bb, bb2)
    end
    bb
end

function project_widths(matrix, vec)
    pr = project(matrix, vec)
    zero = project(matrix, zeros(typeof(vec)))
    pr - zero
end

function boundingbox(x::Text, text::String)
    position = to_value(x[:position])
    @get_attribute x (textsize, font, align, rotation)
    bb = boundingbox(text, position, textsize, font, align, rotation, modelmatrix(x))
    pm = inv(transformationmatrix(parent(x))[])
    wh = widths(bb)
    whp = project_widths(pm, wh)
    aoffset = whp .* to_ndim(Vec3f0, align, 0f0)
    return FRect3D(minimum(bb) .- aoffset, whp)
end

boundingbox(x::Text) = boundingbox(x, to_value(x[1]))

function boundingbox(
        text::String, position, textsize;
        font = "default", align = (:left, :bottom), rotation = 0.0
    )
    boundingbox(
        text, position, textsize,
        to_font(font), to_align(align), to_rotation(rotation)
    )

end

# function boundingbox(
#         text::String, position, textsize, fonts,
#         align, rotation, model = Mat4f0(I)
#     )
#     isempty(text) && return FRect3D()
#     pos_per_char = !isa(position, VecTypes)
#
#     start_pos = Vec(pos_per_char ? first(position) : position)
#     start_pos3d = project(model, to_ndim(Vec3f0, start_pos, 0.0))
#     bb = FRect3D(start_pos3d, Vec3f0(0))
#
#     if pos_per_char
#         broadcast_foreach(position, textsize, fonts, collect(text)) do pos, scale, font, char
#             rect, extent = FreeTypeAbstraction.metrics_bb(char, font, scale)
#             bb = union(FRect3D(rect) + to_ndim(Vec3f0, pos, 0.0), bb)
#             @show pos scale
#         end
#     else
#         y_advance = 0.0
#         line_advance = FreeTypeAbstraction.get_extent(fonts, 'x').advance[2]
#         for line in split(text, r"(\r\n|\r|\n)")
#             rectangles = FreeTypeAbstraction.glyph_rects(line, fonts, textsize)
#             bb2d = reduce(union, rectangles)
#             bb2d = bb2d + Vec2f0(0, y_advance)
#             bb = union(bb, FRect3D(bb2d))
#             y_advance += line_advance
#             @show y_advance
#         end
#     end
#     return bb
# end

function boundingbox(
        text::String, position, textsize, font,
        align, rotation, model = Mat4f0(I)
    )
    atlas = get_texture_atlas()
    N = length(text)
    ctext_state = iterate(text)
    ctext_state === nothing && return FRect3D()
    pos_per_char = !isa(position, VecTypes)
    start_pos = Vec(pos_per_char ? first(position) : position)
    start_pos2D = to_ndim(Point2f0, start_pos, 0.0)
    last_pos = Point2f0(0, 0)
    start_pos3d = project(model, to_ndim(Vec3f0, start_pos, 0.0))
    bb = FRect3D(start_pos3d, Vec3f0(0))
    broadcast_foreach(1:N, rotation, font, textsize) do i, rotation, font, scale
        c, text_state = ctext_state
        ctext_state = iterate(text, text_state)
        # TODO fix center + align + rotation
        if c != '\r'
            pos = if pos_per_char
                to_ndim(Vec3f0, position[i], 0.0)
            else
                last_pos = calc_position(last_pos, Point2f0(0, 0), atlas, c, font, scale)
                start_pos3d .+ (rotation * to_ndim(Vec3f0, last_pos, 0.0))
            end
            s = glyph_scale!(atlas, c, font, scale)
            srot = rotation * to_ndim(Vec3f0, s, 0.0)
            bb = GeometryTypes.update(bb, pos)
            bb = GeometryTypes.update(bb, pos .+ srot)
        end
    end
    return bb
end


text:
function findfst(ind::Array{Int}, idx::Int)
    x = falses(length(ind))
    for k in eachindex(ind)
        if ind[k] == idx
            x[k] = true
            return x
        end
    end
end


text:
"""
    DeviceBuffer{T,A}
A memory address that refers to data of type `T` that is accessible from the GPU. It is the
on-device counterpart of `MtlDAdrv.MtlPtr`, additionally keeping track of the address space
`A` where the data resides (shared, global, constant, etc). This information is used to
provide optimized implementations of operations such as `unsafe_load` and `unsafe_store!.`
"""
DeviceBuffer

# constructors
DeviceBuffer{T,A}(x::Union{Int,UInt,DeviceBuffer,DevicePtr}) where {T,A}            = Base.bitcast(DeviceBuffer{T,A}, x)
DeviceBuffer{T,A}(ptr::MtlBuffer{T})                         where {T,A}            = Base.bitcast(DeviceBuffer{T,A}, ptr.handle)
DeviceBuffer{T}(ptr::MtlBuffer{T})                           where {T}                 = Base.bitcast(DeviceBuffer{T,AS.Device}, ptr.handle)
DeviceBuffer(ptr::MtlBuffer{T})                              where {T}                 = Base.bitcast(DeviceBuffer{T,AS.Device}, ptr.handle)

## conversions
Base.convert(::Type{DeviceBuffer{T,A}}, x::Union{Int,UInt}) where {T,A} = DeviceBuffer{T,A}(x)

# between host and device pointers
Base.convert(::Type{MtlBuffer{T}},      p::DeviceBuffer)  where {T}                   = MtlBuffer{T}(Base.bitcast(MTL.MTLBuffer, p))
Base.convert(::Type{DeviceBuffer{T,A}}, p::MtlBuffer)     where {T,A}                 = Base.bitcast(DeviceBuffer{T,A}, p.handle)
Base.convert(::Type{DeviceBuffer{T}},   p::MtlBuffer)     where {T}                   = Base.bitcast(DeviceBuffer{T,AS.Generic}, p.handle)

# between CPU pointers, for the purpose of working with `ccall`
Base.unsafe_convert(::Type{MTL.MTLBuffer}, x::DeviceBuffer{T}) where {T} = reinterpret(MTL.MTLBuffer, x)
Base.unsafe_convert(::Type{MTL.MTLResource}, x::DeviceBuffer{T}) where {T} = reinterpret(MTL.MTLResource, x)

# between device pointers
Base.convert(::Type{<:DeviceBuffer}, p::DeviceBuffer)                         = throw(ArgumentError("cannot convert between incompatible device pointer types"))
Base.convert(::Type{DeviceBuffer{T,A}}, p::DeviceBuffer{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DeviceBuffer{T,A}}, p::DeviceBuffer) where {T,A}   = Base.bitcast(DeviceBuffer{T,A}, p)
## identical addrspaces
Base.convert(::Type{DeviceBuffer{T,A}}, p::DeviceBuffer{U,A}) where {T,U,A} = Base.unsafe_convert(DeviceBuffer{T,A}, p)
## convert to & from generic
Base.convert(::Type{DeviceBuffer{T,AS.Generic}}, p::DeviceBuffer)               where {T}     = Base.unsafe_convert(DeviceBuffer{T,AS.Generic}, p)
Base.convert(::Type{DeviceBuffer{T,A}},          p::DeviceBuffer{U,AS.Generic}) where {T,U,A} = Base.unsafe_convert(DeviceBuffer{T,A}, p)
Base.convert(::Type{DeviceBuffer{T,AS.Generic}}, p::DeviceBuffer{T,AS.Generic}) where {T}     = p  # avoid ambiguities
## unspecified, preserve source addrspace
Base.convert(::Type{DeviceBuffer{T}}, p::DeviceBuffer{U,A}) where {T,U,A} = Base.unsafe_convert(DeviceBuffer{T,A}, p)

Base.pointer(buf::DeviceBuffer{T,A}) where{T,A} = reinterpret(Core.LLVMPtr{T,A}, buf)
## memory operations

## new set methods
MTL.set_buffer!(cce::MtlArgumentEncoder, buf::DeviceBuffer, offset::Integer, index::Integer) =
    MTL.mtArgumentEncoderSetBufferOffsetAtIndex(cce, buf, offset, index-1)
MTL.set_buffers!(cce::MtlArgumentEncoder, bufs::Vector{<:DeviceBuffer},
             offsets::Vector{Int}, indices::UnitRange{Int}) =
    MTL.mtArgumentSetBuffersOffsetsWithRange(cce, handle_array(bufs), offsets, indices .- 1)

MTL.use!(cce::MtlComputeCommandEncoder, buf::DeviceBuffer, mode::MTL.MtResourceUsage=ReadWriteUsage) =
    MTL.mtComputeCommandEncoderUseResourceUsage(cce, buf, mode)

MTL.use!(cce::MtlComputeCommandEncoder, buf::Vector{DeviceBuffer}, mode::MTL.MtResourceUsage=ReadWriteUsage) =
    MTL.mtComputeCommandEncoderUseResourceCountUsage(cce, handle_array(buf), length(buf), mode)


text:
# 0 -> empty space, 1 -> O, 2 -> X
str = [' ', 'O', 'X']

function playConnect4(;ai::Int=1,method::String="minimaxAB",abDepth::Int=6)
    # initialise board and display
    M = c4_initialise()
    c4_boardDisplay(M)

    # which player(s) are AI?
    if ai == 0
        isAI = [false false]
    elseif ai == 1
        isAI = rand()<0.5 ? [false true] : [true false]
    elseif ai == 2
        isAI = [true true]
    end

    # randomly pick starting player (sets symbols)
    player = rand((1,2))

    while true # take turns until the game ends
        # take a turn, either by human or AI
        if isAI[player]
            M = c4_aiTurn!(M,player,method,abDepth)
        else
            M = c4_takeTurn!(M,player)
        end

        c4_boardDisplay(M) # show the board after the turn

        player = 3-player # swap player
        W = c4_checkWinner(M) # check if someone has won

        if W!=0 # someone has won
            print("$(str[W+1])'s win\n")
            return nothing
        elseif c4_gameOver(M) # board filled but no winner -> draw
            print("It's a draw\n")
            return nothing
        end
    end
end

function c4_initialise()
    # empty board
    return zeros(Int8,6,7)
end

function c4_boardDisplay(M)
    # Displays the board with symbols
    print(" 1 2 3 4 5 6 7\n")
    for row in 1:6
        print("|")
        for column in 1:7
            print("$(str[M[row,column]+1])|")
        end
        print("\n")
    end

    print("---------------\n")
    print("Board eval: $(c4_scoreHelper(M,1))\n")
    print("---------------\n")
end

function c4_takeTurn!(M,player)
    # asks for player input and updates the game state
    col = c4_getMove(M)
    # ids = count(!iszero,M[:,col])
    # M[6-ids,col] = player

    return c4_makeMove!(M,col,player)
end

function c4_getMove(M)
    # gets a move from a human player
    # TODO: More input validation

    while true # repeat until valid location selected
        print("Enter the column for your move: ")
        x = parse(Int,readline())

        if c4_movePossible(M,x) # must be an empty location
            return x
        end
    end
end

function c4_aiTurn!(M,player,method,abDepth)
    # work out the best move and then plays it
    moves = 1:7

    if isequal(method,"minimaxAB")
        col = getBestMoveAB(M,player,c4_score,moves,
            c4_movePossible,c4_makeMove!,c4_gameOver,abDepth)
    else
        col = getBestMove(M,player,c4_score,moves,
            c4_movePossible,c4_makeMove!,c4_gameOver);
    end

    # ids = count(!iszero,M[:,col])
    # M[6-ids,col] = player
    return c4_makeMove!(M,col,player)
end

function c4_score(M,player)
    win = c4_checkWinner(M)
    if win!=0
        return win==player ? 100 : -100
    elseif c4_gameOver(M)
        return 0
    else
        return c4_scoreHelper(M,player)
    end

    # returns 0 if the game isn't over, 1 if the player wins,
    # -1 if the player loses
end

function c4_checkWinner(M)
    # See who has won the game

    for player in (1,2)
        # horizontal lines
        rows = 1:6
        cols = 1:4
        for row in rows, col in cols
            if all(M[row,col:col+3].==player)
                return player
            end
        end

        # vertical lines
        rows = 1:3
        cols = 1:7
        for row in rows, col in cols
            if all(M[row:row+3,col].==player)
                return player
            end
        end

        # diagonals /
        rows = 1:3
        cols = 1:4
        for row in rows, col in cols
            temp = M[row:row+3,col:col+3]
            if all(temp[4:3:13].==player)
                return player
            end
            if all(temp[1:5:16].==player)
                return player
            end
        end
    end
    return 0
end

function c4_scoreHelper(M,player)
    # See who has won the game
    S = [0 0]

    for player in (1,2)
        scr = 0
        # horizontal lines
        rows = 1:6
        cols = 1:4
        for row in rows, col in cols
            if count(!iszero,M[row,col:col+3].==player)>3 && count(!iszero,M[row,col:col+3].==3-player)==0
                scr += 4
            elseif count(!iszero,M[row,col:col+3].==player)>2 && count(!iszero,M[row,col:col+3].==3-player)==0
                scr += 1
            end
        end

        # vertical lines
        rows = 1:3
        cols = 1:7
        for row in rows, col in cols
            if count(!iszero,M[row:row+3,col].==player)>3 && count(!iszero,M[row:row+3,col].==3-player)==0
                scr += 4
            elseif count(!iszero,M[row:row+3,col].==player)>2 && count(!iszero,M[row:row+3,col].==3-player)==0
                scr += 1
            end
        end

        # diagonals /
        rows = 1:3
        cols = 1:4
        for row in rows, col in cols
            temp = M[row:row+3,col:col+3]
            if count(!iszero,temp[4:3:13].==player)>3 && count(!iszero,temp[4:3:13].==3-player)==0
                scr += 4
            elseif count(!iszero,temp[4:3:13].==player)>2 && count(!iszero,temp[4:3:13].==3-player)==0
                scr += 1
            end
            if count(!iszero,temp[1:5:16].==player)>3 && count(!iszero,temp[1:5:16].==3-player)==0
                scr += 4
            elseif count(!iszero,temp[1:5:16].==player)>2 && count(!iszero,temp[1:5:16].==3-player)==0
                scr += 1
            end
        end
        S[player] = scr
    end
    return S[player]-S[3-player]
end

function c4_makeMove!(N,i,player)
    # does not check whether move is legal
    N[6-count(!iszero,N[:,i]),i] = player
    return N
end

function c4_movePossible(M,i)
    return count(!iszero,M[:,i])>=6 ? false : true
end

function c4_gameOver(M)
    # if all the squares are taken (filled with 1s/2s) the game ends
    return all(M.>0) || c4_checkWinner(M)!=0
end


text:
function Base.sort!(df::DataFrame; cols=Any[], alg=nothing,
                    lt=isless, by=identity, rev=false, order=Forward)
    if !(isa(by, Function) || eltype(by) <: Function)
        msg = "'by' must be a Function or a vector of Functions. Perhaps you wanted 'cols'."
        throw(ArgumentError(msg))
    end
    ord = ordering(df, cols, lt, by, rev, order)
    _alg = Sort.defalg(df, ord; alg=alg, cols=cols)
    sort!(df, _alg, ord)
end

function Base.sort!(df::DataFrame, a::Base.Sort.Algorithm, o::Base.Sort.Ordering)
    p = sortperm(df, a, o)
    pp = similar(p)
    c = columns(df)

    for (i,col) in enumerate(c)
        # Check if this column has been sorted already
        if any(j -> c[j]===col, 1:i-1)
            continue
        end

        copy!(pp,p)
        Base.permute!!(col, pp)
    end
    df
end


text:
using Oolong
using Test
using Base.Threads

@testset "Oolong.jl" begin
    include("core.jl")
    include("serve.jl")
end


text:
## Test lag function
@testset "Test lag" begin
    x = collect(1.0:5.0)
    y = collect(10.0:10.0:50.0)

    @test NarmaxLasso.lag(x, 2) == [1, 1, 1, 2, 3]
    @test NarmaxLasso.lag(y, 1) == [10, 10, 20, 30, 40]
    @test NarmaxLasso.lag(x, [1, 2]) == [1 1 2 3 4
                             1 1 1 2 3]'
    @test NarmaxLasso.lag([x y], [[1, 2], [2]]) == [1 1 2 3 4
                                        1 1 1 2 3
                                        10 10 10 20 30]'
end

## Test matrix2subviews
@testset "Test matrix2subviews" begin
    x, N, p = NarmaxLasso.matrix2subviews(ones(10, 10))
    @test x[1] == ones(10)
    @test length(x) == 10
end

## Test subviews2matrix
@testset "Test subviews2matrix" begin
    x, N, p = NarmaxLasso.matrix2subviews(ones(10, 10))
    X = NarmaxLasso.subviews2matrix(x, N, p)
    @test X == ones(10, 10)
end

## Test random input generation
@testset "Test random input generation" begin
    u = NarmaxLasso.generate_random_input(100, 5)
    @test length(u) == 100
    @test u[1] == u[2] == u[3] == u[4] == u[5]
    @test u[6] == u[7] == u[8] == u[9] == u[10]
    @test u[1] != u[6]
end