diff --git a/cutgeneratingfunctionology/igp/parametric.sage b/cutgeneratingfunctionology/igp/parametric.sage
index 908b10e9a..b1a6bf2a8 100644
--- a/cutgeneratingfunctionology/igp/parametric.sage
+++ b/cutgeneratingfunctionology/igp/parametric.sage
@@ -21,6 +21,7 @@ from cutgeneratingfunctionology.spam.semialgebraic_mathematica import BasicSemia
 from cutgeneratingfunctionology.spam.basic_semialgebraic_groebner_basis import BasicSemialgebraicSet_groebner_basis
 from cutgeneratingfunctionology.spam.polyhedral_complex import PolyhedralComplex
 from .parametric_family import Classcall, ParametricFamily_base, ParametricFamily
+from cutgeneratingfunctionology.spam.exceptions import *
 
 debug_new_factors = False
 debug_cell_exceptions = False
@@ -36,31 +37,19 @@ def bigcellify_igp():
     import cutgeneratingfunctionology.spam.big_cells as big_cells
     big_cells.bigcellify_module(igp)
 
+
 ###############################
 # Parametric Real Number Field
 ###############################
 
 from cutgeneratingfunctionology.spam.parametric_real_field_element import ParametricRealFieldElement, is_parametric_element
-
 from sage.rings.ring import Field
 import sage.rings.number_field.number_field_base as number_field_base
 from sage.structure.coerce_maps import CallableConvertMap
 
 
-class ParametricRealFieldFrozenError(ValueError):
-    pass
-
-class ParametricRealFieldInconsistencyError(ValueError):
-    pass
-
-class ParametricRealFieldRefinementError(ValueError):
-    pass
-
 from contextlib import contextmanager
 
-class FactorUndetermined(Exception):
-    pass
-
 allow_refinement_default = True
 big_cells_default = 'if_not_allow_refinement'
 mutable_values_default = False
@@ -158,26 +147,23 @@ class ParametricRealField(Field):
         sage: f[0]*f[1] <= 4
         True
 
-    Test-point free mode (limited functionality and MUCH slower because of many more polynoial
-    evaluations via libsingular)::
+    Test-point free descriptions can be written which every comparison is assumed to be true.
+    MUCH slower because of many more polynoial evaluations via libsingular::
 
         sage: K.<a, b> = ParametricRealField(None, mutable_values=True)
         sage: a <= 2
-        Traceback (most recent call last):
-        ...
-        FactorUndetermined: a cannot be evaluated because the test point is not complete
-        sage: K.assume_comparison(a.sym(), operator.le, 3)
+        True
+        sage: K.<a, b> = ParametricRealField(None, mutable_values=True)
+        sage: K.assume_comparison(a.sym(), operator.le, 2)
 
-    Partial test point mode::
+    Comparisons with test-points that are partially defined are supported. Comparisons made in
+    unspecified variables are assumed to be true::
 
         sage: K.<a, b> = ParametricRealField([None, 1], mutable_values=True)
         sage: a <= 2
-        Traceback (most recent call last):
-        ...
-        FactorUndetermined: a cannot be evaluated because the test point is not complete
+        True
         sage: b <= 11
         True
-
     """
     Element = ParametricRealFieldElement
 
@@ -210,7 +196,7 @@ class ParametricRealField(Field):
         self._big_cells = big_cells
 
         self._zero_element = ParametricRealFieldElement(self, 0)
-        self._one_element =  ParametricRealFieldElement(self, 1)
+        self._one_element = ParametricRealFieldElement(self, 1)
         ## REFACTOR: Maybe replace this by an instance of BasicSemialgebraicSet_eq_lt_le_sets - but careful - this class right now assumes polynomials
         self._eq = set([])
         self._lt = set([])
@@ -312,19 +298,19 @@ class ParametricRealField(Field):
             TypeError: unsupported operand parent(s)...
 
         Test that real number field elements can be upgraded to ``ParametricRealFieldElement``s.
-        Note that this requires setting up the ParametricRealField with a specific base ring, 
+        Note that this requires setting up the ParametricRealField with a specific base ring,
         because there is no common parent of QQ(x) and a RealNumberField``::
 
             sage: sqrt2, = nice_field_values([sqrt(2)])
             sage: K.<f> = ParametricRealField([0], base_ring=sqrt2.parent())
             sage: f + sqrt2
-            (f + (a))~
+            (f + a)~
 
         This currently does not work for Sage's built-in embedded number field elements...
         """
         if isinstance(S, ParametricRealField) and self is not S:
             return None
-        if  S is sage.interfaces.mathematica.MathematicaElement or isinstance(S, RealNumberField_absolute) or isinstance(S, RealNumberField_quadratic) or AA.has_coerce_map_from(S):
+        if S is sage.interfaces.mathematica.MathematicaElement or isinstance(S, RealNumberField_absolute) or isinstance(S, RealNumberField_quadratic) or AA.has_coerce_map_from(S):
             # Does the test with MathematicaElement actually work?
             # We test whether S coerces into AA. This rules out inexact fields such as RDF.
             return True
@@ -410,9 +396,7 @@ class ParametricRealField(Field):
             ....:     with K.temporary_assumptions():
             ....:         K.assume_comparison(a.sym(), operator.le, 3)
             ....:         a <= 4
-            Traceback (most recent call last):
-            ...
-            FactorUndetermined: a cannot be evaluated because the test point is not complete...
+            True
         """
         self._values = [ None for n in self._names ]
 
@@ -563,6 +547,15 @@ class ParametricRealField(Field):
                 pass
         raise FactorUndetermined("{} cannot be evaluated because the test point is not complete".format(fac))
 
+    def _partial_eval_factor(self, fac):
+        """
+        Partially evaluate ``fac`` on the test point.
+
+        This function is only intended to be called after ``FactorUndetermined`` is raised from ``_eval_factor``.
+        """
+        val_dict = {sym:val for sym, val in zip(fac.parent().gens() , self._values) if val is not None}
+        return fac.subs(val_dict)
+
     def _factor_sign(self, fac):
         """
         Determine the sign of ``fac`` evaluated on the test point.
