diff --git a/PEPit/__init__.py b/PEPit/__init__.py index e1ce342b..c22ad4a4 100644 --- a/PEPit/__init__.py +++ b/PEPit/__init__.py @@ -7,6 +7,7 @@ from .wrapper import Wrapper from .pep import PEP from .point import Point, null_point +from .tools.symbolic_scalar import SymbolicScalar __all__ = ['block_partition', 'BlockPartition', 'examples', @@ -24,4 +25,5 @@ 'pep', 'PEP', 'point', 'Point', 'null_point', 'wrapper', 'Wrapper', + 'SymbolicScalar', ] diff --git a/PEPit/expression.py b/PEPit/expression.py index 9eead3b8..863d6d78 100644 --- a/PEPit/expression.py +++ b/PEPit/expression.py @@ -4,6 +4,7 @@ from PEPit.constraint import Constraint from PEPit.tools.dict_operations import merge_dict, prune_dict +from PEPit.tools.symbolic_scalar import is_scalar, evaluate_scalar class Expression(object): @@ -151,7 +152,7 @@ def __add__(self, other): if isinstance(other, Expression): merged_decomposition_dict = merge_dict(self.decomposition_dict, other.decomposition_dict) # If other is a scalar constant, add it to the decomposition_dict of self - elif isinstance(other, int) or isinstance(other, float): + elif is_scalar(other): merged_decomposition_dict = merge_dict(self.decomposition_dict, {1: other}) # Raise an Exception in any other scenario else: @@ -249,7 +250,7 @@ def __rmul__(self, other): """ # Verify other is a scalar constant - assert isinstance(other, int) or isinstance(other, float) + assert is_scalar(other) # Multiply uniformly self's decomposition_dict by other new_decomposition_dict = dict() @@ -414,14 +415,14 @@ def eval(self): # Distinguish 3 cases: function values, inner products, and constant values if type(key) == Expression: assert key.get_is_leaf() - value += weight * key.eval() + value += evaluate_scalar(weight) * key.eval() elif type(key) == tuple: point1, point2 = key assert point1.get_is_leaf() assert point2.get_is_leaf() - value += weight * np.dot(point1.eval(), point2.eval()) + value += evaluate_scalar(weight) * np.dot(point1.eval(), point2.eval()) elif key == 1: - value += weight + value += evaluate_scalar(weight) # Raise Exception out of those 3 cases else: raise TypeError("Expressions are made of function values, inner products and constants only!" diff --git a/PEPit/function.py b/PEPit/function.py index 4285b7ea..b04d3c3f 100644 --- a/PEPit/function.py +++ b/PEPit/function.py @@ -7,6 +7,7 @@ from PEPit.psd_matrix import PSDMatrix from PEPit.tools.dict_operations import merge_dict, prune_dict +from PEPit.tools.symbolic_scalar import is_scalar class Function(object): @@ -241,7 +242,7 @@ def __rmul__(self, other): """ # Verify other is a scalar constant - assert isinstance(other, float) or isinstance(other, int) + assert is_scalar(other) # Multiply uniformly self's decomposition_dict by other new_decomposition_dict = dict() @@ -493,7 +494,7 @@ def get_class_constraints_duals(self): if isinstance(element, Constraint): new_row.append(element.eval_dual()) # If there is a scalar number, simply return it. - elif isinstance(element, float) or isinstance(element, int): + elif is_scalar(element): new_row.append(element) else: raise TypeError("The elements of table of constraints must either be Constraints objects" diff --git a/PEPit/pep.py b/PEPit/pep.py index 73bfa8e9..f305c653 100644 --- a/PEPit/pep.py +++ b/PEPit/pep.py @@ -11,6 +11,7 @@ from PEPit.function import Function from PEPit.psd_matrix import PSDMatrix from PEPit.block_partition import BlockPartition +from PEPit.tools.symbolic_scalar import evaluate_scalar class PEP(object): @@ -832,11 +833,11 @@ def check_feasibility(self, wc_value, verbose=1): ) # Get the actual dual_objective from its dict if 1 in dual_objective_expression_decomposition_dict.keys(): - dual_objective = dual_objective_expression_decomposition_dict[1] + dual_objective = evaluate_scalar(dual_objective_expression_decomposition_dict[1]) else: dual_objective = 0. # Compute the remaining terms, that should be small and only due to numerical stability errors - remaining_terms = np.sum(np.abs([value for key, value in dual_objective_expression_decomposition_dict.items() + remaining_terms = np.sum(np.abs([evaluate_scalar(value) for key, value in dual_objective_expression_decomposition_dict.items() if key != 1])) if verbose: message = "(PEPit) The worst-case guarantee proof is perfectly reconstituted" diff --git a/PEPit/point.py b/PEPit/point.py index 503f2b17..7371e0f5 100644 --- a/PEPit/point.py +++ b/PEPit/point.py @@ -3,6 +3,7 @@ from PEPit.expression import Expression from PEPit.tools.dict_operations import merge_dict, prune_dict, multiply_dicts +from PEPit.tools.symbolic_scalar import is_scalar, evaluate_scalar class Point(object): @@ -200,7 +201,7 @@ def __rmul__(self, other): """ # Multiplying by a scalar value is applying a homothety - if isinstance(other, int) or isinstance(other, float): + if is_scalar(other): # Build the decomposition of the new point new_decomposition_dict = dict() for key, value in self.decomposition_dict.items(): @@ -293,7 +294,7 @@ def eval(self): else: value = np.zeros(Point.counter) for point, weight in self.decomposition_dict.items(): - value += weight * point.eval() + value += evaluate_scalar(weight) * point.eval() self._value = value return self._value diff --git a/PEPit/psd_matrix.py b/PEPit/psd_matrix.py index c319b372..11539520 100644 --- a/PEPit/psd_matrix.py +++ b/PEPit/psd_matrix.py @@ -1,6 +1,7 @@ import numpy as np from PEPit.expression import Expression +from PEPit.tools.symbolic_scalar import is_scalar class PSDMatrix(object): @@ -141,7 +142,7 @@ def _store(matrix_of_expressions): if isinstance(matrix[i, j], Expression): pass # ... or a python scalar. If so, store it as an Expression - elif isinstance(matrix[i, j], int) or isinstance(matrix[i, j], float): + elif is_scalar(matrix[i, j]): matrix[i, j] = Expression(is_leaf=False, decomposition_dict={1: matrix[i, j]}) # Raise an Exception in any other scenario else: diff --git a/PEPit/tools/dict_operations.py b/PEPit/tools/dict_operations.py index bf55b719..7d1d4775 100644 --- a/PEPit/tools/dict_operations.py +++ b/PEPit/tools/dict_operations.py @@ -1,3 +1,6 @@ +from PEPit.tools.symbolic_scalar import evaluate_scalar + + def merge_dict(dict1, dict2): """ Merge keys of dict1 and dict2. @@ -48,11 +51,11 @@ def prune_dict(my_dict): pruned_dict = dict() # Add all entry of dict1 that does not have a null value - for key in my_dict.keys(): + for key, value in my_dict.items(): - if my_dict[key] != 0: + if evaluate_scalar(value) != 0: - pruned_dict[key] = my_dict[key] + pruned_dict[key] = value # Return pruned dict return pruned_dict diff --git a/PEPit/tools/expressions_to_matrices.py b/PEPit/tools/expressions_to_matrices.py index 60268983..dc76a748 100644 --- a/PEPit/tools/expressions_to_matrices.py +++ b/PEPit/tools/expressions_to_matrices.py @@ -2,6 +2,7 @@ from PEPit.point import Point from PEPit.expression import Expression +from PEPit.tools.symbolic_scalar import evaluate_scalar def expression_to_matrices(expression): @@ -34,16 +35,16 @@ def expression_to_matrices(expression): # Function values are stored in F if type(key) == Expression: assert key.get_is_leaf() - Fweights[key.counter] = weight + Fweights[key.counter] = evaluate_scalar(weight) # Inner products are stored in G elif type(key) == tuple: point1, point2 = key assert point1.get_is_leaf() assert point2.get_is_leaf() - Gweights[point1.counter, point2.counter] = weight + Gweights[point1.counter, point2.counter] = evaluate_scalar(weight) # Constants are simply constants elif key == 1: - cons = weight + cons = evaluate_scalar(weight) # Others don't exist and raise an Exception else: raise TypeError("Expressions are made of function values, inner products and constants only!") @@ -91,7 +92,7 @@ def expression_to_sparse_matrices(expression): if type(key) == Expression: assert key.get_is_leaf() Fweights_ind.append(key.counter) - Fweights_val.append(weight) + Fweights_val.append(evaluate_scalar(weight)) # Inner products are stored in G elif type(key) == tuple: point1, point2 = key @@ -102,16 +103,16 @@ def expression_to_sparse_matrices(expression): if (point2, point1) in expression.decomposition_dict: if point1.counter >= point2.counter: # if both entry and symmetrical entry: only append in one case weight_sym = expression.decomposition_dict[(point2, point1)] - Gweights_val.append((weight + weight_sym) / 2) + Gweights_val.append((evaluate_scalar(weight) + evaluate_scalar(weight_sym)) / 2) Gweights_indi.append(point1.counter) Gweights_indj.append(point2.counter) else: - Gweights_val.append((weight + weight_sym) / 2) + Gweights_val.append((evaluate_scalar(weight) + weight_sym) / 2) Gweights_indi.append(max(point1.counter, point2.counter)) Gweights_indj.append(min(point1.counter, point2.counter)) # Constants are simply constants elif key == 1: - cons_val = weight + cons_val = evaluate_scalar(weight) # Others don't exist and raise an Exception else: raise TypeError("Expressions are made of function values, inner products and constants only!") diff --git a/PEPit/tools/symbolic_scalar.py b/PEPit/tools/symbolic_scalar.py new file mode 100644 index 00000000..39e7c5fa --- /dev/null +++ b/PEPit/tools/symbolic_scalar.py @@ -0,0 +1,163 @@ +from fractions import Fraction +import numbers + +import sympy as sp + + +def _to_sympy_number(value, strict=False): + if isinstance(value, bool): + return sp.Integer(int(value)) + if isinstance(value, sp.Rational): + return value + if isinstance(value, Fraction): + return sp.Rational(value.numerator, value.denominator) + if isinstance(value, numbers.Integral): + return sp.Integer(int(value)) + if strict: + raise TypeError( + "SymbolicScalar only accepts exact values: int, fractions.Fraction, or sympy.Rational." + ) + if isinstance(value, numbers.Real): + return sp.Rational(str(float(value))) + raise TypeError("Unsupported scalar type: {}".format(type(value))) + + +class SymbolicScalar: + """ + Symbolic scalar with deferred numeric evaluation. + + Arithmetic operations only manipulate symbolic expressions. Numeric values are + computed when calling ``evaluate`` (directly or through ``evaluate_scalar``). + """ + + def __init__(self, value, symbol=None, _expr=None, _defaults=None): + if _expr is not None: + self._expr = sp.simplify(_expr) + self._defaults = dict(_defaults or {}) + return + + if symbol is None: + self._expr = _to_sympy_number(value, strict=True) + self._defaults = {} + else: + name = str(symbol) + self._expr = sp.Symbol(name) + self._defaults = {name: _to_sympy_number(value, strict=True)} + + @property + def symbol(self): + return str(sp.simplify(self._expr)) + + @staticmethod + def _is_numeric(value): + return isinstance(value, (SymbolicScalar, numbers.Real)) + + @staticmethod + def _merge_defaults(defaults_1, defaults_2): + merged = dict(defaults_1) + for key, value in defaults_2.items(): + if key in merged and sp.simplify(merged[key] - value) != 0: + raise ValueError("Conflicting default value for symbol '{}'.".format(key)) + merged[key] = value + return merged + + @classmethod + def _coerce(cls, value): + if isinstance(value, SymbolicScalar): + return value._expr, dict(value._defaults) + return _to_sympy_number(value, strict=False), {} + + @classmethod + def _from_op(cls, expr, defaults): + return cls(0, _expr=expr, _defaults=defaults) + + @staticmethod + def _format_substitutions(substitutions): + if substitutions is None: + return {} + formatted = {} + for key, value in substitutions.items(): + name = str(key) + formatted[name] = _to_sympy_number(value, strict=False) + return formatted + + def evaluate(self, substitutions=None): + subs = dict(self._defaults) + subs.update(self._format_substitutions(substitutions)) + + sympy_subs = {sp.Symbol(name): value for name, value in subs.items()} + evaluated = sp.simplify(self._expr.subs(sympy_subs)) + if evaluated.free_symbols: + raise ValueError( + "Missing substitutions for symbols: {}".format( + ", ".join(sorted(str(s) for s in evaluated.free_symbols)) + ) + ) + return float(sp.N(evaluated)) + + def __add__(self, other): + if not self._is_numeric(other): + return NotImplemented + other_expr, other_defaults = self._coerce(other) + defaults = self._merge_defaults(self._defaults, other_defaults) + return self._from_op(self._expr + other_expr, defaults) + + def __radd__(self, other): + return self.__add__(other) + + def __sub__(self, other): + if not self._is_numeric(other): + return NotImplemented + other_expr, other_defaults = self._coerce(other) + defaults = self._merge_defaults(self._defaults, other_defaults) + return self._from_op(self._expr - other_expr, defaults) + + def __rsub__(self, other): + if not self._is_numeric(other): + return NotImplemented + other_expr, other_defaults = self._coerce(other) + defaults = self._merge_defaults(other_defaults, self._defaults) + return self._from_op(other_expr - self._expr, defaults) + + def __mul__(self, other): + if not self._is_numeric(other): + return NotImplemented + other_expr, other_defaults = self._coerce(other) + defaults = self._merge_defaults(self._defaults, other_defaults) + return self._from_op(self._expr * other_expr, defaults) + + def __rmul__(self, other): + return self.__mul__(other) + + def __truediv__(self, other): + if not self._is_numeric(other): + return NotImplemented + other_expr, other_defaults = self._coerce(other) + defaults = self._merge_defaults(self._defaults, other_defaults) + return self._from_op(self._expr / other_expr, defaults) + + def __rtruediv__(self, other): + if not self._is_numeric(other): + return NotImplemented + other_expr, other_defaults = self._coerce(other) + defaults = self._merge_defaults(other_defaults, self._defaults) + return self._from_op(other_expr / self._expr, defaults) + + def __neg__(self): + return self._from_op(-self._expr, self._defaults) + + def __float__(self): + return self.evaluate() + + def __repr__(self): + return "SymbolicScalar({})".format(self.symbol) + + +def is_scalar(value): + return isinstance(value, (SymbolicScalar, numbers.Real)) + + +def evaluate_scalar(value, substitutions=None): + if isinstance(value, SymbolicScalar): + return value.evaluate(substitutions=substitutions) + return float(value) diff --git a/README.md b/README.md index 7216a585..22b4c3f4 100644 --- a/README.md +++ b/README.md @@ -57,6 +57,7 @@ It relies on the following Python modules: - Numpy - Scipy - Cvxpy +- Sympy - Matplotlib (for the demo notebook) diff --git a/requirements.txt b/requirements.txt index 4e650ef3..fa53a937 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,4 @@ numpy pandas scipy matplotlib +sympy diff --git a/setup.py b/setup.py index 0d8491c5..7556ea15 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ "to pep their optimization algorithms as easily as they implement them", long_description=long_description, long_description_content_type="text/markdown", - install_requires=["cvxpy>=1.1.17", "pandas>=1.0.0"], + install_requires=["cvxpy>=1.1.17", "pandas>=1.0.0", "sympy>=1.10"], url="https://github.com/PerformanceEstimation/PEPit", project_urls={ "Documentation": "https://pepit.readthedocs.io/en/{}/".format(version), diff --git a/tests/test_symbolic_substitution.py b/tests/test_symbolic_substitution.py new file mode 100644 index 00000000..08b36683 --- /dev/null +++ b/tests/test_symbolic_substitution.py @@ -0,0 +1,69 @@ +import unittest +from fractions import Fraction + +from PEPit import PEP, SymbolicScalar +from PEPit.functions import SmoothStronglyConvexQuadraticFunction + + +def build_gd_quadratic_problem(mu, L, gamma, n=1): + problem = PEP() + func = problem.declare_function(SmoothStronglyConvexQuadraticFunction, mu=mu, L=L) + + xs = func.stationary_point() + fs = func(xs) + x0 = problem.set_initial_point() + problem.set_initial_condition((x0 - xs) ** 2 <= 1) + + x = x0 + for _ in range(n): + x = x - gamma * func.gradient(x) + + problem.set_performance_metric(func(x) - fs) + return problem + + +class TestSymbolicSubstitution(unittest.TestCase): + + def test_symbolic_scalar_rejects_float_constructor(self): + with self.assertRaises(TypeError): + SymbolicScalar(0.3, "mu") + + def test_symbolic_substitution_matches_numeric_solve(self): + mu_num = 0.3 + L_num = 3.0 + gamma_num = 1 / 3 + + numeric_problem = build_gd_quadratic_problem( + mu=mu_num, L=L_num, gamma=gamma_num, n=1 + ) + numeric_tau = numeric_problem.solve(verbose=0, solver="SCS") + + mu = SymbolicScalar(Fraction(3, 10), "mu") + L = SymbolicScalar(3, "L") + gamma = SymbolicScalar(Fraction(1, 3), "gamma") + symbolic_problem = build_gd_quadratic_problem(mu=mu, L=L, gamma=gamma, n=1) + symbolic_tau = symbolic_problem.solve(verbose=0, solver="SCS") + + self.assertAlmostEqual(numeric_tau, symbolic_tau, places=5) + + def test_symbolic_defaults_are_enough_without_solve_substitution(self): + mu_num = 0.3 + L_num = 3.0 + gamma_num = 1 / 3 + + numeric_problem = build_gd_quadratic_problem( + mu=mu_num, L=L_num, gamma=gamma_num, n=1 + ) + numeric_tau = numeric_problem.solve(verbose=0, solver="SCS") + + mu = SymbolicScalar(Fraction(3, 10), "mu") + L = SymbolicScalar(3, "L") + gamma = SymbolicScalar(Fraction(1, 3), "gamma") + symbolic_problem = build_gd_quadratic_problem(mu=mu, L=L, gamma=gamma, n=1) + symbolic_tau = symbolic_problem.solve(verbose=0, solver="SCS") + + self.assertAlmostEqual(numeric_tau, symbolic_tau, places=5) + + +if __name__ == "__main__": + unittest.main()