Piecewise conditions as ineqalities
32 views
Skip to first unread message
Plamen Dimitrov
Dec 6, 2024, 4:08:03 PM12/6/24
to sympy
Hi all, I think this is the closest to a forum regarding sympy I can think of so I want to RFC here if possible. I need to find general attractors (degenrate, limit cycles, etc.) beyond equilibrium points in a dynamical system that may have some discontinuous regions. For the purpose I needed to solve systems of piecewise function containing equations and wrote this code for the purpose:
```python
import sympy as sp
from itertools import product
from typing import List, Dict, Any, Tuple, Set, Union
class MatrixPiecewiseAnalyzer:
def __init__(self, expr: Union[sp.Matrix, sp.Expr], variables: List[sp.Symbol]):
"""
Initialize analyzer with a matrix/vector of piecewise functions or a single expression.
Args:
expr: SymPy Matrix, Vector, or expression containing piecewise functions
variables: List of SymPy symbols used in the function
"""
self.expr = expr
self.variables = variables
self.is_matrix = isinstance(expr, sp.Matrix)
self.pieces_by_element = self._extract_all_pieces()
self.conditions = self._extract_all_conditions()
def _extract_pieces(self, expr: sp.Expr) -> List[Tuple[sp.Expr, sp.Basic]]:
"""
Recursively extract pieces and conditions from a single expression.
"""
expr = sp.piecewise_fold(expr)
if isinstance(expr, sp.Piecewise):
pieces = []
for arg in expr.args:
if isinstance(arg[0], sp.Piecewise):
nested_pieces = self._extract_pieces(arg[0])
pieces.extend([
(piece[0], sp.And(piece[1], arg[1]))
for piece in nested_pieces
])
else:
pieces.append(arg)
return pieces
else:
# Handle non-piecewise terms
return [(expr, sp.true)]
def _extract_all_pieces(self) -> Dict[Tuple[int, int], List[Tuple[sp.Expr, sp.Basic]]]:
"""
Extract pieces from all matrix/vector elements.
Returns:
Dictionary mapping (i,j) positions to lists of (expression, condition) tuples
"""
pieces_dict = {}
if self.is_matrix:
for i in range(self.expr.rows):
for j in range(self.expr.cols):
pieces_dict[(i,j)] = self._extract_pieces(self.expr[i,j])
else:
# Treat as single expression
pieces_dict[(0,0)] = self._extract_pieces(self.expr)
return pieces_dict
def _extract_all_conditions(self) -> Set[sp.Basic]:
"""
Extract all unique atomic conditions from all matrix/vector elements.
"""
conditions = set()
def extract_atomic_conditions(cond):
if isinstance(cond, (sp.And, sp.Or)):
for arg in cond.args:
extract_atomic_conditions(arg)
else:
conditions.add(cond)
for pieces in self.pieces_by_element.values():
for _, cond in pieces:
if cond != sp.true:
extract_atomic_conditions(cond)
return conditions
def _generate_region_combinations(self) -> List[sp.Basic]:
"""
Generate all possible combinations of conditions and their negations.
"""
atomic_conditions = list(self.conditions)
regions = []
for combo in product([True, False], repeat=len(atomic_conditions)):
region_conditions = []
for cond, use_positive in zip(atomic_conditions, combo):
region_conditions.append(cond if use_positive else sp.Not(cond))
regions.append(sp.And(*region_conditions))
return regions
def _get_expression_for_region(self, pos: Tuple[int, int], region: sp.Basic) -> sp.Expr:
"""
Find the expression that applies in a given region for a specific matrix position.
"""
pieces = self.pieces_by_element[pos]
for expr, cond in pieces:
if sp.simplify(sp.And(cond, region)) != False:
return expr
return None
def _build_matrix_for_region(self, region: sp.Basic) -> sp.Matrix:
"""
Build the complete matrix expression for a given region.
"""
if not self.is_matrix:
expr = self._get_expression_for_region((0,0), region)
return expr if expr is not None else sp.zeros(1, 1)
rows, cols = self.expr.rows, self.expr.cols
matrix_elements = []
for i in range(rows):
row_elements = []
for j in range(cols):
expr = self._get_expression_for_region((i,j), region)
row_elements.append(expr if expr is not None else 0)
matrix_elements.append(row_elements)
return sp.Matrix(matrix_elements)
def _solve_region(self, equations: list[sp.Expr], region: sp.Basic) -> list[sp.Basic]:
"""Helper to solve in specific region using solveset"""
try:
# Try to solve the selected expression
# Note: we need to convert to FiniteSet for intersection with region
solutions = sp.nonlinsolve(equations, self.variables)
# Filter solutions that satisfy the region conditions
valid_solutions = []
for sol in solutions:
substitutions = dict(zip(self.variables, sol))
region_intersection = region.subs(substitutions)
# Check if solution components are still symbolic and thus general
#if all(isinstance(x, (sp.Number, sp.Symbol)) for x in sol):
if not all(s.is_real for s in sol):
raise NotImplementedError("Cannot handle degraded, limit cycle, or other attractors - only equilibrium points")
if region_intersection == True:
valid_solutions.append(sol)
return valid_solutions
except Exception as e:
return f"Could not solve: {e}"
def solve_critical_points(self) -> Dict[str, Any]:
"""
Find critical points in each region by solving the system of equations.
"""
results = {}
regions = self._generate_region_combinations()
for region in regions:
# Build matrix for this region
matrix_expr = self._build_matrix_for_region(region)
# For matrix expressions, solve the system matrix = 0
if self.is_matrix:
equations = [elem for elem in matrix_expr]
else:
equations = [matrix_expr]
try:
results[str(region)] = self._solve_region(equations, region)
except Exception as e:
results[str(region)] = f"Could not solve: {str(e)}"
return results
def analyze(self) -> Dict[str, Any]:
"""
Perform complete analysis of the piecewise matrix/vector function.
"""
return {
'pieces_by_element': self.pieces_by_element,
'atomic_conditions': self.conditions,
'critical_points': self.solve_critical_points()
}
# Example usage with your matrix
if __name__ == "__main__":
# Create the matrix from your example
x1, x2, x3, y1, y2, y3 = sp.symbols('x1 x2 x3 y1 y2 y3')
matrix = sp.Matrix([
[2*x1 + 2*x2 + 2*x3 - 2*y1 - 2*y2 - 2*y3 + sp.Piecewise((2*x1 - 2*y1, x1 < y1), (0, True))],
[2*x1 + 2*x2 + 2*x3 - 2*y1 - 2*y2 - 2*y3 + sp.Piecewise((2*x2 - 2*y2, x2 < y2), (0, True))],
[2*x1 + 2*x2 + 2*x3 - 2*y1 - 2*y2 - 2*y3 + sp.Piecewise((2*x3 - 2*y3, x3 < y3), (0, True))]
])
#variables = x1, x2, x3 = sp.symbols('x1 x2 x3')
#matrix = sp.Matrix([
# [1*x1 + 2*x2 + 3*x3 + sp.Piecewise((2*x1 - 2*x3, x1 >= 0), (0, True))],
# [1*x1 + 2*x2 + 4*x3 + sp.Piecewise((3*x1 - 3*x2, x1 >= 0), (0, True))],
# [1*x1 + 2*x2 + 5*x3 + sp.Piecewise((x1, x1 >= 0), (0, True))],
#])
# Analyze the matrix
analyzer = MatrixPiecewiseAnalyzer(matrix, variables)
results = analyzer.analyze()
print("\nPiecewise vector results:")
print("\nAtomic conditions:")
print(results['atomic_conditions'])
print("\nCritical points by region:")
for region, solutions in results['critical_points'].items():
print(f" Region: {region}")
print(f" Solutions: {solutions}")
```
My questions are twofold:
1. Is there any way to treat the `Relational` instance in the piecewise conditions as inequality part of the solved system? If not how about treating it is some solution set to intersect with `solveset` results? So far I could not find a method within sympy to do either.
2. Does anyone thing any part of this code would be useful upstream for contribution in any way?
```python
import sympy as sp
from itertools import product
from typing import List, Dict, Any, Tuple, Set, Union
class MatrixPiecewiseAnalyzer:
def __init__(self, expr: Union[sp.Matrix, sp.Expr], variables: List[sp.Symbol]):
"""
Initialize analyzer with a matrix/vector of piecewise functions or a single expression.
Args:
expr: SymPy Matrix, Vector, or expression containing piecewise functions
variables: List of SymPy symbols used in the function
"""
self.expr = expr
self.variables = variables
self.is_matrix = isinstance(expr, sp.Matrix)
self.pieces_by_element = self._extract_all_pieces()
self.conditions = self._extract_all_conditions()
def _extract_pieces(self, expr: sp.Expr) -> List[Tuple[sp.Expr, sp.Basic]]:
"""
Recursively extract pieces and conditions from a single expression.
"""
expr = sp.piecewise_fold(expr)
if isinstance(expr, sp.Piecewise):
pieces = []
for arg in expr.args:
if isinstance(arg[0], sp.Piecewise):
nested_pieces = self._extract_pieces(arg[0])
pieces.extend([
(piece[0], sp.And(piece[1], arg[1]))
for piece in nested_pieces
])
else:
pieces.append(arg)
return pieces
else:
# Handle non-piecewise terms
return [(expr, sp.true)]
def _extract_all_pieces(self) -> Dict[Tuple[int, int], List[Tuple[sp.Expr, sp.Basic]]]:
"""
Extract pieces from all matrix/vector elements.
Returns:
Dictionary mapping (i,j) positions to lists of (expression, condition) tuples
"""
pieces_dict = {}
if self.is_matrix:
for i in range(self.expr.rows):
for j in range(self.expr.cols):
pieces_dict[(i,j)] = self._extract_pieces(self.expr[i,j])
else:
# Treat as single expression
pieces_dict[(0,0)] = self._extract_pieces(self.expr)
return pieces_dict
def _extract_all_conditions(self) -> Set[sp.Basic]:
"""
Extract all unique atomic conditions from all matrix/vector elements.
"""
conditions = set()
def extract_atomic_conditions(cond):
if isinstance(cond, (sp.And, sp.Or)):
for arg in cond.args:
extract_atomic_conditions(arg)
else:
conditions.add(cond)
for pieces in self.pieces_by_element.values():
for _, cond in pieces:
if cond != sp.true:
extract_atomic_conditions(cond)
return conditions
def _generate_region_combinations(self) -> List[sp.Basic]:
"""
Generate all possible combinations of conditions and their negations.
"""
atomic_conditions = list(self.conditions)
regions = []
for combo in product([True, False], repeat=len(atomic_conditions)):
region_conditions = []
for cond, use_positive in zip(atomic_conditions, combo):
region_conditions.append(cond if use_positive else sp.Not(cond))
regions.append(sp.And(*region_conditions))
return regions
def _get_expression_for_region(self, pos: Tuple[int, int], region: sp.Basic) -> sp.Expr:
"""
Find the expression that applies in a given region for a specific matrix position.
"""
pieces = self.pieces_by_element[pos]
for expr, cond in pieces:
if sp.simplify(sp.And(cond, region)) != False:
return expr
return None
def _build_matrix_for_region(self, region: sp.Basic) -> sp.Matrix:
"""
Build the complete matrix expression for a given region.
"""
if not self.is_matrix:
expr = self._get_expression_for_region((0,0), region)
return expr if expr is not None else sp.zeros(1, 1)
rows, cols = self.expr.rows, self.expr.cols
matrix_elements = []
for i in range(rows):
row_elements = []
for j in range(cols):
expr = self._get_expression_for_region((i,j), region)
row_elements.append(expr if expr is not None else 0)
matrix_elements.append(row_elements)
return sp.Matrix(matrix_elements)
def _solve_region(self, equations: list[sp.Expr], region: sp.Basic) -> list[sp.Basic]:
"""Helper to solve in specific region using solveset"""
try:
# Try to solve the selected expression
# Note: we need to convert to FiniteSet for intersection with region
solutions = sp.nonlinsolve(equations, self.variables)
# Filter solutions that satisfy the region conditions
valid_solutions = []
for sol in solutions:
substitutions = dict(zip(self.variables, sol))
region_intersection = region.subs(substitutions)
# Check if solution components are still symbolic and thus general
#if all(isinstance(x, (sp.Number, sp.Symbol)) for x in sol):
if not all(s.is_real for s in sol):
raise NotImplementedError("Cannot handle degraded, limit cycle, or other attractors - only equilibrium points")
if region_intersection == True:
valid_solutions.append(sol)
return valid_solutions
except Exception as e:
return f"Could not solve: {e}"
def solve_critical_points(self) -> Dict[str, Any]:
"""
Find critical points in each region by solving the system of equations.
"""
results = {}
regions = self._generate_region_combinations()
for region in regions:
# Build matrix for this region
matrix_expr = self._build_matrix_for_region(region)
# For matrix expressions, solve the system matrix = 0
if self.is_matrix:
equations = [elem for elem in matrix_expr]
else:
equations = [matrix_expr]
try:
results[str(region)] = self._solve_region(equations, region)
except Exception as e:
results[str(region)] = f"Could not solve: {str(e)}"
return results
def analyze(self) -> Dict[str, Any]:
"""
Perform complete analysis of the piecewise matrix/vector function.
"""
return {
'pieces_by_element': self.pieces_by_element,
'atomic_conditions': self.conditions,
'critical_points': self.solve_critical_points()
}
# Example usage with your matrix
if __name__ == "__main__":
# Create the matrix from your example
x1, x2, x3, y1, y2, y3 = sp.symbols('x1 x2 x3 y1 y2 y3')
matrix = sp.Matrix([
[2*x1 + 2*x2 + 2*x3 - 2*y1 - 2*y2 - 2*y3 + sp.Piecewise((2*x1 - 2*y1, x1 < y1), (0, True))],
[2*x1 + 2*x2 + 2*x3 - 2*y1 - 2*y2 - 2*y3 + sp.Piecewise((2*x2 - 2*y2, x2 < y2), (0, True))],
[2*x1 + 2*x2 + 2*x3 - 2*y1 - 2*y2 - 2*y3 + sp.Piecewise((2*x3 - 2*y3, x3 < y3), (0, True))]
])
#variables = x1, x2, x3 = sp.symbols('x1 x2 x3')
#matrix = sp.Matrix([
# [1*x1 + 2*x2 + 3*x3 + sp.Piecewise((2*x1 - 2*x3, x1 >= 0), (0, True))],
# [1*x1 + 2*x2 + 4*x3 + sp.Piecewise((3*x1 - 3*x2, x1 >= 0), (0, True))],
# [1*x1 + 2*x2 + 5*x3 + sp.Piecewise((x1, x1 >= 0), (0, True))],
#])
# Analyze the matrix
analyzer = MatrixPiecewiseAnalyzer(matrix, variables)
results = analyzer.analyze()
print("\nPiecewise vector results:")
print("\nAtomic conditions:")
print(results['atomic_conditions'])
print("\nCritical points by region:")
for region, solutions in results['critical_points'].items():
print(f" Region: {region}")
print(f" Solutions: {solutions}")
```
My questions are twofold:
1. Is there any way to treat the `Relational` instance in the piecewise conditions as inequality part of the solved system? If not how about treating it is some solution set to intersect with `solveset` results? So far I could not find a method within sympy to do either.
2. Does anyone thing any part of this code would be useful upstream for contribution in any way?
Reply all
Reply to author
Forward
0 new messages