@@ -596,15 +589,14 @@ class ParametricRealField(Field):
             ...
             ParametricRealFieldInconsistencyError: New constant constraint...
 
-        TESTS for consistency checks when testpoint values is None::
+        When a test point value is ``None``, consistency is delegated to the underlying BSA::
 
             sage: K.<a,b> = ParametricRealField([2, 1], big_cells=True, mutable_values=True, allow_refinement=False)
             sage: assert b>0
             sage: K.remove_test_point()
             sage: K.assume_comparison(b.sym(), operator.lt, 0)
-            Traceback (most recent call last):
-            ...
-            ParametricRealFieldInconsistencyError...
+            sage: K._bsa
+            BasicSemialgebraicSet_veronese(BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(Constraint_System {-1==0}, names=[x0]), polynomial_map=[b])
 
         User code should not rely on whether assumptions regarding
         the nonvanishing of denominators are recorded::
@@ -854,6 +846,15 @@ class ParametricRealField(Field):
                     raise ParametricRealFieldInconsistencyError("New constant constraint {} {} {} is not satisfied".format(lhs, op, rhs))
                 else:
                     raise ParametricRealFieldInconsistencyError("New constraint {} {}  {} is not satisfied by the test point".format(lhs, op, rhs))
+        else: #A numerical evaluation of the expression has failed. Assume the partial evaluation of the expression holds.
+            comparison_val_or_expr = self._partial_eval_factor(comparison)
+            if comparison_val_or_expr in base_ring:
+                if not op(comparison_val_or_expr, 0):
+                    raise ParametricRealFieldInconsistencyError("New constant constraint {} {} {} is not satisfied".format(lhs, op, rhs))
+            else: # comparision_val_or_expr is symbolic expression, assume the comparison here is comparionsion_val_or_expr.
+                comparison = comparison_val_or_expr
+                self.record_factor(comparison, op)
+                return
         if comparison in base_ring:
             return
         if comparison.denominator() == 1 and comparison.numerator().degree() == 1:
@@ -884,7 +885,8 @@ class ParametricRealField(Field):
                         return
         if len(factors) == 1 and factors[0][1] == 1 and comparison_val is not None:
             the_fac, d = factors[0]
-            the_sign = sign(factors.unit() * comparison_val) 
+            the_sign = sign(factors.unit() * comparison_val)
+
             def factor_sign(fac):
                 if fac == the_fac:
                     return the_sign
@@ -1050,7 +1052,9 @@ class ParametricRealField(Field):
                             eq_factors.append(new_fac)
                     except FactorUndetermined:
                         eq_factors.append(new_fac) #??? or pass?
-                        #This doesn't matter for now because under this branch  self._allow_refinement=True, and testpoint value is not removed, so  the exception FactorUndetermined should not happen.  In the future if the code in big_cells_impl changes and this exception happens, it is still safer to eq_factors.append(new_fac). This could result in smaller cell (which is allowed). With "pass" it may raise error in corner cases, such as when the removed testpoint value was the only eq factor without which there would be a sign contradiction.
+                        #This doesn't matter for now because under this branch  self._allow_refinement=True, and testpoint value is not removed, so  the exception FactorUndetermined should not happen.
+                        #In the future if the code in big_cells_impl changes and this exception happens, it is still safer to eq_factors.append(new_fac). T
+                        #This could result in smaller cell (which is allowed). With "pass" it may raise error in corner cases, such as when the removed testpoint value was the only eq factor without which there would be a sign contradiction.
                 for new_fac in lt_factors:
                     self.record_factor(new_fac, operator.le)
                 if not self._big_cells:
@@ -1058,7 +1062,7 @@ class ParametricRealField(Field):
                         self.record_factor(new_fac, operator.eq)
                 else: #self._big_cells is True, self._allow_refinement is True.
                     undecided_eq = []
-                    for new_fac in eq_factors:
+                    for new_fac in eq_factors: # eq_factors is not guantreed to be non empty
                         if self.is_factor_known(new_fac, operator.le):
                             unit_sign = -unit_sign
                         elif not self.is_factor_known(-new_fac, operator.le):
@@ -1108,7 +1112,7 @@ class ParametricRealField(Field):
     def make_proof_cell(self, **opt):
         r"""
         Make a :class:`SemialgebraicComplexComponent` from a :class:`ParametricRealField`.
-        
+
         In **opt, one can provide: region_type, function, find_region_type, default_var_bound, bddbsa, kwds_dict.
 
         EXAMPLES::
@@ -1148,7 +1152,7 @@ class ParametricRealField(Field):
 ###############################
 def find_polynomial_map(eqs=[], poly_ring=None):
     """
-    BAD FUCNTION! It is used in 'mathematica' approach for non-linear case. Can we avoid it?
+    BAD FUNCTION! It is used in 'mathematica' approach for non-linear case. Can we avoid it?
     Return a polynomial map that eliminates linear variables in eqs, and a dictionary recording which equations were used to eliminate those linear variables.
     Assume that gaussian elimination has been performed by PPL.minimized_constraints() on the input list of equations eqs.
     It is only called in SemialgebraicComplex.add_new_component in the case polynomial_map is not provided but bddbsa has equations.
@@ -1210,8 +1214,10 @@ def find_polynomial_map(eqs=[], poly_ring=None):
 # Functions with ParametricRealField K
 ######################################
 
+
 from sage.misc.sageinspect import sage_getargspec, sage_getvariablename
 
