Class: DataStructuresRMolinari::MaxPrioritySearchTree

Inherits:

Object

• Object
• DataStructuresRMolinari::MaxPrioritySearchTree

Includes:

BinaryTreeArithmetic, Shared

Defined in:

lib/data_structures_rmolinari/max_priority_search_tree.rb

Overview

A priority search tree (PST) stores a set, P, of two-dimensional points (x,y) in a way that allows efficient answers to certain questions about P.

The data structure was introduced in 1985 by Edward McCreight. Later, De, Maheshwari, Nandy, and Smid showed how to construct a PST in-place (using only O(1) extra memory), at the expense of some slightly more complicated code for the various supported operations. It is their approach that we have implemented.

The PST structure is an implicit, balanced binary tree with the following properties:

• The tree is a max-heap in the y coordinate. That is, the point at each node has a y-value no greater than its parent.

• For each node p, the x-values of all the nodes in the left subtree of p are less than the x-values of all the nodes in the right subtree of p. Note that this says nothing about the x-value at the node p itself. The tree is thus almost a binary search tree in the x coordinate.

Given a set of n points, we can answer the following questions quickly:

• smallest_x_in_ne: for x0 and y0, what is the leftmost point (x, y) in P satisfying x >= x0 and y >= y0?

• largest_x_in_nw: for x0 and y0, what is the rightmost point (x, y) in P satisfying x <= x0 and y >= y0?

• largest_y_in_ne: for x0 and y0, what is the highest point (x, y) in P satisfying x >= x0 and y >= y0?

• largest_y_in_nw: for x0 and y0, what is the highest point (x, y) in P satisfying x <= x0 and y >= y0?

• largest_y_in_3_sided: for x0, x1, and y0, what is the highest point (x, y) in P satisfying x >= x0, x <= x1 and y >= y0?

• enumerate_3_sided: for x0, x1, and y0, enumerate all points in P satisfying x >= x0, x <= x1 and y >= y0.

(Here, “leftmost/rightmost” means “minimal/maximal x”, and “highest” means “maximal y”.)

The first 5 operations take O(log n) time.

The final operation (enumerate) takes O(m + log n) time, where m is the number of points that are enumerated.

Each of these methods has a named parameter open: that makes the search region an open set. For example, if we call smallest_x_in_ne with open: true then we consider points satisifying x > x0 and y > y0. The default value for this parameter is always false. See below for limitations in this functionality.

If the MaxPST is constructed to be “dynamic” we also have an operation that deletes the top element.

• delete_top!: remove the top (max-y) element of the tree and return it.

It runs in O(log n) time, where n is the size of the PST when it was initially created.

In the current implementation no two points can share an x-value. This restriction can be relaxed with some more complicated code, but it hasn’t been written yet. See issue #9.

There is a related data structure called a Min-max priority search tree so we have called this a “Max priority search tree”, or MaxPST.

## Open regions: limitations

Calls involving open regions - using the open: argument - are implemented internally using closed regions in which the boundaries have been “nudged” by a tiny amount so as to exclude points on the boundary. Since there are only finitely many points in the PST there are no limit points, and the open region given by x > x0 and y > y0 contains the same PST members as a closed region x >= x0 + e and y >= y0 + e for small enough values of e.

But it is hard to determine e robustly. Indeed, assume for the moment that all x- and y-values are floating-point, i.e., IEEE754 double-precision. We can easily determine a value for e: the smallest difference between any two distinct x-values or any two distinct y-values. But the scaling of floating-point numbers makes this buggy. Consider the case in which we have consecutive x-values

0, 5e-324, 1, and 2.

(5e-324 is the smallest positive Float value in Ruby). Our value for e is thus 5e-324 and, because of the way floating point values are represented, 1 + e = 1.0. Any query on open region with x0 = 1 will be run on a closed region with x0 = 1 + e = 1.0, and we may get the wrong result.

