Merge the six PRs whose conflicts were registry-level

blueprint/src/python/bound_mu.py

@@ -220,7 +220,7 @@ def get_bounds(hypothesis_set):
 def compute_convex_hull(bounds, hypothesis_set):
     conv = ConvexHull(np.array([[b.data.sigma, b.data.mu] for b in bounds]))
     vertices = [bounds[v] for v in conv.vertices]
-    hypothesis_set.data["convex_hull"] = vertices
+    hypothesis_set.data['mu_convex_hull'] = vertices
     hypothesis_set.data_valid = True
 
 
@@ -236,7 +236,7 @@ def best_mu_bound(sigma, hypothesis_set):
 
     # Computed convex hull is stored within `hypothesis_list` the first time
     # to avoid repeatedly computing it for multiple values of sigma.
-    if not hypothesis_set.data_valid or "convex_hull" not in hypothesis_set.data:
+    if not hypothesis_set.data_valid or 'mu_convex_hull' not in hypothesis_set.data:
         # Generate set of bounds on mu implied by hypothesis set
         bounds = get_bounds(hypothesis_set)
         # The convex hull package requires us to provide at least 3 points, while some
@@ -246,7 +246,7 @@ def best_mu_bound(sigma, hypothesis_set):
         compute_convex_hull(pts, hypothesis_set)
 
     # The vertices are guaranteed to be in counterclockwise order for 2D hulls
-    verts = hypothesis_set.data["convex_hull"]
+    verts = hypothesis_set.data['mu_convex_hull']
     for i in range(len(verts)):
         b1 = verts[i]
         b2 = verts[(i + 1) % len(verts)]
@@ -280,33 +280,43 @@ def best_mu_bound_piecewise(
     if sigma0 > sigma1:
         raise ValueError("sigma0 must be <= sigma1")
 
-    # deal with edge pieces
+    if sigma_interval.is_empty():
+        return []
+
+    # Outside the critical strip the trivial bound is exact. Clip each piece
+    # to the requested interval, including requests wholly outside [0, 1].
+    if sigma1 <= 0:
+        return [Affine(-1, frac(1, 2), sigma_interval.deep_copy())]
+    if sigma0 >= 1:
+        return [Affine(0, 0, sigma_interval.deep_copy())]
     if sigma0 < 0:
-        return [("-x + 1/2", sigma0, 0)] + best_mu_bound_piecewise(
-            Interval(0, sigma1, include_lower=True, include_upper=sigma_interval.include_upper),
-            hypothesis_set
-        )
+        left = Interval(sigma0, 0, sigma_interval.include_lower, False)
+        return [Affine(-1, frac(1, 2), left)] + best_mu_bound_piecewise(
+            Interval(0, sigma1, True, sigma_interval.include_upper), hypothesis_set)
     if sigma1 > 1:
+        right = Interval(1, sigma1, False, sigma_interval.include_upper)
         return best_mu_bound_piecewise(
-            Interval(sigma0, 1, include_lower=sigma_interval.include_lower, include_upper=True),
-            hypothesis_set) + [("0", 1, sigma1)]
+            Interval(sigma0, 1, sigma_interval.include_lower, True), hypothesis_set
+        ) + [Affine(0, 0, right)]
 
     # Computed convex hull is stored within `hypothesis_set` the first time
     # to avoid repeatedly computing it for multiple values of sigma.
-    if not hypothesis_set.data_valid or "convex_hull" not in hypothesis_set.data:
+    if not hypothesis_set.data_valid or 'mu_convex_hull' not in hypothesis_set.data:
         bounds = get_bounds(hypothesis_set)
+        # As in best_mu_bound, ensure the trivial two-point case has a hull.
+        bounds.append(conjecture_bound_mu(frac(1, 2), 10, "Placeholder"))
         compute_convex_hull(bounds, hypothesis_set)
 
     # The vertices are guaranteed to be in counterclockwise order for 2D hulls
-    verts = hypothesis_set.data["convex_hull"]
+    verts = hypothesis_set.data['mu_convex_hull']
     mu_bounds = []
     for i in range(len(verts)):
         b1 = verts[i].data
         b2 = verts[(i + 1) % len(verts)].data
-        if b2.sigma < sigma0 or b1.sigma > sigma1:
+        if b1.sigma >= b2.sigma:
             continue