+
 def read_default_args(function, **opt_non_default):
     r"""
     Return the default values of arguments of the function.
@@ -1239,7 +1245,7 @@ def read_default_args(function, **opt_non_default):
         default_args = {}
         if defaults is not None:
             for i in range(len(defaults)):
-                default_args[args[-i-1]]=defaults[-i-1]
+                default_args[args[-i-1]] = defaults[-i-1]
     for (opt_name, opt_value) in opt_non_default.items():
         if opt_name in default_args:
             default_args[opt_name] = opt_value
@@ -1298,7 +1304,7 @@ class SemialgebraicComplexComponent(SageObject):    # FIXME: Rename this to be m
         sage: sorted(component.bsa.lt_poly())
         [-x, 3*x - 4]
 
-    In ProofCell region_type should alreay consider the polynomial_map::
+    In ProofCell region_type should already consider the polynomial_map::
 
         sage: K.<x,y> = ParametricRealField([1,1/2])
         sage: assert(x == 2*y)
@@ -1332,7 +1338,11 @@ class SemialgebraicComplexComponent(SageObject):    # FIXME: Rename this to be m
             polynomial_map = find_polynomial_map(eqs=eqs, poly_ring=poly_ring)
             #self.bsa = K._bsa.section(polynomial_map, bsa_class='veronese', poly_ring=poly_ring)  # this is a bigger_bsa
             self.bsa = BasicSemialgebraicSet_veronese.from_bsa(BasicSemialgebraicSet_local(K._bsa.section(polynomial_map, poly_ring=poly_ring), self.var_value)) # TODO:, polynomial_map=list(poly_ring.gens()))
-            # WHY is this input polynomial_map sometimes not compatible with the variable elimination done in bddbsa? Because upstairs ppl bsa eliminates large x_i in the inequalities, and x_i doesn't necessarily correspond to the i-th variable in poly_ring. Since polynomial_map and v_dict were not given at the initialization of veronese, the variable first encounted in the constraints is considered as x0 by upstairs ppl bsa. # In old code, we fixed the order of upstairs variables by adding initial space dimensions. We don't do that in the current code. Instead, we take the section of bddbsa to eliminate the varibles in the equations. # Is the given bddbsa required to be veronese with upstairs being ppl_bsa? Convert it anyway. # It's the same as BasicSemialgebraicSet_veronese.from_bsa(bddbsa.section(self.polynomial_map), poly_ring=poly_ring)
+            # WHY is this input polynomial_map sometimes not compatible with the variable elimination done in bddbsa?
+            # Because upstairs ppl bsa eliminates large x_i in the inequalities, and x_i doesn't necessarily correspond to the i-th variable in poly_ring.
+            # Since polynomial_map and v_dict were not given at the initialization of veronese, the variable first encountered, in the constraints is considered as x0 by upstairs ppl bsa.
+            # In old code, we fixed the order of upstairs variables by adding initial space dimensions. We don't do that in the current code. Instead, we take the section of bddbsa to eliminate the variables in the equations.
+            # Is the given bddbsa required to be veronese with upstairs being ppl_bsa? Convert it anyway. # It's the same as BasicSemialgebraicSet_veronese.from_bsa(bddbsa.section(self.polynomial_map), poly_ring=poly_ring)
             self.bddbsa = BasicSemialgebraicSet_veronese.from_bsa(BasicSemialgebraicSet_local(bddbsa.section(polynomial_map, poly_ring=poly_ring), self.var_value))
             # Taking section forgets the equations. Then add back the equations  # Finally self.bsa should be the same as K._bsa, but its inequalities don't have variables eliminated by polynomial map, so that heuristic wall crossing can be done later.
             for i in range(len(self.var_name)):
@@ -1400,7 +1410,7 @@ class SemialgebraicComplexComponent(SageObject):    # FIXME: Rename this to be m
             else:
                 ptcolor = 'white'
             if (xmin <= pt[0] <= xmax) and (ymin <= pt[1] <= ymax):
-                g += point(pt, color = ptcolor, size = 2, zorder=10)
+                g += point(pt, color=ptcolor, size=2, zorder=10)
         return g
 
     def find_neighbour_candidates(self, flip_ineq_step, wall_crossing_method='heuristic', goto_lower_dim=False, pos_poly=None):
@@ -1412,13 +1422,13 @@ class SemialgebraicComplexComponent(SageObject):    # FIXME: Rename this to be m
         - if goto_lower_dim=False, the cell is considered as its formal closure, so no recursion into test points in lower dimensional cells.
         - pos_poly is a polynomial. The return test point must satisfy pos_poly(new test point) > 0.
 
-        OUTPUT new_points is a dictionary of dictionaries. The new_points[i] is a dictionay whose keys = candidate neighbour testpoints, values = (bddbsa whose eq_poly has i elements, polynomial_map, no_crossing_l) of the candidate neighbour cell that contains the candidate neighbour testpoint. bddbsa is recorded so that (1) complex.bddbsa is always respected; and (2) can recursively go into lower dimensional cells. polynomial_map is recorded and passed to the constructor of the neighbour cell. no_crossing is passed to the neighour cell for its find_neighbour_candidates method. We no longer update self.bsa by removing (obvious) redundant eq, lt, le constraints from its description at the end, even when 'mathematica' is used.
+        OUTPUT new_points is a dictionary of dictionaries. The new_points[i] is a dictionary whose keys = candidate neighbour testpoints, values = (bddbsa whose eq_poly has i elements, polynomial_map, no_crossing_l) of the candidate neighbour cell that contains the candidate neighbour testpoint. bddbsa is recorded so that (1) complex.bddbsa is always respected; and (2) can recursively go into lower dimensional cells. polynomial_map is recorded and passed to the constructor of the neighbour cell. no_crossing is passed to the neighbour cell for its find_neighbour_candidates method. We no longer update self.bsa by removing (obvious) redundant eq, lt, le constraints from its description at the end, even when 'mathematica' is used.
         """
         bsa_eq_poly = list(self.bsa.eq_poly())
         bsa_le_poly = list(self.bsa.le_poly())
         bsa_lt_poly = list(self.bsa.lt_poly())
         num_eq = len(bsa_eq_poly) #was len(list(self.bddbsa.eq_poly()))
-        new_points = {}  #dictionary with key=num_eq, value=dictionay of pt: (bddbsa, polynomial_map).
+        new_points = {}  #dictionary with key=num_eq, value=dictionary of pt: (bddbsa, polynomial_map).
         #bddbsa = copy(self.bddbsa)
         #for l in bsa_eq_poly: # should be already in bddbsa
         #    bddbsa.add_polynomial_constraint(l, operator.eq)
@@ -1447,7 +1457,7 @@ class SemialgebraicComplexComponent(SageObject):    # FIXME: Rename this to be m
                 pt = find_point_flip_ineq_heuristic(self.var_value, l, list(bsa.lt_poly())+list(bsa.le_poly()), flip_ineq_step)
                 if pt is not None:
                     # Find a new point, use polynomial map to recover the values of those eliminated variables.
-                    pt_across_wall = tuple(p(pt) for p in self.polynomial_map)                 
+                    pt_across_wall = tuple(p(pt) for p in self.polynomial_map)
             if wall_crossing_method == 'mathematica' or wall_crossing_method == 'heuristic_then_mathematica' and (pt_across_wall is None):
                 bsa_mathematica = bsa.formal_relint(bsa_class='mathematica') # was BasicSemialgebraicSet_mathematica.from_bsa(bsa)
                 bsa_mathematica.add_polynomial_constraint(l, operator.gt)
@@ -1534,7 +1544,9 @@ class SemialgebraicComplexComponent(SageObject):    # FIXME: Rename this to be m
                         if bsa_section.upstairs()._polyhedron.is_empty():
                             has_pt_on_wall = False
                         else:
-                            lts = []; les = []; eqs = []
+                            lts = []
+                            les = []
+                            eqs = []
                             for ll in list(bsa_section.lt_poly()):
                                 factors = ll.factor()
                                 if len(factors) == 1:
@@ -1662,7 +1674,7 @@ class ProofCell(SemialgebraicComplexComponent, Classcall):
         sage: C_copy._init_args == C._init_args
         True
 