The solution here is to replace (x0, y0) with (x0.next_float, y0.next_float) rather than (x0 + e, y0 + e). Float::next_float is an implementation of the IEEE754 operation nextafter which gives, for a floating point value z, the smallest floating point value larger than z. If our PST contains only points with finite, floating point coordinates, then this approach implements open search regions correctly. This is what the implementation currently does.

However, when coordinates are not all Floats there are cases when this approach will fail. Consider the case in which we have the following consecutive x-values:

0, 1e-324, 2e-324, 1.

(Here 1e-324 and 2e-324 are the Ruby values Rational(1, 10**324) and Rational(2, 10**324).) Then, given an argument x0 = 1e-324, x0.to_f == 0.0 and so x0.to_f.next_float == 5e-324 and the region we use internally for our query incorrectly excludes the point with x value 2e-324. This is a bug in the code.

There are also issues with numeric values (Integer or Rational) that are larger than the maximum floating point value, approximately 1.8e308. For such values z, z.to_f == Float::INFINITY, and we will incorrectly exclude any larger-but-finite coordinates from the search region.

Yet more issues arise when the coordinates of points in the PST aren’t numeric at all, but are some other sort of comparable objects, such as arrays.

So, we may say that queries on open regions will work as expected if either

• all coordinates of the points in the PST are finite Ruby Floats, or

• all coordinates of the points are finite Numeric values and for no such pair of x-values s, t (or pair of y-values) is it such that s.to_f.next_float > t.

Otherwise, use this functionality at your own risk, and not at all with coordinates that do not respond reasonably to to_f.

References:

• E.M. McCreight, _Priority search trees_, SIAM J. Comput., 14(2):257-276, 1985.

• 1. De, A. Maheshwari, S. C. Nandy, M. Smid, _An In-Place Priority Search Tree_, 23rd Canadian Conference on Computational Geometry, 2011

Constant Summary

Constants included from Shared

Shared::INFINITY

Instance Method Summary

• #delete_top! ⇒ Point

  Delete the top (max-y) element of the PST.

• #empty? ⇒ Boolean
• #enumerate_3_sided(x0, x1, y0, open: false) ⇒ Object

  Enumerate the points of P in the box bounded by x0, x1, and y0.

• #initialize(data, dynamic: false, verify: false) ⇒ MaxPrioritySearchTree constructor

  Construct a MaxPST from the collection of points in data.

• #largest_x_in_nw(x0, y0, open: false) ⇒ Object

  Return the rightmost (max-x) point in P to the northwest of (x0, y0).

• #largest_y_in_3_sided(x0, x1, y0, open: false) ⇒ Object

  Return the highest point of P in the box bounded by x0, x1, and y0.

• #largest_y_in_ne(x0, y0, open: false) ⇒ Object

  Return the highest point in P to the “northeast” of (x0, y0).

• #largest_y_in_nw(x0, y0, open: false) ⇒ Object

  Return the highest point in P to the “northwest” of (x0, y0).

• #smallest_x_in_ne(x0, y0, open: false) ⇒ Object

  Return the leftmost (min-x) point in P to the northeast of (x0, y0).

Methods included from Shared

#contains_duplicates?

Constructor Details

#initialize(data, dynamic: false, verify: false) ⇒ MaxPrioritySearchTree

Construct a MaxPST from the collection of points in data.

Parameters:

• data (Array) —

  the set P of points presented as an array. The tree is built in the array in-place without cloning.

  • Each element of the array must respond to #x and #y.

    • This is not checked explicitly but a missing method exception will be thrown when we try to call one of them.

  • The x values must be distinct. We raise a Shared::DataError if this isn’t the case.

    • This is a restriction that simplifies some of the algorithm code. It can be removed as the cost of some extra work. Issue #9.

• dynamic (Boolean) (defaults to: false) —

  when truthy the PST is dynamic. This means the root can be deleted, which is useful in certain algorithms than use a PST.

  • a dynamic PST needs more bookwork for some internal work and so slows things down a little.

• verify (Boolean) (defaults to: false) —

  when truthy, check that the properties of a PST are satisified after construction, raising an exception if not.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 179