-        interval = Interval(max(sigma0, b1.sigma), min(sigma1, b2.sigma))
-        if b1.sigma < b2.sigma and not interval.is_empty():
+        interval = sigma_interval.intersect(Interval(b1.sigma, b2.sigma, True, True))
+        if not interval.is_empty():
             mu_bounds.append(
                 Affine(
                     (b2.mu - b1.mu) / (b2.sigma - b1.sigma),

blueprint/src/python/exponent_pair.py

@@ -69,6 +69,20 @@ def derived_exp_pair(k, l, proof, dependencies):
 
 ###############################################################################
 
+def _convex_hull_vertices(pairs):
+    """Return hull vertices, including point and segment hulls unsupported by Qhull."""
+    pairs = list({(p.data.k, p.data.l): p for p in pairs}.values())
+    if len(pairs) <= 2:
+        return pairs
+    ordered = sorted(pairs, key=lambda p: (p.data.k, p.data.l))
+    a, b = ordered[0].data, ordered[-1].data
+    if all((p.data.k - a.k) * (b.l - a.l) ==
+           (p.data.l - a.l) * (b.k - a.k) for p in ordered):
+        return [ordered[0], ordered[-1]]
+    conv = ConvexHull(np.array([[p.data.k, p.data.l] for p in pairs]))
+    return [pairs[v] for v in conv.vertices]
+
+
 def compute_exp_pairs(
         hypothesis_set: Hypothesis_Set,
         search_depth: int = 5,
@@ -112,9 +126,7 @@ def compute_exp_pairs(
         if prune:
             if len(pairs) < 3:
                 continue  # Don't prune if set is too small
-            pairs = [p for p in pairs.values()]
-            conv = ConvexHull(np.array([[p.data.k, p.data.l] for p in pairs]))
-            verts = [pairs[v] for v in conv.vertices]
+            verts = _convex_hull_vertices(pairs.values())
             pairs = {(p.data.k, p.data.l): p for p in verts}
 
     return [p for p in pairs.values()]
@@ -141,16 +153,11 @@ def compute_convex_hull(hypothesis_set: Hypothesis_Set) -> ConvexHull:
 
     # Computed convex hull is stored within `hypothesis_list` the first time
     # to avoid repeatedly computing it for multiple values of sigma.
-    if not hypothesis_set.data_valid or "convex_hull" not in hypothesis_set.data:
-        if len(pairs) < 3:
-            hypothesis_set.data["convex_hull"] = list(pairs)
-        else:
-            conv = ConvexHull(np.array([[p.data.k, p.data.l] for p in pairs]))
-            vertices = [pairs[v] for v in conv.vertices]
-            hypothesis_set.data["convex_hull"] = vertices
-            hypothesis_set.data_valid = True
+    if not hypothesis_set.data_valid or 'exponent_pair_convex_hull' not in hypothesis_set.data:
+        hypothesis_set.data['exponent_pair_convex_hull'] = _convex_hull_vertices(pairs)
+        hypothesis_set.data_valid = True
 
-    return hypothesis_set.data["convex_hull"]
+    return hypothesis_set.data['exponent_pair_convex_hull']
 
 def beta_bounds_to_exponent_pairs(
         hypothesis_set: Hypothesis_Set,

blueprint/src/python/functions.py

@@ -819,68 +819,69 @@ def _compute_optimal(
             track_dependencies=True
         ) -> list:
 
-        # Start with default bound with reference index -1
-        best_bound = [(RationalFunction([default]), domain, -1)]
-
-        for ref1 in range(len(func_and_domains)):
-            (func1, in1) = func_and_domains[ref1]
-
-            new_best_bound = []
-            for (func2, in2, ref2) in best_bound:
-                solns = sympy.solve(func1.num / func1.den - func2.num / func2.den)
-                crits = set(SympyHelper.to_frac(soln) for soln in solns if soln.is_real)
-                crits.update([in1.x0, in1.x1])
-                crits = set(c for c in crits if in2.contains(c))
-                crits.update([in2.x0, in2.x1])
-
-                crits = list(crits)
-                crits.sort()
-                for i in range(1, len(crits)):
-                    inter = Interval(crits[i - 1], crits[i])
-                    mid = inter.midpoint()
-                    if not in1.contains(mid) or comparator(func2.at(mid), func1.at(mid)):
-                        new_best_bound.append((func2, inter, ref2)) # f2 is the better bound
-                    else:
-                        new_best_bound.append((func1, inter, ref1))
-
-            # Simplify
-            best_bound = []
-            i = 0
-            while i < len(new_best_bound):
-                (fi, inti, refi) = new_best_bound[i]
-                left = inti.x0
-                right = inti.x1
-
-                j = i + 1
-                while j < len(new_best_bound):
-                    (fj, intj, refj) = new_best_bound[j]
-                    #print(type(refi), type(refj))
-                    if not (fi == fj and right == intj.x0 and (not track_dependencies or refi == refj)):
-                        break
-                    right = intj.x1
-                    j += 1
-
-                best_bound.append((fi, Interval(left, right), refi))
-                i = j
+        if domain.is_empty():
+            return []
+        x = RationalFunction.x
+
+        def real_roots(expression):
+            expression = sympy.sympify(expression)
+            if expression.is_constant():
+                return []
+            return sympy.real_roots(expression)
+
+        # A winner can change at a domain endpoint, an intersection, or a pole.
+        crits = {domain.x0, domain.x1}
+        for i, (func, interval) in enumerate(func_and_domains):
+            crits.update([interval.x0, interval.x1])
+            crits.update(real_roots(func.den))
+            for other, _ in func_and_domains[:i]:
+                crits.update(real_roots(func.num * other.den - other.num * func.den))
+        crits = sorted(c for c in crits if domain.x0 <= c <= domain.x1)
+        cells = []
+        for i, point in enumerate(crits):
+            if domain.contains(point):
+                cells.append(Interval(point, point, True, True))
+            if i + 1 < len(crits):
+                cells.append(Interval(point, crits[i + 1], False, False))
+
+        fallback = RationalFunction([default])
+        best_bound = []
+        for cell in cells:
+            mid = cell.midpoint()
+            winner, ref = fallback, -1
+            for i, (func, interval) in enumerate(func_and_domains):
+                if not interval.contains(mid) or sympy.sympify(func.den).subs(x, mid) == 0:
+                    continue
+                if ref == -1 or comparator(func.at(mid), winner.at(mid)):
+                    winner, ref = func, i
+            if best_bound:
+                prev, previous_domain, previous_ref = best_bound[-1]
+                same = (previous_ref == ref if track_dependencies else
+                        previous_ref == ref or (previous_ref != -1 and ref != -1 and prev == winner))
+                if (same and previous_domain.x1 == cell.x0 and
+                        (previous_domain.include_upper or cell.include_lower)):
+                    previous_domain.x1 = cell.x1
+                    previous_domain.include_upper = cell.include_upper
+                    continue
+            best_bound.append((winner, cell, ref))
 
         if track_dependencies:
             return best_bound
-        else:
-            return [(b[0], b[1]) for b in best_bound]
+        return [(f, interval) for f, interval, _ in best_bound]
 
     # Given a list of (RationalFunction, Interval) tuples, compute their (piecewise) minimum.
     # The result is returned as a list of (RationalFunction, Interval, Integer) tuples, where
     # the last element indicates the index of the RationalFunction piece in the original
     # list of functions.
-    # If a function is not defined on an interval, it is treated as -infinity.
+    # If a function is not defined on an interval, it is treated as +infinity.
     def min(func_and_domains: list, domain: Interval, track_dependencies=True) -> list:
         if not isinstance(func_and_domains, list):
             raise ValueError("Parameter piecewise_funcs must be of type list")
         if not isinstance(domain, Interval):
             raise ValueError("Parameter domain must be of type Interval")
 
         return RationalFunction._compute_optimal(
-            func_and_domains, domain, lambda x, y: x < y, 1000000, track_dependencies)
+            func_and_domains, domain, lambda x, y: x < y, sympy.oo, track_dependencies)
 
     # Given a list of (RationalFunction, Interval) tuples, compute their maximum. The
     # result is returned as a list of (RationalFunction, Interval, Integer) tuples, where
@@ -893,7 +894,7 @@ def max(func_and_domains: list, domain: Interval, track_dependencies=True) -> li
         if not isinstance(domain, Interval):
             raise ValueError("Parameter domain must be of type Interval")
         return RationalFunction._compute_optimal(
-            func_and_domains, domain, lambda x, y: x > y, -1000000, track_dependencies)
+            func_and_domains, domain, lambda x, y: x > y, -sympy.oo, track_dependencies)
 
     # ------------------------------------------------------------------
 

blueprint/src/python/hypotheses.py

@@ -169,10 +169,10 @@ def as_year(y):
 class Hypothesis_Set:
     def __init__(self, hypotheses=None):
         self.hypotheses = set()
-        if hypotheses is not None:
-            self.add_hypotheses(hypotheses)
         self.data = {}
         self.data_valid = False  # set to false whenever data needs to be recomputed
+        if hypotheses is not None:
+            self.add_hypotheses(hypotheses)
 
     def __repr__(self):
         return f'Set of {len(self.hypotheses)} hypotheses: [{",".join(h.name for h in self.hypotheses)}]'
@@ -207,6 +207,7 @@ def list_proofs(self):
     def add_hypothesis(self, hypothesis, invalidate_data=True):
         self.hypotheses.add(hypothesis)
         if invalidate_data:
+            self.data.clear()
             self.data_valid = False
 
     # Adds a set or list of hypotheses, an individual hypothesis, or a Hypothesis_Set
@@ -220,6 +221,7 @@ def add_hypotheses(self, new_hypotheses, invalidate_data=True):
         else:
             self.hypotheses.update(new_hypotheses)
             if invalidate_data:
+                self.data.clear()
                 self.data_valid = False
 
     # return all hypotheses of a given type, and (optionally) up to a given year.  Note: is now returning a Hypothesis_Set rather than a list

blueprint/src/python/reference.py

@@ -1,6 +1,16 @@
 # Code to read and interpret bibtex references
 
 
+def _year_matches(stored, wanted):
+    """Compare a stored year to a query, allowing int vs str."""
+    if stored == wanted:
+        return True
+    try:
+        return int(stored) == int(wanted)
+    except (TypeError, ValueError):
+        return False
+
+
 class Reference:
     def __init__(self, label, reftype, fields):
         self.label = label
@@ -60,9 +70,13 @@ def trivial():
     def max_year(refs):
         if len(refs) == 0:
             return -1
-        if any(r.year() == "Unknown date" for r in refs):
-            return "Unknown date"
-        return max(r.year() for r in refs)
+        years = []
+        for r in refs:
+            y = r.year()
+            if y == "Unknown date":
+                return "Unknown date"
+            years.append(int(y))
+        return max(years)
 
 
 class Reference_Manager:
@@ -93,8 +107,8 @@ def find_all(self, author="Any", year="Any"):
         refs = []
         for k in self.refs:
             ref = self.refs[k]
-            if (author == "Any" or author in ref.entries["author"]) and (
-                year == "Any" or ref.entries["year"] == year
+            if (author == "Any" or author in ref.entries.get("author", "")) and (
+                year == "Any" or _year_matches(ref.year(), year)
             ):
                 refs.append(ref)
         return refs

blueprint/src/python/region.py

@@ -272,11 +272,30 @@ def project(self, dims: set[int]) -> 'Region':
             return Region(Region_Type.UNION, ps)
 
         if self.region_type == Region_Type.INTERSECT:
-            ps = [c.project(dims) for c in self.child]
-            ps = [p for p in ps if p is not None]
-            if len(ps) == 0: return None
-            if len(ps) == 1: return ps[0]
-            return Region(Region_Type.INTERSECT, ps)
+            def convex_parts(region):
+                if region.region_type == Region_Type.POLYTOPE:
+                    return [] if region.child.is_empty(True) else [region.child]
+                if region.region_type in {Region_Type.UNION, Region_Type.DISJOINT_UNION}:
+                    return [p for child in region.child for p in convex_parts(child)]
+                if region.region_type == Region_Type.INTERSECT:
+                    if not region.child:
+                        raise ValueError("Cannot infer the dimension of an empty intersection")
+                    parts = convex_parts(region.child[0])
+                    for child in region.child[1:]:
+                        intersections = [Polytope.intersection([p, q])
+                                         for p in parts for q in convex_parts(child)]
+                        parts = [p for p in intersections if not p.is_empty(True)]
+                    return parts
+                raise NotImplementedError("Projection of complements is not supported")
+
+            projected = [p.project(dims) for p in convex_parts(self)]
+            children = [Region(Region_Type.POLYTOPE, p)
+                        for p in projected if p is not None]
+            if not children:
+                return None
+            if len(children) == 1:
+                return children[0]
+            return Region(Region_Type.UNION, children)
 
         # In particular, complements are not yet supported.
         raise NotImplementedError()

blueprint/src/python/tests/test_all.py

@@ -14,6 +14,10 @@
 
 import test_beta_reflection
 
+import test_hypothesis_cache
+
+import test_piecewise_mu
+
 import test_exppair
 
 print("All exponent pair test cases passed.")
@@ -47,4 +51,10 @@
 
 print("All Hypothesis test cases passed.")
 
+import test_reference
+
+print("All Reference test cases passed.")
+
+import test_rational_envelope_regression
+
 print("All test cases passed.")