diff --git a/cutgeneratingfunctionology/igp/__init__.py b/cutgeneratingfunctionology/igp/__init__.py index a5a2da2ba..1f2db68c8 100644 --- a/cutgeneratingfunctionology/igp/__init__.py +++ b/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: diff --git a/cutgeneratingfunctionology/igp/functions.sage b/cutgeneratingfunctionology/igp/functions.sage index 8dfa6cc16..f12bc04be 100644 --- a/cutgeneratingfunctionology/igp/functions.sage +++ b/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""" diff --git a/cutgeneratingfunctionology/igp/minimal_function_cell_description.sage b/cutgeneratingfunctionology/igp/minimal_function_cell_description.sage new file mode 100644 index 000000000..591216fad --- /dev/null +++ b/cutgeneratingfunctionology/igp/minimal_function_cell_description.sage @@ -0,0 +1,872 @@ +from cutgeneratingfunctionology.igp import * +from cutgeneratingfunctionology.spam.basic_semialgebraic import EmptyBSA +from cutgeneratingfunctionology.spam.basic_semialgebraic_pplite import BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron +import csv +import logging +import os + +minimal_funciton_cell_description_logger = logging.getLogger("cutgeneratingfunctionology.igp.minimal_funciton_cell_description") +minimal_funciton_cell_description_logger.setLevel(logging.INFO) + +### Note to future reader: ### +### bkpt is assumed to be a breakpoint sequence of length n>= 2. +### Breakpoint sequence are sorted lists of real numbers in [0,1). +### A breakpoint sequences should always have 0 as an element. +### This is never strictly enforced in this file and it is assumed that +### the user is always provided a breakpoint sequence. + +class RepElemGenFailure(Exception): + pass + + +def nnc_poly_from_bkpt_sequence(bkpt, backend=None, log_paramateric_real_field = False,log_pw_functions = False): + r""" + Defines an NNC polyhedron P such that for all b in P the delta complex of b is isomoprhic to the delta complex of bkpt. + + INPUT: + - ``bkpt`` - sorted list of length 1 or more of sage type in the interval [0,1). + - ```backend`` - ``None``, ``str(pplite)`` + + OUTPUT: class::``BasicSemialgebraicSet_veronese`` + + EXAMPLES:: + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: nnc_poly_from_bkpt_sequence([0, 4/5]) + BasicSemialgebraicSet_veronese(BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(Constraint_System {x0==0, -x1+1>0, 2*x1-1>0}, names=[x0, x1]), polynomial_map=[lambda0, lambda1]) + """ + n = len(bkpt) + coord_names = [] + bkpt_vals = bkpt + vals = bkpt_vals[0:n] + bkpt_extd = list(bkpt)+[1] + if not log_paramateric_real_field: + parametric_logging_level = logging.getLogger("cutgeneratingfunctionology.igp.parametric").getEffectiveLevel() + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(logging.ERROR) + for i in range(0,n): + coord_names.append('lambda'+str(i)) + K = ParametricRealField(names=coord_names, values=vals, mutable_values=True, big_cells=True, default_backend=backend) + K.gens()[0] == 0 + for i in range(n-1): + K.gens()[i] < K.gens()[i+1] + K.gens()[n-1] < 1 + for i in range(n): + for j in range(n): + if bkpt[i]+bkpt[j] >= 1: + w = 1 + else: + w = 0 + for k in range(n): + if bkpt_extd[k] < bkpt[i]+bkpt[j] - w and bkpt[i]+bkpt[j] - w < bkpt_extd[k+1]: + if k != n-1: + K.gens()[k] < K.gens()[i] + K.gens()[j] - w + K.gens()[i] + K.gens()[j] - w < K.gens()[k+1] + else: + K.gens()[k] < K.gens()[i] + K.gens()[j] - w + K.gens()[i] + K.gens()[j] - w < 1 + elif bkpt_extd[k] == bkpt[i]+bkpt[j] - w: + if k != n-1: + K.gens()[k] == K.gens()[i] + K.gens()[j] - w + K.gens()[i] + K.gens()[j] - w < K.gens()[k+1] + else: + K.gens()[k] == K.gens()[i] + K.gens()[j] - w + K.gens()[i] + K.gens()[j] - w < 1 + if not log_paramateric_real_field: + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(parametric_logging_level) + return K._bsa + +def add_breakpoints_and_find_equiv_classes(bkpt_poly): + r""" + Takes dim k-1 breakpoint NNC polyhedron (as a :class:`BasicSemialgebraicSet_base`), adds a dimension, + and finds all possible representative elements for equivalence classes of polyhedral complexes. + + INPUT: :class:`BasicSemialgebraicSet_Polyhedral` + + OUTPUT: unique list of breakpoint sequnes of length k (as tuples). + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: add_breakpoints_and_find_equiv_classes(nnc_poly_from_bkpt_sequence([0,4/5]).upstairs()) + [(0, 9/14, 83/112), (0, 13/20, 33/40), (0, 9/14, 101/112), (0, 7/10, 17/20)] + """ + # BSAs are highly mutable, work only with copies. + B_cap_N_b = copy(bkpt_poly) + B_cap_N_b.add_space_dimensions_and_embed(1) + # get new number of breakpoints. + k = B_cap_N_b.ambient_dim() + if k < 1: + raise ValueError("bkpt_poly should have space dim at least 1.") + model_bound_bkpts = [0]*k + model_bound_bkpts[k-1] = 1 + # 0 < lambda_k < 1 + # we assume the bsa is "polyhedra type" + B_cap_N_b.add_linear_constraint(model_bound_bkpts, QQ(-1), operator.lt) # model bounds + B_cap_N_b.add_linear_constraint(model_bound_bkpts, QQ(0), operator.gt) # model bounds + bkpt_order = [0]*k + bkpt_order[k-2] = 1 + bkpt_order[k-1] = -1 + B_cap_N_b.add_linear_constraint(bkpt_order, QQ(0), operator.lt) # order on bkpts + rep_elems = [] + for j in range(k-1): + for i in range(k): + for interval_w in [0,1]: + for line_w in [0,1]: + # which interval is (lambda_k,lambda_k) located in? + # modeled lambda_i op 2lambda_k - w < lambda_{i+1} + for interval_op in [operator.lt, operator.eq, operator.gt]: + for line_op in [operator.lt, operator.eq, operator.gt]: + # highly mutable objects, operate on the copy. + B_cap_N_b_copy = copy(B_cap_N_b) + lhs_i = [0]*k + lhs_i[k-1] = -2 + lhs_i[i] = 1 + B_cap_N_b_copy.add_linear_constraint(lhs_i, QQ(interval_w), interval_op) + lhs_i_plus_1 = [0]*k + lhs_i_plus_1[k-1] = -2 + if i < k-1: + lhs_i_plus_1[i+1] = 1 + B_cap_N_b_copy.add_linear_constraint(lhs_i_plus_1, QQ(interval_w), operator.gt) + else: + B_cap_N_b_copy.add_linear_constraint(lhs_i_plus_1, QQ(interval_w + 1), operator.gt) + if not B_cap_N_b_copy.is_empty(): + # does the line x+y equiv lambda_k mod 1 lie on/above/below (lambda_j,lambda_j)? + # modeled by 2lambda_j op lambda_k + w + lhs_j = [0]*k + lhs_j[j] = 2 + lhs_j[k-1] = -1 + B_cap_N_b_copy.add_linear_constraint(lhs_j, QQ(-line_w), line_op) + try: + rep_elem = B_cap_N_b_copy.find_point() + rep_elems.append(tuple(rep_elem)) + except EmptyBSA: + pass + return unique_list(rep_elems) + +def make_rep_bkpts_with_len_n(n, k=1, bkpts=None, backend=None): + r""" + Produce representative elements of every isomorphism class of breakpoints complexes for breakpoint sequences of length n. + + INPUT: + - n, integer, maximum length of breakpoint sequence. + - k, assumed length of breakpoint sequences in ``bkpts``. + - bkpts, list of breakpoint sequences of length k. + + OUTPUT: A list of representative elements of every isomorphism class of breakpoints complexes for breakpoint sequences of length n extrapolated from bkpts. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: make_rep_bkpts_with_len_n(2) + [(0, 1/2), (0, 13/18), (0, 5/18)] + + The number of representative elements grows quickly:: + + sage: bkpts_rep_with_len_3 = make_rep_bkpts_with_len_n(3) + sage: len(bkpts_rep_with_len_3) + 34 + + Previous computations can be reused:: + + sage: bkpts_rep_with_len_4 = make_rep_bkpts_with_len_n(4, 3, bkpts_rep_with_len_3) + sage: len(bkpts_rep_with_len_4) + 329 + """ + # Matthias has suggested looking at a directed tree. + # An alternative approach would be to look into using a (graded) lattice as a data structure. + # We have bkpt \leq bkpt' if and only if dim(NNC(bkpt)) \leq dim(NNC(bkpt')) + # and embed(NNC(bkpt), dim(NNC(bkpt')) \leq NNC(bkpt') in the poset of NNC polyhedra. + # This might speed up/less the load of verifying unquiness of cells which is the time bounding task here. + new_bkpts = [] + if n < 2: + raise ValueError("n>=2") + if k == n and bkpts is not None: + minimal_funciton_cell_description_logger.warning(f"Initial inputs suggest that the bkpts provided are already correct. Returning the initial bkpts.") + return bkpts + if k == n and bkpts is None: + raise ValueError("k= 2) + assert(f_index >= 1) + assert(f_index <= n - 1) + if not isinstance(bkpt, list): + bkpt = list(bkpt) + coord_names = [] + val = [None]*(n) + for i in range(n): + coord_names.append('gamma'+str(i)) + K = ParametricRealField(names=coord_names, values=val, mutable_values=True, big_cells=True, allow_refinement=False, default_backend=backend) + K.gens()[0] == 0 + for i in range(1, n): + K.gens()[i] <= 1 + K.gens()[i] > 0 + h = piecewise_function_from_breakpoints_and_values(bkpt + [1], K.gens() + [0], merge=False) + # Assumes minimality for the partially defined function. + for vert in generate_type_1_vertices_continuous(h, operator.ge, bkpt + [1]): + vert + for vert in generate_type_2_vertices_continuous(h, operator.ge, bkpt + [1]): + vert + for vert in generate_assumed_symmetric_vertices_continuous(h, bkpt[f_index], bkpt + [1]): + vert + if not log_paramateric_real_field: + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(parametric_logging_level) + if not log_pw_functions: + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(pw_logging_level) + return K._bsa + +def value_nnc_polyhedron(bkpt, f_index, backend=None, log_paramateric_real_field = False,log_pw_functions = False): + r""" + For a given breakpoints sequence and f index, write the value polyhedron as a basic semialgebraic set in the full space of parameters. + + INPUTS: + - ``bkpt`` - sorted list of length 2 or more of sage type in the interval [0,1). + - ``f_index`` - integer between 1 and length of ``len(bkpt) -1``. + - ``backend`` - ``None``, ``str(pplite)`` + + OUTPUT: + - class::``BasicSemialgebraicSet_veronese`` - The value polyehdron in only the breakpoint and value parameter space. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: value_nnc_polyhedron([0,4/5], 1) # gmic with f=4/5, a trivial value polyhedron + BasicSemialgebraicSet_veronese(BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(Constraint_System {x3-1==0, x2==0, 5*x1-4==0, x0==0}, names=[x0, x1, x2, x3]), polynomial_map=[lambda0, lambda1, gamma0, gamma1]) + sage: P = value_nnc_polyhedron([0,2/15,2/3, 4/5], 3) # A non trivial value polyhedron. + sage: P.upstairs()._polyhedron + A 1-dimensional polyhedron in QQ^8 defined as the convex hull of 2 points + """ + n = len(bkpt) + if not log_paramateric_real_field: + parametric_logging_level = logging.getLogger("cutgeneratingfunctionology.igp.parametric").getEffectiveLevel() + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(logging.ERROR) + if not log_pw_functions: + pw_logging_level = logging.getLogger("cutgeneratingfunctionology.igp.functions").getEffectiveLevel() + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(logging.ERROR) + assert(n >= 2) + assert(f_index >= 1) + assert(f_index <= n) + coord_names = [] + bkpt_vals = list(bkpt) + vals = bkpt_vals + [None]*(n) + for i in range(n): + coord_names.append('lambda'+str(i)) + for i in range(n): + coord_names.append('gamma'+str(i)) + K = ParametricRealField(names=coord_names, values=vals, mutable_values=True, big_cells=True, allow_refinement=False, default_backend=backend) + # breakpoint parameters are the measured breakpoint values. + for i in range(n): + K.gens()[i] == bkpt[i] + # Necessary conditions on value parameters. + K.gens()[n] == 0 + for i in range(1, n): + K.gens()[i+n] <= 1 + K.gens()[i+n] > 0 + h = piecewise_function_from_breakpoints_and_values(K.gens()[0:n] + [1], K.gens()[n:2*n] + [0], merge=False) + # Assumes minimality for the partially defined function. + for vert in generate_type_1_vertices_continuous(h, operator.ge, K.gens()[0:n] + [1]): + vert + for vert in generate_type_2_vertices_continuous(h, operator.ge, K.gens()[0:n] + [1]): + vert + for vert in generate_assumed_symmetric_vertices_continuous(h, K.gens()[f_index], [0] + K.gens()[0:n] + [1]): + vert + if not log_paramateric_real_field: + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(parametric_logging_level) + if not log_pw_functions: + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(pw_logging_level) + return K._bsa + +def bsa_of_rep_element(bkpt, vals, backend=None, log_paramateric_real_field = False,log_pw_functions = False): + r""" + Given pi_(bkpt, vals) is {minimal, not minimal}, find BSA subset of R^(2n) such that (bkpt, vals) in BSA and for all p + in BSA, pi_p is {minimal, not minimal}. + + INPUT: + - ``bkpt`` - a breakpoint sequence + - ``vals`` - list like of sage numerical types corresponding values for the breakpoint sequence. + - ``backend`` - None, ``str(pplite)`` + + OUTPUT: A basic semialgebraic set. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: bsa_of_rep_element([0,4/5], [0,1]) # bsa for GMIC + BasicSemialgebraicSet_veronese(BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(Constraint_System {x5-1==0, x4==0, x0==0, 2*x3-1>0, -2*x1+x2-x3+1>=0, -x3+1>0}, names=[x0, x1, x2, x3, x4, x5]), polynomial_map=[gamma0, lambda1^2*gamma0, lambda1*gamma0, lambda1, lambda0, gamma1]) + """ + n = len(bkpt) + if not log_paramateric_real_field: + parametric_logging_level = logging.getLogger("cutgeneratingfunctionology.igp.parametric").getEffectiveLevel() + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(logging.ERROR) + if not log_pw_functions: + pw_logging_level = logging.getLogger("cutgeneratingfunctionology.igp.functions").getEffectiveLevel() + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(logging.ERROR) + assert(n >= 2) + coord_names = [] + for i in range(n): + coord_names.append('lambda'+str(i)) + for i in range(n): + coord_names.append('gamma'+str(i)) + K = ParametricRealField(names=coord_names, values=bkpt+vals, big_cells=True, default_backend=backend) + h = piecewise_function_from_breakpoints_and_values(K.gens()[0:n] + [1], K.gens()[n:2*n] + [0], merge=False) + minimality_test(h) + if not log_paramateric_real_field: + logging.getLogger("cutgeneratingfunctionology.igp.parametric").setLevel(parametric_logging_level) + if not log_pw_functions: + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(pw_logging_level) + return K.make_proof_cell().bsa + +def find_minimal_function_reps_from_bkpts(bkpts, prove_minimality=True, backend=None): + r""" + Finds representative elements of minimal functions from a given breakpoint sequence. + + INPUT: + - ``bkpts`` - an iterable of breakpoint sequences. + - ``prove_minimality`` - bool, proves minimality of parameterized function. + - ``backend`` - None, ``str(pplite)`` + + OUTPUT: + - List of (breakpoint, value) pairs. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: bkpts = make_rep_bkpts_with_len_n(2) + sage: find_minimal_function_reps_from_bkpts(bkpts) + [([0, 1/2], [0, 1]), ([0, 13/18], [0, 1]), ([0, 5/18], [0, 1])] + + """ + rep_elems = [] + if not log_pw_functions: + pw_logging_level = logging.getLogger("cutgeneratingfunctionology.igp.functions").getEffectiveLevel() + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(logging.ERROR) + for bkpt in bkpts: + n = len(bkpt) + for f_index in range(1, n): + poly_bsa = value_nnc_polyhedron_value_cords(list(bkpt), f_index, backend) + gammas = poly_bsa.polynomial_map()[0].parent().gens() + try: + test_point = poly_bsa.upstairs().find_point() + except EmptyBSA: + raise RepElemGenFailure("The value polyhedron {} is empty. This should not be empty. Double check inputs".format(poly_bsa)) + test_val = [] + for gamma_i in gammas: + test_val.append(test_point[poly_bsa.v_dict()[gamma_i]]) + if prove_minimality: + h = piecewise_function_from_breakpoints_and_values(list(bkpt)+[1], test_val+[0]) + if not minimality_test(h): # The following error should never be raised when this function is used as intended. + raise ValueError(f"({bkpt}, {test_val}) parameterized by breakpoints and values is not a minimal function but assuming a breakpoint sequence is input, this should be minimal.") + rep_elems.append((list(bkpt), test_val)) + if not log_pw_functions: + logging.getLogger("cutgeneratingfunctionology.igp.functions").setLevel(pw_logging_level) + return rep_elems + +class BreakpointComplexClassContainer: + """ + A container for the family of breakpoint complexes for piecewise linear functions + with at most n breakpoints. + + The container will attempt to load from the minimalFunctionCache module or will generate + data if the minimalFunctionCache is not available. + + Loading files generated from the container's `write_data()` method is supported. + When manually loading data, it is assume that all loaded data is correct + for the Initialized number of breakpoints. + + Initialization: + n - integer 2 or larger. + backend - (optional) `None` or `str(pplite)` + manually_load_breakpoint_cache - (optional) bool + loading_kwrds - file_or_folder, path_to_file_or_folder + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: bkpt_of_len_2 = BreakpointComplexClassContainer(2) + sage: bkpt_of_len_2.num_rep_elems() + 3 + sage: [elem for elem in bkpt_of_len_2.get_rep_elems()] + [(0, 1/2), (0, 13/18), (0, 5/18)] + sage: make_rep_bkpts_with_len_n(2) + [(0, 1/2), (0, 13/18), (0, 5/18)] + + A cell description of semialgebraic sets can be accessed:: + + sage: all([isinstance(cell, BasicSemialgebraicSet_base) for cell in bkpt_of_len_2.get_nnc_poly_from_bkpt()]) + True + + (Advanced use) Data can be written and reused:: + + sage: bkpt_of_len_2.write_data() # not tested + sage: bkpt_of_len_2_loaded = BreakpointComplexClassContainer(2, manually_load_breakpoint_cache=True, file_or_folder="file", path_to_file_or_folder="bkpts_of_len_2.csv") # not tested + """ + def __init__(self, n, backend=None, manually_load_breakpoint_cache=False, **loading_kwrds): + self._n = n + assert(self._n >= 2) + self._backend = backend + if not manually_load_breakpoint_cache: + try: + # Load the minimal function cache. + from minimalFunctionCache.utils import minimal_function_cache_info, minimal_function_cache_loader + cache_info = minimal_function_cache_info() + # minimal_function_cache_loader throws a value error if a cache for n is not found. + try: + self._data = minimal_function_cache_loader(self._n, "breakpoints") + except ValueError: + minimal_funciton_cell_description_logger.info(f"The cache for {n} breakpoints has not been generated or could not be found.") + minimal_funciton_cell_description_logger.info("Generating representative breakpoints.") + avail_bkpts = [ i for i in cache_info["avail_bkpts"] if i < self._n] + if len(avail_bkpts) > 0: + k = max(avail_bkpts) + prev_gen_bkpts = self._data = minimal_function_cache_loader(k, "breakpoints") + if self._n > 4: + minimal_funciton_cell_description_logger.warning(f"This may take a while.") + self._data = make_rep_bkpts_with_len_n(self._n, k, prev_gen_bkpts, backend=self._backend) + else: + if self._n > 4: + minimal_funciton_cell_description_logger.warning(f"This may take a while.") + self._data = make_rep_bkpts_with_len_n(self._n, backend=self._backend) + except ImportError: + if self._n > 4: + minimal_funciton_cell_description_logger.warning(f"This may take a while. Try installing the minimalfunctioncache.") + self._data = make_rep_bkpts_with_len_n(self._n, backend=self._backend) + minimal_funciton_cell_description_logger.info("Finished generating breakpoints.") + # For generating the cache and advanced use. + elif manually_load_breakpoint_cache: + try: + if loading_kwrds["file_or_folder"].strip(" ").lower() == "folder": + loaded_data = [] + path = loading_kwrds["path_to_file_or_folder"] + import os + files_in_folder = os.listdir(path) + for file in files_in_folder: + with open(os.path.join(path, file)) as csvfile: + file_reader = csv.reader(csvfile) + for row in file_reader: + loaded_data.append([QQ(data) for data in row]) + k = len(loaded_data[0]) # assume all the data is correct and of the same length. + self._data = make_rep_bkpts_with_len_n(self._n, k, loaded_data, backend=self._backend) + elif loading_kwrds["file_or_folder"].strip(" ").lower() == "file": + loaded_data = [] + with open(loading_kwrds["path_to_file_or_folder"]) as csvfile: + file_reader = csv.reader(csvfile) + for row in file_reader: + loaded_data.append([QQ(data) for data in row]) + k = len(loaded_data[0]) # assume all the data is correct and of the same length. + self._data = make_rep_bkpts_with_len_n(self._n, k, loaded_data, backend=self._backend) + else: + raise ValueError("Check spelling of folder or file.") + except KeyError: + raise ValueError("No valid inputs provided. Maybe check your spelling?") + else: + raise ValueError("No elements have been loaded. Check inputs") + + def __repr__(self): + return f"Container for the space breakpoint sequences of length {self._n} under equivalence of polyhedral complexes." + + def get_rep_elems(self): + """ + Iterator yielding a breakpint sequence. + """ + for bkpt in self._data: + yield bkpt + + def get_nnc_poly_from_bkpt(self): + """ + Iterator yielding a cell such that all elements in the cell correspond to 2d polyhedral complexes of the same combinatorial type. + """ + for bkpt in self._data: + yield nnc_poly_from_bkpt_sequence(bkpt, backend=self._backend) + + def num_rep_elems(self): + """ + Number of repersentative elements. + """ + return len(self._data) + + def n(self): + """ + The length of breakpoint sequences. + """ + return self._n + + def covers_space(self): + """ + Not Implemented. Future use is intended to be a proof of correctness that all breakpoint sequences are covered. + """ ### TODO: Impl. + raise NotImplementedError + + def refine_space(self): + """ + Not Implemented. Future use is to prove that the current cell description for breakpoints is, with respect to inclusion, minimal or in the case that it is not minimal, find a minimal description. + """ + raise NotImplementedError + ### TODO: Add a contaimenet checking method to BSA classes and finish impl. + self._data = unique_list(self._data) + cells_found = [] + for bkpt in self._data: + bkpt_contained_in_cell = False + for found_cell in cells_found: + if found_cell.contains(bkpt): + bkpt_contained_in_cell = True + break + if not contained_in_cell: + cells_found.append(nnc_poly_from_bkpt_sequence(bkpt, backend=self._backend)) + new_data = [] + for cell in cells_found: + new_data.append(cell.find_point()) + self._data = new_data + + def write_data(self, output_file_name_style=None, max_rows=None): + """ + Writes representative element data to a `.csv` file with one column and rows of representative elements. + Optionally, write many `.csv` files with at `most max_rows` rows per file. + Files are named output_file_file_name_style_filenumber.csv. + The default output_file_name_style="bkpts_of_len_n". + """ + # TODO: Future, support writing different types of data such as polyhedra data. + if output_file_name_style is None: + file_name_base = "bkpts_of_len_{}".format(self._n) + else: + file_name_base = output_file_name_style + if max_rows is not None: + assert(max_rows >= 1) + num_files = len(self._data)//max_rows + 1 + file_name_base = file_name_base+"_part_0" + if max_rows is None: + max_rows = len(self._data) + num_files = 1 + output_file = file_name_base+".csv" + for file_number in range(num_files): + out_file = open(output_file, "w") + data_writer = csv.writer(out_file, csv.QUOTE_NONE) + for row in range(max_rows): + try: + data_writer.writerow(self._data[max_rows * file_number + row]) + except IndexError: + break + out_file.close() + output_file = file_name_base[:-1]+"{}".format(file_number+1)+".csv" + + +class PiMinContContainer: + """ + A container for the space of continuous piecewise linear minimal functions with at + most n breakpoints parameterized by breakpoints and values using semialgebraic sets. + + The container will attempt to load from the `minimalFunctionCache module` or will generate + data if the minimalFunctionCache is not available. + + Loading files generated from the container's `write_data()` method is supported. + When manually loading data, it is assume that all loaded data is correct + for the Initialized number of breakpoints. + + INPUTS: + - n, an integer + - backend, None or str(pplite) + - loading_kwrds + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.igp import * + sage: logging.disable(logging.INFO) # suppress logging for tests + sage: PiMin_2 = PiMinContContainer(2) + sage: all([minimality_test(pi) for pi in PiMin_2.get_rep_functions()]) + True + sage: PiMin_2 + Space of minimal functions with at most 2 breakpoints parameterized by breakpoints and values using semialgebraic sets. + + A cell description of semialgebraic sets can be accessed:: + + sage: all([isinstance(cell, BasicSemialgebraicSet_base) for cell in PiMin_2.get_semialgebraic_sets()]) + True + + Data is stored as repersenative elements. The number of repersentative elements grows quickly. :: + + sage: PiMin_4 = PiMinContContainer(4) # not tested + sage: len([rep_elem for rep_elem in PiMin_4.get_rep_elems()]) # not tested + 987 + sage: len([rep_elem for rep_elem in PiMin_2.get_rep_elems()]) + 3 + + The container provides methods of writing data.:: + + sage: PiMin_2.write_data() # not tested + + Written data can be reused.:: + + sage: PiMin_2_loaded_data = PiMinContContainer(2, manually_load_function_cache=True, file_or_folder="file", path_to_file_or_folder="Pi_Min_2.csv", breakpoints_or_rep_elems="rep_elems") # not tested + sage: len([rep_elem for rep_elem in PiMin_2_loaded_data.get_rep_elems()]) # not tested + 3 + """ + def __init__(self, n, backend=None, manually_load_function_cache=False, **loading_kwrds): + self._n = n + assert(self._n >= 2) + self._backend = backend + if not manually_load_function_cache: + try: + # Load the minimal function cache. + from minimalFunctionCache.utils import minimal_function_cache_loader + try: + self._data = minimal_function_cache_loader(self._n, "rep_elems") + # cache loader throws a value error if a cache for n is not found. + except ValueError: + minimal_funciton_cell_description_logger.info(f"The cache for {n} breakpoints has not been generated or could not be found.") + minimal_funciton_cell_description_logger.info("Finding or generating representative breakpoints.") + bkpts = BreakpointComplexClassContainer(self._n, backend=self._backend).get_rep_elems() + if self._n > 4: + minimal_funciton_cell_description_logger.warning(f"This may take a while.") + minimal_funciton_cell_description_logger.info("Finished generating breakpoints.") + minimal_funciton_cell_description_logger.info("Computing repersentative elements.") + self._data = find_minimal_function_reps_from_bkpts(bkpts, backend=self._backend) + except ImportError: + if self._n > 4: + minimal_funciton_cell_description_logger.warning(f"This may take a while. Try installing the minimalFunctionCache.") + bkpts = make_rep_bkpts_with_len_n(self._n, backend=self._backend) + minimal_funciton_cell_description_logger.info("Finished generating breakpoints.") + minimal_funciton_cell_description_logger.info("Computing repersentative elements.") + self._data = find_minimal_function_reps_from_bkpts(bkpts, backend=self._backend) + minimal_funciton_cell_description_logger.info("PiMin container, Reportin' for duty.") + elif manually_load_function_cache: + # this is for generating the cache and advanced use. + minimal_funciton_cell_description_logger.info("loading files...") + try: + if loading_kwrds["breakpoints_or_rep_elems"].strip(" ").lower() == "breakpoints": + bkpts = BreakpointComplexClassContainer(self._n, backend=self._backend, manually_load_breakpoint_cache=True, file_or_folder=loading_kwrds["file_or_folder"], path_to_file_or_folder=loading_kwrds["path_to_file_or_folder"]).get_rep_elems() + if self._n > 4: + minimal_funciton_cell_description_logger.warning("This may take a while. Try installing the minimalFunctionCache.") + self._data = find_minimal_function_reps_from_bkpts(bkpts, backend=self._backend) + minimal_funciton_cell_description_logger.info("PiMin container, Reportin' for duty.") + elif loading_kwrds["breakpoints_or_rep_elems"].strip(" ").lower() == "rep_elems": + if loading_kwrds["file_or_folder"].strip(" ").lower() == "folder": + self._data = [] + path = loading_kwrds["path_to_file_or_folder"] + import os + files_in_folder = os.listdir(path) + for file in files_in_folder: + with open(os.path.join(path, file)) as csvfile: + file_reader = csv.reader(csvfile) + for row in file_reader: + bkpt = [QQ(data) for data in row[0].strip("[]").split(",")] + val = [QQ(data) for data in row[1].strip("[]").split(",")] + self._data.append((bkpt, val)) + minimal_funciton_cell_description_logger.info("PiMin container, Reportin' for duty.") + elif loading_kwrds["file_or_folder"].strip(" ").lower() == "file": + self._data = [] + with open(loading_kwrds["path_to_file_or_folder"]) as csvfile: + file_reader = csv.reader(csvfile) + for row in file_reader: + bkpt = [QQ(data) for data in row[0].strip("[]").split(",")] + val = [QQ(data) for data in row[1].strip("[]").split(",")] + self._data.append((bkpt, val)) + minimal_funciton_cell_description_logger.info("PiMin container, Reportin' for duty.") + else: + raise ValueError("check spelling of folder or file.") + except KeyError: + raise ValueError("No valid inputs provided. Maybe check your spelling?") + else: + raise ValueError("No elements have been loaded. Check inputs") + + def __repr__(self): + return "Space of minimal functions with at most {} breakpoints parameterized by breakpoints and values using semialgebraic sets.".format(self._n) + + def get_semialgebraic_sets(self): + """Iterator for semialgebraic set description.""" + for b, v in self._data: + yield bsa_of_rep_element(list(b), list(v), backend=self._backend) + + def get_rep_elems(self): + """Iterator for representative elements.""" + for b, v in self._data: + yield (list(b), list(v)) + + def get_rep_functions(self): + """Iterator for representative functions.""" + for b, v in self._data: + yield piecewise_function_from_breakpoints_and_values(list(b)+[1], list(v)+[0]) + + def n(self): + """The maximum number of proper breakpoints of parameterized functions in the space.""" + return self._n + + def covers_space(self): + """ + Not Implemented. Future use is intended to be a proof of correctness that all breakpoint sequences are covered. + """ + ### TODO: Impl. + raise NotImplementedError + + def refine_space(self): + """ + Not Implemented. Future use is to prove that the current cell description for breakpoints is, with respect to inclusion, minimal or in the case that it is not minimal, find a minimal description. + """ + ### TODO: Impl. + raise NotImplementedError + + def write_data(self, output_file_name_style=None, max_rows=None): + """ + Writes representative element data to a `.csv` file with one column and rows of representative elements. + Optionally, write many `.csv` files with at `most max_rows` rows per file. + Files are named output_file_file_name_style_filenumber.csv. + The default output_file_name_style=f"Pi_Min_{n}". + """ + # TODO: Future, support writing different types of data such as polyhedra data. + if output_file_name_style is None: + file_name_base = "Pi_Min_{}".format(self._n) + else: + file_name_base = output_file_name_style + if max_rows is not None: + assert(max_rows >= 1) + num_files = len(self._data)//max_rows + 1 + file_name_base = file_name_base + "_part_0" + if max_rows is None: + max_rows = len(self._data) + num_files = 1 + output_file = file_name_base + ".csv" + for file_number in range(num_files): + out_file = open(output_file, "w") + data_writer = csv.writer(out_file, csv.QUOTE_NONE) + for row in range(max_rows): + try: + data_writer.writerow(self._data[max_rows * file_number + row]) + except IndexError: + break + out_file.close() + output_file = file_name_base[:-1]+"{}".format(file_number+1)+".csv" + + +def plot_2_d_polyhedral_complex(bkpt, highlight_index=-1, highlight_color='blue', **kwds): + r""" + Returns a plot of 2-d polyehdral complex of the breakpoint sequence of length n, Complex Delta P_bkpt. + Lines generated from a partiuclar breakpoint can be highlighted. + Assumes that the breakpoint sequence is correct and ``highlight_index`` is between 0 and ``len(bkpt)-1``, inclusive. + + EXAMPLES:: + sage: from cutgeneratingfunctiology.igp import * + sage: plot_2_d_polyhedral_complex([0, 4/5]) # not tested + + A breakpoint can be highlighted:: + sage: plot_2_d_polyhedral_complex([0, 1/2, 4/5], highlight_index=1) # not tested + + An the color changed:: + sage: plot_2_d_polyhedral_complex([0, 4/5], highlight_index=1, highlight_color='green') # not tested + """ + bkpt_ext = bkpt + [1] + x = var('x') + p = Graphics() + n = len(bkpt) + if highlight_index is None: + highlight_index = -1 + for i in range(1,n+1): + if i == highlight_index: + p += line([(0, bkpt_ext[i]), (bkpt_ext[i], 0)], color=highlight_color, linestyle='dashed') + else: + p += line([(0, bkpt_ext[i]), (bkpt_ext[i], 0)], color='grey') + for i in range(1,n): + if i == highlight_index: + p += line([(bkpt_ext[i], 1), (1, bkpt_ext[i])], color=highlight_color, linestyle='dashed') + else: + p += line([(bkpt_ext[i], 1), (1, bkpt_ext[i])], color='grey') + for i in range(n+1): + if i == highlight_index: + p += plot(bkpt_ext[i], (0, 1), color=highlight_color, linestyle='dashed') + else: + p += plot(bkpt_ext[i], (0, 1), color='grey') + y = var('y') + for i in range(n): + if i == highlight_index: + p += parametric_plot((bkpt_ext[i],y), (y,0,1), color=highlight_color, linestyle='dashed') + else: + p += parametric_plot((bkpt_ext[i],y), (y,0,1), color='grey') + p += line([(1,0), (1,1)], color='grey') + return p + +def plot_2_d_polyhedral_complex_and_descendants(bkpt, max_row_length=5, max_number_additional_diagrams=15, **kwds): + r""" + Returns a plot of 2-d polyehdral complex of the breakpoint sequence of length n, Complex Delta P_bkpt + together with plots of Complex :math:`Delta P_{bkpt \cup \lambda^*}` where Complex Delta P_bkpt is a + subcomplex of :math:`Delta P_{bkpt \cup \lambda^*}`. + + EXAMPLES:: + sage: from cutgeneratingfunctiology.igp import * + sage: plot_2_d_polyhedral_complex_and_descendants([0, 1/2]) # not tested + """ + n = len(bkpt) + n_plus_one_breakpoints = unique_list(make_rep_bkpts_with_len_n(n+1, k=n, bkpts=[bkpt])) + number_of_new_diagrams = len(n_plus_one_breakpoints) + m = max_row_length + i = 0 + plots = [] + while i < max_number_additional_diagrams and i < number_of_new_diagrams: + plots.append((plot_2_d_polyhedral_complex(list(n_plus_one_breakpoints[i]), n, **kwds), (i % m , int(i/m) ,.75,.75))) + i += 1 + if i < m: + plots.append((plot_2_d_polyhedral_complex(bkpt), (int(i/2), -1,.75,.75))) + else: + plots.append((plot_2_d_polyhedral_complex(bkpt), (int(m/2), -1 ,.75,.75))) + return multi_graphics(plots) + diff --git a/cutgeneratingfunctionology/igp/parametric.sage b/cutgeneratingfunctionology/igp/parametric.sage index 26bf3eb39..a8f7b3b47 100644 --- a/cutgeneratingfunctionology/igp/parametric.sage +++ b/cutgeneratingfunctionology/igp/parametric.sage @@ -20,6 +20,7 @@ from cutgeneratingfunctionology.spam.basic_semialgebraic_local import BasicSemia from cutgeneratingfunctionology.spam.semialgebraic_mathematica import BasicSemialgebraicSet_mathematica, from_mathematica from cutgeneratingfunctionology.spam.basic_semialgebraic_groebner_basis import BasicSemialgebraicSet_groebner_basis from cutgeneratingfunctionology.spam.polyhedral_complex import PolyhedralComplex +from cutgeneratingfunctionology.spam.EvaluationExceptions import FactorUndetermined from .parametric_family import Classcall, ParametricFamily_base, ParametricFamily import logging @@ -40,12 +41,12 @@ 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 @@ -54,16 +55,16 @@ 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 +from contextlib import contextmanager allow_refinement_default = True big_cells_default = 'if_not_allow_refinement' @@ -162,32 +163,29 @@ 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. = 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. = 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. = 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 def __init__(self, values=None, names=None, allow_coercion_to_float=True, mutable_values=None, allow_refinement=None, big_cells=None, - base_ring=None, sym_ring=None, bsa=None): + base_ring=None, sym_ring=None, bsa=None, default_backend=None): Field.__init__(self, self) if mutable_values is None: @@ -214,7 +212,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([]) @@ -257,7 +255,10 @@ class ParametricRealField(Field): # do the computation of the polyhedron incrementally, # rather than first building a huge list and then in a second step processing it. # the upstairs polyhedron defined by all constraints in self._eq/lt_factor - polyhedron = BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(0) + if default_backend == "pplite": + polyhedron = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(0) + else: + polyhedron = BasicSemialgebraicSet_polyhedral_ppl_NNC_Polyhedron(0) # monomial_list records the monomials that appear in self._eq/lt_factor. # v_dict is a dictionary that maps each monomial to the index of its corresponding Variable in polyhedron bsa = BasicSemialgebraicSet_veronese(polyhedron, polynomial_map=[], poly_ring=sym_ring, v_dict={}) @@ -316,19 +317,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. = 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 @@ -414,9 +415,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 ] @@ -567,6 +566,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. @@ -858,6 +866,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 symobolic 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: @@ -888,7 +905,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 @@ -1112,7 +1130,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:: @@ -1152,7 +1170,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. @@ -1214,8 +1232,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. @@ -1243,7 +1263,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 @@ -1302,7 +1322,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. = ParametricRealField([1,1/2]) sage: assert(x == 2*y) @@ -1336,7 +1356,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)): @@ -1404,7 +1428,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): @@ -1416,13 +1440,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) @@ -1451,7 +1475,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) @@ -1538,7 +1562,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: @@ -1666,7 +1692,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.) """ @@ -1776,8 +1802,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. @@ -1811,7 +1839,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 @@ -1906,7 +1934,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'``, @@ -1985,7 +2013,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. @@ -2011,11 +2039,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: @@ -2076,12 +2104,13 @@ 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 * @@ -2110,7 +2139,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 * @@ -2187,18 +2216,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: parametric_logger.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: @@ -2259,7 +2288,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:: @@ -2296,7 +2325,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'] @@ -2315,7 +2344,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. @@ -2337,11 +2366,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())) @@ -2351,9 +2380,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 @@ -2390,7 +2419,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. @@ -2549,9 +2578,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 @@ -2684,8 +2713,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. @@ -2707,8 +2736,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. @@ -2732,7 +2761,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 @@ -2757,7 +2786,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) @@ -2803,7 +2832,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 @@ -2837,7 +2866,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""" @@ -2882,11 +2912,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. @@ -2972,7 +3002,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.=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) @@ -2992,7 +3022,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: @@ -3016,8 +3046,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 @@ -3039,7 +3069,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.=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 @@ -3058,8 +3088,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.=QQ[]; - sage: current_point = vector((71582788/143165577, 4673/377000)) # came from vector((RR(70727/150800), RR(4673/377000))), + sage: P.=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 @@ -3074,7 +3104,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. @@ -3136,7 +3166,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. @@ -3197,13 +3227,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 @@ -3246,6 +3277,7 @@ def embed_function_into_family(given_function, parametric_family, check_completi # plot_cpl_components(complex.components) return {} + """ EXAMPLES:: diff --git a/cutgeneratingfunctionology/spam/EvaluationExceptions.py b/cutgeneratingfunctionology/spam/EvaluationExceptions.py new file mode 100644 index 000000000..7e7dd6e04 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/cutgeneratingfunctionology/spam/basic_semialgebraic.py b/cutgeneratingfunctionology/spam/basic_semialgebraic.py index 796aec7aa..4ca1471a5 100644 --- a/cutgeneratingfunctionology/spam/basic_semialgebraic.py +++ b/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: diff --git a/cutgeneratingfunctionology/spam/basic_semialgebraic_pplite.py b/cutgeneratingfunctionology/spam/basic_semialgebraic_pplite.py new file mode 100644 index 000000000..559be7403 --- /dev/null +++ b/cutgeneratingfunctionology/spam/basic_semialgebraic_pplite.py @@ -0,0 +1,367 @@ +from cutgeneratingfunctionology.spam.basic_semialgebraic import * + +from pplite import Variable as pplite_Var, Constraint as pplite_Con, Linear_Expression as pplite_Lin_Expr, Affine_Expression as pplite_Aff_expr, NNC_Polyhedron as pplite_NNC_Polyhedron, PPliteGenerator, Polyhedron_Constraint_Rel, Polyhedron_Generator_Rel + +poly_is_included_pplite = Polyhedron_Constraint_Rel.is_included() +point_is_included_pplite = Polyhedron_Generator_Rel.subsumes() + +class BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(BasicSemialgebraicSet_polyhedral): + + r""" + A (possibly half-open) polyhedral basic semialgebraic set, + represented by a PPLite ``NNC_Polyhedron`` + + """ + + def __init__(self, ambient_dim=None, polyhedron=None, base_ring=None, poly_ring=None, **options): + r""" + Initialize a basic semialgebraic set as the universe in + ``ambient_dim``, or, if ``polyhedron`` (an ``NNC_Polyhedron``, + which after that belongs to this object) is provided, as + that. + + TEST:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: P = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(2) + sage: P.add_linear_constraint([0,1],0,operator.ge) + sage: P.add_linear_constraint([1,0],0,operator.ge) + sage: P.add_linear_constraint([2,3],-6,operator.lt) + sage: P + BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron([x1>=0, x0>=0, -2*x0-3*x1+6>0], names=[x0, x1]) + sage: sorted(P.eq_poly()) + [] + sage: sorted(P.lt_poly()) + [2*x0 + 3*x1 - 6] + sage: sorted(P.le_poly()) + [-x1, -x0] + """ + if ambient_dim is None and polyhedron is not None: + ambient_dim = polyhedron.space_dimension() + if base_ring is None and poly_ring is None: + base_ring = QQ + poly_ring, base_ring, ambient_dim, names = self._poly_ring_from_options( + ambient_dim=ambient_dim, base_ring=base_ring, poly_ring=poly_ring, **options) + if base_ring is not QQ: + raise ValueError("only base_ring=QQ is supported") + super(BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron, self).__init__(poly_ring=poly_ring) + if polyhedron is None: + self._polyhedron = pplite_NNC_Polyhedron(dim_type=int(ambient_dim), spec_elem='universe', topology="nnc") # To work with pplite, ambient_dim is required to be type int + else: + self._polyhedron = polyhedron + + @staticmethod + def _pplite_constraint(lhs, cst, op): + r""" + Make a PPL ``Constraint`` ``lhs`` * x + cst ``op`` 0, + where ``lhs`` is be a vector of length ambient_dim. + """ + lcd = lcm(lcm(x.denominator() for x in lhs), cst.denominator()) + lin_expr = sum([int((lcd*lhs[i]))*pplite_Var(i) for i in range(len(lhs))]) + aff_expr = pplite_Aff_expr(lin_expr, int(lcd*cst)) + #linexpr = pplite_Lin_Expr(lhs * lcd, cst * lcd) + if op == operator.lt: + return (aff_expr < 0) + elif op == operator.gt: + return (aff_expr > 0) + elif op == operator.eq: + return (aff_expr == 0) + elif op == operator.le: + return (aff_expr <= 0) + elif op == operator.ge: + return (aff_expr >= 0) + else: + raise ValueError("{} is not a supported operator".format(op)) + + def __copy__(self): + r""" + Make a copy of ``self``. + + TESTS: + + Test that it is actually making a copy of the (mutable!) NNC_Polyhedron:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: P = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(2) + sage: P._polyhedron is copy(P)._polyhedron + False + """ + return self.__class__(polyhedron=pplite_NNC_Polyhedron(nnc_poly=self._polyhedron), poly_ring=self.poly_ring()) + + def _repr_(self): + constraints = self._polyhedron.constraints() + names = list(self.poly_ring().gens()) + return 'BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron({}, names={})'.format( + constraints, names) + + def closure(self, bsa_class='formal_closure'): + r""" + Return the basic semialgebraic set that is the topological closure + of ``self``. + """ + # Because our description consists of minimized constraints, the closure is + # just the formal closure. + return self.formal_closure(bsa_class=bsa_class) + + def relint(self, bsa_class='formal_relint'): + r""" + Return the basic semialgebraic set that is the topological relative interior + of ``self``. + """ + # Because our description consists of minimized constraints, the relint is + # just the formal relint. + return self.formal_relint(bsa_class=bsa_class) + + def eq_poly(self): + r""" + Return a list of the polynomials `f` in equations `f(x) = 0` + in the description of ``self``. + + Together, ``eq_poly``, ``lt_poly``, and ``le_poly`` describe ``self``. + """ + # add tests + for c in self._polyhedron.constraints(): + if c.is_equality(): + coeff = [c.coefficient(pplite_Var(i)) for i in range(c.space_dimension())] + # observe: coeffients in a constraint of NNC_Polyhedron could have gcd != 1. + gcd_c = gcd(gcd(coeff), c.inhomogeneous_term()) + t = sum(QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) + QQ(c.inhomogeneous_term())/gcd_c # not type stable, make it type stable + yield self.poly_ring()(t) + + def lt_poly(self): + r""" + Return a list of the polynomials `f` in strict inequalities `f(x) < 0` + in the description of ``self``. + + Together, ``eq_poly``, ``lt_poly``, and ``le_poly`` describe ``self``. + """ + for c in self._polyhedron.constraints(): + if c.is_strict_inequality(): + coeff = [c.coefficient(pplite_Var(i)) for i in range(c.space_dimension())] + gcd_c = gcd(gcd(coeff), c.inhomogeneous_term()) + # constraint is written with '>', while lt_poly records '<' relation + t = sum(-QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) - QQ(c.inhomogeneous_term())/gcd_c + yield self.poly_ring()(t) + + def le_poly(self): + r""" + Return a list of the polynomials `f` in inequalities `f(x) \leq 0` + in the description of ``self``. + + Together, ``eq_poly``, ``lt_poly``, and ``le_poly`` describe ``self``. + """ + for c in self._polyhedron.constraints(): + if c.is_nonstrict_inequality(): + coeff = [c.coefficient(pplite_Var(i)) for i in range(c.space_dimension())] + gcd_c = gcd(gcd(coeff), c.inhomogeneous_term()) + # constraint is written with '>=', while lt_poly records '<=' relation + t = sum(-QQ(x)/gcd_c*y for x, y in zip(coeff, self.poly_ring().gens())) - QQ(c.inhomogeneous_term())/gcd_c + yield self.poly_ring()(t) + + # override the default implementation + def __contains__(self, point): + r""" + Whether the set contains the ``point`` (vector). + """ + rational_list = [ QQ(x) for x in point ] + num_list = [x.numerator() for x in rational_list] + den_list = [x.denominator() for x in rational_list] + common_den = lcm(den_list) + coef = [common_den // den_list[i] * num_list[i] for i in range(len(rational_list))] + pt = ppl_point(Linear_Expression(coef, 0), common_den) + return self._polyhedron.relation_with(pt).implies(point_is_included_pplite) + + # override the abstract methods + def find_point(self): + r""" + Find a point in ``self``. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: P = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(2) + sage: P.add_linear_constraint([0,1],0,operator.ge) + sage: P.add_linear_constraint([1,0],0,operator.ge) + sage: P.add_linear_constraint([2,3],-6,operator.lt) + sage: P.find_point() + (1, 2/3) + """ + # pplite has a different representation points, closure points, of NNC polys compared to ppl + # so the find_point method yields different results + + def to_point(g, ambient_dim): + den = g.divisor() + # g.set_space_dimension(ambient_dim) # PPlite generators have space dim of largest dimension of variables in point expression. + # To sum points as vectors in sagemath vectors need to have the same dimension. + # To fix, update the points space dim to be a defined ambient dimension. + # TODO: update after this gets fixed in pplite. + return vector(QQ, (QQ(x)/den for x in [g.coefficient(v) for v in range(ambient_dim)])) # based on email this should in theory works + + def to_vector(g, ambient_dim): + den = g.divisor() + # g.set_space_dimension(ambient_dim) + return vector(QQ, (QQ(x)/den for x in [g.coefficient(v) for v in range(ambient_dim)])) + points = [to_point(g, self._polyhedron.space_dimension()) for g in self._polyhedron.generators() + if g.is_point() or g.is_closure_point()] + rays = [to_vector(g, self._polyhedron.space_dimension()) 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: + p += sum(rays) / len(rays) + return p + raise NotImplementedError("find_test_point implementation cannot handle this case") + + def add_space_dimensions_and_embed(self, space_dim_to_add): + r""" + Mutate ``self`` by injecting it into a higher dimensional space. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: P = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(2) + sage: P.add_linear_constraint([0,1],0,operator.ge) + sage: P.add_linear_constraint([1,0],0,operator.ge) + sage: P.add_linear_constraint([2,3],-6,operator.lt) + sage: P.add_space_dimensions_and_embed(2) + sage: P + BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron([x1>=0, x0>=0, -2*x0-3*x1+6>0], names=[x0, x1, x2, x3]) + sage: P.ambient_dim() + 4 + """ + super(BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron, self).add_space_dimensions_and_embed(space_dim_to_add) + self._polyhedron.add_space_dimensions(int(space_dim_to_add), False) + + @staticmethod + def _pplite_constraint(lhs, cst, op): + r""" + Make a PPLite ``Constraint`` ``lhs`` * x + cst ``op`` 0, + where ``lhs`` is be a vector of length ambient_dim. + + TESTS:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: from pplite import Constraint as pplite_Con + sage: test_constraint = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron._pplite_constraint([0,1], 0, operator.ge) + sage: test_constraint + x1>=0 + sage: isinstance(test_constraint, pplite_Con) + True + + """ + lcd = lcm(lcm(x.denominator() for x in lhs), cst.denominator()) + aff_expr = sum([int((lcd*lhs[i]))*pplite_Var(i) for i in range(len(lhs))]) + int(lcd*cst) + if op == operator.lt: + return (aff_expr < 0) + elif op == operator.gt: + return (aff_expr > 0) + elif op == operator.eq: + return (aff_expr == 0) + elif op == operator.le: + return (aff_expr <= 0) + elif op == operator.ge: + return (aff_expr >= 0) + else: + raise ValueError("{} is not a supported operator".format(op)) + + def linear_function_upper_bound(self, form): + r""" + Find an upper bound for ``form`` (a vector) on ``self``. + This upper bound is the supremum. + + If ``self`` is empty, it returns -oo + """ + + def to_point(g): + den = g.divisor() + return vector(QQ, (QQ(x)/den for x in [g.coefficient(v) for v in range(g.space_dimision())])) + + def to_vector(g): + return vector(QQ, (QQ(x) for x in [g.coefficient(v) for v in range(g.space_dimision())])) + if self._polyhedron.is_empty(): + return -Infinity + form = vector(form) + for g in self._polyhedron.generators(): + if g.is_line(): + if to_vector(g) * form != 0: + return +Infinity + if g.is_ray(): + if to_vector(g) * form > 0: + return +Infinity + points = [to_point(g) for g in self._polyhedron.generators() + if g.is_point() or g.is_closure_point()] + return max(p * form for p in points) + + def is_linear_constraint_valid(self, lhs, cst, op): + r""" + Whether the constraint ``lhs`` * x + cst ``op`` 0 + is satisfied for all points of ``self``, + where ``lhs`` is be a vector of length ambient_dim. + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: P = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(2) + sage: P.add_linear_constraint([0,1],0,operator.ge) + sage: P.add_linear_constraint([1,0],0,operator.ge) + sage: P.add_linear_constraint([2,3],-6,operator.lt) + sage: P.is_linear_constraint_valid([1,1],-3,operator.lt) + True + sage: P.is_linear_constraint_valid([0,1],0,operator.gt) + False + """ + lhs = vector(lhs) + constraint = self._pplite_constraint(lhs, cst, op) + return self._polyhedron.relation_with(constraint).implies(poly_is_included_pplite) + + def add_linear_constraint(self, lhs, cst, op): + r""" + Add the constraint ``lhs`` * x + cst ``op`` 0, + where ``lhs`` is a vector of length ambient_dim, and + ``op`` is one of ``operator.lt``, ``operator.gt``, ``operator.eq``, + ``operator.le``, ``operator.ge`` + + EXAMPLES:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: P = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(2) + sage: P.add_linear_constraint([2,3],-6,operator.gt) + sage: sorted(P.lt_poly()) + [-2*x0 - 3*x1 + 6] + """ + lhs = vector(lhs) + constraint = self._pplite_constraint(lhs, cst, op) + self._polyhedron.add_constraint(constraint) + + def is_empty(self): + """ + EXAMPLES:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: S = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(1) + sage: S.add_linear_constraint([1], -1, operator.ge) + sage: S.is_empty() + False + sage: S.add_linear_constraint([1], +1, operator.le) + sage: S.is_empty() + True + """ + return self._polyhedron.is_empty() + + def is_universe(self): + """ + EXAMPLES:: + + sage: from cutgeneratingfunctionology.spam.basic_semialgebraic import * + sage: S = BasicSemialgebraicSet_polyhedral_pplite_NNC_Polyhedron(1) + sage: S.add_linear_constraint([0], 0, operator.eq) + sage: S.is_universe() + True + sage: S.add_linear_constraint([1], 1, operator.le) + sage: S.is_universe() + False + """ + self._polyhedron.minimize() + return self._polyhedron.is_universe() diff --git a/cutgeneratingfunctionology/spam/parametric_real_field_element.py b/cutgeneratingfunctionology/spam/parametric_real_field_element.py index c6688eeeb..c6ba3bfa7 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.EvaluationExceptions 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) - else: - left.parent().assume_comparison(right, operator.lt, left) - return richcmp(left.val(), right.val(), op) + 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: + 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. diff --git a/requirements.txt b/requirements.txt index a54e3a07a..aba99df62 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,3 +4,4 @@ sphinx sphinxcontrib-bibtex sphinxcontrib-websupport pynormaliz +pplitepy diff --git a/setup.py b/setup.py index a0713411a..cfc589fbd 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,14 @@ +## -*- encoding: utf-8 -*- +import os +import sys from setuptools import setup +from codecs import open # To open the README file with proper encoding + + +# Get information from separate files (README, VERSION) +def readfile(filename): + with open(filename, encoding='utf-8') as f: + return f.read() setup( packages = ['cutgeneratingfunctionology', 'cutgeneratingfunctionology.igp', 'cutgeneratingfunctionology.multirow', 'cutgeneratingfunctionology.dff', 'cutgeneratingfunctionology.spam', 'cutgeneratingfunctionology.igp.subadditivity_slack_diagrams', 'cutgeneratingfunctionology.igp.procedures'],