From 3b509c000e2ece8a68953f3cdc289099f42c4868 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Wed, 13 Mar 2024 14:50:54 +0100 Subject: [PATCH] allowing symbolic powers --- qbee/polynomialization.py | 25 ++++++++++++++++++++++--- tests/test_polynomialization.py | 24 ++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 3 deletions(-) diff --git a/qbee/polynomialization.py b/qbee/polynomialization.py index 546c132..7ff38d1 100644 --- a/qbee/polynomialization.py +++ b/qbee/polynomialization.py @@ -230,10 +230,11 @@ def _try_convert_to_polynomial(self, expr: sp.Expr) -> sp.Expr | None: replaced_non_pow_lazy = expr.subs( {v: k for k, v in self._substitution_equations.items() if not v.is_Pow}) replaced_pow_lazy = replaced_non_pow_lazy \ + .replace(sp.Pow, self._replace_symbolic_pow) \ .replace(sp.Pow, self._replace_irrational_pow) \ .replace(sp.Pow, self._replace_negative_integer_pow) except Exception as e: - warnings.warn("Substituting new variables failed to produce the expected calculations," + earnings.warn("Substituting new variables failed to produce the expected calculations," " so the calculation was done in an alternative mode. " "If you see this message, please let us know in the Issues: " "https://github.com/AndreyBychkov/QBee/issues", RuntimeWarning) @@ -246,6 +247,7 @@ def _try_convert_to_polynomial(self, expr: sp.Expr) -> sp.Expr | None: # eager evaluation replaced_non_pow = expr.subs({v: k for k, v in self._substitution_equations.items() if not v.is_Pow}) replaced_pow = replaced_non_pow \ + .replace(sp.Pow, self._replace_symbolic_pow) \ .replace(sp.Pow, self._replace_irrational_pow) \ .replace(sp.Pow, self._replace_negative_integer_pow) return replaced_pow if self._is_expr_polynomial(replaced_pow) else None @@ -263,6 +265,14 @@ def _replace_irrational_pow(self, base, exp): return new_vars[0] if len(new_vars) != 0 else base ** exp return base ** exp + def _replace_symbolic_pow(self, base, exp): + if not (exp.is_Float or exp.is_Integer): + for k, v in self._substitution_equations.items(): + if v.is_Pow and v.base == base: + if sp.simplify(exp / v.exp).is_Integer: + return k ** sp.simplify(exp / v.exp) + return base ** exp + def _is_expr_polynomial(self, expr: sp.Expr): if not self.variables.laurent: # performance optimization return expr.is_polynomial(*self.variables.state, *self.variables.input) @@ -457,8 +467,17 @@ def available_substitutions(system: EquationSystem) -> set[sp.Expr]: return filter_laurent_monoms(system, unused_subs) -def find_pow(values: Collection[sp.Expr], pow: sp.Pow) -> sp.Pow | None: - all_matches = [v for v in values if v.is_Pow and v.base == pow.base and math.isclose(v.exp, pow.exp)] +def find_pow(values: Collection[sp.Expr], pw: sp.Pow) -> sp.Pow | None: + all_matches = [] + numtypes = (sp.core.numbers.Float, sp.core.numbers.Integer, float, int) + for v in values: + if v.is_Pow and v.base == pw.base: + if type(v.exp) in numtypes and type(pw.exp) in numtypes: + if math.isclose(v.exp, pw.exp): + all_matches.append(v) + else: + if sp.simplify(v.exp - pw.exp) == 0: + all_matches.append(v) if all_matches: return all_matches[0] return None diff --git a/tests/test_polynomialization.py b/tests/test_polynomialization.py index 7d0c769..0ecdb60 100644 --- a/tests/test_polynomialization.py +++ b/tests/test_polynomialization.py @@ -94,6 +94,30 @@ def test_nested_functions(): res = polynomialize(system) assert len(res) == 4 +def test_parametric_pow(): + x = functions("x") + a = parameters("a") + system = [(x, x**a)] + res = polynomialize(system) + assert len(res) == 2 + + x = functions("x", laurent=False) + res = polynomialize([(x, x**a)]) + assert len(res) == 3 + + res = polynomialize([(x, x**(-a))]) + assert len(res) == 3 + +def test_hill(): + x = functions("x") + a, b = parameters("a b") + system = [(x, 1 / (a + x**b))] + res = polynomialize(system) + assert len(res) == 3 + + x = functions("x", laurent=False) + res = polynomialize([(x, 1 / (a + x**b))]) + assert len(res) == 4 def test_parameter(): x, y = functions("x, y", laurent=False)