Cell Descriptions of Minimal Functions - Inital version

cutgeneratingfunctionology/igp/__init__.py

@@ -76,9 +76,9 @@ def igp_load(fpath):
 igp_load(igp_dir + "plot_options.sage")
 igp_load(igp_dir + "faster_subadditivity_test.sage")
 igp_load(igp_dir + "faster_subadditivity_test_discontinuous.sage")
+igp_load(igp_dir + "minimal_function_cell_description.sage")
 
 from . import extreme_functions, procedures
-
 try:
     igp_load(igp_dir + "config.sage")
 except IOError:

cutgeneratingfunctionology/igp/functions.sage

@@ -1517,9 +1517,15 @@ def piecewise_function_from_breakpoints_slopes_and_values(bkpt, slopes, values,
     if field is None:
         field = default_field
     # global symb_values
-    symb_values = bkpt + slopes + values
-    field_values = nice_field_values(symb_values, field)
-    bkpt, slopes, values = field_values[0:len(bkpt)], field_values[len(bkpt):len(bkpt)+len(slopes)], field_values[-len(values):]
+    if slopes is None: # make order assumptions in these functions, remove != when possible.
+        symb_values = bkpt + values
+        field_values = nice_field_values(symb_values, field)
+        bkpt, values = field_values[0:len(bkpt)], field_values[-len(values):]
+        slopes = [(values[i+1]-values[i])/(bkpt[i+1]-bkpt[i]) if bkpt[i+1] > bkpt[i] else 0 for i in range(len(bkpt)-1)] 
+    else:
+        symb_values = bkpt + slopes + values
+        field_values = nice_field_values(symb_values, field)
+        bkpt, slopes, values = field_values[0:len(bkpt)], field_values[len(bkpt):len(bkpt)+len(slopes)], field_values[-len(values):]
     intercepts = [ values[i] - slopes[i]*bkpt[i] for i in range(len(slopes)) ]
     # Make numbers nice
     ## slopes = [ canonicalize_number(slope) for slope in slopes ]
@@ -1550,8 +1556,7 @@ def piecewise_function_from_breakpoints_and_values(bkpt, values, field=None, mer
     """
     if len(bkpt)!=len(values):
         raise ValueError("Need to have the same number of breakpoints and values.")
-    slopes = [ (values[i+1]-values[i])/(bkpt[i+1]-bkpt[i]) if bkpt[i+1] != bkpt[i] else 0 for i in range(len(bkpt)-1) ]
-    return piecewise_function_from_breakpoints_slopes_and_values(bkpt, slopes, values, field, merge=merge)
+    return piecewise_function_from_breakpoints_slopes_and_values(bkpt, None, values, field, merge=merge)
 
 def piecewise_function_from_breakpoints_and_slopes(bkpt, slopes, field=None, merge=True):
     r"""

cutgeneratingfunctionology/spam/EvaluationExceptions.py

@@ -0,0 +1,2 @@
+class FactorUndetermined(Exception): # FactorUndetermined is raised when an expression can not be evaluated with a test point.
+    pass

cutgeneratingfunctionology/spam/basic_semialgebraic.py

@@ -35,6 +35,9 @@
 from sage.structure.sage_object import SageObject
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 
+class EmptyBSA(Exception):
+    pass
+
 def _bsa_class(bsa_class):
     r"""
     Translate a class nickname to a class.
@@ -1070,8 +1073,8 @@ def __contains__(self, point):
     # override the abstract methods
     def find_point(self):
         r"""
-        Find a point in ``self``.
-
+            Find a point in ``self``.
+            
         EXAMPLES::
 
             sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import *
@@ -1081,6 +1084,12 @@ def find_point(self):
             sage: P.add_linear_constraint([2,3],-6,operator.lt)
             sage: P.find_point()
             (11/10, 11/15)
+            sage: P = BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(1)
+            sage: P.add_linear_constraint([1],0,operator.ge)
+            sage: P.add_linear_constraint([1],-1,operator.lt)
+            sage: P.find_point()
+            EmptyBSA
+
         """
         def to_point(g):
             den = g.divisor()
@@ -1091,6 +1100,8 @@ def to_vector(g):
                    if g.is_point() or g.is_closure_point() ]
         rays = [ to_vector(g) for g in self._polyhedron.generators()
                  if g.is_ray() ]
+        if self.is_empty():
+            raise EmptyBSA("The underlying bsa is empty set.")
         if points:
             p = sum(points) / len(points)
             if rays: