Class: DataStructuresRMolinari::MaxPrioritySearchTree

Inherits:
Object
  • Object
show all
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.

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 the Min-max priority search tree so we have called this a “Max priority search tree”, or MaxPST.

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 collapse

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.



66
67
68
69
70
71
72
73
74
75
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 66

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:



1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 1045

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)


77
78
79
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 77

def empty?
  @member_count.zero?
end

#enumerate_3_sided(x0, x1, y0) ⇒ Object

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

Let Q = [x0, x1] X [y0, infty) be the “three-sided” box bounded by x0, x1, and y0, and let P be the set of points in the MaxPST. (Note that Q is empty if x1 < x0.) We find an 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.



582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 582

def enumerate_3_sided(x0, x1, y0)
  # 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) ⇒ Object

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

Let Q = (-infty, x0] X [y0, infty) be the northwest quadrant defined by the point (x0, y0) and 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.



223
224
225
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 223

def largest_x_in_nw(x0, y0)
  extremal_in_x_dimension(x0, y0, :nw)
end

#largest_y_in_3_sided(x0, x1, y0) ⇒ Object

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

Let Q = [x0, x1] X [y0, infty) be the “three-sided” box bounded by x0, x1, and y0, and let P be the set of points in the MaxPST. (Note that Q is empty if x1 < x0.) 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.



379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 379

def largest_y_in_3_sided(x0, x1, y0)
  # 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) ⇒ Object

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

Let Q = [x0, infty) X [y0, infty) be the northeast quadrant defined by the point (x0, y0) and 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.



93
94
95
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 93

def largest_y_in_ne(x0, y0)
  largest_y_in_quadrant(x0, y0, :ne)
end

#largest_y_in_nw(x0, y0) ⇒ Object

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

Let Q = (-infty, x0] X [y0, infty) be the northwest quadrant defined by the point (x0, y0) and 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.



106
107
108
# File 'lib/data_structures_rmolinari/max_priority_search_tree.rb', line 106

def largest_y_in_nw(x0, y0)
  largest_y_in_quadrant(x0, y0, :nw)
end

#smallest_x_in_ne(x0, y0) ⇒ Object

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

Let Q = [x0, infty) X [y0, infty) be the northeast quadrant defined by the point (x0, y0) and 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.



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

def smallest_x_in_ne(x0, y0)
  extremal_in_x_dimension(x0, y0, :ne)
end