def initialize(data, dynamic: false, verify: false)
  @data = data
  @size = @data.size
  @member_count = @size # these can diverge for dynamic PSTs
  @dynamic = dynamic

  construct_pst

  verify_properties if verify
end

Instance Method Details

#delete_top! ⇒ Point

Delete the top (max-y) element of the PST. This is possible only for dynamic PSTs

It runs in guaranteed O(log n) time, where n is the size of the PST when it was intially constructed. As elements are deleted the internal tree structure is no longer guaranteed to be balanced and so we cannot guarantee operation in O(log n’) time, where n’ is the current size. In practice, “random” deletion is likely to leave the tree almost balanced.

Returns:

• (Point) —

  the top element that was deleted

Raises:

• (LogicError)

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 1222

def delete_top!
  raise LogicError, 'delete_top! not supported for PSTs that are not dynamic' unless dynamic?
  raise DataError, 'delete_top! not possible for empty PSTs' unless @member_count.positive?

  i = root
  while !leaf?(i)
    if (child = one_child?(i))
      next_node = child
    else
      next_node = left(i)

      if better_y?(right(i), next_node)
        next_node = right(i)
      end
    end
    swap(i, next_node)
    i = next_node
  end
  @member_count -= 1
  @data[i]
end

#empty? ⇒ Boolean

Returns:

• (Boolean)

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 190

def empty?
  @member_count.zero?
end

#enumerate_3_sided(x0, x1, y0, open: false) ⇒ Object

Enumerate the points of P in the box bounded by x0, x1, and y0.

Let Q be the “three-sided” box bounded by x0, x1, and y0:

• [x0, x1] X [y0, infty) if open is false and

• (x0, x1) X (y0, infty) if open is true.

Note that Q is empty if x1 < x0 or if open is true and x1 <= x0.

Let P be the set of points in the MaxPST.

We find and enumerate all the points in Q intersect P.

If the calling code provides a block then we yield each point to it. Otherwise we return a set containing all the points in the intersection.

This method runs in O(m + log n) time and O(1) extra space, where m is the number of points found.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 753

