diff options
| author | Michael Schneeberger <michael.schneeberger@fhnw.ch> | 2022-07-30 11:34:55 +0200 |
|---|---|---|
| committer | Michael Schneeberger <michael.schneeberger@fhnw.ch> | 2022-07-30 11:34:55 +0200 |
| commit | 05d372c6a05e0184fb02ad6fc4b612de7edd91a8 (patch) | |
| tree | a720ee1e82fc10e6c2fc100343836b3b2d8fdb6e | |
| parent | add eye, sum and symmetric operation (diff) | |
| download | polymatrix-05d372c6a05e0184fb02ad6fc4b612de7edd91a8.tar.gz polymatrix-05d372c6a05e0184fb02ad6fc4b612de7edd91a8.zip | |
add polynomial operations for sos optimization
94 files changed, 1955 insertions, 1303 deletions
diff --git a/polymatrix/__init__.py b/polymatrix/__init__.py index 85def18..aa96d21 100644 --- a/polymatrix/__init__.py +++ b/polymatrix/__init__.py @@ -8,7 +8,6 @@ import scipy.sparse from polymatrix.expression.expression import Expression from polymatrix.expression.expressionstate import ExpressionState -from polymatrix.expression.init.initaccumulateexpr import init_accumulate_expr from polymatrix.expression.init.initblockdiagexpr import init_block_diag_expr from polymatrix.expression.init.initexpression import init_expression from polymatrix.expression.init.initeyeexpr import init_eye_expr @@ -38,31 +37,31 @@ def from_polymatrix( ) -def accumulate( - expr, - func, - initial = None, -): - def lifted_func(acc, polymat: PolyMatrix): +# def accumulate( +# expr, +# func, +# initial = None, +# ): +# def lifted_func(acc, polymat: PolyMatrix): - # # print(f'{terms=}') - # print(f'{terms=}') - - lifeted_val = init_expression( - underlying=init_from_terms_expr( - terms=polymat.terms, - shape=polymat.shape, - ), - ) - return func(acc, lifeted_val) +# # # print(f'{terms=}') +# # print(f'{terms=}') + +# lifeted_val = init_expression( +# underlying=init_from_terms_expr( +# terms=polymat.terms, +# shape=polymat.shape, +# ), +# ) +# return func(acc, lifeted_val) - return init_expression( - underlying=init_accumulate_expr( - underlying=expr.underlying, - acc_func=lifted_func, - initial=initial, - ), - ) +# return init_expression( +# underlying=init_accumulate_expr( +# underlying=expr.underlying, +# acc_func=lifted_func, +# initial=initial, +# ), +# ) def v_stack( expressions: tuple[Expression], @@ -96,20 +95,14 @@ def eye( ) -# def kkt( -# cost: Expression, -# variables: Expression, -# equality: Expression = None, -# inequality: Expression = None, +# def sos_matrix( +# underlying: Expression, # ): -# return init_expression( -# init_kkt_expr( -# cost=cost, -# equality=equality, -# variables=variables, -# inequality=inequality, -# ) -# ) +# def func(state: ExpressionState): +# underlying = + +# return init_state_monad(func) + def kkt_equality( variables: Expression, @@ -293,37 +286,6 @@ def kkt_inequality( return init_state_monad(func) -# def to_polymatrix( -# expr: Expression, -# ): -# def func(state: ExpressionState): -# state, polymatrix = expr.apply(state) - -# return state, polymatrix - -# return init_state_monad(func) - -# def to_linear_matrix( -# expr: Expression, -# variables: Expression, -# ) -> StateMonad[ExpressionState, tuple[Expression, ...]]: -# def func(state: ExpressionState): -# state, variable_indices = get_variable_indices(state, variables) - -# def gen_matrices(): -# for variable_index in variable_indices: -# yield init_linear_matrix_in_expr( -# underlying=expr, -# variable=variable_index, -# ) - -# matrices = tuple(gen_matrices()) - -# return state, matrices - -# return StateMonad.init(func) - - def shape( expr: Expression, ) -> StateMonadMixin[ExpressionState, tuple[int, ...]]: @@ -337,8 +299,8 @@ def shape( @dataclasses.dataclass class MatrixEquations: - matrix_equations: tuple[tuple[np.ndarray, ...], ...] - auxillary_matrix_equations: typing.Optional[tuple[np.ndarray, ...]] + matrix_equations: tuple[dict[int, np.ndarray], ...] + auxillary_matrix_equations: typing.Optional[dict[int, np.ndarray]] variable_index: tuple[int, ...] state: ExpressionState @@ -351,15 +313,33 @@ class MatrixEquations: if index < len(self.auxillary_matrix_equations): yield self.auxillary_matrix_equations[index] - matrix_1 = np.vstack(tuple(gen_matrices(0))) - matrix_2 = np.vstack(tuple(gen_matrices(1))) - matrix_3 = scipy.sparse.vstack(tuple(gen_matrices(2))) + indices = set(key for equations in self.matrix_equations + (self.auxillary_matrix_equations,) for key in equations.keys()) - return (matrix_1, matrix_2, matrix_3) + def gen_matrices(): + for index in indices: + if index <= 1: + yield index, np.vstack(tuple(gen_matrices(index))) + else: + yield index, scipy.sparse.vstack(tuple(gen_matrices(index))) + + # matrix_1 = np.vstack(tuple(gen_matrices(0))) + # matrix_2 = np.vstack(tuple(gen_matrices(1))) + # matrix_3 = scipy.sparse.vstack(tuple(gen_matrices(2))) + + return dict(gen_matrices()) def get_value(self, variable, value): variable_indices = get_variable_indices(self.state, variable)[1] - value_index = list(self.variable_index.index(variable_index) for variable_index in variable_indices) + + def gen_value_index(): + for variable_index in variable_indices: + try: + yield self.variable_index.index(variable_index) + except ValueError: + raise ValueError(f'{variable_index} not found in {self.variable_index}') + + value_index = list(gen_value_index()) + return value[value_index] def set_value(self, variable, value): @@ -371,112 +351,122 @@ class MatrixEquations: def to_matrix_equations( - expr: tuple[Expression], + expressions: tuple[Expression], + variables: Expression = None, ) -> StateMonadMixin[ExpressionState, tuple[tuple[tuple[np.ndarray, ...], ...], tuple[int, ...]]]: - if isinstance(expr, Expression): - expr = (expr,) + if isinstance(expressions, Expression): + expressions = (expressions,) def func(state: ExpressionState): - def acc_underlying(acc, v): + def acc_underlying_application(acc, v): state, underlying_list = acc state, underlying = v.apply(state) - assert underlying.shape[1] == 1, f'{underlying.shape[1]} is not 1' + assert underlying.shape[1] == 1, f'{underlying.shape[1]=} is not 1' return state, underlying_list + (underlying,) *_, (state, underlying_list) = tuple(itertools.accumulate( - expr, - acc_underlying, + expressions, + acc_underlying_application, initial=(state, tuple()), )) - # state, underlying = expr.apply(state) + if variables is None: - def gen_used_variables(): - def gen_used_auxillary_variables(considered): - monomial_terms = state.auxillary_equations[considered[-1]] - for monomial in monomial_terms.keys(): - for variable in monomial: - yield variable + def gen_used_variables(): + def gen_used_auxillary_variables(considered): + monomial_terms = state.auxillary_equations[considered[-1]] + for monomial in monomial_terms.keys(): + for variable in monomial: + yield variable - if variable not in considered and variable in state.auxillary_equations: - yield from gen_used_auxillary_variables(considered + (variable,)) + if variable not in considered and variable in state.auxillary_equations: + yield from gen_used_auxillary_variables(considered + (variable,)) - for underlying in underlying_list: - for row in range(underlying.shape[0]): - for col in range(underlying.shape[1]): + for underlying in underlying_list: + for row in range(underlying.shape[0]): + for col in range(underlying.shape[1]): - try: - underlying_terms = underlying.get_poly(row, col) - except KeyError: - continue + try: + underlying_terms = underlying.get_poly(row, col) + except KeyError: + continue - for monomial in underlying_terms.keys(): - for variable in monomial: - yield variable + for monomial in underlying_terms.keys(): + for variable, _ in monomial: + yield variable - if variable in state.auxillary_equations: - yield from gen_used_auxillary_variables((variable,)) + if variable in state.auxillary_equations: + yield from gen_used_auxillary_variables((variable,)) + + ordered_variable_index = tuple(sorted(set(gen_used_variables()))) + + else: + state, ordered_variable_index = get_variable_indices(state, variables) - used_variables = set(gen_used_variables()) - ordered_variable_index = tuple(sorted(used_variables)) variable_index_map = {old: new for new, old in enumerate(ordered_variable_index)} n_param = len(ordered_variable_index) + @dataclasses.dataclass + class MatrixBuffer: + data: dict + n_row: int + n_param: int + + def add_buffer(self, index: int): + if index <= 1: + buffer = np.zeros((self.n_row, self.n_param**index), dtype=np.double) + + else: + buffer = scipy.sparse.dok_array((self.n_row, self.n_param**index), dtype=np.double) + + self.data[index] = buffer + + def add(self, row: int, col: int, index: int, value: float): + if index not in self.data: + self.add_buffer(index) + + self.data[index][row, col] = value + + def __getitem__(self, key): + if key not in self.data: + self.add_buffer(key) + + return self.data[key] + def gen_underlying_matrices(): for underlying in underlying_list: + n_row = underlying.shape[0] - A = np.zeros((n_row, 1), dtype=np.double) - B = np.zeros((n_row, n_param), dtype=np.double) - C = scipy.sparse.dok_array((n_row, n_param**2), dtype=np.double) - - # def populate_matrices(monomial_terms, row): - # for monomial, value in monomial_terms.items(): - # new_monomial = tuple(variable_index_map[var] for var in monomial) - # col = monomial_to_index(n_param, new_monomial) - - # match len(new_monomial): - # case 0: - # A[row, col] = value - # case 1: - # B[row, col] = value - # case 2: - # C[row, col] = value - # case _: - # raise Exception(f'illegal case {new_monomial=}') - - for row in range(underlying.shape[0]): + buffer = MatrixBuffer( + data={}, + n_row=n_row, + n_param=n_param, + ) + + for row in range(n_row): try: underlying_terms = underlying.get_poly(row, 0) except KeyError: continue - - # populate_matrices( - # monomial_terms=underlying_terms, - # row=row, - # ) for monomial, value in underlying_terms.items(): - new_monomial = tuple(variable_index_map[var] for var in monomial) + try: + new_monomial = tuple(variable_index_map[var] for var, count in monomial for _ in range(count)) + except KeyError: + raise KeyError(f'{monomial=} is incompatible with {variable_index_map=}') col = monomial_to_index(n_param, new_monomial) - match len(new_monomial): - case 0: - A[row, col] = value - case 1: - B[row, col] = value - case 2: - C[row, col] = value - case _: - raise Exception(f'illegal case {new_monomial=}') + buffer.add(row, col, sum(count for _, count in monomial), value) - yield A, B, C + # yield A, B, C + yield buffer underlying_matrices = tuple(gen_underlying_matrices()) @@ -495,26 +485,21 @@ def to_matrix_equations( auxillary_matrix_equations = None else: - A = np.zeros((n_row, 1), dtype=np.double) - B = np.zeros((n_row, n_param), dtype=np.double) - C = scipy.sparse.dok_array((n_row, n_param**2), dtype=np.double) + buffer = MatrixBuffer( + data={}, + n_row=n_row, + n_param=n_param, + ) for row, (key, monomial_terms) in enumerate(auxillary_equations): for monomial, value in monomial_terms.items(): - new_monomial = tuple(variable_index_map[var] for var in monomial) + # new_monomial = tuple(variable_index_map[var] for var, _ in monomial) + new_monomial = tuple(variable_index_map[var] for var, count in monomial for _ in range(count)) col = monomial_to_index(n_param, new_monomial) - match len(new_monomial): - case 0: - A[row, col] = value - case 1: - B[row, col] = value - case 2: - C[row, col] = value - case _: - raise Exception(f'illegal case {new_monomial=}') - - auxillary_matrix_equations = (A, B, C) + buffer.add(row, col, sum(count for _, count in monomial), value) + + auxillary_matrix_equations = buffer result = MatrixEquations( matrix_equations=underlying_matrices, diff --git a/polymatrix/expression/accumulateexpr.py b/polymatrix/expression/accumulateexpr.py deleted file mode 100644 index efcbc39..0000000 --- a/polymatrix/expression/accumulateexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.accumulateexprmixin import AccumulateExprMixin - -class AccumulateExpr(AccumulateExprMixin): - pass diff --git a/polymatrix/expression/divergenceexpr.py b/polymatrix/expression/divergenceexpr.py new file mode 100644 index 0000000..9207bd2 --- /dev/null +++ b/polymatrix/expression/divergenceexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.divergenceexprmixin import DivergenceExprMixin + +class DivergenceExpr(DivergenceExprMixin): + pass diff --git a/polymatrix/expression/impl/cacheexprimpl.py b/polymatrix/expression/impl/cacheexprimpl.py index a6a74a1..f6e45d3 100644 --- a/polymatrix/expression/impl/cacheexprimpl.py +++ b/polymatrix/expression/impl/cacheexprimpl.py @@ -3,6 +3,6 @@ from polymatrix.expression.cacheexpr import CacheExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -@dataclass_abc.dataclass_abc(frozen=True) +@dataclass_abc.dataclass_abc(frozen=True, eq=False) class CacheExprImpl(CacheExpr): underlying: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/combinationsexprimpl.py b/polymatrix/expression/impl/combinationsexprimpl.py index 2e809a2..2438e01 100644 --- a/polymatrix/expression/impl/combinationsexprimpl.py +++ b/polymatrix/expression/impl/combinationsexprimpl.py @@ -5,5 +5,5 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) class CombinationsExprImpl(CombinationsExpr): - underlying: ExpressionBaseMixin + monomials: ExpressionBaseMixin number: int diff --git a/polymatrix/expression/impl/divergenceexprimpl.py b/polymatrix/expression/impl/divergenceexprimpl.py new file mode 100644 index 0000000..0ea9cdb --- /dev/null +++ b/polymatrix/expression/impl/divergenceexprimpl.py @@ -0,0 +1,10 @@ +import dataclass_abc +from polymatrix.expression.divergenceexpr import DivergenceExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from typing import Union + +@dataclass_abc.dataclass_abc(frozen=True) +class DivergenceExprImpl(DivergenceExpr): + underlying: ExpressionBaseMixin + variables: tuple diff --git a/polymatrix/expression/impl/divisionexprimpl.py b/polymatrix/expression/impl/divisionexprimpl.py index c2c8966..639210e 100644 --- a/polymatrix/expression/impl/divisionexprimpl.py +++ b/polymatrix/expression/impl/divisionexprimpl.py @@ -3,7 +3,7 @@ from polymatrix.expression.divisionexpr import DivisionExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -@dataclass_abc.dataclass_abc(frozen=True) +@dataclass_abc.dataclass_abc(frozen=True, eq=False) class DivisionExprImpl(DivisionExpr): left: ExpressionBaseMixin right: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/kktexprimpl.py b/polymatrix/expression/impl/kktexprimpl.py deleted file mode 100644 index 415b206..0000000 --- a/polymatrix/expression/impl/kktexprimpl.py +++ /dev/null @@ -1,12 +0,0 @@ -import dataclass_abc -from polymatrix.expression.kktexpr import KKTExpr - -from polymatrix.expression.mixins.expressionmixin import ExpressionMixin -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin - -@dataclass_abc.dataclass_abc(frozen=True) -class KKTExprImpl(KKTExpr): - cost: ExpressionMixin - equality: ExpressionMixin - inequality: ExpressionMixin - variables: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/linearinexprimpl.py b/polymatrix/expression/impl/linearinexprimpl.py index d97051f..346b19b 100644 --- a/polymatrix/expression/impl/linearinexprimpl.py +++ b/polymatrix/expression/impl/linearinexprimpl.py @@ -6,4 +6,6 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) class LinearInExprImpl(LinearInExpr): underlying: ExpressionBaseMixin - variables: tuple + monomials: ExpressionBaseMixin + variables: ExpressionBaseMixin + ignore_unmatched: bool diff --git a/polymatrix/expression/impl/linearinmonomialsinexprimpl.py b/polymatrix/expression/impl/linearinmonomialsinexprimpl.py new file mode 100644 index 0000000..537f164 --- /dev/null +++ b/polymatrix/expression/impl/linearinmonomialsinexprimpl.py @@ -0,0 +1,9 @@ +import dataclass_abc +from polymatrix.expression.linearinmonomialsinexpr import LinearInMonomialsInExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class LinearInMonomialsInExprImpl(LinearInMonomialsInExpr): + underlying: ExpressionBaseMixin + variables: tuple diff --git a/polymatrix/expression/impl/linearmonomialsexprimpl.py b/polymatrix/expression/impl/linearmonomialsexprimpl.py new file mode 100644 index 0000000..905ec77 --- /dev/null +++ b/polymatrix/expression/impl/linearmonomialsexprimpl.py @@ -0,0 +1,9 @@ +import dataclass_abc +from polymatrix.expression.linearmonomialsexpr import LinearMonomialsExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class LinearMonomialsExprImpl(LinearMonomialsExpr): + underlying: ExpressionBaseMixin + variables: tuple diff --git a/polymatrix/expression/impl/parametrizeexprimpl.py b/polymatrix/expression/impl/parametrizeexprimpl.py new file mode 100644 index 0000000..a1fca6e --- /dev/null +++ b/polymatrix/expression/impl/parametrizeexprimpl.py @@ -0,0 +1,9 @@ +import dataclass_abc +from polymatrix.expression.parametrizeexpr import ParametrizeExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True, eq=False) +class ParametrizeExprImpl(ParametrizeExpr): + underlying: ExpressionBaseMixin + name: str diff --git a/polymatrix/expression/impl/parametrizetermsexprimpl.py b/polymatrix/expression/impl/parametrizetermsexprimpl.py index b851306..32beac2 100644 --- a/polymatrix/expression/impl/parametrizetermsexprimpl.py +++ b/polymatrix/expression/impl/parametrizetermsexprimpl.py @@ -1,10 +1,10 @@ import dataclass_abc -from polymatrix.expression.parametrizesymbolsexpr import ParametrizeSymbolsExpr +from polymatrix.expression.parametrizetermsexpr import ParametrizeTermsExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) -class ParametrizeTermsExprImpl(ParametrizeSymbolsExpr): - name: str +class ParametrizeTermsExprImpl(ParametrizeTermsExpr): + # name: str underlying: ExpressionBaseMixin - variables: tuple + # variables: tuple diff --git a/polymatrix/expression/impl/quadraticinexprimpl.py b/polymatrix/expression/impl/quadraticinexprimpl.py index 8c7ac17..da8fa21 100644 --- a/polymatrix/expression/impl/quadraticinexprimpl.py +++ b/polymatrix/expression/impl/quadraticinexprimpl.py @@ -6,4 +6,5 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) class QuadraticInExprImpl(QuadraticInExpr): underlying: ExpressionBaseMixin + monomials: ExpressionBaseMixin variables: tuple diff --git a/polymatrix/expression/impl/setelementatexprimpl.py b/polymatrix/expression/impl/setelementatexprimpl.py new file mode 100644 index 0000000..6d070ac --- /dev/null +++ b/polymatrix/expression/impl/setelementatexprimpl.py @@ -0,0 +1,10 @@ +import dataclass_abc +from polymatrix.expression.setelementatexpr import SetElementAtExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class SetElementAtExprImpl(SetElementAtExpr): + underlying: ExpressionBaseMixin + index: tuple + value: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/linearin2exprimpl.py b/polymatrix/expression/impl/sosmonomialsexprimpl.py index 34b42bf..9e43b23 100644 --- a/polymatrix/expression/impl/linearin2exprimpl.py +++ b/polymatrix/expression/impl/sosmonomialsexprimpl.py @@ -1,9 +1,9 @@ import dataclass_abc -from polymatrix.expression.linearin2expr import LinearIn2Expr +from polymatrix.expression.sosmonomialsexpr import SOSMonomialsExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) -class LinearIn2ExprImpl(LinearIn2Expr): +class SOSMonomialsExprImpl(SOSMonomialsExpr): underlying: ExpressionBaseMixin - expression: ExpressionBaseMixin + variables: tuple diff --git a/polymatrix/expression/impl/substituteexprimpl.py b/polymatrix/expression/impl/substituteexprimpl.py new file mode 100644 index 0000000..1eae940 --- /dev/null +++ b/polymatrix/expression/impl/substituteexprimpl.py @@ -0,0 +1,10 @@ +import dataclass_abc +from polymatrix.expression.substituteexpr import SubstituteExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class SubstituteExprImpl(SubstituteExpr): + underlying: ExpressionBaseMixin + variables: tuple + substitutions: tuple diff --git a/polymatrix/expression/impl/subtractmonomialsexprimpl.py b/polymatrix/expression/impl/subtractmonomialsexprimpl.py new file mode 100644 index 0000000..211f7b1 --- /dev/null +++ b/polymatrix/expression/impl/subtractmonomialsexprimpl.py @@ -0,0 +1,9 @@ +import dataclass_abc +from polymatrix.expression.subtractmonomialsexpr import SubtractMonomialsExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class SubtractMonomialsExprImpl(SubtractMonomialsExpr): + underlying: ExpressionBaseMixin + monomials: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/truncateexprimpl.py b/polymatrix/expression/impl/truncateexprimpl.py new file mode 100644 index 0000000..01c47d0 --- /dev/null +++ b/polymatrix/expression/impl/truncateexprimpl.py @@ -0,0 +1,10 @@ +import dataclass_abc +from polymatrix.expression.truncateexpr import TruncateExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class TruncateExprImpl(TruncateExpr): + underlying: ExpressionBaseMixin + variables: ExpressionBaseMixin + degree: int diff --git a/polymatrix/expression/init/initaccumulateexpr.py b/polymatrix/expression/init/initaccumulateexpr.py deleted file mode 100644 index 30297bf..0000000 --- a/polymatrix/expression/init/initaccumulateexpr.py +++ /dev/null @@ -1,25 +0,0 @@ -import dataclass_abc -import typing - -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.accumulateexpr import AccumulateExpr - - -def init_accumulate_expr( - underlying: ExpressionBaseMixin, - acc_func: typing.Callable, - initial: ExpressionBaseMixin = None, -): - - @dataclass_abc.dataclass_abc(frozen=True) - class AccumulateExprImpl(AccumulateExpr): - underlying: ExpressionBaseMixin - initial: ExpressionBaseMixin - - def acc_func(self, acc, v): - return acc_func(acc, v) - - return AccumulateExprImpl( - underlying=underlying, - initial=initial, - ) diff --git a/polymatrix/expression/init/initcombinationsexpr.py b/polymatrix/expression/init/initcombinationsexpr.py index 7a1df2f..80bf4d0 100644 --- a/polymatrix/expression/init/initcombinationsexpr.py +++ b/polymatrix/expression/init/initcombinationsexpr.py @@ -3,10 +3,10 @@ from polymatrix.expression.impl.combinationsexprimpl import CombinationsExprImpl def init_combinations_expr( - underlying: ExpressionBaseMixin, + monomials: ExpressionBaseMixin, number: int, ): return CombinationsExprImpl( - underlying=underlying, + monomials=monomials, number=number, ) diff --git a/polymatrix/expression/init/initdivergenceexpr.py b/polymatrix/expression/init/initdivergenceexpr.py new file mode 100644 index 0000000..d1aa02c --- /dev/null +++ b/polymatrix/expression/init/initdivergenceexpr.py @@ -0,0 +1,13 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from typing import Union +from polymatrix.expression.impl.divergenceexprimpl import DivergenceExprImpl + + +def init_divergence_expr( + underlying: ExpressionBaseMixin, + variables: tuple, +): + return DivergenceExprImpl( + underlying=underlying, + variables=variables, +) diff --git a/polymatrix/expression/init/initkktexpr.py b/polymatrix/expression/init/initkktexpr.py deleted file mode 100644 index f634d97..0000000 --- a/polymatrix/expression/init/initkktexpr.py +++ /dev/null @@ -1,17 +0,0 @@ -from polymatrix.expression.mixins.expressionmixin import ExpressionMixin -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.impl.kktexprimpl import KKTExprImpl - - -def init_kkt_expr( - cost: ExpressionMixin, - variables: ExpressionBaseMixin, - equality: ExpressionMixin = None, - inequality: ExpressionMixin = None, -): - return KKTExprImpl( - cost=cost, - equality=equality, - inequality=inequality, - variables=variables, -) diff --git a/polymatrix/expression/init/initlinearin2expr.py b/polymatrix/expression/init/initlinearin2expr.py deleted file mode 100644 index 7225510..0000000 --- a/polymatrix/expression/init/initlinearin2expr.py +++ /dev/null @@ -1,12 +0,0 @@ -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.impl.linearin2exprimpl import LinearIn2ExprImpl - - -def init_linear_in2_expr( - underlying: ExpressionBaseMixin, - expression: ExpressionBaseMixin, -): - return LinearIn2ExprImpl( - underlying=underlying, - expression=expression, -) diff --git a/polymatrix/expression/init/initlinearinexpr.py b/polymatrix/expression/init/initlinearinexpr.py index f7f76e4..5cd172c 100644 --- a/polymatrix/expression/init/initlinearinexpr.py +++ b/polymatrix/expression/init/initlinearinexpr.py @@ -4,9 +4,13 @@ from polymatrix.expression.impl.linearinexprimpl import LinearInExprImpl def init_linear_in_expr( underlying: ExpressionBaseMixin, - variables: tuple, + monomials: ExpressionBaseMixin, + variables: ExpressionBaseMixin, + ignore_unmatched: bool = None, ): return LinearInExprImpl( underlying=underlying, + monomials=monomials, variables=variables, + ignore_unmatched = ignore_unmatched, ) diff --git a/polymatrix/expression/init/initlinearinmonomialsinexpr.py b/polymatrix/expression/init/initlinearinmonomialsinexpr.py new file mode 100644 index 0000000..77208a7 --- /dev/null +++ b/polymatrix/expression/init/initlinearinmonomialsinexpr.py @@ -0,0 +1,12 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.linearinmonomialsinexprimpl import LinearInMonomialsInExprImpl + + +def init_linear_in_monomials_in( + underlying: ExpressionBaseMixin, + variables: tuple, +): + return LinearInMonomialsInExprImpl( + underlying=underlying, + variables=variables, +) diff --git a/polymatrix/expression/init/initlinearmonomialsexpr.py b/polymatrix/expression/init/initlinearmonomialsexpr.py new file mode 100644 index 0000000..f116562 --- /dev/null +++ b/polymatrix/expression/init/initlinearmonomialsexpr.py @@ -0,0 +1,12 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.linearmonomialsexprimpl import LinearMonomialsExprImpl + + +def init_linear_monomials_expr( + underlying: ExpressionBaseMixin, + variables: tuple, +): + return LinearMonomialsExprImpl( + underlying=underlying, + variables=variables, +) diff --git a/polymatrix/expression/init/initparametrizeexpr.py b/polymatrix/expression/init/initparametrizeexpr.py new file mode 100644 index 0000000..c5b77a3 --- /dev/null +++ b/polymatrix/expression/init/initparametrizeexpr.py @@ -0,0 +1,15 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.parametrizeexprimpl import ParametrizeExprImpl + + +def init_parametrize_expr( + underlying: ExpressionBaseMixin, + name: str = None, +): + if name is None: + name = 'undefined' + + return ParametrizeExprImpl( + underlying=underlying, + name=name, +) diff --git a/polymatrix/expression/init/initparametrizetermsexpr.py b/polymatrix/expression/init/initparametrizetermsexpr.py index b932825..6e2a446 100644 --- a/polymatrix/expression/init/initparametrizetermsexpr.py +++ b/polymatrix/expression/init/initparametrizetermsexpr.py @@ -3,12 +3,12 @@ from polymatrix.expression.impl.parametrizetermsexprimpl import ParametrizeTerms def init_parametrize_terms_expr( - name: str, + # name: str, underlying: ExpressionBaseMixin, - variables: tuple, + # variables: tuple, ): return ParametrizeTermsExprImpl( - name=name, + # name=name, underlying=underlying, - variables=variables, + # variables=variables, ) diff --git a/polymatrix/expression/init/initquadraticinexpr.py b/polymatrix/expression/init/initquadraticinexpr.py index 1e73745..5aa40a5 100644 --- a/polymatrix/expression/init/initquadraticinexpr.py +++ b/polymatrix/expression/init/initquadraticinexpr.py @@ -4,9 +4,11 @@ from polymatrix.expression.impl.quadraticinexprimpl import QuadraticInExprImpl def init_quadratic_in_expr( underlying: ExpressionBaseMixin, + monomials: ExpressionBaseMixin, variables: tuple, ): return QuadraticInExprImpl( underlying=underlying, + monomials=monomials, variables=variables, ) diff --git a/polymatrix/expression/init/initsetelementatexpr.py b/polymatrix/expression/init/initsetelementatexpr.py new file mode 100644 index 0000000..86df5aa --- /dev/null +++ b/polymatrix/expression/init/initsetelementatexpr.py @@ -0,0 +1,14 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.setelementatexprimpl import SetElementAtExprImpl + + +def init_set_element_at_expr( + underlying: ExpressionBaseMixin, + index: tuple, + value: ExpressionBaseMixin, +): + return SetElementAtExprImpl( + underlying=underlying, + index=index, + value=value, +) diff --git a/polymatrix/expression/init/initsosmonomialsexpr.py b/polymatrix/expression/init/initsosmonomialsexpr.py new file mode 100644 index 0000000..0f5d034 --- /dev/null +++ b/polymatrix/expression/init/initsosmonomialsexpr.py @@ -0,0 +1,12 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.sosmonomialsexprimpl import SOSMonomialsExprImpl + + +def init_sos_monomials_expr( + underlying: ExpressionBaseMixin, + variables: tuple, +): + return SOSMonomialsExprImpl( + underlying=underlying, + variables=variables, +) diff --git a/polymatrix/expression/init/initsubstituteexpr.py b/polymatrix/expression/init/initsubstituteexpr.py new file mode 100644 index 0000000..63e1123 --- /dev/null +++ b/polymatrix/expression/init/initsubstituteexpr.py @@ -0,0 +1,27 @@ +import numpy as np + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.substituteexprimpl import SubstituteExprImpl + + +def init_substitute_expr( + underlying: ExpressionBaseMixin, + variables: tuple, + substitutions: tuple = None, +): + if substitutions is None: + assert isinstance(variables, tuple) + + variables, substitutions = tuple(zip(*variables)) + + elif isinstance(substitutions, np.ndarray): + substitutions = tuple(substitutions.reshape(-1)) + + elif not isinstance(substitutions, tuple): + substitutions = (substitutions,) + + return SubstituteExprImpl( + underlying=underlying, + variables=variables, + substitutions=substitutions, +) diff --git a/polymatrix/expression/init/initsubtractmonomialsexpr.py b/polymatrix/expression/init/initsubtractmonomialsexpr.py new file mode 100644 index 0000000..131401b --- /dev/null +++ b/polymatrix/expression/init/initsubtractmonomialsexpr.py @@ -0,0 +1,12 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.subtractmonomialsexprimpl import SubtractMonomialsExprImpl + + +def init_subtract_monomials_expr( + underlying: ExpressionBaseMixin, + monomials: ExpressionBaseMixin, +): + return SubtractMonomialsExprImpl( + underlying=underlying, + monomials=monomials, +) diff --git a/polymatrix/expression/init/inittruncateexpr.py b/polymatrix/expression/init/inittruncateexpr.py new file mode 100644 index 0000000..818981e --- /dev/null +++ b/polymatrix/expression/init/inittruncateexpr.py @@ -0,0 +1,14 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.truncateexprimpl import TruncateExprImpl + + +def init_truncate_expr( + underlying: ExpressionBaseMixin, + variables: ExpressionBaseMixin, + degree: int, +): + return TruncateExprImpl( + underlying=underlying, + variables=variables, + degree=degree, +) diff --git a/polymatrix/expression/init/initvstackexpr.py b/polymatrix/expression/init/initvstackexpr.py index 1d7834e..70764ef 100644 --- a/polymatrix/expression/init/initvstackexpr.py +++ b/polymatrix/expression/init/initvstackexpr.py @@ -1,9 +1,20 @@ from polymatrix.expression.impl.vstackexprimpl import VStackExprImpl +from polymatrix.expression.init.initfromarrayexpr import init_from_array_expr +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin def init_v_stack_expr( underlying: tuple, ): + + def gen_underlying(): + + for e in underlying: + if isinstance(e, ExpressionBaseMixin): + yield e + else: + yield init_from_array_expr(e) + return VStackExprImpl( - underlying=underlying, + underlying=tuple(gen_underlying()), ) diff --git a/polymatrix/expression/kktexpr.py b/polymatrix/expression/kktexpr.py deleted file mode 100644 index aa5be52..0000000 --- a/polymatrix/expression/kktexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.kktexprmixin import KKTExprMixin - -class KKTExpr(KKTExprMixin): - pass diff --git a/polymatrix/expression/linearin2expr.py b/polymatrix/expression/linearin2expr.py deleted file mode 100644 index 2c40dc8..0000000 --- a/polymatrix/expression/linearin2expr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.linearinexprmixin import LinearInExprMixin - -class LinearIn2Expr(LinearInExprMixin): - pass diff --git a/polymatrix/expression/linearinexpr.py b/polymatrix/expression/linearinexpr.py index 9054937..4edf8b3 100644 --- a/polymatrix/expression/linearinexpr.py +++ b/polymatrix/expression/linearinexpr.py @@ -1,4 +1,4 @@ -from polymatrix.expression.mixins.oldlinearexprmixin import OldLinearExprMixin +from polymatrix.expression.mixins.linearinexprmixin import LinearInExprMixin -class LinearInExpr(OldLinearExprMixin): +class LinearInExpr(LinearInExprMixin): pass diff --git a/polymatrix/expression/linearinmonomialsinexpr.py b/polymatrix/expression/linearinmonomialsinexpr.py new file mode 100644 index 0000000..a78fb40 --- /dev/null +++ b/polymatrix/expression/linearinmonomialsinexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.linearinmutipleexprmixin import LinearInMultipleExprMixin + +class LinearInMonomialsInExpr(LinearInMultipleExprMixin): + pass diff --git a/polymatrix/expression/linearmatrixinexpr.py b/polymatrix/expression/linearmatrixinexpr.py index 296e902..2bce2e7 100644 --- a/polymatrix/expression/linearmatrixinexpr.py +++ b/polymatrix/expression/linearmatrixinexpr.py @@ -1,4 +1,4 @@ -from polymatrix.expression.mixins.filterlinearpartexprmixin import FilterLinearPartExprMixin +from polymatrix.expression.mixins.linearmatrixinexprmixin import LinearMatrixInExprMixin -class LinearMatrixInExpr(FilterLinearPartExprMixin): +class LinearMatrixInExpr(LinearMatrixInExprMixin): pass diff --git a/polymatrix/expression/linearmonomialsexpr.py b/polymatrix/expression/linearmonomialsexpr.py new file mode 100644 index 0000000..9ff6337 --- /dev/null +++ b/polymatrix/expression/linearmonomialsexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.linearmonomialsexprmixin import LinearMonomialsExprMixin + +class LinearMonomialsExpr(LinearMonomialsExprMixin): + pass diff --git a/polymatrix/expression/mixins/addauxequationsexprmixin.py b/polymatrix/expression/mixins/addauxequationsexprmixin.py deleted file mode 100644 index e5353fe..0000000 --- a/polymatrix/expression/mixins/addauxequationsexprmixin.py +++ /dev/null @@ -1,58 +0,0 @@ - -import abc -import dataclasses -import itertools -import dataclass_abc -from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin - -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.expression.expressionstate import ExpressionState - - -# todo: is this really needed? -class AddAuxEquationsExprMixin(ExpressionBaseMixin): - @property - @abc.abstractmethod - def underlying(self) -> ExpressionBaseMixin: - ... - - # overwrites abstract method of `ExpressionBaseMixin` - def apply( - self, - state: ExpressionState, - ) -> tuple[ExpressionState, PolyMatrix]: - - state, underlying = self.underlying.apply(state=state) - - assert underlying.shape[1] == 1 - - @dataclass_abc.dataclass_abc(frozen=True) - class AddAuxEquationsPolyMatrix(PolyMatrixMixin): - underlying: tuple[PolyMatrixMixin] - shape: tuple[int, int] - n_row: int - auxillary_equations: tuple[dict[tuple[int], float]] - - def get_poly(self, row: int, col: int) -> dict[tuple[int, ...], float]: - if row < self.n_row: - return self.underlying.get_poly(row, col) - - elif row < self.shape[0]: - return self.auxillary_equations[row - self.n_row] - - else: - raise Exception(f'row {row} is out of bounds') - - auxillary_equations = tuple(state.auxillary_equations.values()) - - polymat = AddAuxEquationsPolyMatrix( - underlying=underlying, - shape=(underlying.shape[0] + len(auxillary_equations), 1), - n_row=underlying.shape[0], - auxillary_equations=auxillary_equations, - ) - - state = dataclasses.replace(state, auxillary_equations={}) - - return state, polymat diff --git a/polymatrix/expression/mixins/assertdegreeexprmixin.py b/polymatrix/expression/mixins/assertdegreeexprmixin.py new file mode 100644 index 0000000..af0a3e6 --- /dev/null +++ b/polymatrix/expression/mixins/assertdegreeexprmixin.py @@ -0,0 +1,49 @@ + +import abc +import itertools +import dataclass_abc +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState + + +class AssertDegreeExprMixin(ExpressionBaseMixin): + @property + @abc.abstractmethod + def underlying(self) -> tuple[ExpressionBaseMixin, ...]: + ... + + @property + @abc.abstractmethod + def degrees(self) -> tuple[int, ...]: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + + state, underlying = self.underlying.apply(state) + + @dataclass_abc.dataclass_abc(frozen=True) + class AssertDegreePolyMatrix(PolyMatrixMixin): + underlying: PolyMatrixMixin + degrees: tuple[int, ...] + + def get_poly(self, row: int, col: int) -> dict[tuple[int, ...], float]: + terms = self.underlying.get_poly(row, col) + + for monomial in terms.keys(): + pass + + return terms + + polymatrix = AssertDegreePolyMatrix( + underlying=underlying, + degrees=self.degrees, + ) + + return state, polymatrix diff --git a/polymatrix/expression/mixins/blockdiagexprmixin.py b/polymatrix/expression/mixins/blockdiagexprmixin.py index df4bd80..b0006e3 100644 --- a/polymatrix/expression/mixins/blockdiagexprmixin.py +++ b/polymatrix/expression/mixins/blockdiagexprmixin.py @@ -35,11 +35,11 @@ class BlockDiagExprMixin(ExpressionBaseMixin): shape: tuple[int, int] def get_poly(self, row: int, col: int) -> dict[tuple[int, ...], float]: - for polymat, ((row_start, col_start), (row_end, col_end)) in zip(self.all_underlying, self.underlying_row_col_range): + for polymatrix, ((row_start, col_start), (row_end, col_end)) in zip(self.all_underlying, self.underlying_row_col_range): if row_start <= row < row_end: if col_start <= col < col_end: - return polymat.get_poly( + return polymatrix.get_poly( row=row-row_start, col=col-col_start, ) @@ -58,10 +58,10 @@ class BlockDiagExprMixin(ExpressionBaseMixin): shape = underlying_row_col_range[-1][1] - polymat = BlockDiagPolyMatrix( + polymatrix = BlockDiagPolyMatrix( all_underlying=all_underlying, shape=shape, underlying_row_col_range=underlying_row_col_range, ) - return state, polymat + return state, polymatrix diff --git a/polymatrix/expression/mixins/combinationsexprmixin.py b/polymatrix/expression/mixins/combinationsexprmixin.py index 6d0f6eb..824ebf1 100644 --- a/polymatrix/expression/mixins/combinationsexprmixin.py +++ b/polymatrix/expression/mixins/combinationsexprmixin.py @@ -7,12 +7,15 @@ from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getmonomialindices import get_monomial_indices +from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices +from polymatrix.expression.utils.sortmonomials import sort_monomials class CombinationsExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod - def underlying(self) -> ExpressionBaseMixin: + def monomials(self) -> ExpressionBaseMixin: ... @property @@ -25,51 +28,34 @@ class CombinationsExprMixin(ExpressionBaseMixin): self, state: ExpressionState, ) -> tuple[ExpressionState, PolyMatrix]: - state, underlying = self.underlying.apply(state=state) + if self.number == 0: + terms = {(0, 0): {tuple(): 1.0}} - assert underlying.shape[1] == 1 - - def gen_monomials(): - for _, monomial_terms in underlying.get_terms(): - assert len(monomial_terms) == 1, f'{monomial_terms} has more than 1 element' - - for monomial, values in monomial_terms.items(): - yield monomial, values - - monomial_terms = tuple(gen_monomials()) + poly_matrix = init_poly_matrix( + terms=terms, + shape=(1, 1), + ) - assert underlying.shape[0] == len(monomial_terms) + elif self.number == 1: + state, monomials = self.monomials.apply(state=state) + poly_matrix = monomials - combinations = tuple(itertools.combinations_with_replacement(range(underlying.shape[0]), self.number)) + else: - # print(combinations) - # print(math.comb(underlying.shape[0]+self.number-1, self.number)) + state, monomials = get_monomial_indices(state, self.monomials) - terms = {} + combinations = tuple(itertools.combinations_with_replacement(monomials, self.number)) - for row, combination in enumerate(combinations): - def gen_combination_monomials(): - for index in combination: - monomial, _ = monomial_terms[index] - yield from monomial + terms = {} - combination_monomial = tuple(gen_combination_monomials()) + for row, combination in enumerate(combinations): + combination_monomial = merge_monomial_indices(combination) - def acc_combination_value(acc, index): - _, value = monomial_terms[index] - return acc * value + terms[row, 0] = {combination_monomial: 1.0} - *_, combination_value = itertools.accumulate( - combination, - acc_combination_value, - initial=1.0, + poly_matrix = init_poly_matrix( + terms=terms, + shape=(math.comb(len(monomials) + self.number - 1, self.number), 1), ) - terms[row, 0] = {combination_monomial: combination_value} - - poly_matrix = init_poly_matrix( - terms=terms, - shape=(math.comb(underlying.shape[0]+self.number-1, self.number), 1), - ) - return state, poly_matrix diff --git a/polymatrix/expression/mixins/derivativeexprmixin.py b/polymatrix/expression/mixins/derivativeexprmixin.py index b403873..4fc48a5 100644 --- a/polymatrix/expression/mixins/derivativeexprmixin.py +++ b/polymatrix/expression/mixins/derivativeexprmixin.py @@ -10,6 +10,7 @@ from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getderivativemonomials import get_derivative_monomials from polymatrix.expression.utils.getvariableindices import get_variable_indices @@ -36,112 +37,119 @@ class DerivativeExprMixin(ExpressionBaseMixin): ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) + + assert underlying.shape[1] == 1, f'{underlying.shape=}' state, diff_wrt_variables = get_variable_indices(state, self.variables) - def get_derivative_terms( - monomial_terms, - diff_wrt_variable: int, - state: ExpressionState, - considered_variables: set, - ): + # def get_derivative_terms( + # monomial_terms, + # diff_wrt_variable: int, + # state: ExpressionState, + # considered_variables: set, + # ): + + # if self.introduce_derivatives: - if self.introduce_derivatives: + # raise Exception('not implemented') - def gen_new_variables(): - for monomial in monomial_terms.keys(): - for var in monomial: - if var not in diff_wrt_variables and var not in considered_variables: - yield var + # def gen_new_variables(): + # for monomial in monomial_terms.keys(): + # for var in monomial: + # if var not in diff_wrt_variables and var not in considered_variables: + # yield var - new_variables = set(gen_new_variables()) + # new_variables = set(gen_new_variables()) - new_considered_variables = considered_variables | new_variables + # new_considered_variables = considered_variables | new_variables - def acc_state_candidates(acc, new_variable): - state, candidates = acc + # def acc_state_candidates(acc, new_variable): + # state, candidates = acc - key = init_derivative_key( - variable=new_variable, - with_respect_to=diff_wrt_variable, - ) - state = state.register(key=key, n_param=1) + # key = init_derivative_key( + # variable=new_variable, + # with_respect_to=diff_wrt_variable, + # ) + # state = state.register(key=key, n_param=1) - state, auxillary_derivation_terms = get_derivative_terms( - monomial_terms=state.auxillary_equations[new_variable], - diff_wrt_variable=diff_wrt_variable, - state=state, - considered_variables=new_considered_variables, - ) + # state, auxillary_derivation_terms = get_derivative_terms( + # monomial_terms=state.auxillary_equations[new_variable], + # diff_wrt_variable=diff_wrt_variable, + # state=state, + # considered_variables=new_considered_variables, + # ) - if 1 < len(auxillary_derivation_terms): - derivation_variable = state.offset_dict[key][0] + # if 1 < len(auxillary_derivation_terms): + # derivation_variable = state.offset_dict[key][0] - state = dataclasses.replace( - state, - auxillary_equations=state.auxillary_equations | {derivation_variable: auxillary_derivation_terms}, - ) + # state = dataclasses.replace( + # state, + # auxillary_equations=state.auxillary_equations | {derivation_variable: auxillary_derivation_terms}, + # ) - return state, candidates + (new_variable,) + # return state, candidates + (new_variable,) - else: - return state, candidates + # else: + # return state, candidates - *_, (state, confirmed_variables) = itertools.accumulate( - new_variables, - acc_state_candidates, - initial=(state, tuple()), - ) + # *_, (state, confirmed_variables) = itertools.accumulate( + # new_variables, + # acc_state_candidates, + # initial=(state, tuple()), + # ) - else: - confirmed_variables = tuple() + # else: + # confirmed_variables = tuple() - derivation_terms = collections.defaultdict(float) + # derivation_terms = collections.defaultdict(float) - for monomial, value in monomial_terms.items(): + # for monomial, value in monomial_terms.items(): - # count powers for each variable - monomial_cnt = dict(collections.Counter(monomial)) + # # # count powers for each variable + # # monomial_cnt = dict(collections.Counter(monomial)) + # monomial_cnt = dict(monomial) - def differentiate_monomial(dependent_variable, derivation_variable=None): - def gen_diff_monomial(): - for current_variable, current_count in monomial_cnt.items(): + # def differentiate_monomial(dependent_variable, derivation_variable=None): + # def gen_diff_monomial(): + # for current_variable, current_count in monomial: - if current_variable is dependent_variable: - sel_counter = current_count - 1 + # if current_variable is dependent_variable: + # sel_counter = current_count - 1 - else: - sel_counter = current_count + # else: + # sel_counter = current_count - for _ in range(sel_counter): - yield current_variable + # # for _ in range(sel_counter): + # # yield current_variable + # yield current_variable, sel_counter - if derivation_variable is not None: - yield derivation_variable + # if derivation_variable is not None: + # yield derivation_variable - diff_monomial = tuple(sorted(gen_diff_monomial())) + # diff_monomial = tuple(sorted(gen_diff_monomial())) - return diff_monomial, value * monomial_cnt[dependent_variable] + # return diff_monomial, value * monomial_cnt[dependent_variable] - if diff_wrt_variable in monomial_cnt: - diff_monomial, value = differentiate_monomial(diff_wrt_variable) - derivation_terms[diff_monomial] += value + # if diff_wrt_variable in monomial_cnt: + # diff_monomial, value = differentiate_monomial(diff_wrt_variable) + # derivation_terms[diff_monomial] += value - for candidate_variable in monomial_cnt.keys(): - if candidate_variable in considered_variables or candidate_variable in confirmed_variables: - key = init_derivative_key( - variable=candidate_variable, - with_respect_to=diff_wrt_variable, - ) - derivation_variable = state.offset_dict[key][0] + # # only used if introduce_derivatives == True + # for candidate_variable in monomial_cnt.keys(): + # if candidate_variable in considered_variables or candidate_variable in confirmed_variables: + # key = init_derivative_key( + # variable=candidate_variable, + # with_respect_to=diff_wrt_variable, + # ) + # derivation_variable = state.offset_dict[key][0] - diff_monomial, value = differentiate_monomial( - dependent_variable=candidate_variable, - derivation_variable=derivation_variable, - ) - derivation_terms[diff_monomial] += value + # diff_monomial, value = differentiate_monomial( + # dependent_variable=candidate_variable, + # derivation_variable=derivation_variable, + # ) + # derivation_terms[diff_monomial] += value - return state, dict(derivation_terms) + # return state, dict(derivation_terms) terms = {} @@ -155,11 +163,12 @@ class DerivativeExprMixin(ExpressionBaseMixin): # derivate each variable and map result to the corresponding column for col, diff_wrt_variable in enumerate(diff_wrt_variables): - state, derivation_terms = get_derivative_terms( + state, derivation_terms = get_derivative_monomials( monomial_terms=underlying_terms, diff_wrt_variable=diff_wrt_variable, state=state, considered_variables=set(), + introduce_derivatives=self.introduce_derivatives, ) if 0 < len(derivation_terms): @@ -171,299 +180,3 @@ class DerivativeExprMixin(ExpressionBaseMixin): ) return state, poly_matrix - - # def get_derivative_terms( - # monomial_terms, - # diff_wrt_variable: int, - # state: PolyMatrixExprState, - # considered_variables: set, - # # implement_derivation: bool, - # ): - # derivation_terms = collections.defaultdict(float) - - # other_independent_variables = tuple(var for var in diff_wrt_variables if var is not diff_wrt_variable) - - # # print(other_independent_variables) - # # print(tuple(variable for monomial in monomial_terms.keys() for variable in monomial)) - - # if sum(variable not in other_independent_variables for monomial in monomial_terms.keys() for variable in monomial) < 2: - # return {}, state - - # # if not implement_derivation: - # # implement_derivation = any(diff_wrt_variable in monomial for monomial in monomial_terms.keys()) - - # for monomial, value in monomial_terms.items(): - - # # count powers for each variable - # monomial_cnt = dict(collections.Counter(monomial)) - - # def differentiate_monomial(dependent_variable, derivation_variable=None): - # def gen_diff_monomial(): - # for current_variable, current_count in monomial_cnt.items(): - - # if current_variable is dependent_variable: - # sel_counter = current_count - 1 - - # else: - # sel_counter = current_count - - # for _ in range(sel_counter): - # yield current_variable - - # if derivation_variable is not None: - # yield derivation_variable - - # diff_monomial = tuple(gen_diff_monomial()) - - # return diff_monomial, value * monomial_cnt[dependent_variable] - - # if diff_wrt_variable in monomial_cnt: - # diff_monomial, value = differentiate_monomial(diff_wrt_variable) - # derivation_terms[diff_monomial] += value - - # if self.introduce_derivatives: - - # def gen_derivation_keys(): - # for variable in monomial_cnt.keys(): - # if variable not in diff_wrt_variables: - # yield variable - - # candidate_variables = tuple(gen_derivation_keys()) - - # new_considered_derivations = considered_variables | set(candidate_variables) - - # for candidate_variable in candidate_variables: - - # # introduce new auxillary equation - # if candidate_variable not in considered_variables: - # auxillary_derivation_terms, state = get_derivative_terms( - # monomial_terms=state.auxillary_equations[candidate_variable], - # diff_wrt_variable=diff_wrt_variable, - # state=state, - # considered_variables=new_considered_derivations, - # # implement_derivation=implement_derivation, - # ) - - # if 0 < len(auxillary_derivation_terms): - # key = init_derivative_key( - # variable=candidate_variable, - # with_respect_to=diff_wrt_variable, - # ) - # state = state.register(key=key, n_param=1) - # derivation_variable = state.offset_dict[key][0] - - # state = dataclasses.replace( - # state, - # auxillary_equations=state.auxillary_equations | {derivation_variable: auxillary_derivation_terms}, - # ) - - # else: - - # key = init_derivative_key( - # variable=candidate_variable, - # with_respect_to=diff_wrt_variable, - # ) - # state = state.register(key=key, n_param=1) - # derivation_variable = state.offset_dict[key][0] - - # diff_monomial, value = differentiate_monomial( - # dependent_variable=candidate_variable, - # derivation_variable=derivation_variable, - # ) - # derivation_terms[diff_monomial] += value - - # return dict(derivation_terms), state - - # terms = {} - - # for row in range(self.shape[0]): - - # try: - # underlying_terms = underlying.get_poly(row, 0) - # except KeyError: - # continue - - # # derivate each variable and map result to the corresponding column - # for col, diff_wrt_variable in enumerate(diff_wrt_variables): - - # derivation_terms, state = get_derivative_terms( - # monomial_terms=underlying_terms, - # diff_wrt_variable=diff_wrt_variable, - # state=state, - # considered_variables=set(), - # # implement_derivation=False, - # ) - - # if 0 < len(derivation_terms): - # terms[row, col] = derivation_terms - - # poly_matrix = init_poly_matrix( - # terms=terms, - # shape=self.shape, - # ) - - # return state, poly_matrix - - - - # state = [state] - - # state, underlying = self.underlying.apply(state=state) - - # match self.variables: - # case ExpressionBaseMixin(): - # assert self.variables.shape[1] == 1 - - # state, dependent_variables = self.variables.apply(state) - - # def gen_indices(): - # for row in range(dependent_variables.shape[0]): - # for monomial in dependent_variables.get_poly(row, 0).keys(): - # yield monomial[0] - - # variable_indices = tuple(sorted(gen_indices())) - - # case _: - # def gen_indices(): - # for variable in self.variables: - # if variable in state.offset_dict: - # yield state.offset_dict[variable][0] - - # variable_indices = tuple(sorted(gen_indices())) - - # terms = {} - # derivations_keys = set() - - # # derivate each variable and map result to the corresponding column - # for col, derivation_variable in enumerate(variable_indices): - - # def get_derivative_terms(monomial_terms): - - # terms_row_col = collections.defaultdict(float) - - # for monomial, value in monomial_terms.items(): - - # # count powers for each variable - # monomial_cnt = dict(collections.Counter(monomial)) - - # variable_candidates = tuple() - - # if derivation_variable in monomial_cnt: - # variable_candidates += ((derivation_variable, None),) - - # if self.introduce_derivatives: - # def gen_dependent_variables(): - # for dependent_variable in monomial_cnt.keys(): - # if dependent_variable not in variable_indices: - # derivation_key = init_derivative_key( - # variable=dependent_variable, - # with_respect_to=derivation_variable, - # ) - # derivations_keys.add(derivation_key) - # state = state.register(key=derivation_key, n_param=1) - # yield dependent_variable, derivation_key - - # variable_candidates += tuple(gen_dependent_variables()) - - # for variable_candidate, derivation_key in variable_candidates: - - # def generate_monomial(): - # for current_variable, current_count in monomial_cnt.items(): - - # if current_variable is variable_candidate: - # sel_counter = current_count - 1 - - # else: - # sel_counter = current_count - - # for _ in range(sel_counter): - # yield current_variable - - # if derivation_key is not None: - # yield state.offset_dict[derivation_key][0] - - # col_monomial = tuple(generate_monomial()) - - # terms_row_col[col_monomial] += value * monomial_cnt[variable_candidate] - - # return dict(terms_row_col) - - # for row in range(self.shape[0]): - - # try: - # underlying_terms = underlying.get_poly(row, 0) - # except KeyError: - # continue - - # derivative_terms = get_derivative_terms(underlying_terms) - - # if 0 < len(derivative_terms): - # terms[row, col] = derivative_terms - - # derivation_variables = collections.defaultdict(list) - # for derivation_key in derivations_keys: - # derivation_variables[derivation_key.with_respect_to].append(derivation_key) - - # aux_der_terms = [] - - # for derivation_variable, derivation_keys in derivation_variables.items(): - - # dependent_variables = tuple(derivation_key.variable for derivation_key in derivation_keys) - - - # for aux_terms in state.auxillary_equations: - - # # only intoduce a new auxillary equation if there is a monomial containing at least one dependent variable - # if any(variable in dependent_variables for monomial in aux_terms.keys() for variable in monomial): - - # terms_row_col = collections.defaultdict(float) - - # # for each monomial - # for aux_monomial, value in aux_terms.items(): - - # # count powers for each variable - # monomial_cnt = dict(collections.Counter(aux_monomial)) - - # variable_candidates = tuple() - - # if derivation_variable in monomial_cnt: - # variable_candidates += ((derivation_variable, None),) - - # # add dependent variables - # variable_candidates += tuple((derivation_key.variable, derivation_key) for derivation_key in derivation_keys if derivation_key.variable in monomial_cnt) - - # for variable_candidate, derivative_key in variable_candidates: - - # def generate_monomial(): - # for current_variable, current_count in monomial_cnt.items(): - - # if current_variable is variable_candidate: - # sel_counter = current_count - 1 - - # else: - # sel_counter = current_count - - # for _ in range(sel_counter): - # yield current_variable - - # if derivative_key is not None: - # yield state.offset_dict[derivative_key][0] - - # col_monomial = tuple(generate_monomial()) - - # terms_row_col[col_monomial] += value * monomial_cnt[variable_candidate] - - # if 0 < len(terms_row_col): - # aux_der_terms.append(dict(terms_row_col)) - - # state = dataclasses.replace( - # state, - # auxillary_equations=state.auxillary_equations + tuple(aux_der_terms), - # ) - - # poly_matrix = init_poly_matrix( - # terms=terms, - # shape=self.shape, - # ) - - # return state, poly_matrix diff --git a/polymatrix/expression/mixins/derivativekeymixin.py b/polymatrix/expression/mixins/derivativekeymixin.py index e0adddc..db4951a 100644 --- a/polymatrix/expression/mixins/derivativekeymixin.py +++ b/polymatrix/expression/mixins/derivativekeymixin.py @@ -1,5 +1,6 @@ import abc +# move this class class DerivativeKeyMixin(abc.ABC): @property @abc.abstractmethod diff --git a/polymatrix/expression/mixins/determinantexprmixin.py b/polymatrix/expression/mixins/determinantexprmixin.py index f002f0b..8d41e3d 100644 --- a/polymatrix/expression/mixins/determinantexprmixin.py +++ b/polymatrix/expression/mixins/determinantexprmixin.py @@ -26,6 +26,8 @@ class DeterminantExprMixin(ExpressionBaseMixin): self, state: ExpressionState, ) -> tuple[ExpressionState, PolyMatrix]: + raise Exception('not implemented') + if self in state.cache: return state, state.cache[self] diff --git a/polymatrix/expression/mixins/divergenceexprmixin.py b/polymatrix/expression/mixins/divergenceexprmixin.py new file mode 100644 index 0000000..b523b2d --- /dev/null +++ b/polymatrix/expression/mixins/divergenceexprmixin.py @@ -0,0 +1,65 @@ + +import abc +import collections +import dataclasses +import itertools +import typing +from polymatrix.expression.init.initderivativekey import init_derivative_key + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getderivativemonomials import get_derivative_monomials +from polymatrix.expression.utils.getvariableindices import get_variable_indices + + +class DivergenceExprMixin(ExpressionBaseMixin): + @property + @abc.abstractmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> typing.Union[tuple, ExpressionBaseMixin]: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + + state, underlying = self.underlying.apply(state=state) + state, variables = get_variable_indices(state, self.variables) + + assert underlying.shape[1] == 1, f'{underlying.shape=}' + assert len(variables) == underlying.shape[0], f'{variables=}, {underlying.shape=}' + + monomial_terms = collections.defaultdict(float) + + for row, variable in enumerate(variables): + + try: + underlying_terms = underlying.get_poly(row, 0) + except KeyError: + continue + + state, derivation_terms = get_derivative_monomials( + monomial_terms=underlying_terms, + diff_wrt_variable=variable, + state=state, + considered_variables=set(), + introduce_derivatives=False, + ) + + for monomial, value in derivation_terms.items(): + monomial_terms[monomial] += value + + poly_matrix = init_poly_matrix( + terms={(0, 0): dict(monomial_terms)}, + shape=(1, 1), + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/divisionexprmixin.py b/polymatrix/expression/mixins/divisionexprmixin.py index dc5f616..7551cda 100644 --- a/polymatrix/expression/mixins/divisionexprmixin.py +++ b/polymatrix/expression/mixins/divisionexprmixin.py @@ -19,17 +19,13 @@ class DivisionExprMixin(ExpressionBaseMixin): def right(self) -> ExpressionBaseMixin: ... - # # overwrites abstract method of `ExpressionBaseMixin` - # @property - # def shape(self) -> tuple[int, int]: - # return self.left.shape - # overwrites abstract method of `ExpressionBaseMixin` def apply( self, state: ExpressionState, ) -> tuple[ExpressionState, PolyMatrix]: + # add an auxillary equation and, therefore, needs to be cached if self in state.cache: return state, state.cache[self] diff --git a/polymatrix/expression/mixins/elemmultexprmixin.py b/polymatrix/expression/mixins/elemmultexprmixin.py index 931e022..7e0bc8c 100644 --- a/polymatrix/expression/mixins/elemmultexprmixin.py +++ b/polymatrix/expression/mixins/elemmultexprmixin.py @@ -1,11 +1,15 @@ import abc import itertools +import typing +import dataclass_abc from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices class ElemMultExprMixin(ExpressionBaseMixin): @@ -32,9 +36,59 @@ class ElemMultExprMixin(ExpressionBaseMixin): state, left = self.left.apply(state=state) state, right = self.right.apply(state=state) - assert right.shape == (1, 1) + if right.shape == (1, 1): + right_poly = right.get_poly(0, 0) - right_terms = right.get_poly(0, 0) + @dataclass_abc.dataclass_abc(frozen=True) + class BroadCastedPolyMatrix(PolyMatrixMixin): + underlying: tuple[tuple[int], float] + shape: tuple[int, int] + + def get_poly(self, row: int, col: int) -> typing.Optional[dict[tuple[int, ...], float]]: + return self.underlying + + right = BroadCastedPolyMatrix( + underlying=right_poly, + shape=left.shape, + ) + + # if right.shape == (1, 1): + + # right_terms = right.get_poly(0, 0) + + # terms = {} + + # for poly_row in range(left.shape[0]): + # for poly_col in range(left.shape[1]): + + # terms_row_col = {} + + # try: + # left_terms = left.get_poly(poly_row, poly_col) + # except KeyError: + # continue + + # for (left_monomial, left_value), (right_monomial, right_value) \ + # in itertools.product(left_terms.items(), right_terms.items()): + + # value = left_value * right_value + + # # if value == 0: + # # continue + + # # monomial = tuple(sorted(left_monomial + right_monomial)) + + # new_monomial = merge_monomial_indices((left_monomial, right_monomial)) + + # if new_monomial not in terms_row_col: + # terms_row_col[new_monomial] = 0 + + # terms_row_col[new_monomial] += value + + # if 0 < len(terms_row_col): + # terms[poly_row, poly_col] = terms_row_col + + # elif left.shape == right.shape: terms = {} @@ -48,20 +102,27 @@ class ElemMultExprMixin(ExpressionBaseMixin): except KeyError: continue + try: + right_terms = right.get_poly(poly_row, poly_col) + except KeyError: + continue + for (left_monomial, left_value), (right_monomial, right_value) \ in itertools.product(left_terms.items(), right_terms.items()): value = left_value * right_value - if value == 0: - continue + # if value == 0: + # continue + + # monomial = tuple(sorted(left_monomial + right_monomial)) - monomial = tuple(sorted(left_monomial + right_monomial)) + new_monomial = merge_monomial_indices((left_monomial, right_monomial)) - if monomial not in terms_row_col: - terms_row_col[monomial] = 0 + if new_monomial not in terms_row_col: + terms_row_col[new_monomial] = 0 - terms_row_col[monomial] += value + terms_row_col[new_monomial] += value if 0 < len(terms_row_col): terms[poly_row, poly_col] = terms_row_col diff --git a/polymatrix/expression/mixins/evalexprmixin.py b/polymatrix/expression/mixins/evalexprmixin.py index 3089bbb..7fa149e 100644 --- a/polymatrix/expression/mixins/evalexprmixin.py +++ b/polymatrix/expression/mixins/evalexprmixin.py @@ -1,6 +1,7 @@ import abc import itertools +import math from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @@ -56,16 +57,17 @@ class EvalExprMixin(ExpressionBaseMixin): for monomial, value in underlying_terms.items(): - def acc_monomial(acc, variable): - monomial, value = acc + def acc_monomial(acc, next): + new_monomial, value = acc + variable, count = next if variable in variable_indices: index = variable_indices.index(variable) - new_value = value * values[index] - return monomial, new_value + new_value = value * values[index] * count + return new_monomial, new_value else: - return monomial + (variable,), value + return new_monomial + (next,), value *_, (new_monomial, new_value) = tuple(itertools.accumulate( monomial, @@ -79,10 +81,11 @@ class EvalExprMixin(ExpressionBaseMixin): terms_row_col[new_monomial] += new_value # delete zero entries - if terms_row_col[new_monomial] == 0: + if math.isclose(terms_row_col[new_monomial], 0, abs_tol=1e-12): del terms_row_col[new_monomial] - terms[row, col] = terms_row_col + if 0 < len(terms_row_col): + terms[row, col] = terms_row_col poly_matrix = init_poly_matrix( terms=terms, diff --git a/polymatrix/expression/mixins/expressionmixin.py b/polymatrix/expression/mixins/expressionmixin.py index f8b5ef7..e75973e 100644 --- a/polymatrix/expression/mixins/expressionmixin.py +++ b/polymatrix/expression/mixins/expressionmixin.py @@ -3,31 +3,36 @@ import dataclasses import typing import numpy as np from sympy import re -from polymatrix.expression.init.initaccumulateexpr import init_accumulate_expr + from polymatrix.expression.init.initadditionexpr import init_addition_expr from polymatrix.expression.init.initcacheexpr import init_cache_expr from polymatrix.expression.init.initcombinationsexpr import init_combinations_expr from polymatrix.expression.init.initderivativeexpr import init_derivative_expr from polymatrix.expression.init.initdeterminantexpr import init_determinant_expr +from polymatrix.expression.init.initdivergenceexpr import init_divergence_expr from polymatrix.expression.init.initdivisionexpr import init_division_expr from polymatrix.expression.init.initelemmultexpr import init_elem_mult_expr from polymatrix.expression.init.initevalexpr import init_eval_expr -from polymatrix.expression.init.initlinearin2expr import init_linear_in2_expr from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr +from polymatrix.expression.init.initlinearinmonomialsinexpr import init_linear_in_monomials_in from polymatrix.expression.init.initfromarrayexpr import init_from_array_expr from polymatrix.expression.init.initgetitemexpr import init_get_item_expr from polymatrix.expression.init.initlinearmatrixinexpr import init_linear_matrix_in_expr +from polymatrix.expression.init.initlinearmonomialsexpr import init_linear_monomials_expr from polymatrix.expression.init.initmatrixmultexpr import init_matrix_mult_expr -from polymatrix.expression.init.initparametrizetermsexpr import init_parametrize_terms_expr +from polymatrix.expression.init.initparametrizeexpr import init_parametrize_expr from polymatrix.expression.init.initquadraticinexpr import init_quadratic_in_expr from polymatrix.expression.init.initrepmatexpr import init_rep_mat_expr from polymatrix.expression.init.initreshapeexpr import init_reshape_expr -from polymatrix.expression.init.initsqueezeexpr import init_squeeze_expr +from polymatrix.expression.init.initsetelementatexpr import init_set_element_at_expr +from polymatrix.expression.init.initsosmonomialsexpr import init_sos_monomials_expr +from polymatrix.expression.init.initsubtractmonomialsexpr import init_subtract_monomials_expr from polymatrix.expression.init.initsumexpr import init_sum_expr from polymatrix.expression.init.initsymmetricexpr import init_symmetric_expr from polymatrix.expression.init.inittoconstantexpr import init_to_constant_expr from polymatrix.expression.init.inittoquadraticexpr import init_to_quadratic_expr from polymatrix.expression.init.inittransposeexpr import init_transpose_expr +from polymatrix.expression.init.inittruncateexpr import init_truncate_expr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix @@ -48,15 +53,6 @@ class ExpressionMixin( def apply(self, state: ExpressionState) -> tuple[ExpressionState, PolyMatrix]: return self.underlying.apply(state) - # # overwrites abstract method of `PolyMatrixExprBaseMixin` - # @property - # def degrees(self) -> set[int]: - # return self.underlying.degrees - - # def __iter__(self): - # for row in range(self.shape[0]): - # yield self[row, 0] - def __add__(self, other: ExpressionBaseMixin) -> 'ExpressionMixin': if other is None: return self @@ -181,7 +177,7 @@ class ExpressionMixin( return dataclasses.replace( self, underlying=init_combinations_expr( - underlying=self.underlying, + monomials=self.underlying, number=number, ), ) @@ -208,6 +204,18 @@ class ExpressionMixin( ), ) + def divergence( + self, + variables: tuple, + ) -> 'ExpressionMixin': + return dataclasses.replace( + self, + underlying=init_divergence_expr( + underlying=self.underlying, + variables=variables, + ), + ) + def eval( self, variable: tuple, @@ -222,7 +230,7 @@ class ExpressionMixin( ), ) - def filter_linear_part(self, variable) -> 'ExpressionMixin': + def linear_matrix_in(self, variable: 'ExpressionMixin') -> 'ExpressionMixin': return dataclasses.replace( self, underlying=init_linear_matrix_in_expr( @@ -231,30 +239,73 @@ class ExpressionMixin( ), ) - def linear_in(self, variables: tuple) -> 'ExpressionMixin': + def linear_monomials(self, variables: 'ExpressionMixin') -> 'ExpressionMixin': return dataclasses.replace( self, - underlying=init_linear_in2_expr( + underlying=init_linear_monomials_expr( underlying=self.underlying, - expression=variables, + variables=variables, ), ) - def parametrize(self, name: str, variables: tuple) -> 'ExpressionMixin': + def linear_in( + self, + variables: 'ExpressionMixin', + monomials: 'ExpressionMixin' = None, + ignore_unmatched: bool = None, + ) -> 'ExpressionMixin': + if monomials is None: + monomials = self.linear_monomials(variables) + return dataclasses.replace( self, - underlying=init_parametrize_terms_expr( - name=name, + underlying=init_linear_in_expr( underlying=self.underlying, + monomials=monomials, variables=variables, + ignore_unmatched=ignore_unmatched, ), ) - def quadratic_in(self, variables: tuple) -> 'ExpressionMixin': + def expressions_in(self, variables: tuple): + return dataclasses.replace( + self, + underlying=init_linear_in_monomials_in( + underlying=self.underlying, + variables=variables, + ), + ) + + def parametrize(self, name: str = None) -> 'ExpressionMixin': + return dataclasses.replace( + self, + underlying=init_parametrize_expr( + underlying=self.underlying, + name=name, + ), + ) + + def quadratic_in(self, variables: 'ExpressionMixin', monomials: 'ExpressionMixin' = None) -> 'ExpressionMixin': + if monomials is None: + monomials = self.quadratic_monomials(variables) + return dataclasses.replace( self, underlying=init_quadratic_in_expr( underlying=self.underlying, + monomials=monomials, + variables=variables, + ), + ) + + def quadratic_monomials( + self, + variables: tuple, + ) -> 'ExpressionMixin': + return dataclasses.replace( + self, + underlying=init_sos_monomials_expr( + underlying=self.underlying, variables=variables, ), ) @@ -277,11 +328,57 @@ class ExpressionMixin( ), ) - def squeeze(self): + # def squeeze(self): + # return dataclasses.replace( + # self, + # underlying=init_squeeze_expr( + # underlying=self.underlying, + # ), + # ) + + def set_element_at( + self, + row: int, + col: int, + value: 'ExpressionMixin', + ) -> 'ExpressionMixin': + if isinstance(value, ExpressionMixin): + value = value.underlying + else: + value = init_from_array_expr(value) + + return dataclasses.replace( + self, + underlying=init_set_element_at_expr( + underlying=self.underlying, + index=(row, col), + value=value, + ), + ) + + def subtract_monomials( + self, + monomials: 'ExpressionMixin', + ) -> 'ExpressionMixin': return dataclasses.replace( self, - underlying=init_squeeze_expr( + underlying=init_subtract_monomials_expr( underlying=self.underlying, + monomials=monomials, + ), + ) + + def substitute( + self, + variable: tuple, + substitutions: tuple['ExpressionMixin', ...] = None, + ) -> 'ExpressionMixin': + return dataclasses.replace( + self, + underlying=init_eval_expr( + underlying=self.underlying, + variables=variable, + substitutions=substitutions, ), ) @@ -316,3 +413,13 @@ class ExpressionMixin( underlying=self.underlying, ), ) + + def truncate(self, variables: tuple, degree: int): + return dataclasses.replace( + self, + underlying=init_truncate_expr( + underlying=self.underlying, + variables=variables, + degree=degree, + ), + ) diff --git a/polymatrix/expression/mixins/fromarrayexprmixin.py b/polymatrix/expression/mixins/fromarrayexprmixin.py index 179cec4..2b066d4 100644 --- a/polymatrix/expression/mixins/fromarrayexprmixin.py +++ b/polymatrix/expression/mixins/fromarrayexprmixin.py @@ -1,5 +1,6 @@ import abc +import math from numpy import poly import sympy @@ -29,7 +30,9 @@ class FromArrayExprMixin(ExpressionBaseMixin): try: poly = sympy.poly(poly_data) except sympy.polys.polyerrors.GeneratorsNeeded: - terms[poly_row, poly_col] = {tuple(): poly_data} + + if not math.isclose(poly_data, 0): + terms[poly_row, poly_col] = {tuple(): poly_data} continue for symbol in poly.gens: @@ -41,17 +44,15 @@ class FromArrayExprMixin(ExpressionBaseMixin): # a5 x1 x3**2 -> c=a5, m_cnt=(1, 0, 2) for value, monomial_count in zip(poly.coeffs(), poly.monoms()): - if value == 0.0: + if math.isclose(value, 0): continue # m_cnt=(1, 0, 2) -> m=(0, 2, 2) def gen_monomial(): for rel_idx, p in enumerate(monomial_count): - - idx, _ = state.offset_dict[poly.gens[rel_idx]] - - for _ in range(p): - yield idx + if 0 < p: + idx, _ = state.offset_dict[poly.gens[rel_idx]] + yield idx, p monomial = tuple(gen_monomial()) diff --git a/polymatrix/expression/mixins/kktexprmixin.py b/polymatrix/expression/mixins/kktexprmixin.py deleted file mode 100644 index d4a3d97..0000000 --- a/polymatrix/expression/mixins/kktexprmixin.py +++ /dev/null @@ -1,191 +0,0 @@ - -import abc -import itertools -import typing -import dataclass_abc -from polymatrix.expression.derivativekey import DerivativeKey - -from polymatrix.expression.init.initpolymatrix import init_poly_matrix -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.mixins.expressionmixin import ExpressionMixin -from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.expression.expressionstate import ExpressionState - - -# todo: remove this? -class KKTExprMixin(ExpressionBaseMixin): - @property - @abc.abstractmethod - def cost(self) -> ExpressionMixin: - ... - - @property - @abc.abstractmethod - def equality(self) -> ExpressionMixin: - ... - - @property - @abc.abstractmethod - def equality(self) -> typing.Optional[ExpressionMixin]: - ... - - @property - @abc.abstractmethod - def inequality(self) -> typing.Optional[ExpressionMixin]: - ... - - @property - @abc.abstractmethod - def variables(self) -> ExpressionBaseMixin: - ... - - # overwrites abstract method of `ExpressionBaseMixin` - def apply( - self, - state: ExpressionState, - ) -> tuple[ExpressionState, PolyMatrix]: - - state, cost = self.cost.apply(state=state) - - assert cost.shape[1] == 1 - - if self.equality is not None: - - state, equality = self.equality.apply(state=state) - - state, equality_der = self.equality.diff( - self.variables, - introduce_derivatives=True, - ).apply(state) - - assert cost.shape[0] == equality_der.shape[1] - - def acc_nu_variables(acc, v): - state, nu_variables = acc - - nu_variable = state.n_param - state = state.register(n_param=1) - - return state, nu_variables + [nu_variable] - - *_, (state, nu_variables) = tuple(itertools.accumulate( - range(equality.shape[0]), - acc_nu_variables, - initial=(state, []), - )) - - else: - nu_variables = tuple() - - if self.inequality is not None: - - state, inequality = self.inequality.apply(state=state) - - state, inequality_der = self.inequality.diff( - self.variables, - introduce_derivatives=True, - ).apply(state) - - assert cost.shape[0] == inequality_der.shape[1] - - def acc_lambda_variables(acc, v): - state, lambda_variables = acc - - lambda_variable = state.n_param - state = state.register(n_param=3) - - return state, lambda_variables + [lambda_variable] - - *_, (state, lambda_variables) = tuple(itertools.accumulate( - range(inequality.shape[0]), - acc_lambda_variables, - initial=(state, []), - )) - - else: - lambda_variables = tuple() - - idx_start = 0 - - terms = {} - - for row in range(cost.shape[0]): - try: - monomial_terms = cost.get_poly(row, 0) - except KeyError: - monomial_terms = {} - - for eq_idx, nu_variable in enumerate(nu_variables): - - try: - underlying_terms = equality_der.get_poly(eq_idx, row) - except KeyError: - continue - - for monomial, value in underlying_terms.items(): - new_monomial = monomial + (nu_variable,) - - if new_monomial not in monomial_terms: - monomial_terms[new_monomial] = 0 - - monomial_terms[new_monomial] += value - - for inequality_idx, lambda_variable in enumerate(lambda_variables): - try: - underlying_terms = inequality_der.get_poly(inequality_idx, row) - except KeyError: - continue - - for monomial, value in underlying_terms.items(): - new_monomial = monomial + (lambda_variable,) - - if new_monomial not in monomial_terms: - monomial_terms[new_monomial] = 0 - - monomial_terms[new_monomial] += value - - terms[idx_start, 0] = monomial_terms - idx_start += 1 - - # equality constraints - # -------------------- - - for eq_idx, nu_variable in enumerate(nu_variables): - try: - underlying_terms = equality.get_poly(eq_idx, 0) - except KeyError: - continue - - terms[idx_start, 0] = underlying_terms - idx_start += 1 - - # inequality constraints - # -------------------- - - for inequality_idx, lambda_variable in enumerate(lambda_variables): - r_lambda = lambda_variable + 1 - r_inequality = lambda_variable + 2 - - try: - underlying_terms = inequality.get_poly(inequality_idx, 0) - except KeyError: - continue - - # f(x) <= -0.01 - terms[idx_start, 0] = underlying_terms | {(r_inequality, r_inequality): 1} - idx_start += 1 - - # dual feasibility, lambda >= 0 - terms[idx_start, 0] = {(lambda_variable,): 1, (r_lambda, r_lambda): -1} - idx_start += 1 - - # complementary slackness - terms[idx_start, 0] = {(r_lambda, r_inequality): 1} - idx_start += 1 - - poly_matrix = init_poly_matrix( - terms=terms, - shape=(idx_start, 1), - ) - - return state, poly_matrix diff --git a/polymatrix/expression/mixins/linearinexprmixin.py b/polymatrix/expression/mixins/linearinexprmixin.py index b690891..369467a 100644 --- a/polymatrix/expression/mixins/linearinexprmixin.py +++ b/polymatrix/expression/mixins/linearinexprmixin.py @@ -7,6 +7,7 @@ from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getmonomialindices import get_monomial_indices from polymatrix.expression.utils.getvariableindices import get_variable_indices @@ -18,7 +19,17 @@ class LinearInExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod - def expression(self) -> ExpressionBaseMixin: + def monomials(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> tuple: + ... + + @property + @abc.abstractmethod + def ignore_unmatched(self) -> bool: ... # overwrites abstract method of `ExpressionBaseMixin` @@ -27,59 +38,38 @@ class LinearInExprMixin(ExpressionBaseMixin): state: ExpressionState, ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) + state, monomials = get_monomial_indices(state, self.monomials) + state, variable_indices = get_variable_indices(state, self.variables) - state, variables = self.expression.apply(state=state) - - def gen_variable_monomials(): - for _, term in variables.get_terms(): - assert len(term) == 1, f'{term} should have only a single monomial' - - for monomial in term.keys(): - yield monomial - - variable_monomials = tuple(gen_variable_monomials()) - all_variables = set(variable for monomial in variable_monomials for variable in monomial) - - # print(variables) - # print(variable_monomials) + assert underlying.shape[1] == 1 terms = collections.defaultdict(dict) for row in range(underlying.shape[0]): - for col in range(underlying.shape[1]): - - try: - underlying_terms = underlying.get_poly(row, col) - except KeyError: - continue - - for monomial, value in underlying_terms.items(): - for variable_index, variable_monomial in enumerate(variable_monomials): + try: + underlying_terms = underlying.get_poly(row, 0) + except KeyError: + continue - remainder = list(monomial) + for monomial, value in underlying_terms.items(): - try: - for variable in variable_monomial: - remainder.remove(variable) + x_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx in variable_indices) + p_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx not in variable_indices) - except ValueError: - continue - - # print(remainder) - # print(variable_monomial) - # print(all(variable not in remainder for variable in all_variables)) - - # take the first that matches - if all(variable not in remainder for variable in all_variables): - - terms[row, variable_index][tuple(remainder)] = value - - break + try: + col = monomials.index(x_monomial) + except ValueError: + if self.ignore_unmatched: + continue + else: + raise Exception(f'{x_monomial} not in {monomials}') + + terms[row, col][p_monomial] = value poly_matrix = init_poly_matrix( terms=dict(terms), - shape=(underlying.shape[0], len(variable_monomials)), + shape=(underlying.shape[0], len(monomials)), ) return state, poly_matrix diff --git a/polymatrix/expression/mixins/linearinmutipleexprmixin.py b/polymatrix/expression/mixins/linearinmutipleexprmixin.py new file mode 100644 index 0000000..7d4f00e --- /dev/null +++ b/polymatrix/expression/mixins/linearinmutipleexprmixin.py @@ -0,0 +1,144 @@ + +import abc +import collections + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getvariableindices import get_variable_indices + + +class LinearInMultipleExprMixin(ExpressionBaseMixin): + @property + @abc.abstractmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> ExpressionBaseMixin: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + state, underlying = self.underlying.apply(state=state) + + state, variables = self.variables.apply(state=state) + + def gen_variable_monomials(): + for _, term in variables.get_terms(): + assert len(term) == 1, f'{term} should have only a single monomial' + + for monomial in term.keys(): + yield monomial + + variable_monomials = tuple(gen_variable_monomials()) + all_variables = set(variable for monomial in variable_monomials for variable in monomial) + + terms = {} + idx_row = 0 + + for row in range(underlying.shape[0]): + for col in range(underlying.shape[1]): + + try: + underlying_terms = underlying.get_poly(row, col) + except KeyError: + continue + + x_monomial_terms = collections.defaultdict(lambda: collections.defaultdict(float)) + + for monomial, value in underlying_terms.items(): + + assert tuple(sorted(monomial)) == monomial, f'{monomial} is not sorted' + + x_monomial = tuple((var_idx for var_idx in monomial if var_idx in all_variables)) + p_monomial = tuple(var_idx for var_idx in monomial if var_idx not in all_variables) + + x_monomial_terms[x_monomial][p_monomial] += value + + for data in x_monomial_terms.values(): + terms[idx_row, 0] = dict(data) + idx_row += 1 + + poly_matrix = init_poly_matrix( + terms=terms, + shape=(idx_row, 1), + ) + + return state, poly_matrix + + # state, variable_indices = get_variable_indices(state, self.variables) + + # underlying_terms = underlying.get_terms() + + # def gen_variable_terms(): + # for _, monomial_term in underlying_terms: + # for monomial in monomial_term.keys(): + + # x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) + # yield x_monomial + + # variable_terms = tuple(set(gen_variable_terms())) + + # terms = {} + + # for (row, _), monomial_term in underlying.get_terms(): + + # x_monomial_terms = collections.defaultdict(lambda: collections.defaultdict(float)) + + # for monomial, value in monomial_term.items(): + + # x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) + # p_monomial = tuple(var_idx for var_idx in monomial if var_idx not in variable_indices) + + # assert tuple(sorted(x_monomial)) == x_monomial, f'{x_monomial} is not sorted' + + # x_monomial_terms[x_monomial][p_monomial] += value + + # for x_monomial, data in x_monomial_terms.items(): + # terms[row, variable_terms.index(x_monomial)] = dict(data) + + # poly_matrix = init_poly_matrix( + # terms=terms, + # shape=(underlying.shape[0], len(variable_terms)), + # ) + + # return state, poly_matrix + + # terms = {} + # idx_row = 0 + + # for row in range(underlying.shape[0]): + # for col in range(underlying.shape[1]): + + # try: + # underlying_terms = underlying.get_poly(row, col) + # except KeyError: + # continue + + # x_monomial_terms = collections.defaultdict(lambda: collections.defaultdict(float)) + + # for monomial, value in underlying_terms.items(): + + # x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) + # p_monomial = tuple(var_idx for var_idx in monomial if var_idx not in variable_indices) + + # assert tuple(sorted(x_monomial)) == x_monomial, f'{x_monomial} is not sorted' + + # x_monomial_terms[x_monomial][p_monomial] += value + + # for data in x_monomial_terms.values(): + # terms[idx_row, 0] = dict(data) + # idx_row += 1 + + # poly_matrix = init_poly_matrix( + # terms=terms, + # shape=(idx_row, 1), + # ) + + # return state, poly_matrix diff --git a/polymatrix/expression/mixins/linearmatrixinexprmixin.py b/polymatrix/expression/mixins/linearmatrixinexprmixin.py new file mode 100644 index 0000000..956b01a --- /dev/null +++ b/polymatrix/expression/mixins/linearmatrixinexprmixin.py @@ -0,0 +1,61 @@ + +import abc +import collections +from numpy import var + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getvariableindices import get_variable_indices + + +class LinearMatrixInExprMixin(ExpressionBaseMixin): + @property + @abc.abstractmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> tuple: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + state, underlying = self.underlying.apply(state=state) + state, variable_index = get_variable_indices(state, variables=self.variable) + + assert len(variable_index) == 1 + + terms = collections.defaultdict(dict) + + for row in range(underlying.shape[0]): + for col in range(underlying.shape[1]): + + try: + underlying_terms = underlying.get_poly(row, col) + except KeyError: + continue + + for monomial, value in underlying_terms.items(): + + if len(monomial) == 1: + + variable, _ = monomial[0] + + if variable == variable_index: + + terms[row, col][tuple()] = value + + break + + poly_matrix = init_poly_matrix( + terms=dict(terms), + shape=underlying.shape, + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/linearmonomialsexprmixin.py b/polymatrix/expression/mixins/linearmonomialsexprmixin.py new file mode 100644 index 0000000..42a5b95 --- /dev/null +++ b/polymatrix/expression/mixins/linearmonomialsexprmixin.py @@ -0,0 +1,67 @@ + +import abc +import dataclasses + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin +from polymatrix.expression.mixins.quadraticinexprmixin import QuadraticInExprMixin +from polymatrix.expression.utils.getvariableindices import get_variable_indices +from polymatrix.expression.utils.sortmonomials import sort_monomials +from polymatrix.expression.utils.splitmonomialindices import split_monomial_indices + + +class LinearMonomialsExprMixin(ExpressionBaseMixin): + @property + @abc.abstractclassmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> tuple: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionStateMixin, + ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: + + state, underlying = self.underlying.apply(state=state) + state, variable_indices = get_variable_indices(state, self.variables) + + def gen_linear_monomials(): + for row in range(underlying.shape[0]): + for col in range(underlying.shape[1]): + + try: + polynomial = underlying.get_poly(row, col) + except KeyError: + continue + + for monomial in polynomial.keys(): + x_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx in variable_indices) + + yield x_monomial + + linear_monomials = sort_monomials(set(gen_linear_monomials())) + + def gen_terms(): + for index, monomial in enumerate(linear_monomials): + yield (index, 0), {monomial: 1.0} + + terms = dict(gen_terms()) + + poly_matrix = init_poly_matrix( + terms=terms, + shape=(len(linear_monomials), 1), + ) + + state = dataclasses.replace( + state, + cache=state.cache | {self: poly_matrix}, + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/matrixmultexprmixin.py b/polymatrix/expression/mixins/matrixmultexprmixin.py index 1f09fdc..029f700 100644 --- a/polymatrix/expression/mixins/matrixmultexprmixin.py +++ b/polymatrix/expression/mixins/matrixmultexprmixin.py @@ -1,11 +1,14 @@ import abc import itertools +import math from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices +from polymatrix.expression.utils.multiplypolynomial import multiply_polynomial class MatrixMultExprMixin(ExpressionBaseMixin): @@ -44,23 +47,27 @@ class MatrixMultExprMixin(ExpressionBaseMixin): except KeyError: continue - for (left_monomial, left_value), (right_monomial, right_value) \ - in itertools.product(left_terms.items(), right_terms.items()): + # for (left_monomial, left_value), (right_monomial, right_value) \ + # in itertools.product(left_terms.items(), right_terms.items()): - value = left_value * right_value + # value = left_value * right_value - if value == 0: - continue + # if value == 0: + # continue - monomial = tuple(sorted(left_monomial + right_monomial)) + # # monomial = tuple(sorted(left_monomial + right_monomial)) + # monomial = merge_monomial_indices((left_monomial, right_monomial)) - if monomial not in terms_row_col: - terms_row_col[monomial] = 0 + # if monomial not in terms_row_col: + # terms_row_col[monomial] = 0 - terms_row_col[monomial] += value + # terms_row_col[monomial] += value + + multiply_polynomial(left_terms, right_terms, terms_row_col) if 0 < len(terms_row_col): terms[poly_row, poly_col] = terms_row_col + # terms[poly_row, poly_col] = {key: val for key, val in terms_row_col.items() if not math.isclose(val, 0, abs_tol=1e-12)} poly_matrix = init_poly_matrix( terms=terms, diff --git a/polymatrix/expression/mixins/oldlinearexprmixin.py b/polymatrix/expression/mixins/oldlinearexprmixin.py deleted file mode 100644 index 1bdebd9..0000000 --- a/polymatrix/expression/mixins/oldlinearexprmixin.py +++ /dev/null @@ -1,131 +0,0 @@ - -import abc -import collections -from numpy import var - -from polymatrix.expression.init.initpolymatrix import init_poly_matrix -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.expression.expressionstate import ExpressionState -from polymatrix.expression.utils.getvariableindices import get_variable_indices - - -class OldLinearExprMixin(ExpressionBaseMixin): - @property - @abc.abstractmethod - def underlying(self) -> ExpressionBaseMixin: - ... - - @property - @abc.abstractmethod - def variables(self) -> ExpressionBaseMixin: - ... - - # overwrites abstract method of `ExpressionBaseMixin` - def apply( - self, - state: ExpressionState, - ) -> tuple[ExpressionState, PolyMatrix]: - state, underlying = self.underlying.apply(state=state) - - state, variable_indices = get_variable_indices(state, self.variables) - - underlying_terms = underlying.get_terms() - - def gen_variable_terms(): - for _, monomial_term in underlying_terms: - for monomial in monomial_term.keys(): - - x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) - yield x_monomial - - variable_terms = tuple(set(gen_variable_terms())) - - terms = {} - - for (row, _), monomial_term in underlying_terms: - - x_monomial_terms = collections.defaultdict(lambda: collections.defaultdict(float)) - - for monomial, value in monomial_term.items(): - - x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) - p_monomial = tuple(var_idx for var_idx in monomial if var_idx not in variable_indices) - - assert tuple(sorted(x_monomial)) == x_monomial, f'{x_monomial} is not sorted' - - x_monomial_terms[x_monomial][p_monomial] += value - - for x_monomial, data in x_monomial_terms.items(): - terms[row, variable_terms.index(x_monomial)] = dict(data) - - poly_matrix = init_poly_matrix( - terms=terms, - shape=(underlying.shape[0], len(variable_terms)), - ) - - return state, poly_matrix - - - # for row in range(underlying.shape[0]): - # for col in range(underlying.shape[1]): - - # try: - # underlying_terms = underlying.get_poly(row, col) - # except KeyError: - # continue - - # x_monomial_terms = collections.defaultdict(lambda: collections.defaultdict(float)) - - # for monomial, value in underlying_terms.items(): - - # x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) - # p_monomial = tuple(var_idx for var_idx in monomial if var_idx not in variable_indices) - - # assert tuple(sorted(x_monomial)) == x_monomial, f'{x_monomial} is not sorted' - - # x_monomial_terms[x_monomial][p_monomial] += value - - # for data in x_monomial_terms.values(): - # terms[idx_row, 0] = dict(data) - # idx_row += 1 - - # poly_matrix = init_poly_matrix( - # terms=terms, - # shape=(idx_row, 1), - # ) - - # return state, poly_matrix - - # terms = {} - # idx_row = 0 - - # for row in range(underlying.shape[0]): - # for col in range(underlying.shape[1]): - - # try: - # underlying_terms = underlying.get_poly(row, col) - # except KeyError: - # continue - - # x_monomial_terms = collections.defaultdict(lambda: collections.defaultdict(float)) - - # for monomial, value in underlying_terms.items(): - - # x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) - # p_monomial = tuple(var_idx for var_idx in monomial if var_idx not in variable_indices) - - # assert tuple(sorted(x_monomial)) == x_monomial, f'{x_monomial} is not sorted' - - # x_monomial_terms[x_monomial][p_monomial] += value - - # for data in x_monomial_terms.values(): - # terms[idx_row, 0] = dict(data) - # idx_row += 1 - - # poly_matrix = init_poly_matrix( - # terms=terms, - # shape=(idx_row, 1), - # ) - - # return state, poly_matrix diff --git a/polymatrix/expression/mixins/parametrizeexprmixin.py b/polymatrix/expression/mixins/parametrizeexprmixin.py new file mode 100644 index 0000000..11a767f --- /dev/null +++ b/polymatrix/expression/mixins/parametrizeexprmixin.py @@ -0,0 +1,58 @@ + +import abc +import dataclasses + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin + + +class ParametrizeExprMixin(ExpressionBaseMixin): + @property + @abc.abstractclassmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractclassmethod + def name(self) -> str: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionStateMixin, + ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: + + # cache polymatrix to not re-parametrize at every apply call + if self in state.cache: + return state, state.cache[self] + + state, underlying = self.underlying.apply(state) + + assert underlying.shape[1] == 1 + + start_index = state.n_param + terms = {} + + for row in range(underlying.shape[0]): + var_index = start_index + row + terms[row, 0] = {((var_index, 1),): 1.0} + + state = state.register( + key=self, + n_param=underlying.shape[0], + ) + + poly_matrix = init_poly_matrix( + terms=terms, + shape=underlying.shape, + ) + + state = dataclasses.replace( + state, + cache=state.cache | {self: poly_matrix}, + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/parametrizetermsexprmixin.py b/polymatrix/expression/mixins/parametrizetermsexprmixin.py index d678db5..51244d4 100644 --- a/polymatrix/expression/mixins/parametrizetermsexprmixin.py +++ b/polymatrix/expression/mixins/parametrizetermsexprmixin.py @@ -7,23 +7,24 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin from polymatrix.expression.utils.getvariableindices import get_variable_indices +from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices class ParametrizeTermsExprMixin(ExpressionBaseMixin): - @property - @abc.abstractclassmethod - def name(self) -> str: - ... + # @property + # @abc.abstractclassmethod + # def name(self) -> str: + # ... @property @abc.abstractclassmethod def underlying(self) -> ExpressionBaseMixin: ... - @property - @abc.abstractclassmethod - def variables(self) -> tuple: - ... + # @property + # @abc.abstractclassmethod + # def variables(self) -> tuple: + # ... # overwrites abstract method of `ExpressionBaseMixin` def apply( @@ -31,17 +32,13 @@ class ParametrizeTermsExprMixin(ExpressionBaseMixin): state: ExpressionStateMixin, ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: + # cache polymatrix to not re-parametrize at every apply call if self in state.cache: return state, state.cache[self] state, underlying = self.underlying.apply(state) - # for variable in self.variables: - # state = state.register(key=variable, n_param=1) - - # variable_indices = tuple(state.offset_dict[variable][0] for variable in self.variables if variable in state.offset_dict) - - state, variable_indices = get_variable_indices(state, self.variables) + # state, variable_indices = get_variable_indices(state, self.variables) start_index = state.n_param terms = {} @@ -54,22 +51,33 @@ class ParametrizeTermsExprMixin(ExpressionBaseMixin): except KeyError: continue - def gen_x_monomial_terms(): - for monomial, value in underlying_terms.items(): - x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) - yield monomial, x_monomial, value + # def gen_x_monomial_terms(): + # for monomial, value in underlying_terms.items(): + # x_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx in variable_indices) + # yield monomial, x_monomial, value - x_monomial_terms = tuple(gen_x_monomial_terms()) + # x_monomial_terms = tuple(gen_x_monomial_terms()) - collected_terms = tuple(sorted(set((x_monomial for _, x_monomial, _ in x_monomial_terms)))) + # # use sort to make the parametrization predictable + # # x_monomials = (0, 1), then the parametrization is (2, 3) + # collected_terms = tuple(sorted( + # set((x_monomial for _, x_monomial, _ in x_monomial_terms)), + # key=lambda m: (sum(count for _, count in m), m)), + # ) + collected_terms = tuple(underlying_terms.keys()) + + # print(tuple(x_monomial for _, x_monomial, _ in x_monomial_terms)) + # print(x_monomial_terms) terms_row_col = {} - for monomial, x_monomial, value in x_monomial_terms: + for monomial, value in underlying_terms.items(): + + param_index = start_index + collected_terms.index(monomial) - new_monomial = monomial + (start_index + collected_terms.index(x_monomial),) + monomial_with_param = merge_monomial_indices(monomial, ((param_index, 1),)) - terms_row_col[new_monomial] = value + terms_row_col[monomial_with_param] = value terms[row, col] = terms_row_col diff --git a/polymatrix/expression/mixins/quadraticinexprmixin.py b/polymatrix/expression/mixins/quadraticinexprmixin.py new file mode 100644 index 0000000..dbfb76d --- /dev/null +++ b/polymatrix/expression/mixins/quadraticinexprmixin.py @@ -0,0 +1,67 @@ + +import abc +import collections + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getmonomialindices import get_monomial_indices +from polymatrix.expression.utils.getvariableindices import get_variable_indices +from polymatrix.expression.utils.splitmonomialindices import split_monomial_indices + + +class QuadraticInExprMixin(ExpressionBaseMixin): + @property + @abc.abstractmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def monomials(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> tuple: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + state, underlying = self.underlying.apply(state=state) + state, sos_monomials = get_monomial_indices(state, self.monomials) + state, variable_indices = get_variable_indices(state, self.variables) + + assert underlying.shape == (1, 1), f'underlying shape is {underlying.shape}' + + underlying_terms = underlying.get_poly(0, 0) + + terms = collections.defaultdict(dict) + + for monomial, value in underlying_terms.items(): + + x_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx in variable_indices) + p_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx not in variable_indices) + + left, right = split_monomial_indices(x_monomial) + + row = sos_monomials.index(left) + col = sos_monomials.index(right) + + monomial_terms = terms[row, col] + + if p_monomial not in monomial_terms: + monomial_terms[p_monomial] = 0 + + monomial_terms[p_monomial] += value + + poly_matrix = init_poly_matrix( + terms=dict(terms), + shape=2*(len(sos_monomials),), + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/quadraticmonomialsexprmixin.py b/polymatrix/expression/mixins/quadraticmonomialsexprmixin.py new file mode 100644 index 0000000..b953d17 --- /dev/null +++ b/polymatrix/expression/mixins/quadraticmonomialsexprmixin.py @@ -0,0 +1,69 @@ + +import abc +import dataclasses + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin +from polymatrix.expression.mixins.quadraticinexprmixin import QuadraticInExprMixin +from polymatrix.expression.utils.getvariableindices import get_variable_indices +from polymatrix.expression.utils.splitmonomialindices import split_monomial_indices + + +class QuadraticMonomialsExprMixin(ExpressionBaseMixin): + @property + @abc.abstractclassmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> tuple: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionStateMixin, + ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: + + state, underlying = self.underlying.apply(state=state) + state, variable_indices = get_variable_indices(state, self.variables) + + def gen_sos_monomials(): + for row in range(underlying.shape[0]): + for col in range(underlying.shape[1]): + + try: + polynomial = underlying.get_poly(row, col) + except KeyError: + continue + + for monomial in polynomial.keys(): + x_monomial = tuple((var_idx, count) for var_idx, count in monomial if var_idx in variable_indices) + + left, right = split_monomial_indices(x_monomial) + + yield left + yield right + + sos_monomials = tuple(sorted(set(gen_sos_monomials()), key=lambda m: (len(m), m))) + + def gen_terms(): + for index, monomial in enumerate(sos_monomials): + yield (index, 0), {monomial: 1.0} + + terms = dict(gen_terms()) + + poly_matrix = init_poly_matrix( + terms=terms, + shape=(len(sos_monomials), 1), + ) + + state = dataclasses.replace( + state, + cache=state.cache | {self: poly_matrix}, + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/quadraticreprexprmixin.py b/polymatrix/expression/mixins/quadraticreprexprmixin.py deleted file mode 100644 index 29ef950..0000000 --- a/polymatrix/expression/mixins/quadraticreprexprmixin.py +++ /dev/null @@ -1,87 +0,0 @@ - -import abc -import collections - -from polymatrix.expression.init.initpolymatrix import init_poly_matrix -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.expression.expressionstate import ExpressionState -from polymatrix.expression.utils.getvariableindices import get_variable_indices - - -class QuadraticReprExprMixin(ExpressionBaseMixin): - @property - @abc.abstractmethod - def underlying(self) -> ExpressionBaseMixin: - ... - - @property - @abc.abstractmethod - def variables(self) -> tuple: - ... - - # overwrites abstract method of `ExpressionBaseMixin` - def apply( - self, - state: ExpressionState, - ) -> tuple[ExpressionState, PolyMatrix]: - state, underlying = self.underlying.apply(state=state) - - assert underlying.shape == (1, 1) - - state, variable_indices = get_variable_indices(state, self.variables) - - # state, variables = self.expression.apply(state=state) - - def gen_variable_monomials(): - _, monomials = underlying.get_terms()[0] - - for monomial, value in monomials.items(): - x_monomial = tuple(variable for variable in monomial if variable in variable_indices) - p_monomial = tuple(variable for variable in monomial if variable not in variable_indices) - - assert len(x_monomial) % 2 == 0, f'length of {x_monomial} should be a multiple of 2' - - x_monomial_sorted = tuple(sorted(x_monomial)) - - idx_half = int(len(x_monomial)/2) - - yield x_monomial_sorted[:idx_half], x_monomial_sorted[idx_half:], p_monomial, value - - underlying_monomials = tuple(gen_variable_monomials()) - - def gen_sos_monomials(): - for left_monomial, right_monomial, _, _ in underlying_monomials: - yield left_monomial - yield right_monomial - - sos_monomials = tuple(sorted(set(gen_sos_monomials()), key=lambda m: (len(m), m))) - - # print(f'{underlying_monomials=}') - # print(f'{sos_monomials=}') - - terms = collections.defaultdict(dict) - - for left_monomial, right_monomial, p_monomial, value in underlying_monomials: - - col = sos_monomials.index(left_monomial) - row = sos_monomials.index(right_monomial) - key = (row, col) - - # print(key) - - monomial_terms = terms[key] - - if p_monomial not in monomial_terms: - monomial_terms[p_monomial] = 0 - - monomial_terms[p_monomial] += value - - # print(terms) - - poly_matrix = init_poly_matrix( - terms=dict(terms), - shape=2*(len(sos_monomials),), - ) - - return state, poly_matrix diff --git a/polymatrix/expression/mixins/setelementatexprmixin.py b/polymatrix/expression/mixins/setelementatexprmixin.py new file mode 100644 index 0000000..5714631 --- /dev/null +++ b/polymatrix/expression/mixins/setelementatexprmixin.py @@ -0,0 +1,69 @@ + +import abc +import dataclasses +import typing +import dataclass_abc +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState + + +class SetElementAtExprMixin(ExpressionBaseMixin): + @dataclasses.dataclass(frozen=True) + class Slice: + start: int + step: int + stop: int + + @property + @abc.abstractmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def index(self) -> tuple[tuple[int, ...], tuple[int, ...]]: + ... + + @property + @abc.abstractmethod + def value(self) -> ExpressionBaseMixin: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + state, underlying = self.underlying.apply(state=state) + state, value = self.value.apply(state=state) + + assert value.shape == (1, 1) + + @dataclass_abc.dataclass_abc(frozen=True) + class SetElementAtPolyMatrix(PolyMatrixMixin): + underlying: PolyMatrixMixin + shape: tuple[int, int] + index: tuple[int, int] + value: dict + + # @property + # def shape(self) -> tuple[int, int]: + # return self.shape + + def get_poly(self, row: int, col: int) -> dict[tuple[int, ...], float]: + if (row, col) == self.index: + return self.value + else: + return self.underlying.get_poly(row, col) + + return state, SetElementAtPolyMatrix( + underlying=underlying, + index=self.index, + shape=underlying.shape, + value=value.get_poly(0, 0) + ) +
\ No newline at end of file diff --git a/polymatrix/expression/mixins/squeezeexprmixin.py b/polymatrix/expression/mixins/squeezeexprmixin.py index 47f2fcf..62fe7de 100644 --- a/polymatrix/expression/mixins/squeezeexprmixin.py +++ b/polymatrix/expression/mixins/squeezeexprmixin.py @@ -10,6 +10,7 @@ from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +# remove? class SqueezeExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod diff --git a/polymatrix/expression/mixins/substituteexprmixin.py b/polymatrix/expression/mixins/substituteexprmixin.py new file mode 100644 index 0000000..fc71f73 --- /dev/null +++ b/polymatrix/expression/mixins/substituteexprmixin.py @@ -0,0 +1,107 @@ + +import abc +import collections +import itertools +import math + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.polymatrix import PolyMatrix +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getvariableindices import get_variable_indices +from polymatrix.expression.utils.multiplypolynomial import multiply_polynomial + + +class SubstituteExprMixin(ExpressionBaseMixin): + @property + @abc.abstractmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> tuple: + ... + + @property + @abc.abstractmethod + def substitutions(self) -> tuple[ExpressionBaseMixin, ...]: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: + state, underlying = self.underlying.apply(state=state) + state, variable_indices = get_variable_indices(state, self.variables) + + def acc_substitutions(acc, substitution_expr): + state, result = acc + + # for expr in self.expressions: + if isinstance(substitution_expr, ExpressionBaseMixin): + state, substitution = substitution_expr.apply(state) + + assert substitution.shape == (1, 1), f'{substitution=}' + + polynomial = substitution.get_poly(0, 0) + + elif isinstance(substitution_expr, int) or isinstance(substitution_expr, float): + polynomial = {tuple(): substitution_expr} + + return state, result + (polynomial,) + + *_, (state, substitutions) = tuple(itertools.accumulate( + self.substitutions, + acc_substitutions, + initial=(state, tuple()), + )) + + terms = {} + + for row in range(underlying.shape[0]): + for col in range(underlying.shape[1]): + + try: + underlying_terms = underlying.get_poly(row, col) + except KeyError: + continue + + terms_row_col = collections.defaultdict(float) + + for monomial, value in underlying_terms.items(): + + terms_monomial = {tuple(): value} + + for variable, count in monomial: + if variable in variable_indices: + index = variable_indices.index(variable) + substitution = substitutions[index] + + + for _ in range(count): + + next = {} + multiply_polynomial(terms_monomial, substitution, next) + terms_monomial = next + + else: + next = {} + multiply_polynomial(terms_monomial, {((variable, count),): 1.0}, next) + terms_monomial = next + + for monomial, value in terms_monomial.items(): + terms_row_col[monomial] += value + + terms_row_col = {key: val for key, val in terms_row_col.items() if not math.isclose(val, 0, abs_tol=1e-12)} + + if 0 < len(terms_row_col): + terms[row, col] = terms_row_col + + poly_matrix = init_poly_matrix( + terms=terms, + shape=underlying.shape, + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/subtractmonomialsexprmixin.py b/polymatrix/expression/mixins/subtractmonomialsexprmixin.py new file mode 100644 index 0000000..45ef444 --- /dev/null +++ b/polymatrix/expression/mixins/subtractmonomialsexprmixin.py @@ -0,0 +1,60 @@ + +import abc +import dataclasses + +from polymatrix.expression.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin +from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin +from polymatrix.expression.mixins.quadraticinexprmixin import QuadraticInExprMixin +from polymatrix.expression.utils.getmonomialindices import get_monomial_indices +from polymatrix.expression.utils.getvariableindices import get_variable_indices +from polymatrix.expression.utils.sortmonomials import sort_monomials +from polymatrix.expression.utils.subtractmonomialindices import subtract_monomial_indices + + +class SubtractMonomialsExprMixin(ExpressionBaseMixin): + @property + @abc.abstractclassmethod + def underlying(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def monomials(self) -> ExpressionBaseMixin: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionStateMixin, + ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: + + state, underlying = get_monomial_indices(state, self.underlying) + state, sub_monomials = get_monomial_indices(state, self.monomials) + + def gen_remainders(): + for m1 in underlying: + for m2 in sub_monomials: + try: + remainder = subtract_monomial_indices(m1, m2) + except KeyError: + continue + + yield remainder + + remainders = sort_monomials(set(gen_remainders())) + + terms = {(row, 0): {remainder: 1.0} for row, remainder in enumerate(remainders)} + + poly_matrix = init_poly_matrix( + terms=terms, + shape=(len(remainders), 1), + ) + + state = dataclasses.replace( + state, + cache=state.cache | {self: poly_matrix}, + ) + + return state, poly_matrix diff --git a/polymatrix/expression/mixins/sumexprmixin.py b/polymatrix/expression/mixins/sumexprmixin.py index 96bf3da..b3afd51 100644 --- a/polymatrix/expression/mixins/sumexprmixin.py +++ b/polymatrix/expression/mixins/sumexprmixin.py @@ -35,11 +35,11 @@ class SumExprMixin(ExpressionBaseMixin): term_monomials = terms[row, 0] - for monomial, coeff in underlying_terms.items(): + for monomial, value in underlying_terms.items(): if monomial in term_monomials: - term_monomials[monomial] += coeff + term_monomials[monomial] += value else: - term_monomials[monomial] = coeff + term_monomials[monomial] = value poly_matrix = init_poly_matrix( terms=terms, diff --git a/polymatrix/expression/mixins/symmetricexprmixin.py b/polymatrix/expression/mixins/symmetricexprmixin.py index e914895..da07715 100644 --- a/polymatrix/expression/mixins/symmetricexprmixin.py +++ b/polymatrix/expression/mixins/symmetricexprmixin.py @@ -1,5 +1,6 @@ import abc +import collections import itertools import dataclass_abc from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin @@ -50,19 +51,17 @@ class SymmetricExprMixin(ExpressionBaseMixin): raise KeyError() else: - terms = {} + terms = collections.defaultdict(float) # merge monomials for monomials in all_monomials: for monomial, value in monomials.items(): - if monomial in terms: - terms[monomial] = (terms[monomial] + value) / 2 - else: - terms[monomial] = value - return terms + terms[monomial] += value / 2 - polymat = SymmetricPolyMatrix( + return dict(terms) + + polymatrix = SymmetricPolyMatrix( underlying=underlying, ) - return state, polymat + return state, polymatrix diff --git a/polymatrix/expression/mixins/toquadraticexprmixin.py b/polymatrix/expression/mixins/toquadraticexprmixin.py index 2b62a89..02e366d 100644 --- a/polymatrix/expression/mixins/toquadraticexprmixin.py +++ b/polymatrix/expression/mixins/toquadraticexprmixin.py @@ -9,6 +9,7 @@ from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +# to be deleted? class ToQuadraticExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod diff --git a/polymatrix/expression/mixins/accumulateexprmixin.py b/polymatrix/expression/mixins/truncateexprmixin.py index 47df335..7d0927a 100644 --- a/polymatrix/expression/mixins/accumulateexprmixin.py +++ b/polymatrix/expression/mixins/truncateexprmixin.py @@ -1,27 +1,29 @@ import abc -import typing +import collections +from numpy import var from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getvariableindices import get_variable_indices -# todo: is this needed? -class AccumulateExprMixin(ExpressionBaseMixin): +class TruncateExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod def underlying(self) -> ExpressionBaseMixin: ... + @property @abc.abstractmethod - def acc_func(self) -> typing.Callable: + def variables(self) -> ExpressionBaseMixin: ... @property @abc.abstractmethod - def initial(self) -> ExpressionBaseMixin: + def degree(self) -> int: ... # overwrites abstract method of `ExpressionBaseMixin` @@ -31,27 +33,32 @@ class AccumulateExprMixin(ExpressionBaseMixin): ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) - acc = self.initial + state, variable_indices = get_variable_indices(state, self.variables) - for row in range(underlying.shape[0]): + terms = {} - col_terms = {} - + for row in range(underlying.shape[0]): for col in range(underlying.shape[1]): + try: underlying_terms = underlying.get_poly(row, col) except KeyError: continue - col_terms[0, col] = underlying_terms + terms_row_col = {} - poly_matrix = init_poly_matrix( - terms=col_terms, - shape=(1, underlying.shape[1]), - ) + for monomial, value in underlying_terms.items(): - acc = self.acc_func(acc, poly_matrix) + degree = sum((count for var_idx, count in monomial if var_idx in variable_indices)) - state, poly_matrix = acc.apply(state=state) + if degree <= self.degree: + terms_row_col[monomial] = value + + terms[row, col] = terms_row_col + + poly_matrix = init_poly_matrix( + terms=dict(terms), + shape=underlying.shape, + ) return state, poly_matrix diff --git a/polymatrix/expression/mixins/vstackexprmixin.py b/polymatrix/expression/mixins/vstackexprmixin.py index b27828e..2929f7a 100644 --- a/polymatrix/expression/mixins/vstackexprmixin.py +++ b/polymatrix/expression/mixins/vstackexprmixin.py @@ -23,8 +23,8 @@ class VStackExprMixin(ExpressionBaseMixin): all_underlying = [] for expr in self.underlying: - state, polymat = expr.apply(state=state) - all_underlying.append(polymat) + state, polymatrix = expr.apply(state=state) + all_underlying.append(polymatrix) for underlying in all_underlying: assert underlying.shape[1] == all_underlying[0].shape[1], f'{underlying.shape[1]} not equal {all_underlying[0].shape[1]}' @@ -36,9 +36,9 @@ class VStackExprMixin(ExpressionBaseMixin): shape: tuple[int, int] def get_poly(self, row: int, col: int) -> dict[tuple[int, ...], float]: - for polymat, (row_start, row_end) in zip(self.all_underlying, self.underlying_row_range): + for polymatrix, (row_start, row_end) in zip(self.all_underlying, self.underlying_row_range): if row_start <= row < row_end: - return polymat.get_poly( + return polymatrix.get_poly( row=row - row_start, col=col, ) @@ -51,12 +51,12 @@ class VStackExprMixin(ExpressionBaseMixin): initial=0) )) - n_row = sum(polymat.shape[0] for polymat in all_underlying) + n_row = sum(polymatrix.shape[0] for polymatrix in all_underlying) - polymat = VStackPolyMatrix( + polymatrix = VStackPolyMatrix( all_underlying=all_underlying, shape=(n_row, all_underlying[0].shape[1]), underlying_row_range=underlying_row_range, ) - return state, polymat + return state, polymatrix diff --git a/polymatrix/expression/parametrizeexpr.py b/polymatrix/expression/parametrizeexpr.py new file mode 100644 index 0000000..f98fdcd --- /dev/null +++ b/polymatrix/expression/parametrizeexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.parametrizeexprmixin import ParametrizeExprMixin + +class ParametrizeExpr(ParametrizeExprMixin): + pass diff --git a/polymatrix/expression/parametrizesymbolsexpr.py b/polymatrix/expression/parametrizetermsexpr.py index a81ec23..18b4de2 100644 --- a/polymatrix/expression/parametrizesymbolsexpr.py +++ b/polymatrix/expression/parametrizetermsexpr.py @@ -1,4 +1,4 @@ from polymatrix.expression.mixins.parametrizetermsexprmixin import ParametrizeTermsExprMixin -class ParametrizeSymbolsExpr(ParametrizeTermsExprMixin): +class ParametrizeTermsExpr(ParametrizeTermsExprMixin): pass diff --git a/polymatrix/expression/quadraticinexpr.py b/polymatrix/expression/quadraticinexpr.py index 59f452b..bab76f6 100644 --- a/polymatrix/expression/quadraticinexpr.py +++ b/polymatrix/expression/quadraticinexpr.py @@ -1,4 +1,4 @@ -from polymatrix.expression.mixins.quadraticreprexprmixin import QuadraticReprExprMixin +from polymatrix.expression.mixins.quadraticinexprmixin import QuadraticInExprMixin -class QuadraticInExpr(QuadraticReprExprMixin): +class QuadraticInExpr(QuadraticInExprMixin): pass diff --git a/polymatrix/expression/setelementatexpr.py b/polymatrix/expression/setelementatexpr.py new file mode 100644 index 0000000..147d1f8 --- /dev/null +++ b/polymatrix/expression/setelementatexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.setelementatexprmixin import SetElementAtExprMixin + +class SetElementAtExpr(SetElementAtExprMixin): + pass diff --git a/polymatrix/expression/sosmonomialsexpr.py b/polymatrix/expression/sosmonomialsexpr.py new file mode 100644 index 0000000..63dfcef --- /dev/null +++ b/polymatrix/expression/sosmonomialsexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.quadraticmonomialsexprmixin import QuadraticMonomialsExprMixin + +class SOSMonomialsExpr(QuadraticMonomialsExprMixin): + pass diff --git a/polymatrix/expression/substituteexpr.py b/polymatrix/expression/substituteexpr.py new file mode 100644 index 0000000..3fd4916 --- /dev/null +++ b/polymatrix/expression/substituteexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.substituteexprmixin import SubstituteExprMixin + +class SubstituteExpr(SubstituteExprMixin): + pass diff --git a/polymatrix/expression/subtractmonomialsexpr.py b/polymatrix/expression/subtractmonomialsexpr.py new file mode 100644 index 0000000..9bc92ea --- /dev/null +++ b/polymatrix/expression/subtractmonomialsexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.subtractmonomialsexprmixin import SubtractMonomialsExprMixin + +class SubtractMonomialsExpr(SubtractMonomialsExprMixin): + pass diff --git a/polymatrix/expression/truncateexpr.py b/polymatrix/expression/truncateexpr.py new file mode 100644 index 0000000..b555480 --- /dev/null +++ b/polymatrix/expression/truncateexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.truncateexprmixin import TruncateExprMixin + +class TruncateExpr(TruncateExprMixin): + pass diff --git a/polymatrix/expression/utils/getderivativemonomials.py b/polymatrix/expression/utils/getderivativemonomials.py new file mode 100644 index 0000000..9c5ac98 --- /dev/null +++ b/polymatrix/expression/utils/getderivativemonomials.py @@ -0,0 +1,114 @@ +import collections +from polymatrix.expression.expressionstate import ExpressionState + + +def get_derivative_monomials( + monomial_terms, + diff_wrt_variable: int, + state: ExpressionState, + considered_variables: set, + introduce_derivatives: bool, +): + + # if self.introduce_derivatives: + + # raise Exception('not implemented') + + # def gen_new_variables(): + # for monomial in monomial_terms.keys(): + # for var in monomial: + # if var not in diff_wrt_variables and var not in considered_variables: + # yield var + + # new_variables = set(gen_new_variables()) + + # new_considered_variables = considered_variables | new_variables + + # def acc_state_candidates(acc, new_variable): + # state, candidates = acc + + # key = init_derivative_key( + # variable=new_variable, + # with_respect_to=diff_wrt_variable, + # ) + # state = state.register(key=key, n_param=1) + + # state, auxillary_derivation_terms = get_derivative_terms( + # monomial_terms=state.auxillary_equations[new_variable], + # diff_wrt_variable=diff_wrt_variable, + # state=state, + # considered_variables=new_considered_variables, + # ) + + # if 1 < len(auxillary_derivation_terms): + # derivation_variable = state.offset_dict[key][0] + + # state = dataclasses.replace( + # state, + # auxillary_equations=state.auxillary_equations | {derivation_variable: auxillary_derivation_terms}, + # ) + + # return state, candidates + (new_variable,) + + # else: + # return state, candidates + + # *_, (state, confirmed_variables) = itertools.accumulate( + # new_variables, + # acc_state_candidates, + # initial=(state, tuple()), + # ) + + # else: + # confirmed_variables = tuple() + + derivation_terms = collections.defaultdict(float) + + for monomial, value in monomial_terms.items(): + + # # count powers for each variable + # monomial_cnt = dict(collections.Counter(monomial)) + monomial_cnt = dict(monomial) + + def differentiate_monomial(dependent_variable, derivation_variable=None): + def gen_diff_monomial(): + for current_variable, current_count in monomial: + + if current_variable is dependent_variable: + sel_counter = current_count - 1 + + else: + sel_counter = current_count + + # for _ in range(sel_counter): + # yield current_variable + if 0 < sel_counter: + yield current_variable, sel_counter + + if derivation_variable is not None: + yield derivation_variable + + diff_monomial = tuple(sorted(gen_diff_monomial())) + + return diff_monomial, value * monomial_cnt[dependent_variable] + + if diff_wrt_variable in monomial_cnt: + diff_monomial, value = differentiate_monomial(diff_wrt_variable) + derivation_terms[diff_monomial] += value + + # # only used if introduce_derivatives == True + # for candidate_variable in monomial_cnt.keys(): + # if candidate_variable in considered_variables or candidate_variable in confirmed_variables: + # key = init_derivative_key( + # variable=candidate_variable, + # with_respect_to=diff_wrt_variable, + # ) + # derivation_variable = state.offset_dict[key][0] + + # diff_monomial, value = differentiate_monomial( + # dependent_variable=candidate_variable, + # derivation_variable=derivation_variable, + # ) + # derivation_terms[diff_monomial] += value + + return state, dict(derivation_terms)
\ No newline at end of file diff --git a/polymatrix/expression/utils/getmonomialindices.py b/polymatrix/expression/utils/getmonomialindices.py new file mode 100644 index 0000000..6d6c599 --- /dev/null +++ b/polymatrix/expression/utils/getmonomialindices.py @@ -0,0 +1,19 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + + +def get_monomial_indices(state, monomials) -> tuple: + + state, monomials_obj = monomials.apply(state) + + assert monomials_obj.shape[1] == 1, f'{monomials_obj.shape=}' + + def gen_indices(): + for row in range(monomials_obj.shape[0]): + row_terms = monomials_obj.get_poly(row, 0) + + assert len(row_terms) == 1, f'{row_terms} contains more than one term' + + yield from row_terms.keys() + + return state, tuple(gen_indices()) + diff --git a/polymatrix/expression/utils/gettupleremainder.py b/polymatrix/expression/utils/gettupleremainder.py new file mode 100644 index 0000000..2a22ead --- /dev/null +++ b/polymatrix/expression/utils/gettupleremainder.py @@ -0,0 +1,9 @@ +def get_tuple_remainder(l1, l2): + # create a copy of l1 + l1_copy = list(l1) + + # remove from the copy all elements of l2 + for e in l2: + l1_copy.remove(e) + + return tuple(l1_copy)
\ No newline at end of file diff --git a/polymatrix/expression/utils/getvariableindices.py b/polymatrix/expression/utils/getvariableindices.py index 71a91f6..9aa3619 100644 --- a/polymatrix/expression/utils/getvariableindices.py +++ b/polymatrix/expression/utils/getvariableindices.py @@ -17,14 +17,16 @@ def get_variable_indices(state, variables): assert len(row_terms) == 1, f'{row_terms} contains more than one term' for monomial in row_terms.keys(): - assert len(monomial) == 1, f'{monomial} contains more than one variable' - yield monomial[0] + assert len(monomial) == 1, f'{monomial=} contains more than one variable' + assert monomial[0][1] == 1, f'{monomial[0]=}' + + yield monomial[0][0] return state, tuple(gen_indices()) else: - raise Exception('not supported anymore') + # raise Exception('not supported anymore') if not isinstance(variables, tuple): variables = (variables,) diff --git a/polymatrix/expression/utils/mergemonomialindices.py b/polymatrix/expression/utils/mergemonomialindices.py new file mode 100644 index 0000000..7f8ba15 --- /dev/null +++ b/polymatrix/expression/utils/mergemonomialindices.py @@ -0,0 +1,27 @@ +from polymatrix.expression.utils.sortmonomialindices import sort_monomial_indices + + +def merge_monomial_indices(monomials): + if len(monomials) == 0: + return tuple() + + elif len(monomials) == 1: + return monomials[0] + + else: + + m1_dict = dict(monomials[0]) + + for other in monomials[1:]: + for index, count in other: + if index in m1_dict: + m1_dict[index] += count + else: + m1_dict[index] = count + + # return tuple(sorted( + # m1_dict.items(), + # key=lambda m: m[0], + # )) + + return sort_monomial_indices(m1_dict.items()) diff --git a/polymatrix/expression/utils/multiplypolynomial.py b/polymatrix/expression/utils/multiplypolynomial.py new file mode 100644 index 0000000..0b6f5b6 --- /dev/null +++ b/polymatrix/expression/utils/multiplypolynomial.py @@ -0,0 +1,24 @@ +import itertools +import math + +from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices + + +def multiply_polynomial(left, right, terms): + for (left_monomial, left_value), (right_monomial, right_value) \ + in itertools.product(left.items(), right.items()): + + value = left_value * right_value + + if math.isclose(value, 0, abs_tol=1e-12): + continue + + monomial = merge_monomial_indices((left_monomial, right_monomial)) + + if monomial not in terms: + terms[monomial] = 0 + + terms[monomial] += value + + if math.isclose(terms[monomial], 0, abs_tol=1e-12): + del terms[monomial] diff --git a/polymatrix/expression/utils/sortmonomialindices.py b/polymatrix/expression/utils/sortmonomialindices.py new file mode 100644 index 0000000..3d4734e --- /dev/null +++ b/polymatrix/expression/utils/sortmonomialindices.py @@ -0,0 +1,5 @@ +def sort_monomial_indices(monomial): + return tuple(sorted( + monomial, + key=lambda m: m[0], + ))
\ No newline at end of file diff --git a/polymatrix/expression/utils/sortmonomials.py b/polymatrix/expression/utils/sortmonomials.py new file mode 100644 index 0000000..46e7118 --- /dev/null +++ b/polymatrix/expression/utils/sortmonomials.py @@ -0,0 +1,5 @@ +def sort_monomials(monomials): + return tuple(sorted( + monomials, + key=lambda m: (sum(count for _, count in m), len(m), m), + ))
\ No newline at end of file diff --git a/polymatrix/expression/utils/splitmonomialindices.py b/polymatrix/expression/utils/splitmonomialindices.py new file mode 100644 index 0000000..82e5454 --- /dev/null +++ b/polymatrix/expression/utils/splitmonomialindices.py @@ -0,0 +1,27 @@ +def split_monomial_indices(monomial): + left = [] + right = [] + + is_left = True + + for idx, count in monomial: + count_left = count // 2 + + is_uneven = count % 2 + if is_uneven: + if is_left: + count_left = count_left + 1 + + is_left = not is_left + + count_right = count - count_left + + if 0 < count_left: + left.append((idx, count_left)) + + if 0 < count_right: + right.append((idx, count - count_left)) + + # print((monomial, tuple(left), tuple(right))) + + return tuple(left), tuple(right)
\ No newline at end of file diff --git a/polymatrix/expression/utils/subtractmonomialindices.py b/polymatrix/expression/utils/subtractmonomialindices.py new file mode 100644 index 0000000..a8aef1c --- /dev/null +++ b/polymatrix/expression/utils/subtractmonomialindices.py @@ -0,0 +1,18 @@ +from polymatrix.expression.utils.sortmonomialindices import sort_monomial_indices + + +def subtract_monomial_indices(m1, m2): + m1_dict = dict(m1) + + for index, count in m2: + m1_dict[index] -= count + + if m1_dict[index] == 0: + del m1_dict[index] + + # return tuple(sorted( + # m1_dict.items(), + # key=lambda m: m[0], + # )) + + return sort_monomial_indices(m1_dict.items()) |