-    (We do not test for equality C_copy == C -- we have not even decided yet what the semantics 
+    (We do not test for equality C_copy == C -- we have not even decided yet what the semantics
     of equality of bsa is.)
     """
 
@@ -1772,8 +1784,10 @@ class ProofCell(SemialgebraicComplexComponent, Classcall):
         super(ProofCell, self).__init__(K, region_type, bddbsa, polynomial_map)
         self.family = family
 
+
 from collections import OrderedDict
 
+
 class SemialgebraicComplex(SageObject):
     r"""
     A proof complex for parameter space analysis.
@@ -1807,7 +1821,7 @@ class SemialgebraicComplex(SageObject):
 
         sage: complex.is_complete()                           # optional - mathematica
         True
-        
+
     Example with non-linear wall::
 
         sage: complex = SemialgebraicComplex(lambda x,y: max(x,y^2), ['x','y'], find_region_type=result_symbolic_expression, default_var_bound=(-3,3))  # optional - mathematica
@@ -1902,7 +1916,7 @@ class SemialgebraicComplex(SageObject):
           to a "type" of the parameter region, for example:
 
           - :func:`find_region_type_igp` (the default). The result of the computation is a Gomory-Johnson
-            function `h`; it is passed to :func:`find_region_type_igp` as 2nd arg, 
+            function `h`; it is passed to :func:`find_region_type_igp` as 2nd arg,
             and then :func:`find_region_type_igp`which classifies the region of the
             parameter by returning one of the strings
             ``'is_constructible'``, ``'not_constructible'``,
@@ -1981,7 +1995,7 @@ class SemialgebraicComplex(SageObject):
         r"""
         Return a random point that satisfies var_bounds and is in self.bsa.
 
-        - If var_bounds is not specified, self.default_var_bound is taken. 
+        - If var_bounds is not specified, self.default_var_bound is taken.
         - var_bounds can be a list of 2-tuples whose length equals to the number of parameters, or lambda functions.
         - It is used in random shooting method for functions like ``dg_2_step_mir``, which involve floor/ceil operations. We try to plot one layer for each n = floor(...) and superimpose the layers at the end to get the whole picture.
 
@@ -2007,11 +2021,11 @@ class SemialgebraicComplex(SageObject):
                     x = QQ(uniform(self.default_var_bound[0], self.default_var_bound[1]))
                 else:
                     if hasattr(var_bounds[i][0], '__call__'):
-                        l =  var_bounds[i][0](*var_value)
+                        l = var_bounds[i][0](*var_value)
                     else:
                         l = var_bounds[i][0]
                     if hasattr(var_bounds[i][1], '__call__'):
-                        u =  var_bounds[i][1](*var_value)
+                        u = var_bounds[i][1](*var_value)
                     else:
                         u = var_bounds[i][1]
                     if l > u:
@@ -2072,12 +2086,12 @@ class SemialgebraicComplex(SageObject):
         for c in self.components:
             if var_value in c.bsa:
                 yield c
-        
+
     def find_uncovered_random_point(self, var_bounds=None, max_failings=10000):
         r"""
         Return a random point that satisfies the bounds and is uncovered by any cells in the complex.
-        Return ``None`` if the number of attemps > max_failings.
- 
+        Return ``None`` if the number of attempts > max_failings.
+
         EXAMPLES::
 
             sage: from cutgeneratingfunctionology.igp import *
@@ -2106,7 +2120,7 @@ class SemialgebraicComplex(SageObject):
 
         The argument formal_closure whether inequalities are treated as <= 0 or as < 0.
         If such point does not exist, return ``None``.
- 
+
         EXAMPLES::
 
             sage: from cutgeneratingfunctionology.igp import *
@@ -2183,18 +2197,18 @@ class SemialgebraicComplex(SageObject):
             self.points_to_test[num_eq] = OrderedDict()
             if not num_eq in self.tested_points:
                 self.tested_points[num_eq] = set([])
-        new_component = self._cell_class(self.family, var_value, 
+        new_component = self._cell_class(self.family, var_value,
                                          find_region_type=self.find_region_type, bddbsa=bddbsa, polynomial_map=polynomial_map)
 
         new_num_eq = len(list(new_component.bsa.eq_poly()))
-        if  new_num_eq > num_eq:
-            logging.warning("The cell around %s defined by %s has more equations than boundary %s" %(new_component.var_value, new_component.bsa, bddbsa))
+        if new_num_eq > num_eq:
+            logging.warning("The cell around %s defined by %s has more equations than boundary %s" % (new_component.var_value, new_component.bsa, bddbsa))
             #import pdb; pdb.set_trace()
-            # bsa is lower dimensional as it has more equations than bddbsa, 
+            # bsa is lower dimensional as it has more equations than bddbsa,
             # so we try to perturb the testpoint to obtain a
             # new testpoint in bddbsa that does not fall into a lower dimensional cell.
             # Heuristic code using gradient desecent.  #FIXME.
-            for l in (set(new_component.bsa.eq_poly())- set(bddbsa.eq_poly())):
+            for l in (set(new_component.bsa.eq_poly()) - set(bddbsa.eq_poly())):
                 ineqs = list(new_component.bddbsa.lt_poly())+list(new_component.bddbsa.le_poly())
                 pts = [find_point_flip_ineq_heuristic(var_value, l, ineqs, 1/2017), find_point_flip_ineq_heuristic(var_value, -l, ineqs, 1/2017)]
                 for pt in pts:
@@ -2255,7 +2269,7 @@ class SemialgebraicComplex(SageObject):
         Plot the complex and store the graph.
 
         - If restart is ``False``, plot the newly added cells on top of the last graph; otherwise, start a new graph.
-        - If slice_value is given, it is either a polynomial_map that defines a section, or a list of fixed parameter values with two of them being None. Plot the section. 
+        - If slice_value is given, it is either a polynomial_map that defines a section, or a list of fixed parameter values with two of them being None. Plot the section.
         - plot_points controls the quality of the plotting.
 
         EXAMPLES::
@@ -2292,7 +2306,7 @@ class SemialgebraicComplex(SageObject):
             self.graph.xmin(kwds['xmin'])
             xmin = kwds['xmin']
         else:
-            xmin = self.default_var_bound[0]   # special treatement in the case goto_lower_dim which uses bsa.plot() instead of component.plot() because zorder is broken in region_plot/ContourPlot.
+            xmin = self.default_var_bound[0]   # special treatment in the case goto_lower_dim which uses bsa.plot() instead of component.plot() because zorder is broken in region_plot/ContourPlot.
         if 'xmax' in kwds:
             self.graph.xmax(kwds['xmax'])
             xmax = kwds['xmax']
@@ -2311,7 +2325,7 @@ class SemialgebraicComplex(SageObject):
         # # FIXME: zorder is broken in region_plot/ContourPlot.
         # for c in self.components[self.num_plotted_components::]:
         #     num_eq = len(list(c.bsa.eq_poly()))
-        #     gc = c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=goto_lower_dim, zorder=num_eq, **kwds)  
+        #     gc = c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=goto_lower_dim, zorder=num_eq, **kwds)
         #     if gc: # need this because (empty g + empty gc) forgets about xmin xmax ymin ymax.
         #         self.graph += gc
         # Workaround.
@@ -2333,11 +2347,11 @@ class SemialgebraicComplex(SageObject):
                         new_bsa.add_polynomial_constraint(l, operator.eq)
                         self.graph += new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, color=color, fill_color=color, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
         for c in components:
-            if len(list(c.bsa.eq_poly()))==1:
+            if len(list(c.bsa.eq_poly())) == 1:
                 self.graph += c.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, default_var_bound=self.default_var_bound, goto_lower_dim=False, zorder=0, **kwds)
         if goto_lower_dim:
             for c in components:
-                if len(list(c.bsa.eq_poly()))==1:
+                if len(list(c.bsa.eq_poly())) == 1:
                     color = find_region_color(c.region_type)
                     for l in c.bsa.lt_poly():
                         new_bsa = BasicSemialgebraicSet_eq_lt_le_sets(eq=list(c.bsa.eq_poly())+[l], lt=[ll for ll in c.bsa.lt_poly() if ll != l], le=list(c.bsa.le_poly()))
@@ -2347,9 +2361,9 @@ class SemialgebraicComplex(SageObject):
                         new_bsa.add_polynomial_constraint(l, operator.eq)
                         self.graph += new_bsa.plot(alpha=alpha, plot_points=plot_points, slice_value=slice_value, color=color, fill_color=color, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
         for c in components:
-            if len(list(c.bsa.eq_poly()))==2:
+            if len(list(c.bsa.eq_poly())) == 2:
                 ptcolor = find_region_color(c.region_type)
-                self.graph += point(c.var_value, color = ptcolor, zorder=10)
+                self.graph += point(c.var_value, color=ptcolor, zorder=10)
         self.num_plotted_components = len(self.components)
         return self.graph
 
@@ -2386,7 +2400,7 @@ class SemialgebraicComplex(SageObject):
         - var_value: the starting point. If not given, start with a random point.
         - flip_ineq_step: a small positive number that controls the step length in wall-crossing.
         - check_completion: When check_completion is ``True``, after bfs terminates, check whether the entire parameter space is covered by cells, using Mathematica's ``FindInstance`` (if max_failings=0) or random shooting (if max_failings>0). If an uncovered point has been found, restart the bfs from this point.
-        - wall_crossing_method: 'mathematica' or 'heuristic' or 'heuristic_then_mathematica' 
+        - wall_crossing_method: 'mathematica' or 'heuristic' or 'heuristic_then_mathematica'
               - wall_crossing_method='heuristic_then_mathematica': try heuristic method first. If it does not find a new testpoint, then use Mathematica.
         - goto_lower_dim: whether lower dimensional cells are considered. If goto_lower_dim is ``False`` or if goto_lower_dim is ``True`` and wall_crossing method is 'heuristic' but the wall is non-linear, then find new testpoint across the wall only.
 
@@ -2440,10 +2454,10 @@ class SemialgebraicComplex(SageObject):
                 uncovered_pt = self.find_uncovered_random_point(max_failings=max_failings)
             if uncovered_pt is not None:
                 logging.warning("After bfs, the complex has uncovered point %s." % (uncovered_pt,))
-                self.bfs_completion(var_value=uncovered_pt, \
-                                    flip_ineq_step=flip_ineq_step, \
-                                    check_completion=check_completion, \
-                                    wall_crossing_method=wall_crossing_method, \
+                self.bfs_completion(var_value=uncovered_pt,
+                                    flip_ineq_step=flip_ineq_step,
+                                    check_completion=check_completion,
+                                    wall_crossing_method=wall_crossing_method,
                                     goto_lower_dim=goto_lower_dim)
 
     def is_complete(self, formal_closure=False):
@@ -2545,9 +2559,9 @@ def gradient(ineq):
         [3*z^2 + 6*z]
     """
     if hasattr(ineq, 'gradient'):
-       return ineq.gradient()
+        return ineq.gradient()
     else:
-       return [ineq.derivative()]
+        return [ineq.derivative()]
 
 ####################################
 # Find region type and region color
@@ -2680,8 +2694,8 @@ def return_result(field, result):
 def result_concrete_value(field, result):
     r"""
     Return the concrete values in result as a tuple. See also ``result_symbolic_expression()``.
- 
-    This function can provided to ``find_region_type`` when setting up a :class:`SemialgebraicComplex`. 
+
+    This function can provided to ``find_region_type`` when setting up a :class:`SemialgebraicComplex`.
     In this way, one can compare result of type :class:`ParametricRealFieldElement` or list of :class:`ParametricRealFieldElement`
     with the previous elements in ``region_type_color_map`` which do not necessarily have the same parent.
 
@@ -2703,8 +2717,8 @@ def result_concrete_value(field, result):
 def result_symbolic_expression(field, result):
     r"""
     Return the symbolic expressions in result as a tuple
- 
-    This function can provided to ``find_region_type`` when setting up a :class:`SemialgebraicComplex`. 
+
+    This function can provided to ``find_region_type`` when setting up a :class:`SemialgebraicComplex`.
     In this way, one can compare result of type :class:`ParametricRealFieldElement` or list of :class:`ParametricRealFieldElement`
     with the previous elements in ``region_type_color_map`` which do not necessarily have the same parent.
 
@@ -2728,7 +2742,7 @@ def result_symbolic_expression(field, result):
         sage: vol1 == vol2
         False
         sage: sym_exp1 == sym_exp2
-        True    
+        True
     """
     symbolic_expression = tuple(elt._sym if hasattr(elt, '_sym') else elt for elt in flatten([result]))
     return symbolic_expression
@@ -2753,7 +2767,7 @@ def find_region_type_igp_extreme_big_cells(K, h):
         h = copy(hcopy)
         for x in h.values_at_end_points():
             if (x < 0) or (x > 1):
-                is_extreme  = False
+                is_extreme = False
                 break
     if not is_extreme:
         assert (x < 0) or (x > 1)
@@ -2799,7 +2813,7 @@ def find_region_type_igp_extreme_big_cells(K, h):
         ucs = generate_uncovered_components(h)
         f = find_f(h)
         for uncovered_pt in [f/2, (f+1)/2]:
-            if any((i[0] == uncovered_pt or i[1] == uncovered_pt) for uc in ucs for i in uc if len(uc)==2):
+            if any((i[0] == uncovered_pt or i[1] == uncovered_pt) for uc in ucs for i in uc if len(uc) == 2):
                 uncovered_pts = [uncovered_pt]
                 is_extreme = False
                 break
@@ -2833,7 +2847,8 @@ def find_region_type_igp_extreme_big_cells(K, h):
         return False
     return True
 
-region_type_color_map = [('not_constructible', 'lightgrey'), ('is_constructible', 'black'), ('not_minimal', 'orange'), ('is_minimal', 'darkgrey'),('not_extreme', 'green'), ('is_extreme', 'blue'), ('stop', 'grey'), (True, 'blue'), (False, 'red'), ('constructible', 'darkgrey'), ('extreme', 'red')]
+
+region_type_color_map = [('not_constructible', 'lightgrey'), ('is_constructible', 'black'), ('not_minimal', 'orange'), ('is_minimal', 'darkgrey'), ('not_extreme', 'green'), ('is_extreme', 'blue'), ('stop', 'grey'), (True, 'blue'), (False, 'red'), ('constructible', 'darkgrey'), ('extreme', 'red')]
 
 def find_region_color(region_type):
     r"""
@@ -2878,11 +2893,11 @@ def color_of_ith_region_type(i):
 
 def find_point_flip_ineq_heuristic(current_var_value, ineq, ineqs, flip_ineq_step):
     r"""
-    The current_var_value satisfies that l(current_var_value) <= 0 for l=ineq and for every l in ineqs, 
+    The current_var_value satisfies that l(current_var_value) <= 0 for l=ineq and for every l in ineqs,
     where ineq is a polynomial and ineqs is a list of polynomials.
 
     Use heuristic method (gradient descent method with given small positive step length flip_ineq_step)
-    to find a new_point (type is tuple) such that 
+    to find a new_point (type is tuple) such that
     ineq(new_point) > 0, l(new_point) < 0 for all l in ineqs
     Return new_point, or ``None`` if it fails to find one.
 
@@ -2968,7 +2983,7 @@ def find_point_flip_ineq_heuristic(current_var_value, ineq, ineqs, flip_ineq_ste
         sage: all(l(pt) < 0 for l in ineqs) and ineq(pt)>0
         True
 
-    Bug examle from positive definite matrix [a, b; b, 1/4], where ineq is very negative at the test point. Make big moves first, then small moves::
+    Bug example from positive definite matrix [a, b; b, 1/4], where ineq is very negative at the test point. Make big moves first, then small moves::
 
         sage: P.<a,b>=QQ[]; current_var_value = (5, 4); ineq = -4*b^2 + a; ineqs = [-a]; flip_ineq_step=1/100
         sage: pt = find_point_flip_ineq_heuristic(current_var_value, ineq, ineqs, flip_ineq_step); # got pt (30943/6018, 17803/15716)
@@ -2988,7 +3003,7 @@ def find_point_flip_ineq_heuristic(current_var_value, ineq, ineqs, flip_ineq_ste
             ineq_value = ineq(*current_point)
             try_before_fail -= 1
             # print (current_point, RR(ineq_value))
-    try_before_fail =  max(ceil(2/flip_ineq_step), 2000)  # define maximum number of walks. Considered ceil(-2 * ineq_value /flip_ineq_step) but it is too slow in the impossible cases. Added a loop with 2000 times step length when ineq_value is very negative.
+    try_before_fail = max(ceil(2/flip_ineq_step), 2000)  # define maximum number of walks. Considered ceil(-2 * ineq_value /flip_ineq_step) but it is too slow in the impossible cases. Added a loop with 2000 times step length when ineq_value is very negative.
     while (ineq_value <= 1e-10) and (try_before_fail > 0):
         ineq_direction = vector(g(*current_point) for g in ineq_gradient)
         if ineq.degree() == 1:
@@ -3012,8 +3027,8 @@ def find_point_flip_ineq_heuristic(current_var_value, ineq, ineqs, flip_ineq_ste
 
 def adjust_pt_to_satisfy_ineqs(current_point, ineq, strict_ineqs, nonstrict_ineqs, flip_ineq_step):
     r"""
-    Walk from current_point (type=vector) in the direction perpendicular to 
-    the gradient of ineq with small positive step length flip_ineq_step, 
+    Walk from current_point (type=vector) in the direction perpendicular to
+    the gradient of ineq with small positive step length flip_ineq_step,
     while maintaining the value of ineq(*current_point) which is >= 0.
     until get a new point such that l(new point)<0 for all l in strict_ineqs and l(new point)<=0 for all l in nonstrict_ineqs
 
@@ -3035,7 +3050,7 @@ def adjust_pt_to_satisfy_ineqs(current_point, ineq, strict_ineqs, nonstrict_ineq
         sage: pt = adjust_pt_to_satisfy_ineqs(vector([11/8, 7/8]), a+b-2, [-a+b^2], [a-1], 1/4); pt is None
         True
 
-    Bug example in cpl Cell 9 with test point (499/1250, 488072439572/4866126017667). Without converting input to QQ, output was (0.333000000000000, 0.111333333333333) with  -2*f - 3*z + 1 = -5.55111512312578e-17 , the QQ of the output=(333/1000, 167/1500) has -2*f - 3*z + 1 == 0. Revise the code to take QQ input current point:: 
+    Bug example in cpl Cell 9 with test point (499/1250, 488072439572/4866126017667). Without converting input to QQ, output was (0.333000000000000, 0.111333333333333) with  -2*f - 3*z + 1 = -5.55111512312578e-17 , the QQ of the output=(333/1000, 167/1500) has -2*f - 3*z + 1 == 0. Revise the code to take QQ input current point::
 
         sage: P.<f,z>=QQ[]
         sage: current_point = vector([0.333000000000000, 0.100300000000000]); ineq = -3*f + 1; strict_ineqs = [2*f + 2*z - 1, f + 5*z - 1, -f - 6*z + 1, -2*f - 3*z + 1]; nonstrict_ineqs = []; flip_ineq_step = 1/1000
@@ -3054,8 +3069,8 @@ def adjust_pt_to_satisfy_ineqs(current_point, ineq, strict_ineqs, nonstrict_ineq
 
     Bug example in cpl bigcell 16 with test point (12219/26000, 24/1625). Redo with QQ had infinite loop. Bug comes from find_neighbour_point where it calls bsa_section.upstairs()._polyhedron.is_empty(), which is not strong enough. If we could test bsa_section is empty (perhaps by tighten_upstairs_by_mccormick), then this example should not appear::
 
-        sage: P.<f,z>=QQ[]; 
-        sage: current_point = vector((71582788/143165577, 4673/377000)) # came from vector((RR(70727/150800), RR(4673/377000))), 
+        sage: P.<f,z>=QQ[];
+        sage: current_point = vector((71582788/143165577, 4673/377000)) # came from vector((RR(70727/150800), RR(4673/377000))),
         sage: ineq=None; strict_ineqs=[2*f - 1, -9*f + 2]; nonstrict_ineqs=[4*f^2 - 4*f + 1]; flip_ineq_step=1/1000
         sage: pt = adjust_pt_to_satisfy_ineqs(current_point, None, strict_ineqs, nonstrict_ineqs, flip_ineq_step=1/1000); pt is None #long time
         True
@@ -3070,7 +3085,7 @@ def adjust_pt_to_satisfy_ineqs(current_point, ineq, strict_ineqs, nonstrict_ineq
     #current_point is a vector
     if ineq is not None:
         ineq_gradient = gradient(ineq)
-    if all(x.parent()==QQ for x in current_point):
+    if all(x.parent() == QQ for x in current_point):
         max_walks = min(ceil(2/flip_ineq_step), 20)
     else:
         max_walks = min(ceil(2/flip_ineq_step), 200) #1000? # define maximum number of walks.
@@ -3132,7 +3147,7 @@ def adjust_pt_to_satisfy_ineqs(current_point, ineq, strict_ineqs, nonstrict_ineq
             return None
     if ineq is not None and ineq(*current_point) < 0:
         return None
-    if all(x.parent()==QQ for x in current_point):
+    if all(x.parent() == QQ for x in current_point):
         return tuple(current_point)
     else:
         prec = 30  # We hope to have small denominator for the new point, so we set precision in bits = 30 is about 8 digits.
@@ -3193,13 +3208,14 @@ def embed_function_into_family(given_function, parametric_family, check_completi
         var_name = []
         var_value = []
         for (name, value) in default_args.items():
-            if not isinstance(value, bool) and not value is None:
+            if not isinstance(value, bool) and value is not None:
                 try:
                     RR(value)
                     var_name.append(name)
                     var_value.append(value)
                 except:
                     pass
+
     def frt(K, h):
         if h is None:
             return False
@@ -3242,6 +3258,7 @@ def embed_function_into_family(given_function, parametric_family, check_completi
     # plot_cpl_components(complex.components)
     return {}
 
+
 """
 EXAMPLES::
 
diff --git a/cutgeneratingfunctionology/spam/exceptions.py b/cutgeneratingfunctionology/spam/exceptions.py
new file mode 100644
index 000000000..997328c34
--- /dev/null
+++ b/cutgeneratingfunctionology/spam/exceptions.py
@@ -0,0 +1,14 @@
+class FactorUndetermined(Exception): # FactorUndetermined is raised when an expression can not be evaluated with a test point.
+    pass
+
+
+class ParametricRealFieldFrozenError(ValueError):
+    pass
+
+
+class ParametricRealFieldInconsistencyError(ValueError):
+    pass
+
+
+class ParametricRealFieldRefinementError(ValueError):
+    pass
diff --git a/cutgeneratingfunctionology/spam/parametric_real_field_element.py b/cutgeneratingfunctionology/spam/parametric_real_field_element.py
index c6688eeeb..0d9bb2e86 100644
--- a/cutgeneratingfunctionology/spam/parametric_real_field_element.py
+++ b/cutgeneratingfunctionology/spam/parametric_real_field_element.py
@@ -9,8 +9,10 @@
 from sage.rings.real_mpfr import RR
 from sage.functions.other import ceil, floor
 from sage.functions.generalized import sign
+from cutgeneratingfunctionology.spam.exceptions import FactorUndetermined
 import operator
 
+
 def richcmp_op_negation(op):
     if op == op_LT:
         return op_GE
@@ -27,6 +29,7 @@ def richcmp_op_negation(op):
     else:
         raise ValueError("{} is not a valid richcmp operator".format(op))
 
+
 def format_richcmp_op(op):
     if op == op_LT:
         return '<'
@@ -43,6 +46,7 @@ def format_richcmp_op(op):
     else:
         raise ValueError("{} is not a valid richcmp operator".format(op))
 
+
 class ParametricRealFieldElement(FieldElement):
     r"""
     A :class:`ParametricRealFieldElement` stores a symbolic expression of the parameters in the problem and a concrete value, which is the evaluation of this expression on the given parameter tuple.
@@ -77,7 +81,14 @@ def val(self):
         try:
             return self._val
         except AttributeError:
-            return self.parent()._eval_factor(self._sym)
+            try:
+                return self.parent()._eval_factor(self._sym)
+            except FactorUndetermined:
+                possible_val = self.parent()._partial_eval_factor(self._sym)
+                if possible_val in possible_val.base_ring():
+                    return possible_val
+                else:
+                    raise FactorUndetermined("{} cannot be evaluated because the test point is not complete".format(self.sym()))
 
     def _richcmp_(left, right, op):
         r"""
@@ -126,7 +137,10 @@ def _richcmp_(left, right, op):
             # shouldn't really happen, within coercion
             raise TypeError("comparing elements from different fields")
         if left.parent()._big_cells:
-            result = richcmp(left.val(), right.val(), op)
+            try:
+                result = richcmp((left-right).val(), 0, op)
+            except FactorUndetermined:  # Partial evaluation has happen, assume the result is True.
+                result = True
             if result:
                 true_op = op
             else:
@@ -150,13 +164,33 @@ def _richcmp_(left, right, op):
             return result
         else:
             # Traditional cmp semantics.
-            if (left.val() == right.val()):
-                left.parent().assume_comparison(right, operator.eq, left)
-            elif (left.val() < right.val()):
-                left.parent().assume_comparison(left, operator.lt, right)
+            try:
+                expr_val = (left-right).val()
+                if( expr_val == 0):
+                    left.parent().assume_comparison(right, operator.eq, left)
+                elif (expr_val < 0):
+                    left.parent().assume_comparison(left, operator.lt, right)
+                else:
+                    left.parent().assume_comparison(right, operator.lt, left)
+                return richcmp(left.val(), right.val(), op)
+            except FactorUndetermined:
+                # With a partial evaluation, assume the written inequality is true.
+                true_op = op
+            if true_op == op_LT:
+                left.parent().assume_comparison(left - right, operator.lt)
+            elif true_op == op_GT:
+                left.parent().assume_comparison(left - right, operator.gt)
+            elif true_op == op_EQ:
+                left.parent().assume_comparison(right - left, operator.eq)
+            elif true_op == op_LE:
+                left.parent().assume_comparison(left - right, operator.le)
+            elif true_op == op_GE:
+                left.parent().assume_comparison(left - right, operator.ge)
+            elif true_op == op_NE:
+                left.parent().assume_comparison(right - left, operator.ne)
             else:
-                left.parent().assume_comparison(right, operator.lt, left)
-            return richcmp(left.val(), right.val(), op)
+                raise ValueError("{} is not a valid richcmp operator".format(op))
+            return True
 
     def __abs__(self):
         """
@@ -399,6 +433,7 @@ def __hash__(self):
         """
         return hash(self.val())
 
+
 def is_parametric_element(x):
     # We avoid using isinstance here so that this is robust even if parametric.sage is reloaded.
     # For example, apparently in the test suite.