def enumerate_3_sided(x0, x1, y0, open: false)
  if open
    x0 = slightly_bigger(x0)
    x1 = slightly_smaller(x1)
    y0 = slightly_bigger(y0)
  end

  # From the paper
  #
  #     Given three real numbers x0, x1, and y0 define the three sided range Q = [x0, x1] X [y0, infty). Algorithm
  #     Enumerage3Sided(x0, x1,y0) returns all elements of Q \intersect P. The algorithm uses the same approach as algorithm
  #     Highest3Sided. Besides the two bits L and R it uses two additional bits L' and R'. Each of these four bits ... corresponds
  #     to a subtree of T rooted at the points p, p', q, and q', respectively; if the bit is equal to one, then the subtree may
  #     contain points that are in the query range Q.
  #
  #     The following variant will be maintained:
  #
  #     - If L = 1 then x(p) < x0.
  #     - If L' = 1 then x0 <= x(p') <= x1.
  #     - If R = 1 then x(q) > x1.
  #     - If R' = 1 then x0 <= x(q') <= x1.
  #     - If L' = 1 and R' = 1 then x(p') <= x(q').
  #     - All points in Q \intersect P [other than those in the subtrees of the currently active search nodes] have been reported.
  #
  #
  # My high-level understanding of the algorithm
  # --------------------------------------------
  #
  # We need to find all elements of Q \intersect P, so it isn't enough, as it was in largest_y_in_3_sided simply to keep track of p and
  # q. We need to track four nodes, p, p', q', and q which are (with a little handwaving) respectively
  #
  # - the rightmost node to the left of Q' = [x0, x1] X [-infinity, infinity],
  # - the leftmost node inside Q',
  # - the rightmost node inside Q', and
  # - the leftmost node to the right of Q'.
  #
  # Tracking these is enough. Subtrees of things to the left of p can't have anything in Q by the x-value properties of the PST,
  # and likewise with things to the right of q.
  #
  # And we don't need to track any more nodes inside Q'. If we had r with p' <~ r <~ q' (where s <~ t represents "t is to the
  # right of s"), then all of the subtree rooted at r lies inside Q', and we can visit all of its elements of Q \intersect P via
  # the routine Explore(), which is what we do whenever we need to. The node r is thus exhausted, and we can forget about it.
  #
  # So the algorithm is actually quite simple. There is a large amount of code here because of the many cases that need to be
  # handled at each update.
  #
  # If a block is given, yield each found point to it. Otherwise return all the found points in an enumerable (currently Set).
  x_range = x0..x1
  # Instead of using primes we use "_in"
  left = left_in = right_in = right = false
  p = p_in = q_in = q = nil

  result = Set.new
  if empty?
    return if block_given?

    return result
  end

  report = lambda do |node|
    if block_given?
      yield @data[node]
    else
      result << @data[node]
    end
  end

  # "reports all points in T_t whose y-coordinates are at least y0"
  #
  # We follow the logic from the min-max paper, leaving out the need to worry about the parity of the leval and the min- or max-
  # switching.
  explore = lambda do |t|
    current = t
    state = 0
    while current != t || state != 2
      case state
      when 0
        # State 0: we have arrived at this node for the first time
        # look at current and perhaps descend to left child
        # The paper describes this algorithm as in-order, but isn't this pre-order?
        if @data[current].y >= y0
          report.call(current)
        end
        if !leaf?(current) && in_tree?(left(current)) && @data[left(current)].y >= y0
          current = left(current)
        else
          state = 1
        end
      when 1
        # State 1: we've already handled this node and its left subtree. Should we descend to the right subtree?
        if in_tree?(right(current)) && @data[right(current)].y >= y0
          current = right(current)
          state = 0
        else
          state = 2
        end
      when 2
        # State 2: we're done with this node and its subtrees. Go back up a level, having set state correctly for the logic at the
        # parent node.
        if left_child?(current)
          state = 1
        end
        current = parent(current)
      else
        raise InternalLogicError, "Explore(t) state is somehow #{state} rather than 0, 1, or 2."
      end
    end
  end

  # Helpers for the helpers
  #
  # Invariant: if q_in is active then p_in is active. In other words, if only one "inside" node is active then it is p_in.

  # Mark p_in as inactive. Then, if q_in is active, it becomes p_in.
  deactivate_p_in = lambda do
    left_in = false
    return unless right_in

    p_in = q_in
    left_in = true
    right_in = false
  end

  # Add a new leftmost "in" point. This becomes p_in. We handle existing "inside" points appropriately
  add_leftmost_inner_node = lambda do |node|
    if left_in && right_in
      # the old p_in is squeezed between node and q_in
      explore.call(p_in)
    elsif left_in
      q_in = p_in
      right_in = true
    else
      left_in = true
    end
    p_in = node
  end

  add_rightmost_inner_node = lambda do |node|
    if left_in && right_in
      # the old q_in is squeezed between p_in and node
      explore.call(q_in)
      q_in = node
    elsif left_in
      q_in = node
      right_in = true
    else
      p_in = node
      left_in = true
    end
  end

  ########################################
  # The four key helpers described in the paper

  # Handle the next step of the subtree at p
  enumerate_left = lambda do
    if leaf?(p)
      left = false
      return
    end

    if (only_child = one_child?(p))
      child_val = @data[only_child]
      if x_range.cover? child_val.x
        add_leftmost_inner_node.call(only_child)
        left = false
      elsif child_val.x < x0
        p = only_child
      else
        q = only_child
        right = true
        left = false
      end
      return
    end

    # p has two children
    if @data[left(p)].x < x0
      if @data[right(p)].x < x0
        p = right(p)
      elsif @data[right(p)].x <= x1
        add_leftmost_inner_node.call(right(p))
        p = left(p)
      else
        q = right(p)
        p = left(p)
        right = true
      end
    elsif @data[left(p)].x <= x1
      if @data[right(p)].x > x1
        q = right(p)
        p_in = left(p)
        left = false
        left_in = right = true
      else
        # p_l and p_r both lie inside [x0, x1]
        add_leftmost_inner_node.call(right(p))
        add_leftmost_inner_node.call(left(p))
        left = false
      end
    else
      q = left(p)
      left = false
      right = true
    end
  end

  # Given: p' satisfied x0 <= x(p') <= x1. (Our p_in is the paper's p')
  enumerate_left_in = lambda do
    if @data[p_in].y >= y0
      report.call(p_in)
    end

    if leaf?(p_in) # nothing more to do
      deactivate_p_in.call
      return
    end

    if (only_child = one_child?(p_in))
      child_val = @data[only_child]
      if x_range.cover? child_val.x
        p_in = only_child
      elsif child_val.x < x0
        # We aren't in the [x0, x1] zone any more and have moved out to the left
        p = only_child
        deactivate_p_in.call
        left = true
      else
        # similar, but we've moved out to the right. Note that left(p_in) is the leftmost node to the right of Q.
        raise 'q_in should not be active (by the val of left(p_in))' if right_in

        q = only_child
        deactivate_p_in.call
        right = true
      end
    else
      # p' has two children
      left_val = @data[left(p_in)]
      right_val = @data[right(p_in)]
      if left_val.x < x0
        if right_val.x < x0
          p = right(p_in)
          left = true
          deactivate_p_in.call
        elsif right_val.x <= x1
          p = left(p_in)
          p_in = right(p_in)
          left = true
        else
          raise InternalLogicError, 'q_in cannot be active, by the value in the right child of p_in!' if right_in

          p = left(p_in)
          q = right(p_in)
          deactivate_p_in.call
          left = true
          right = true
        end
      elsif left_val.x <= x1
        if right_val.x > x1
          raise InternalLogicError, 'q_in cannot be active, by the value in the right child of p_in!' if right_in

          q = right(p_in)
          p_in = left(p_in)
          right = true
        elsif right_in
          explore.call(right(p_in))
          p_in = left(p_in)
        else
          q_in = right(p_in)
          p_in = left(p_in)
          right_in = true
        end
      else
        raise InternalLogicError, 'q_in cannot be active, by the value in the right child of p_in!' if right_in

        q = left(p_in)
        deactivate_p_in.call
        right = true
      end
    end
  end

  # This is "just like" enumerate left, but handles q instead of p.
  #
  # The paper doesn't given an implementation, but it should be pretty symmetric. Can we share any logic with enumerate_left?
  #
  # Q: why is my implementation more complicated than enumerate_left? I must be missing something.
  enumerate_right = lambda do
    if leaf?(q)
      right = false
      return
    end

    if (only_child = one_child?(q))
      child_val = @data[only_child]
      if x_range.cover? child_val.x
        add_rightmost_inner_node.call(only_child)
        right = false
      elsif child_val.x < x0
        p = only_child
        left = true
        right = false
      else
        q = only_child
      end
      return
    end

    # q has two children. Cases!
    if @data[left(q)].x < x0
      raise InternalLogicError, 'p_in should not be active, based on the value at left(q)' if left_in
      raise InternalLogicError, 'q_in should not be active, based on the value at left(q)' if right_in

      left = true
      if @data[right(q)].x < x0
        p = right(q)
        right = false
      elsif @data[right(q)].x <= x1
        p_in = right(q)
        p = left(q)
        left_in = true
        right = false
      else
        p = left(q)
        q = right(q)
      end
    elsif @data[left(q)].x <= x1
      add_rightmost_inner_node.call(left(q))
      if @data[right(q)].x > x1
        q = right(q)
      else
        add_rightmost_inner_node.call(right(q))
        right = false
      end
    else
      # x(q_l) > x1
      q = left(q)
    end
  end

  # Given: q' is active and satisfied x0 <= x(q') <= x1
  enumerate_right_in = lambda do
    raise InternalLogicError, 'right_in should be true if we call enumerate_right_in' unless right_in

    if @data[q_in].y >= y0
      report.call(q_in)
    end

    if leaf?(q_in)
      right_in = false
      return
    end

    if (only_child = one_child?(q_in))
      child_val = @data[only_child]
      if x_range.cover? child_val.x
        q_in = only_child
      elsif child_val.x < x0
        # We have moved out to the left
        p = only_child
        right_in = false
        left = true
      else
        # We have moved out to the right
        q = only_child
        right_in = false
        right = true
      end
      return
    end

    # q' has two children
    left_val = @data[left(q_in)]
    right_val = @data[right(q_in)]
    if left_val.x < x0
      raise InternalLogicError, 'p_in cannot be active, by the value in the left child of q_in' if left_in

      if right_val.x < x0
        p = right(q_in)
      elsif right_val.x <= x1
        p = left(q_in)
        p_in = right(q_in) # should this be q_in = right(q_in) ??
        left_in = true
      else
        p = left(q_in)
        q = right(q_in)
        right = true
      end
      right_in = false
      left = true
    elsif left_val.x <= x1
      if right_val.x > x1
        q = right(q_in)
        right = true
        if left_in
          q_in = left(q_in)
        else
          p_in = left(q_in)
          left_in = true
          right_in = false
        end
      else
        if left_in
          explore.call(left(q_in))
        else
          p_in = left(q_in)
          left_in = true
        end
        q_in = right(q_in)
      end
    else
      q = left(q_in)
      right_in = false
      right = true
    end
  end

  val = ->(sym) { { left: p, left_in: p_in, right_in: q_in, right: q }[sym] }

  root_val = @data[root]
  if root_val.y < y0
    # no hope, no op
  elsif root_val.x < x0
    p = root
    left = true
  elsif root_val.x <= x1 # Possible bug in paper, which tests "< x1"
    p_in = root
    left_in = true
  else
    q = root
    right = 1
  end

  while left || left_in || right_in || right
    raise InternalLogicError, 'It should not be that q_in is active but p_in is not' if right_in && !left_in

    set_i = []
    set_i << :left if left
    set_i << :left_in if left_in
    set_i << :right_in if right_in
    set_i << :right if right
    z = set_i.min_by { |sym| level(val.call(sym)) }
    case z
    when :left
      enumerate_left.call
    when :left_in
      enumerate_left_in.call
    when :right_in
      enumerate_right_in.call
    when :right
      enumerate_right.call
    else
      raise InternalLogicError, "bad symbol #{z}"
    end
  end
  return result unless block_given?
end

#largest_x_in_nw(x0, y0, open: false) ⇒ Object

Return the rightmost (max-x) point in P to the northwest of (x0, y0).

Let Q = be the northwest quadrant defined by the point (x0, y0):

• (infty, x0] X [y0, infty) if open is false and

• (infty, x0) X (y0, infty) if open is true.

Let P be the points in this data structure.

Define p* as

• (-infty, infty) if Q intersect P is empty and

• the leftmost (min-x) point in Q intersect P otherwise.

This method returns p* in O(log n) time and O(1) extra space.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 371

def largest_x_in_nw(x0, y0, open: false)
  if open
    extremal_in_x_dimension(slightly_smaller(x0), slightly_bigger(y0), :nw)
  else
    extremal_in_x_dimension(x0, y0, :nw)
  end
end

#largest_y_in_3_sided(x0, x1, y0, open: false) ⇒ Object

Return the highest point of P in the box bounded by x0, x1, and y0.

Let Q be the “three-sided” box bounded by x0, x1, and y0:

• [x0, x1] X [y0, infty) if open is false and

• (x0, x1) X (y0, infty) if open is true.

Note that Q is empty if x1 < x0 or if open is true and x1 <= x0.

Let P be the set of points in the MaxPST.

Define p* as

• (infty, -infty) if Q intersect P is empty and

• the highest (max-y) point in Q intersect P otherwise, breaking ties by preferring smaller x values.

This method returns p* in O(log n) time and O(1) extra space.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 538

def largest_y_in_3_sided(x0, x1, y0, open: false)
  if open
    x0 = slightly_bigger(x0)
    x1 = slightly_smaller(x1)
    y0 = slightly_bigger(y0)
  end
  # From the paper:
  #
  #    The three real numbers x0, x1, and y0 define the three-sided range Q = [x0,x1] X [y0,∞). If Q \intersect P̸ is not \empty,
  #    define p* to be the highest point of P in Q. If Q \intersect P = \empty, define p∗ to be the point (infty, -infty).
  #    Algorithm Highest3Sided(x0,x1,y0) returns the point p∗.
  #
  #    The algorithm uses two bits L and R, and three variables best, p, and q. As before, best stores the highest point in Q
  #    found so far. The bit L indicates whether or not p∗ may be in the subtree of p; if L=1, then p is to the left of
  #    Q. Similarly, the bit R indicates whether or not p∗ may be in the subtree of q; if R=1, then q is to the right of Q.
  #
  # Although there are a lot of lines and cases the overall idea is simple. We maintain in p the rightmost node at its level that
  # is to the left of the area Q. Likewise, q is the leftmost node that is the right of Q. The logic just updates this data at
  # each step. The helper check_left updates p and check_right updates q.
  #
  # A couple of simple observations that show why maintaining just these two points is enough.
  #
  # - We know that x(p) < x0. This tells us nothing about the x values in the subtrees of p (which is why we need to check various
  #   cases), but it does tell us that everything to the left of p has values of x that are too small to bother with.
  # - We don't need to maintain any state inside the region Q because the max-heap property means that if we ever find a node r in
  #   Q we check it for best and then ignore its subtree (which cannot beat r on y-value).
  #
  # Sometimes we don't have a relevant node to the left or right of Q. The booleans L and R (which we call left and right) track
  # whether p and q are defined at the moment.
  best = Point.new(INFINITY, -INFINITY)
  return best if empty?

  p = q = left = right = nil

  x_range = (x0..x1)

  in_q = lambda do |pair|
    x_range.cover?(pair.x) && pair.y >= y0
  end

  # From the paper:
  #
  #   takes as input a point t and does the following: if t \in Q and x(t) < x(best) then it assignes best = t
  #
  # Note that the paper identifies a node in the tree with its value. We need to grab the correct node.
  update_highest = lambda do |node|
    t = @data[node]
    if in_q.call(t) && (t.y > best.y || (t.y == best.y && t.x < best.x))
      best = t
    end
  end

  # "Input: a node p such that x(p) < x0""
  #
  # Step-by-step it is pretty straightforward. As the paper says
  #
  #   [E]ither p moves one level down in the tree T or the bit L is set to 0. In addition, the point q either stays the same or it
  #   become a child of (the original) p.
  check_left = lambda do
    if leaf?(p)
      left = false
    elsif (only_child = one_child?(p))
      if x_range.cover? @data[only_child].x
        update_highest.call(only_child)
        left = false # can't do y-better in the subtree
      elsif @data[only_child].x < x0
        p = only_child
      else
        q = only_child
        right = true
        left = false
      end
    else
      # p has two children
      if @data[left(p)].x < x0
        if @data[right(p)].x < x0
          p = right(p)
        elsif @data[right(p)].x <= x1
          update_highest.call(right(p))
          p = left(p)
        else
          # x(p_r) > x1, so q needs to take it
          q = right(p)
          p = left(p)
          right = true
        end
      elsif @data[left(p)].x <= x1
        update_highest.call(left(p))
        left = false # we won't do better in T(p_l)
        if @data[right(p)].x > x1
          q = right(p)
          right = true
        else
          update_highest.call(right(p))
        end
      else
        q = left(p)
        left = false
        right = true
      end
    end
  end

  # Do "on the right" with q what check_left does on the left with p
  #
  # We know that x(q) > x1
  #
  # TODO: can we share logic between check_left and check_right? At first glance they are too different to parameterize but maybe
  # the bones can be shared.
  #
  # We either push q further down the tree or make right = false. We might also make p a child of (original) q. We never change
  # left from true to false
  check_right = lambda do
    if leaf?(q)
      right = false
    elsif (only_child = one_child?(q))
      if x_range.cover? @data[only_child].x
        update_highest.call(only_child)
        right = false # can't do y-better in the subtree
      elsif @data[only_child].x < x0
        p = only_child
        left = true
        right = false
      else
        q = only_child
      end
    else
      # q has two children
      if @data[left(q)].x < x0
        left = true
        if @data[right(q)].x < x0
          p = right(q)
          right = false
        elsif @data[right(q)].x <= x1
          update_highest.call(right(q))
          p = left(q)
          right = false
        else
          # x(q_r) > x1
          p = left(q)
          q = right(q)
          # left = true
        end
      elsif @data[left(q)].x <= x1
        update_highest.call(left(q))
        if @data[right(q)].x > x1
          q = right(q)
        else
          update_highest.call(right(q))
          right = false
        end
      else
        q = left(q)
      end
    end
  end

  return best if empty?

  root_val = @data[root]

  # If the root value is in the region Q, the max-heap property on y means we can't do better
  if x_range.cover? root_val.x
    # If y(root) is large enough then the root is the winner because of the max heap property in y. And if it isn't large enough
    # then no other point in the tree can be high enough either
    left = right = false
    best = root_val if root_val.y >= y0
  end

  if root_val.x < x0
    p = root
    left = true
    right = false
  else
    q = root
    left = false
    right = true
  end

  val = ->(sym) { sym == :left ? p : q }

  while left || right
    set_i = []
    set_i << :left if left
    set_i << :right if right
    z = set_i.min_by { |s| level(val.call(s)) }
    if z == :left
      check_left.call
    else
      check_right.call
    end
  end

  best
end

#largest_y_in_ne(x0, y0, open: false) ⇒ Object

Return the highest point in P to the “northeast” of (x0, y0).

Let Q = be the northeast quadrant defined by the point (x0, y0):

• [x0, infty) X [y0, infty) if open is false and

• (x0, infty) X (y0, infty) if open is true.

Let P be the points in this data structure.

Define p* as

• (infty, -infty) if Q intersect P is empty and

• the highest (max-y) point in Q intersect P otherwise, breaking ties by preferring smaller values of x

This method returns p* in O(log n) time and O(1) extra space.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 211

def largest_y_in_ne(x0, y0, open: false)
  if open
    largest_y_in_quadrant(slightly_bigger(x0), slightly_bigger(y0), :ne)
  else
    largest_y_in_quadrant(x0, y0, :ne)
  end
end

#largest_y_in_nw(x0, y0, open: false) ⇒ Object

Return the highest point in P to the “northwest” of (x0, y0).

Let Q = be the northwest quadrant defined by the point (x0, y0):

• (infty, x0] X [y0, infty) if open is false and

• (infity, x0) X (y0, infty) if open is true.

Let P be the points in this data structure.

Define p* as

• (-infty, -infty) if Q intersect P is empty and

• the highest (max-y) point in Q intersect P otherwise, breaking ties by preferring smaller values of x

This method returns p* in O(log n) time and O(1) extra space.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 233

def largest_y_in_nw(x0, y0, open: false)
  if open
    largest_y_in_quadrant(slightly_smaller(x0), slightly_bigger(y0), :nw)
  else
    largest_y_in_quadrant(x0, y0, :nw)
  end
end

#smallest_x_in_ne(x0, y0, open: false) ⇒ Object

Return the leftmost (min-x) point in P to the northeast of (x0, y0).

Let Q = be the northeast quadrant defined by the point (x0, y0):

• [x0, infty) X [y0, infty) if open is false and

• (x0, infty) X (y0, infty) if open is true.

Let P be the points in this data structure.

Define p* as

• (infty, infty) if Q intersect P is empty and

• the leftmost (min-x) point in Q intersect P otherwise.

This method returns p* in O(log n) time and O(1) extra space.

# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 349

def smallest_x_in_ne(x0, y0, open: false)
  if open
    extremal_in_x_dimension(slightly_bigger(x0), slightly_bigger(y0), :ne)
  else
    extremal_in_x_dimension(x0, y0, :ne)
  end
end