diff options
| author | Michael Schneeberger <michael.schneeberger@fhnw.ch> | 2022-04-12 09:30:35 +0200 |
|---|---|---|
| committer | Michael Schneeberger <michael.schneeberger@fhnw.ch> | 2022-04-12 09:30:35 +0200 |
| commit | f9ec423683462800c757119c08947b742458639e (patch) | |
| tree | c75a4686696368b160a46ed4f67601be8581492f | |
| parent | remove shape property, introduce accumulate function, reimplement derivative ... (diff) | |
| download | polymatrix-f9ec423683462800c757119c08947b742458639e.tar.gz polymatrix-f9ec423683462800c757119c08947b742458639e.zip | |
bug fixes and clean ups
86 files changed, 539 insertions, 1817 deletions
diff --git a/main.py b/main.py deleted file mode 100644 index 4079ac5..0000000 --- a/main.py +++ /dev/null @@ -1,159 +0,0 @@ -import sympy -import scipy.special -import scipy.sparse -import numpy as np -from polymatrix.polystruct import init_equation, init_poly_matrix, init_poly_vector -from polymatrix.sympyutils import poly_to_data_coord, poly_to_matrix -from polymatrix.utils import variable_to_index - -binom = scipy.special.binom - -def skew_symmetric(degree, p_row, p_col): - if p_col < p_row: - return p_col, p_row, -1 - -def zero_diagonal_for_uneven_degree(degree, p_row, p_col): - if degree % 2 == 1 and p_row == p_col: - return 0 - -def gradient(degree, v_row, monom): - if degree == 1: - factor = sum(v_row==e for e in monom) + 1 - - if monom[-1] < v_row: - n_v_row = monom[-1] - n_monom = sorted(monom + (v_row,), reverse=True) - - if v_row <= monom[-1]: - n_v_row = v_row - n_monom = monom - - return n_v_row, n_monom, factor - -x = sympy.symbols('v_dc, i_d, i_q') -vdc, id, iq = x - -g_dc = 3 -wl = 106 - -jr_poly_list = [ - g_dc, -1, g_dc*wl, - 1, 0, wl, - -g_dc*wl, -wl, 0 -] -jr_terms = poly_to_data_coord(list(sympy.poly(p, x) for p in jr_poly_list)) -# print(jr_terms) - -h = vdc**2/2 + id**2/2 + iq**2/2 -dh_poly_list = [ - sympy.diff(h, vdc), sympy.diff(h, id), sympy.diff(h, iq) -] -dh_terms = poly_to_data_coord(list(sympy.poly(p, x) for p in dh_poly_list)) -# print(dh_terms) - -mg_poly_list = [ - g_dc-id, -iq, - vdc+1, 0, - 0, vdc+1 -] -mg_terms = poly_to_data_coord(list(sympy.poly(p, x) for p in mg_poly_list)) -# print(mg_terms) - -nabla_h = init_poly_vector(subs=dh_terms) -JR = init_poly_matrix(subs=jr_terms) -mg = init_poly_matrix(subs=mg_terms) - -nabla_ha = init_poly_vector( - degrees=(1,), - re_index_func_2=gradient, -) -JRa = init_poly_matrix( - degrees=(0,1,2), - re_index_func=skew_symmetric, - subs_func=zero_diagonal_for_uneven_degree, -) -u = init_poly_vector(degrees=(1,2)) - -eq = init_equation( - terms = [(JR, nabla_ha), (JRa, nabla_ha), (JRa, nabla_h), (mg, u)], - n_var = 2, - # terms = [(JR, nabla_ha)], - # n_var = 2, -) - -# n_var = 2 -# n_param_nabla_ha = sum(n_var*binom(n_var+p-1, p) for p in nabla_ha.degrees) -# n_param_JRa = sum(n_var**2*binom(n_var+p-1, p) for p in JRa.non_zero_degrees) -# n_param_u = sum(n_var*binom(n_var+p-1, p) for p in u.non_zero_degrees) -# total = n_param_nabla_ha+n_param_JRa+n_param_u - -# print(f'{n_param_nabla_ha=}') -# print(f'{n_param_JRa=}') -# print(f'{n_param_u=}') -# print(f'{total=}') - -# print(binom(total+2-1, 2)) - -# mat = init_poly_matrix( -# degrees=(1,), -# re_index_func=skew_symmetric, -# subs_func=zero_diagonal_for_uneven_degree, -# ) -# vec = init_poly_vector( -# subs={0: {(0, 0): 1, (1, 0): 1}}, -# ) - -# # mat = init_poly_matrix( -# # subs={0: {(0, 0): 1, (1, 0): 1, (0, 1): 1, (1, 1): 1}}, -# # ) - -# # vec = init_poly_vector( -# # degrees=(1,), -# # re_index_func_2=gradient, -# # ) - -# eq = init_equation( -# terms = [(mat, vec)], -# n_var = 2, -# ) - -eq_tuples, offset_dict, n_param = eq.create() - -print(f'{list(offset_dict.values())=}') -print(f'{n_param=}') - -# mapping from row index to entry in eq_tuples -rows_to_eq = list(set(key for eq_tuple_degree in eq_tuples.values() for key in eq_tuple_degree.keys())) -eq_to_rows = {eq: idx for idx, eq in enumerate(rows_to_eq)} - -print(f'n_eq = {len(rows_to_eq)}') - -# # mapping from col index to entry in eq_tuples -# cols_to_var = list(set(key for eq in eq_tuples[degree].values() for key in eq.keys())) -# var_to_cols = {var: idx for idx, var in enumerate(cols_to_var)} - -def gen_matrix_per_degree(): - for degree, eq_tuple_degree in eq_tuples.items(): - - if 0 < degree: - def gen_coords(eq_tuple_degree=eq_tuple_degree): - for eq, var_dict in eq_tuple_degree.items(): - zero_degree_val = eq_tuples[0][eq][0] - row = eq_to_rows[eq] - for var, value in var_dict.items(): - # col = variable_to_index(n_param, var) - # col = var_to_cols[var] - col = var - yield row, col, value / zero_degree_val - - row, col, data = list(zip(*gen_coords())) - - n_col = int(binom(n_param+degree-1, degree)) - n_row = int(max(row) + 1) - sparse_matrix = scipy.sparse.coo_matrix((data, (row, col)), dtype=np.float64, shape=(n_row, n_col)).toarray() - - yield degree, sparse_matrix - -result = dict(gen_matrix_per_degree()) - -print(result[1].shape) diff --git a/polymatrix/__init__.py b/polymatrix/__init__.py index 3cf9f35..78b4afa 100644 --- a/polymatrix/__init__.py +++ b/polymatrix/__init__.py @@ -10,10 +10,6 @@ from polymatrix.expression.polymatrix import PolyMatrix def from_( data: tuple[tuple[float]], - # name: str, - # shape: tuple, - # degrees: tuple, - # re_index: typing.Callable[[int, int, int, tuple[int, ...]], tuple[int, int, int, tuple[int, ...], float]] = None, ): return init_expression( init_from_array_expr(data) @@ -57,14 +53,16 @@ def v_stack( def kkt( cost: Expression, - equality: Expression, - variables: Expression + variables: Expression, + equality: Expression = None, + inequality: Expression = None, ): return init_expression( init_kkt_expr( cost=cost, equality=equality, variables=variables, + inequality=inequality, ) ) diff --git a/polymatrix/addexpr.py b/polymatrix/addexpr.py deleted file mode 100644 index 983645e..0000000 --- a/polymatrix/addexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.mixins.addexprmixin import AddExprMixin - -class AddExpr(AddExprMixin): - pass
\ No newline at end of file diff --git a/polymatrix/exprcontainer.py b/polymatrix/exprcontainer.py deleted file mode 100644 index e5766f1..0000000 --- a/polymatrix/exprcontainer.py +++ /dev/null @@ -1,10 +0,0 @@ -import abc -import typing -from polymatrix.mixins.exprcontainermixin import ExprContainerMixin, ExprType - -class ExprContainer( - ExprContainerMixin[ExprType], - typing.Generic[ExprType], - abc.ABC, -): - pass diff --git a/polymatrix/expression/additionexpr.py b/polymatrix/expression/additionexpr.py new file mode 100644 index 0000000..6a43ba3 --- /dev/null +++ b/polymatrix/expression/additionexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.additionexprmixin import AdditionExprMixin + +class AdditionExpr(AdditionExprMixin): + pass diff --git a/polymatrix/impl/polymatrixmultexprimpl.py b/polymatrix/expression/impl/additionexprimpl.py index 7f8375e..2fd7220 100644 --- a/polymatrix/impl/polymatrixmultexprimpl.py +++ b/polymatrix/expression/impl/additionexprimpl.py @@ -1,9 +1,9 @@ import dataclass_abc -from polymatrix.polymatrixmultexpr import PolyMatrixMultExpr +from polymatrix.expression.additionexpr import AdditionExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixMultExprImpl(PolyMatrixMultExpr): +class AdditionExprImpl(AdditionExpr): left: ExpressionBaseMixin right: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/evalexprimpl.py b/polymatrix/expression/impl/evalexprimpl.py index 1730170..cd97155 100644 --- a/polymatrix/expression/impl/evalexprimpl.py +++ b/polymatrix/expression/impl/evalexprimpl.py @@ -7,4 +7,4 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin class EvalExprImpl(EvalExpr): underlying: ExpressionBaseMixin variables: tuple - eval_values: tuple + values: tuple diff --git a/polymatrix/expression/impl/kktexprimpl.py b/polymatrix/expression/impl/kktexprimpl.py index e7df9bf..415b206 100644 --- a/polymatrix/expression/impl/kktexprimpl.py +++ b/polymatrix/expression/impl/kktexprimpl.py @@ -8,4 +8,5 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin class KKTExprImpl(KKTExpr): cost: ExpressionMixin equality: ExpressionMixin + inequality: ExpressionMixin variables: ExpressionBaseMixin diff --git a/polymatrix/impl/polymatrixaddexprimpl.py b/polymatrix/expression/impl/matrixmultexprimpl.py index f30495b..0b093ab 100644 --- a/polymatrix/impl/polymatrixaddexprimpl.py +++ b/polymatrix/expression/impl/matrixmultexprimpl.py @@ -1,9 +1,9 @@ import dataclass_abc -from polymatrix.polymatrixaddexpr import PolyMatrixAddExpr +from polymatrix.expression.matrixmultexpr import MatrixMultExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixAddExprImpl(PolyMatrixAddExpr): +class MatrixMultExprImpl(MatrixMultExpr): left: ExpressionBaseMixin right: ExpressionBaseMixin diff --git a/polymatrix/expression/impl/parametrizesymbolsexprimpl.py b/polymatrix/expression/impl/parametrizetermsexprimpl.py index af1fe3e..b851306 100644 --- a/polymatrix/expression/impl/parametrizesymbolsexprimpl.py +++ b/polymatrix/expression/impl/parametrizetermsexprimpl.py @@ -4,7 +4,7 @@ from polymatrix.expression.parametrizesymbolsexpr import ParametrizeSymbolsExpr from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin @dataclass_abc.dataclass_abc(frozen=True) -class ParametrizeSymbolsExprImpl(ParametrizeSymbolsExpr): +class ParametrizeTermsExprImpl(ParametrizeSymbolsExpr): name: str underlying: ExpressionBaseMixin variables: tuple diff --git a/polymatrix/expression/impl/polymatriximpl.py b/polymatrix/expression/impl/polymatriximpl.py index 8dd6fca..63f577a 100644 --- a/polymatrix/expression/impl/polymatriximpl.py +++ b/polymatrix/expression/impl/polymatriximpl.py @@ -6,4 +6,3 @@ from polymatrix.expression.polymatrix import PolyMatrix class PolyMatrixImpl(PolyMatrix): terms: dict shape: tuple[int, ...] - # aux_terms: tuple[dict[tuple[int, ...], float]] diff --git a/polymatrix/expression/init/__init__.py b/polymatrix/expression/init/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/polymatrix/expression/init/__init__.py diff --git a/polymatrix/init/initpolymatrixmultexpr.py b/polymatrix/expression/init/initadditionexpr.py index 2c423c9..e5b31c8 100644 --- a/polymatrix/init/initpolymatrixmultexpr.py +++ b/polymatrix/expression/init/initadditionexpr.py @@ -1,12 +1,12 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.impl.polymatrixmultexprimpl import PolyMatrixMultExprImpl +from polymatrix.expression.impl.additionexprimpl import AdditionExprImpl -def init_poly_matrix_mult_expr( +def init_addition_expr( left: ExpressionBaseMixin, right: ExpressionBaseMixin, ): - return PolyMatrixMultExprImpl( + return AdditionExprImpl( left=left, right=right, ) diff --git a/polymatrix/expression/init/initevalexpr.py b/polymatrix/expression/init/initevalexpr.py index deeb4e9..9d0fc0d 100644 --- a/polymatrix/expression/init/initevalexpr.py +++ b/polymatrix/expression/init/initevalexpr.py @@ -1,21 +1,30 @@ -from tkinter import Variable +import numpy as np from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.impl.evalexprimpl import EvalExprImpl def init_eval_expr( underlying: ExpressionBaseMixin, - variables: tuple, - eval_values: tuple, + values: tuple, + variables: tuple = None, ): + + if variables is None: + assert isinstance(values, tuple) + + variables, values = tuple(zip(*values)) + + elif isinstance(values, np.ndarray): + values = tuple(values) + + elif not isinstance(values, tuple): + values = (values,) + if not isinstance(variables, tuple): variables = (variables,) - if not isinstance(eval_values, tuple): - eval_values = (eval_values,) - return EvalExprImpl( underlying=underlying, variables=variables, - eval_values=eval_values, + values=values, ) diff --git a/polymatrix/expression/init/initkktexpr.py b/polymatrix/expression/init/initkktexpr.py index d8bfba8..f634d97 100644 --- a/polymatrix/expression/init/initkktexpr.py +++ b/polymatrix/expression/init/initkktexpr.py @@ -5,11 +5,13 @@ from polymatrix.expression.impl.kktexprimpl import KKTExprImpl def init_kkt_expr( cost: ExpressionMixin, - equality: ExpressionMixin, variables: ExpressionBaseMixin, + equality: ExpressionMixin = None, + inequality: ExpressionMixin = None, ): return KKTExprImpl( cost=cost, equality=equality, + inequality=inequality, variables=variables, ) diff --git a/polymatrix/init/initpolymatrixaddexpr.py b/polymatrix/expression/init/initmatrixmultexpr.py index 43c2d1f..9081368 100644 --- a/polymatrix/init/initpolymatrixaddexpr.py +++ b/polymatrix/expression/init/initmatrixmultexpr.py @@ -1,12 +1,12 @@ -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.impl.polymatrixaddexprimpl import PolyMatrixAddExprImpl +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.matrixmultexprimpl import MatrixMultExprImpl -def init_poly_matrix_add_expr( +def init_matrix_mult_expr( left: ExpressionBaseMixin, right: ExpressionBaseMixin, ): - return PolyMatrixAddExprImpl( + return MatrixMultExprImpl( left=left, right=right, ) diff --git a/polymatrix/expression/init/initparametrizesymbolsexpr.py b/polymatrix/expression/init/initparametrizetermsexpr.py index 8d7db55..b932825 100644 --- a/polymatrix/expression/init/initparametrizesymbolsexpr.py +++ b/polymatrix/expression/init/initparametrizetermsexpr.py @@ -1,13 +1,13 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.expression.impl.parametrizesymbolsexprimpl import ParametrizeSymbolsExprImpl +from polymatrix.expression.impl.parametrizetermsexprimpl import ParametrizeTermsExprImpl -def init_parametrize_symbols_expr( +def init_parametrize_terms_expr( name: str, underlying: ExpressionBaseMixin, variables: tuple, ): - return ParametrizeSymbolsExprImpl( + return ParametrizeTermsExprImpl( name=name, underlying=underlying, variables=variables, diff --git a/polymatrix/polymatrixmultexpr.py b/polymatrix/expression/matrixmultexpr.py index 6b7f997..62f0cf3 100644 --- a/polymatrix/polymatrixmultexpr.py +++ b/polymatrix/expression/matrixmultexpr.py @@ -1,4 +1,4 @@ from polymatrix.expression.mixins.matrixmultexprmixin import MatrixMultExprMixin -class PolyMatrixMultExpr(MatrixMultExprMixin): +class MatrixMultExpr(MatrixMultExprMixin): pass diff --git a/polymatrix/expression/mixins/__init__.py b/polymatrix/expression/mixins/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/polymatrix/expression/mixins/__init__.py diff --git a/polymatrix/expression/mixins/accumulateexprmixin.py b/polymatrix/expression/mixins/accumulateexprmixin.py index c92051a..1f834bf 100644 --- a/polymatrix/expression/mixins/accumulateexprmixin.py +++ b/polymatrix/expression/mixins/accumulateexprmixin.py @@ -5,7 +5,7 @@ import typing from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class AccumulateExprMixin(ExpressionBaseMixin): @@ -26,8 +26,8 @@ class AccumulateExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) acc = self.initial diff --git a/polymatrix/expression/mixins/additionexprmixin.py b/polymatrix/expression/mixins/additionexprmixin.py index cab8161..e068ba1 100644 --- a/polymatrix/expression/mixins/additionexprmixin.py +++ b/polymatrix/expression/mixins/additionexprmixin.py @@ -6,10 +6,10 @@ import dataclass_abc from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState -class AddExprMixin(ExpressionBaseMixin): +class AdditionExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod def left(self) -> ExpressionBaseMixin: @@ -28,8 +28,8 @@ class AddExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, left = self.left.apply(state=state) state, right = self.right.apply(state=state) diff --git a/polymatrix/expression/mixins/derivativeexprmixin.py b/polymatrix/expression/mixins/derivativeexprmixin.py index 6098aa9..43c0efe 100644 --- a/polymatrix/expression/mixins/derivativeexprmixin.py +++ b/polymatrix/expression/mixins/derivativeexprmixin.py @@ -11,7 +11,8 @@ 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.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getvariableindices import get_variable_indices class DerivativeExprMixin(ExpressionBaseMixin): @@ -45,39 +46,41 @@ class DerivativeExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: - match self.variables: - case ExpressionBaseMixin(): - state, variables = self.variables.apply(state) + # match self.variables: + # case ExpressionBaseMixin(): + # state, variables = self.variables.apply(state) - assert variables.shape[1] == 1 + # assert variables.shape[1] == 1 - def gen_indices(): - for row in range(variables.shape[0]): - row_terms = variables.get_poly(row, 0) + # def gen_indices(): + # for row in range(variables.shape[0]): + # row_terms = variables.get_poly(row, 0) - assert len(row_terms) == 1 + # assert len(row_terms) == 1 - for monomial in row_terms.keys(): - assert len(monomial) == 1 - yield monomial[0] + # for monomial in row_terms.keys(): + # assert len(monomial) == 1 + # yield monomial[0] - diff_wrt_variables = tuple(gen_indices()) + # diff_wrt_variables = tuple(gen_indices()) - case _: - def gen_indices(): - for variable in self.variables: - if variable in state.offset_dict: - yield state.offset_dict[variable][0] + # case _: + # def gen_indices(): + # for variable in self.variables: + # if variable in state.offset_dict: + # yield state.offset_dict[variable][0] - diff_wrt_variables = tuple(gen_indices()) + # diff_wrt_variables = tuple(gen_indices()) + + state, diff_wrt_variables = get_variable_indices(self.variables, state) def get_derivative_terms( monomial_terms, diff_wrt_variable: int, - state: PolyMatrixExprState, + state: ExpressionState, considered_variables: set, ): diff --git a/polymatrix/expression/mixins/determinantexprmixin.py b/polymatrix/expression/mixins/determinantexprmixin.py index 048b263..f7f76eb 100644 --- a/polymatrix/expression/mixins/determinantexprmixin.py +++ b/polymatrix/expression/mixins/determinantexprmixin.py @@ -7,7 +7,7 @@ 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.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class DeterminantExprMixin(ExpressionBaseMixin): @@ -24,8 +24,8 @@ class DeterminantExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: if self in state.cached_polymatrix: return state, state.cached_polymatrix[self] diff --git a/polymatrix/expression/mixins/divisionexprmixin.py b/polymatrix/expression/mixins/divisionexprmixin.py index de5f685..7b5e0c5 100644 --- a/polymatrix/expression/mixins/divisionexprmixin.py +++ b/polymatrix/expression/mixins/divisionexprmixin.py @@ -5,7 +5,7 @@ import dataclasses from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class DivisionExprMixin(ExpressionBaseMixin): @@ -27,8 +27,8 @@ class DivisionExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: if self in state.cached_polymatrix: return state, state.cached_polymatrix[self] diff --git a/polymatrix/expression/mixins/elemmultexprmixin.py b/polymatrix/expression/mixins/elemmultexprmixin.py index 1f6172e..931e022 100644 --- a/polymatrix/expression/mixins/elemmultexprmixin.py +++ b/polymatrix/expression/mixins/elemmultexprmixin.py @@ -5,7 +5,7 @@ import itertools from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class ElemMultExprMixin(ExpressionBaseMixin): @@ -27,8 +27,8 @@ class ElemMultExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, left = self.left.apply(state=state) state, right = self.right.apply(state=state) diff --git a/polymatrix/expression/mixins/evalexprmixin.py b/polymatrix/expression/mixins/evalexprmixin.py index 88e0174..e61cebd 100644 --- a/polymatrix/expression/mixins/evalexprmixin.py +++ b/polymatrix/expression/mixins/evalexprmixin.py @@ -5,7 +5,8 @@ import itertools from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState +from polymatrix.expression.utils.getvariableindices import get_variable_indices class EvalExprMixin(ExpressionBaseMixin): @@ -21,29 +22,19 @@ class EvalExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod - def eval_values(self) -> tuple[float, ...]: + def values(self) -> tuple[float, ...]: ... - # # overwrites abstract method of `ExpressionBaseMixin` - # @property - # def shape(self) -> tuple[int, int]: - # return self.underlying.shape - # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: - assert len(self.variables) == len(self.eval_values) - + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) - def gen_indices(): - for variable in self.variables: - if variable in state.offset_dict: - yield state.offset_dict[variable][0] + state, variable_indices = get_variable_indices(self.variables, state) - variable_indices = tuple(sorted(gen_indices())) + assert len(variable_indices) == len(self.values) terms = {} @@ -64,7 +55,7 @@ class EvalExprMixin(ExpressionBaseMixin): if variable in variable_indices: index = variable_indices.index(variable) - new_value = value * self.eval_values[index] + new_value = value * self.values[index] return monomial, new_value else: diff --git a/polymatrix/expression/mixins/expressionmixin.py b/polymatrix/expression/mixins/expressionmixin.py index 6806293..b43b3be 100644 --- a/polymatrix/expression/mixins/expressionmixin.py +++ b/polymatrix/expression/mixins/expressionmixin.py @@ -4,6 +4,7 @@ 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.initderivativeexpr import init_derivative_expr from polymatrix.expression.init.initdeterminantexpr import init_determinant_expr from polymatrix.expression.init.initdivisionexpr import init_division_expr @@ -12,18 +13,16 @@ from polymatrix.expression.init.initevalexpr import init_eval_expr from polymatrix.expression.init.initforallexpr import init_for_all_expr from polymatrix.expression.init.initfromarrayexpr import init_from_array_expr from polymatrix.expression.init.initgetitemexpr import init_get_item_expr -from polymatrix.expression.init.initparametrizesymbolsexpr import init_parametrize_symbols_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.initquadraticinexpr import init_quadratic_in_expr from polymatrix.expression.init.initrepmatexpr import init_rep_mat_expr from polymatrix.expression.init.inittoquadraticexpr import init_to_quadratic_expr from polymatrix.expression.init.inittransposeexpr import init_transpose_expr -from polymatrix.init.initpolymatrixaddexpr import init_poly_matrix_add_expr -from polymatrix.init.initpolymatrixmultexpr import init_poly_matrix_mult_expr - from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class ExpressionMixin( @@ -37,7 +36,7 @@ class ExpressionMixin( ... # overwrites abstract method of `PolyMatrixExprBaseMixin` - def apply(self, state: PolyMatrixExprState) -> tuple[PolyMatrixExprState, PolyMatrix]: + def apply(self, state: ExpressionState) -> tuple[ExpressionState, PolyMatrix]: return self.underlying.apply(state) # # overwrites abstract method of `PolyMatrixExprBaseMixin` @@ -55,22 +54,29 @@ class ExpressionMixin( # yield self[row, 0] def __add__(self, other: ExpressionBaseMixin) -> 'ExpressionMixin': - # assert self.underlying.shape == other.shape, f'shapes {(self.shape, other.shape)} of polynomial matrix do not match!' - if other is None: return self + match other: + case ExpressionBaseMixin(): + right = other.underlying + case _: + right = init_from_array_expr(other) + return dataclasses.replace( self, - underlying=init_poly_matrix_add_expr( + underlying=init_addition_expr( left=self.underlying, - right=other.underlying, + right=right, ), ) def __radd__(self, other): return self + other + def __sub__(self, other): + return self + other * (-1) + def __getattr__(self, name): attr = getattr(self.underlying, name) @@ -108,12 +114,17 @@ class ExpressionMixin( return dataclasses.replace( self, - underlying=init_poly_matrix_mult_expr( + underlying=init_matrix_mult_expr( left=self.underlying, right=right, ), ) + def __rmatmul__(self, other): + other = init_from_array_expr(other) + + return other @ self + def __truediv__(self, other: ExpressionBaseMixin): return dataclasses.replace( self, @@ -144,7 +155,7 @@ class ExpressionMixin( def parametrize(self, name: str, variables: tuple) -> 'ExpressionMixin': return dataclasses.replace( self, - underlying=init_parametrize_symbols_expr( + underlying=init_parametrize_terms_expr( name=name, underlying=self.underlying, variables=variables, @@ -201,16 +212,16 @@ class ExpressionMixin( ) def eval( - self, - variable: tuple, - value: tuple[float, ...], + self, + values: tuple[float, ...], + variables: tuple = None, ) -> 'ExpressionMixin': return dataclasses.replace( self, underlying=init_eval_expr( underlying=self.underlying, - variables=variable, - eval_values=value, + variables=variables, + values=values, ), ) diff --git a/polymatrix/expression/mixins/forallexprmixin.py b/polymatrix/expression/mixins/forallexprmixin.py index 9a28283..a988b24 100644 --- a/polymatrix/expression/mixins/forallexprmixin.py +++ b/polymatrix/expression/mixins/forallexprmixin.py @@ -6,7 +6,7 @@ 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.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class ForAllExprMixin(ExpressionBaseMixin): @@ -28,8 +28,8 @@ class ForAllExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) variable_indices = tuple(state.offset_dict[variable][0] for variable in self.variables if variable in state.offset_dict) diff --git a/polymatrix/expression/mixins/getitemexprmixin.py b/polymatrix/expression/mixins/getitemexprmixin.py index 57ffe38..26ed728 100644 --- a/polymatrix/expression/mixins/getitemexprmixin.py +++ b/polymatrix/expression/mixins/getitemexprmixin.py @@ -8,7 +8,7 @@ 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.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class GetItemExprMixin(ExpressionBaseMixin): @@ -30,8 +30,8 @@ class GetItemExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) @dataclass_abc.dataclass_abc(frozen=True) diff --git a/polymatrix/expression/mixins/kktexprmixin.py b/polymatrix/expression/mixins/kktexprmixin.py index ecc69ed..b4c9e72 100644 --- a/polymatrix/expression/mixins/kktexprmixin.py +++ b/polymatrix/expression/mixins/kktexprmixin.py @@ -9,7 +9,7 @@ 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.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class KKTExprMixin(ExpressionBaseMixin): @@ -25,57 +25,90 @@ class KKTExprMixin(ExpressionBaseMixin): @property @abc.abstractmethod - def equality(self) -> ExpressionMixin: + def equality(self) -> typing.Optional[ExpressionMixin]: ... @property @abc.abstractmethod - def variables(self) -> ExpressionBaseMixin: + def inequality(self) -> typing.Optional[ExpressionMixin]: ... - # # overwrites abstract method of `ExpressionBaseMixin` - # @property - # def shape(self) -> tuple[int, int]: - # return self.cost.shape[0] + self.equality.shape[0], 1 + @property + @abc.abstractmethod + def variables(self) -> ExpressionBaseMixin: + ... # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, cost = self.cost.apply(state=state) - state, equality = self.equality.apply(state=state) - - # todo: check shape - # assert cost.shape == self.variables.shape assert cost.shape[1] == 1 - state, equality_der = self.equality.diff( - self.variables, - introduce_derivatives=True, - ).apply(state) + 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 + def acc_nu_variables(acc, v): + state, nu_variables = acc - nu_variable = state.n_param - state = state.register(n_param=1) + 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, []), + )) - return state, nu_variables + [nu_variable] + else: + lambda_variables = tuple() - *_, (state, nu_variables) = tuple(itertools.accumulate( - range(equality.shape[0]), - acc_nu_variables, - initial=(state, []), - )) + idx_start = 0 terms = {} - n_row = cost.shape[0] - - for row in range(n_row): + for row in range(cost.shape[0]): try: monomial_terms = cost.get_poly(row, 0) except KeyError: @@ -96,20 +129,58 @@ class KKTExprMixin(ExpressionBaseMixin): monomial_terms[new_monomial] += value - terms[row, 0] = monomial_terms + 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 - idx_start = n_row - - for row in range(equality.shape[0]): + 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(row, 0) + underlying_terms = equality.get_poly(eq_idx, 0) except KeyError: continue - terms[idx_start + row, 0] = underlying_terms + 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 + 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 - idx_start += equality.shape[0] + # complementary slackness + terms[idx_start, 0] = {(r_lambda, r_inequality): 1} + idx_start += 1 poly_matrix = init_poly_matrix( terms=terms, diff --git a/polymatrix/expression/mixins/matrixmultexprmixin.py b/polymatrix/expression/mixins/matrixmultexprmixin.py index 1fc4977..d1c96d4 100644 --- a/polymatrix/expression/mixins/matrixmultexprmixin.py +++ b/polymatrix/expression/mixins/matrixmultexprmixin.py @@ -5,7 +5,7 @@ import itertools from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class MatrixMultExprMixin(ExpressionBaseMixin): @@ -27,8 +27,8 @@ class MatrixMultExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, left = self.left.apply(state=state) state, right = self.right.apply(state=state) diff --git a/polymatrix/expression/mixins/parametrizetermsexprmixin.py b/polymatrix/expression/mixins/parametrizetermsexprmixin.py index 9c9795d..4d46ece 100644 --- a/polymatrix/expression/mixins/parametrizetermsexprmixin.py +++ b/polymatrix/expression/mixins/parametrizetermsexprmixin.py @@ -1,5 +1,6 @@ import abc +import dataclasses import functools import dataclass_abc from polymatrix.expression.init.initexpressionstate import init_expression_state @@ -7,6 +8,7 @@ from polymatrix.expression.init.initexpressionstate import init_expression_state 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.polymatrixasdictmixin import PolyMatrixAsDictMixin from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin @@ -30,54 +32,12 @@ class ParametrizeTermsExprMixin(ExpressionBaseMixin): # def shape(self) -> tuple[int, int]: # return self.underlying.shape - def _internal_apply(self, state: ExpressionStateMixin): - if not hasattr(self, '_terms'): - 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) - - idx_start = state.n_param - n_param = 0 - 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 = {} - collected_terms = [] - - for monomial, value in underlying_terms.items(): - - x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) - - if x_monomial not in collected_terms: - collected_terms.append(x_monomial) - - idx = idx_start + n_param + collected_terms.index(x_monomial) - - new_monomial = monomial + (idx,) - - terms_row_col[new_monomial] = value - - n_param += len(collected_terms) - terms[row, col] = terms_row_col - - state = state.register(key=self, n_param=n_param) - - self._terms = terms - self._start_index = idx_start - self._shape = underlying.shape - self._n_param = n_param - - return state, self._terms + @dataclass_abc.dataclass_abc + class ParametrizeTermsPolyMatrix(PolyMatrixAsDictMixin): + shape: tuple[int, int] + terms: dict + start_index: int + n_param: int @property def param(self) -> tuple[int, int]: @@ -85,22 +45,15 @@ class ParametrizeTermsExprMixin(ExpressionBaseMixin): @dataclass_abc.dataclass_abc(frozen=True) class ParameterExpr(ExpressionBaseMixin): - # n_param: int - # start_index: int - - # @property - # def shape(self) -> tuple[int, int]: - # return self.n_param, 1 - def apply( self, state: ExpressionStateMixin, ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: - state, _ = outer_self._internal_apply(state) + state, poly_matrix = outer_self.apply(state) - n_param = outer_self._n_param - start_index = outer_self._start_index + n_param = poly_matrix.n_param + start_index = poly_matrix.start_index def gen_monomials(): for rel_index in range(n_param): @@ -123,11 +76,60 @@ class ParametrizeTermsExprMixin(ExpressionBaseMixin): state: ExpressionStateMixin, ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: - state, terms = self._internal_apply(state) + if self in state.cached_polymatrix: + return state, state.cached_polymatrix[self] - poly_matrix = init_poly_matrix( + # if not hasattr(self, '_terms'): + 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) + + idx_start = state.n_param + n_param = 0 + 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 = {} + collected_terms = [] + + for monomial, value in underlying_terms.items(): + + x_monomial = tuple(var_idx for var_idx in monomial if var_idx in variable_indices) + + if x_monomial not in collected_terms: + collected_terms.append(x_monomial) + + idx = idx_start + n_param + collected_terms.index(x_monomial) + + new_monomial = monomial + (idx,) + + terms_row_col[new_monomial] = value + + n_param += len(collected_terms) + terms[row, col] = terms_row_col + + state = state.register(key=self, n_param=n_param) + + poly_matrix = ParametrizeTermsExprMixin.ParametrizeTermsPolyMatrix( terms=terms, - shape=self._shape, + shape=underlying.shape, + start_index=idx_start, + n_param=n_param, + ) + + state = dataclasses.replace( + state, + cached_polymatrix=state.cached_polymatrix | {self: poly_matrix}, ) return state, poly_matrix diff --git a/polymatrix/expression/mixins/polymatrixmixin.py b/polymatrix/expression/mixins/polymatrixmixin.py index e5b803e..c0dcac2 100644 --- a/polymatrix/expression/mixins/polymatrixmixin.py +++ b/polymatrix/expression/mixins/polymatrixmixin.py @@ -1,9 +1,12 @@ import abc import dataclasses -import itertools -import collections +import numpy as np +import scipy.sparse import typing +from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin +from polymatrix.utils import monomial_to_index + class PolyMatrixMixin(abc.ABC): @property @@ -15,8 +18,22 @@ class PolyMatrixMixin(abc.ABC): def get_poly(self, row: int, col: int) -> typing.Optional[dict[tuple[int, ...], float]]: ... - def reduce(self): + def get_equations( + self, + state: ExpressionStateMixin, + ): + assert self.shape[1] == 1 + 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,)) + for row in range(self.shape[0]): for col in range(self.shape[1]): @@ -26,84 +43,59 @@ class PolyMatrixMixin(abc.ABC): continue for monomial in underlying_terms.keys(): - yield from monomial + for variable in monomial: + yield variable + + if variable in state.auxillary_equations: + yield from gen_used_auxillary_variables((variable,)) used_variables = set(gen_used_variables()) - variable_index_reverse_map = tuple(sorted(used_variables)) - variable_index_map = {old: new for new, old in enumerate(variable_index_reverse_map)} - + 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) + + A = np.zeros((n_param, 1), dtype=np.float32) + B = np.zeros((n_param, n_param), dtype=np.float32) + C = scipy.sparse.dok_array((n_param, n_param**2), dtype=np.float32) - terms = {} + 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(self.shape[0]): - for col in range(self.shape[1]): - - try: - underlying_terms = self.get_poly(row, col) - except KeyError: - continue - - def gen_updated_monomials(): - for monomial, value in underlying_terms.items(): - new_monomial = tuple(variable_index_map[var] for var in monomial) - yield new_monomial, value - - terms[row, col] = dict(gen_updated_monomials()) - - poly_matrix = dataclasses.replace( - self, - terms=terms, - shape=self.shape, - ) - - return poly_matrix, variable_index_reverse_map - - - # @property - # @abc.abstractmethod - # def terms(self) -> dict[tuple[int, int], dict[tuple[int, ...], dict[tuple[int, ...], float]]]: - # ... - - # @property - # @abc.abstractclassmethod - # def aux_terms(self) -> tuple[dict[tuple[int, ...], float]]: - # ... - - # def get_equality_equations( - # self, - # offset: int = None - # ) -> dict[int, list[tuple[list[int], float]]]: - # if offset is None: - # offset = 0 - - # to_eq_index_dict = {} - # equality_constraints = collections.defaultdict(list) - - # # def gen_key_equation_mapping(): - # for (poly_row, poly_col), terms_x_monomial in self.terms.items(): - # for x_monomial, terms_p_monomial in terms_x_monomial.items(): - # key = (poly_row, poly_col, x_monomial) - - # if key not in to_eq_index_dict: - # to_eq_index_dict[key] = offset - # offset += 1 - - # for p_monomial, value in terms_p_monomial.items(): - # equality_constraints[to_eq_index_dict[key]] += ((p_monomial, value),) - - # return equality_constraints - - # @property - # @abc.abstractmethod - # def degrees(self) -> set[int]: - # ... - - # @abc.abstractmethod - # def get_term( - # self, - # poly_col: int, - # poly_row: int, - # x_monomial: tuple[int, ...] - # ) -> PolyMatrixTerm: - # ... + try: + underlying_terms = self.get_poly(row, 0) + except KeyError: + continue + + populate_matrices( + monomial_terms=underlying_terms, + row=row, + ) + + current_row = self.shape[0] + + for key, monomial_terms in state.auxillary_equations.items(): + if key in ordered_variable_index: + populate_matrices( + monomial_terms=monomial_terms, + row=current_row, + ) + current_row += 1 + + # assert current_row == n_param, f'{current_row} is not {n_param}' + + return (A, B, C), ordered_variable_index diff --git a/polymatrix/expression/mixins/quadraticinexprmixin.py b/polymatrix/expression/mixins/quadraticinexprmixin.py index ee0e5a9..5aeef1b 100644 --- a/polymatrix/expression/mixins/quadraticinexprmixin.py +++ b/polymatrix/expression/mixins/quadraticinexprmixin.py @@ -4,7 +4,7 @@ import abc from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class QuadraticInExprMixin(ExpressionBaseMixin): @@ -26,8 +26,8 @@ class QuadraticInExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) assert underlying.shape == (1, 1) diff --git a/polymatrix/expression/mixins/toquadraticexprmixin.py b/polymatrix/expression/mixins/toquadraticexprmixin.py index d1b6fce..3079e0f 100644 --- a/polymatrix/expression/mixins/toquadraticexprmixin.py +++ b/polymatrix/expression/mixins/toquadraticexprmixin.py @@ -6,7 +6,7 @@ import dataclasses from polymatrix.expression.init.initpolymatrix import init_poly_matrix from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class ToQuadraticExprMixin(ExpressionBaseMixin): @@ -23,12 +23,42 @@ class ToQuadraticExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) + state = [state] + terms = {} - auxillary_equations = {} + auxillary_equations_from_quadratic = {} + + def to_quadratic(monomial_terms): + terms_row_col = collections.defaultdict(float) + + for monomial, value in monomial_terms.items(): + + if 2 < len(monomial): + current_aux = state[0].n_param + terms_row_col[(monomial[0], current_aux)] += value + state[0] = state[0].register(n_param=1) + + for variable in monomial[1:-2]: + auxillary_equations_from_quadratic[current_aux] = { + (variable, current_aux + 1): 1, + (current_aux,): -1, + } + state[0] = state[0].register(n_param=1) + current_aux += 1 + + auxillary_equations_from_quadratic[current_aux] = { + (monomial[-2], monomial[-1]): 1, + (current_aux,): -1, + } + + else: + terms_row_col[monomial] += value + + return dict(terms_row_col) for row in range(underlying.shape[0]): for col in range(underlying.shape[1]): @@ -38,36 +68,47 @@ class ToQuadraticExprMixin(ExpressionBaseMixin): except KeyError: continue - terms_row_col = collections.defaultdict(float) + terms_row_col = to_quadratic( + monomial_terms=underlying_terms, + ) - for monomial, value in underlying_terms.items(): - - if 2 < len(monomial): - current_aux = state.n_param - terms_row_col[(monomial[0], current_aux)] += value - state = state.register(n_param=1) - - for variable in monomial[1:-2]: - auxillary_equations[current_aux] = { - (variable, current_aux + 1): 1, - (current_aux,): -1, - } - state = state.register(n_param=1) - current_aux += 1 - - auxillary_equations[current_aux] = { - (monomial[-2], monomial[-1]): 1, - (current_aux,): -1, - } - - else: - terms_row_col[monomial] += value + # terms_row_col = collections.defaultdict(float) - terms[row, col] = dict(terms_row_col) + # for monomial, value in underlying_terms.items(): + + # if 2 < len(monomial): + # current_aux = state.n_param + # terms_row_col[(monomial[0], current_aux)] += value + # state = state.register(n_param=1) + + # for variable in monomial[1:-2]: + # auxillary_equations[current_aux] = { + # (variable, current_aux + 1): 1, + # (current_aux,): -1, + # } + # state = state.register(n_param=1) + # current_aux += 1 + + # auxillary_equations[current_aux] = { + # (monomial[-2], monomial[-1]): 1, + # (current_aux,): -1, + # } + + # else: + # terms_row_col[monomial] += value + + terms[row, col] = terms_row_col + + def gen_auxillary_equations(): + for key, monomial_terms in state[0].auxillary_equations.items(): + terms_row_col = to_quadratic( + monomial_terms=monomial_terms, + ) + yield key, terms_row_col state = dataclasses.replace( - state, - auxillary_equations=state.auxillary_equations | auxillary_equations, + state[0], + auxillary_equations=dict(gen_auxillary_equations()) | auxillary_equations_from_quadratic, ) poly_matrix = init_poly_matrix( diff --git a/polymatrix/expression/mixins/transposeexprmixin.py b/polymatrix/expression/mixins/transposeexprmixin.py index 4d09380..d24cf53 100644 --- a/polymatrix/expression/mixins/transposeexprmixin.py +++ b/polymatrix/expression/mixins/transposeexprmixin.py @@ -8,7 +8,7 @@ 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.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class TransposeExprMixin(ExpressionBaseMixin): @@ -25,8 +25,8 @@ class TransposeExprMixin(ExpressionBaseMixin): # overwrites abstract method of `PolyMatrixExprBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: state, underlying = self.underlying.apply(state=state) @dataclass_abc.dataclass_abc(frozen=True) diff --git a/polymatrix/expression/mixins/vstackexprmixin.py b/polymatrix/expression/mixins/vstackexprmixin.py index 0e95016..c192ae3 100644 --- a/polymatrix/expression/mixins/vstackexprmixin.py +++ b/polymatrix/expression/mixins/vstackexprmixin.py @@ -6,7 +6,7 @@ from polymatrix.expression.mixins.polymatrixmixin import PolyMatrixMixin from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin from polymatrix.expression.polymatrix import PolyMatrix -from polymatrix.polymatrixexprstate import PolyMatrixExprState +from polymatrix.expression.expressionstate import ExpressionState class VStackExprMixin(ExpressionBaseMixin): @@ -24,8 +24,8 @@ class VStackExprMixin(ExpressionBaseMixin): # overwrites abstract method of `ExpressionBaseMixin` def apply( self, - state: PolyMatrixExprState, - ) -> tuple[PolyMatrixExprState, PolyMatrix]: + state: ExpressionState, + ) -> tuple[ExpressionState, PolyMatrix]: all_underlying = [] for expr in self.underlying: diff --git a/polymatrix/expression/utils/__init__.py b/polymatrix/expression/utils/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/polymatrix/expression/utils/__init__.py diff --git a/polymatrix/expression/utils/getvariableindices.py b/polymatrix/expression/utils/getvariableindices.py new file mode 100644 index 0000000..9379e5f --- /dev/null +++ b/polymatrix/expression/utils/getvariableindices.py @@ -0,0 +1,37 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + + +def get_variable_indices(variables, state): + + if isinstance(variables, ExpressionBaseMixin): + state, variables = variables.apply(state) + + assert variables.shape[1] == 1 + + def gen_indices(): + for row in range(variables.shape[0]): + row_terms = variables.get_poly(row, 0) + + assert len(row_terms) == 1 + + for monomial in row_terms.keys(): + assert len(monomial) == 1 + yield monomial[0] + + return state, tuple(gen_indices()) + + if not isinstance(variables, tuple): + variables = (variables,) + + # assert all(isinstance(variable, type(variables[0])) for variable in variables) + + def gen_indices(): + for variable in variables: + + if isinstance(variable, int): + yield variable + + else: + yield state.offset_dict[variable][0] + + return state, tuple(gen_indices()) diff --git a/polymatrix/impl/addexprimpl.py b/polymatrix/impl/addexprimpl.py deleted file mode 100644 index a3d9e1f..0000000 --- a/polymatrix/impl/addexprimpl.py +++ /dev/null @@ -1,8 +0,0 @@ -import dataclass_abc -from polymatrix.addexpr import AddExpr -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin - -@dataclass_abc.dataclass_abc(frozen=True) -class AddExprImpl(AddExpr): - left: OldPolyMatrixMixin - right: OldPolyMatrixMixin diff --git a/polymatrix/impl/exprcontainerimpl.py b/polymatrix/impl/exprcontainerimpl.py deleted file mode 100644 index 0706150..0000000 --- a/polymatrix/impl/exprcontainerimpl.py +++ /dev/null @@ -1,8 +0,0 @@ -import dataclass_abc -from polymatrix.exprcontainer import ExprContainer -from polymatrix.mixins.exprcontainermixin import ExprType -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin - -@dataclass_abc.dataclass_abc(frozen=True) -class ExprContainerImpl(ExprContainer[ExprType]): - expr: ExpressionBaseMixin diff --git a/polymatrix/impl/multexprimpl.py b/polymatrix/impl/multexprimpl.py deleted file mode 100644 index 7401ffa..0000000 --- a/polymatrix/impl/multexprimpl.py +++ /dev/null @@ -1,8 +0,0 @@ -import dataclass_abc -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin -from polymatrix.multexpr import MultExpr - -@dataclass_abc.dataclass_abc(frozen=True) -class MultExprImpl(MultExpr): - left: OldPolyMatrixMixin - right: OldPolyMatrixMixin diff --git a/polymatrix/impl/oldpolymatriximpl.py b/polymatrix/impl/oldpolymatriximpl.py deleted file mode 100644 index 29590db..0000000 --- a/polymatrix/impl/oldpolymatriximpl.py +++ /dev/null @@ -1,15 +0,0 @@ -import typing -import dataclass_abc - -from polymatrix.oldpolymatrix import OldPolyMatrix - - -@dataclass_abc.dataclass_abc(frozen=True, eq=False) -class PolyMatrixImpl(OldPolyMatrix): - name: str - degrees: list[int] - subs: dict[int, dict[tuple[int, int], float]] - re_index: typing.Callable[[int, int, int, tuple[int, ...]], tuple[int, int, int, tuple[int, ...], float]] - is_constant: bool - shape: tuple[int, int] - is_negated: bool diff --git a/polymatrix/impl/optimizationimpl.py b/polymatrix/impl/optimizationimpl.py deleted file mode 100644 index 3243175..0000000 --- a/polymatrix/impl/optimizationimpl.py +++ /dev/null @@ -1,13 +0,0 @@ -import dataclass_abc - -from polymatrix.optimization import Optimization -from polymatrix.oldpolymatrixexprstate import OldPolyMatrixExprState - - -@dataclass_abc.dataclass_abc(frozen=True, eq=False) -class OptimizationImpl(Optimization): - # n_var: int - state: OldPolyMatrixExprState - equality_constraints: dict[int, dict[tuple[int, int], float]] - inequality_constraints: dict[int, dict[tuple[int, int], float]] - auxillary_equality: dict[int, dict[tuple[int, int], float]] diff --git a/polymatrix/impl/optimizationstateimpl.py b/polymatrix/impl/optimizationstateimpl.py deleted file mode 100644 index acfaddb..0000000 --- a/polymatrix/impl/optimizationstateimpl.py +++ /dev/null @@ -1,11 +0,0 @@ -import dataclass_abc -from polymatrix.oldpolymatrixexprstate import OldPolyMatrixExprStateMixin -from polymatrix.oldpolymatrix import OldPolyMatrix - - -@dataclass_abc.dataclass_abc(frozen=True, eq=False) -class OptimizationStateImpl(OldPolyMatrixExprStateMixin): - n_var: int - n_param: int - offset_dict: dict[tuple[OldPolyMatrix, int], int] - diff --git a/polymatrix/impl/polymatexprimpl.py b/polymatrix/impl/polymatexprimpl.py deleted file mode 100644 index c0d6abb..0000000 --- a/polymatrix/impl/polymatexprimpl.py +++ /dev/null @@ -1,7 +0,0 @@ -import dataclass_abc -from polymatrix.expression.polymatrix import PolyMatrix - -@dataclass_abc.dataclass_abc(frozen=True) -class PolyMatExprImpl(PolyMatrix): - terms: dict - shape: tuple diff --git a/polymatrix/impl/polymatrixarrayexprimpl.py b/polymatrix/impl/polymatrixarrayexprimpl.py deleted file mode 100644 index a005129..0000000 --- a/polymatrix/impl/polymatrixarrayexprimpl.py +++ /dev/null @@ -1,8 +0,0 @@ -import dataclass_abc -from polymatrix.polymatrixarrayexpr import PolyMatrixArrayExpr - -from numpy import array - -@dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixArrayExprImpl(PolyMatrixArrayExpr): - data: array diff --git a/polymatrix/impl/polymatrixexprimpl.py b/polymatrix/impl/polymatrixexprimpl.py deleted file mode 100644 index a5c52f0..0000000 --- a/polymatrix/impl/polymatrixexprimpl.py +++ /dev/null @@ -1,8 +0,0 @@ -import dataclass_abc -from polymatrix.polymatrixexpr import PolyMatrixExpr - -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin - -@dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixExprImpl(PolyMatrixExpr): - underlying: ExpressionBaseMixin diff --git a/polymatrix/impl/polymatrixexprstateimpl.py b/polymatrix/impl/polymatrixexprstateimpl.py deleted file mode 100644 index 25f5c52..0000000 --- a/polymatrix/impl/polymatrixexprstateimpl.py +++ /dev/null @@ -1,12 +0,0 @@ -import typing -import dataclass_abc -import sympy -from polymatrix.polymatrixexprstate import PolyMatrixExprState - - -@dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixExprStateImpl(PolyMatrixExprState): - n_var: int - n_param: int - offset_dict: dict - symbols: typing.Optional[sympy.Symbol] diff --git a/polymatrix/impl/polymatrixparamexprimpl.py b/polymatrix/impl/polymatrixparamexprimpl.py deleted file mode 100644 index afb5388..0000000 --- a/polymatrix/impl/polymatrixparamexprimpl.py +++ /dev/null @@ -1,9 +0,0 @@ -import dataclass_abc -from polymatrix.polymatrixparamexpr import PolyMatrixParamExpr - - -@dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixParamExprImpl(PolyMatrixParamExpr): - name: str - shape: tuple - degrees: tuple diff --git a/polymatrix/impl/polymatrixvalueimpl.py b/polymatrix/impl/polymatrixvalueimpl.py deleted file mode 100644 index c0c9921..0000000 --- a/polymatrix/impl/polymatrixvalueimpl.py +++ /dev/null @@ -1,7 +0,0 @@ -import dataclass_abc -from polymatrix.polymatrixvalue import PolyMatrixTerm - -@dataclass_abc.dataclass_abc(frozen=True) -class PolyMatrixValueImpl(PolyMatrixTerm): - p_monomial: tuple - value: float diff --git a/polymatrix/impl/scalarmultexprimpl.py b/polymatrix/impl/scalarmultexprimpl.py deleted file mode 100644 index 42bec55..0000000 --- a/polymatrix/impl/scalarmultexprimpl.py +++ /dev/null @@ -1,8 +0,0 @@ -import dataclass_abc -from polymatrix.scalarmultexpr import ScalarMultExpr -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin - -@dataclass_abc.dataclass_abc(frozen=True) -class ScalarMultExprImpl(ScalarMultExpr): - left: float - right: ExpressionBaseMixin diff --git a/polymatrix/init/initaddexpr.py b/polymatrix/init/initaddexpr.py deleted file mode 100644 index e660ed1..0000000 --- a/polymatrix/init/initaddexpr.py +++ /dev/null @@ -1,11 +0,0 @@ -from polymatrix.impl.addexprimpl import AddExprImpl -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin - -def init_add_expr( - left: OldPolyMatrixMixin, - right: OldPolyMatrixMixin, -): - return AddExprImpl( - left=left, - right=right, -) diff --git a/polymatrix/init/initexprcontainer.py b/polymatrix/init/initexprcontainer.py deleted file mode 100644 index f790705..0000000 --- a/polymatrix/init/initexprcontainer.py +++ /dev/null @@ -1,9 +0,0 @@ -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.impl.exprcontainerimpl import ExprContainerImpl - -def init_expr_container( - expr: ExpressionBaseMixin, -): - return ExprContainerImpl( - expr=expr, -) diff --git a/polymatrix/init/initmultexpr.py b/polymatrix/init/initmultexpr.py deleted file mode 100644 index dfffd97..0000000 --- a/polymatrix/init/initmultexpr.py +++ /dev/null @@ -1,11 +0,0 @@ -from polymatrix.impl.multexprimpl import MultExprImpl -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin - -def init_mult_expr( - left: OldPolyMatrixMixin, - right: OldPolyMatrixMixin, -): - return MultExprImpl( - left=left, - right=right, -) diff --git a/polymatrix/init/initoptimization.py b/polymatrix/init/initoptimization.py deleted file mode 100644 index 2a340eb..0000000 --- a/polymatrix/init/initoptimization.py +++ /dev/null @@ -1,28 +0,0 @@ -from polymatrix.impl.optimizationimpl import OptimizationImpl -from polymatrix.init.initoptimizationstate import init_optimization_state -from polymatrix.oldpolymatrixexprstate import OldPolyMatrixExprState - - -def init_optimization( - n_var: int, - equality_constraints: dict[int, dict[tuple[int, int], float]] = None, - inequality_constraints: dict[int, dict[tuple[int, int], float]] = None, - auxillary_equality: dict[int, dict[tuple[int, int], float]] = None, -): - state = init_optimization_state(n_var=n_var) - - if equality_constraints is None: - equality_constraints = {} - - if inequality_constraints is None: - inequality_constraints = {} - - if auxillary_equality is None: - auxillary_equality = {} - - return OptimizationImpl( - state=state, - equality_constraints=equality_constraints, - inequality_constraints=inequality_constraints, - auxillary_equality=auxillary_equality, - ) diff --git a/polymatrix/init/initoptimizationstate.py b/polymatrix/init/initoptimizationstate.py deleted file mode 100644 index edfa33e..0000000 --- a/polymatrix/init/initoptimizationstate.py +++ /dev/null @@ -1,20 +0,0 @@ -from polymatrix.impl.optimizationstateimpl import OptimizationStateImpl -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin - - -def init_optimization_state( - n_var: int, - n_param: int = None, - offset_dict: dict[tuple[OldPolyMatrixMixin, int], int] = None, -): - if n_param is None: - n_param = 0 - - if offset_dict is None: - offset_dict = {} - - return OptimizationStateImpl( - n_var=n_var, - n_param=n_param, - offset_dict=offset_dict, - ) diff --git a/polymatrix/init/initpolymatexpr.py b/polymatrix/init/initpolymatexpr.py deleted file mode 100644 index 9f6cd20..0000000 --- a/polymatrix/init/initpolymatexpr.py +++ /dev/null @@ -1,10 +0,0 @@ -from polymatrix.impl.polymatexprimpl import PolyMatExprImpl - -def init_poly_mat_expr( - terms: dict, - shape: tuple, -): - return PolyMatExprImpl( - terms=terms, - shape=shape, -) diff --git a/polymatrix/init/initpolymatrixarrayexpr.py b/polymatrix/init/initpolymatrixarrayexpr.py deleted file mode 100644 index 8f7dbed..0000000 --- a/polymatrix/init/initpolymatrixarrayexpr.py +++ /dev/null @@ -1,10 +0,0 @@ -from numpy import array -from polymatrix.impl.polymatrixarrayexprimpl import PolyMatrixArrayExprImpl - - -def init_poly_matrix_array_expr( - data: array, -): - return PolyMatrixArrayExprImpl( - data=data, -) diff --git a/polymatrix/init/initpolymatrixexpr.py b/polymatrix/init/initpolymatrixexpr.py deleted file mode 100644 index 2749bd3..0000000 --- a/polymatrix/init/initpolymatrixexpr.py +++ /dev/null @@ -1,10 +0,0 @@ -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.impl.polymatrixexprimpl import PolyMatrixExprImpl - - -def init_poly_matrix_expr( - underlying: ExpressionBaseMixin, -): - return PolyMatrixExprImpl( - underlying=underlying, -) diff --git a/polymatrix/init/initpolymatrixexprstate.py b/polymatrix/init/initpolymatrixexprstate.py deleted file mode 100644 index 03649c8..0000000 --- a/polymatrix/init/initpolymatrixexprstate.py +++ /dev/null @@ -1,26 +0,0 @@ -import sympy -from polymatrix.impl.polymatrixexprstateimpl import PolyMatrixExprStateImpl - -def init_poly_matrix_expr_state( - n_var: int = None, - n_param: int = None, - offset_dict: dict = None, - symbols: sympy.Symbol = None, -): - if n_var is None: - n_var = len(symbols) - elif symbols is not None: - assert n_var == len(symbols), f'{n_var} is not equal {len(symbols)}' - - if n_param is None: - n_param = 0 - - if offset_dict is None: - offset_dict = {} - - return PolyMatrixExprStateImpl( - n_var=n_var, - n_param=n_param, - offset_dict=offset_dict, - symbols=symbols, -) diff --git a/polymatrix/init/initpolymatrixparamexpr.py b/polymatrix/init/initpolymatrixparamexpr.py deleted file mode 100644 index 7e1b6d7..0000000 --- a/polymatrix/init/initpolymatrixparamexpr.py +++ /dev/null @@ -1,33 +0,0 @@ -import typing -from polymatrix.impl.polymatrixparamexprimpl import PolyMatrixParamExprImpl - -def init_poly_matrix_param_expr( - name: str, - shape: tuple, - degrees: tuple, - re_index: typing.Callable[[int, int, int, tuple[int, ...]], tuple[int, int, int, tuple[int, ...], float]] = None, -): - - if re_index is None: - return PolyMatrixParamExprImpl( - name=name, - shape=shape, - degrees=degrees, - ) - - else: - class ExtendedPolyMatrixParamExprImpl(PolyMatrixParamExprImpl): - def re_index( - self, - degree: int, - poly_col: int, - poly_row: int, - x_monomial: tuple[int, ...], - ) -> tuple[int, int, tuple[int, ...], float]: - return re_index(degree, poly_col, poly_row, x_monomial) - - return ExtendedPolyMatrixParamExprImpl( - name=name, - shape=shape, - degrees=degrees, - ) diff --git a/polymatrix/init/initpolymatrixvalue.py b/polymatrix/init/initpolymatrixvalue.py deleted file mode 100644 index 570f399..0000000 --- a/polymatrix/init/initpolymatrixvalue.py +++ /dev/null @@ -1,10 +0,0 @@ -from polymatrix.impl.polymatrixvalueimpl import PolyMatrixValueImpl - -def init_poly_matrix_value( - p_monomial: tuple, - value: float, -): - return PolyMatrixValueImpl( - p_monomial=p_monomial, - value=value, -) diff --git a/polymatrix/init/initscalarmultexpr.py b/polymatrix/init/initscalarmultexpr.py deleted file mode 100644 index a401758..0000000 --- a/polymatrix/init/initscalarmultexpr.py +++ /dev/null @@ -1,11 +0,0 @@ -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.impl.scalarmultexprimpl import ScalarMultExprImpl - -def init_scalar_mult_expr( - left: float, - right: ExpressionBaseMixin, -): - return ScalarMultExprImpl( - left=left, - right=right, -) diff --git a/polymatrix/init/oldinitpolymatrix.py b/polymatrix/init/oldinitpolymatrix.py deleted file mode 100644 index f1b2efd..0000000 --- a/polymatrix/init/oldinitpolymatrix.py +++ /dev/null @@ -1,39 +0,0 @@ -import typing - -from polymatrix.impl.oldpolymatriximpl import PolyMatrixImpl -from polymatrix.init.initexprcontainer import init_expr_container - - -def init_poly_matrix( - name: str, - shape: tuple[int, int], - degrees: list[int] = None, - is_constant: bool = None, - subs: dict[int, dict[tuple[int, int], float]] = None, - re_index: typing.Callable[[int, int, int, tuple[int, ...]], tuple[int, int, int, tuple[int, ...], float]] = None, - is_negated: bool = None, -): - if degrees is None: - assert isinstance(subs, dict) - degrees = list(subs.keys()) - - if is_constant is None: - if subs is None: - is_constant = False - else: - is_constant = True - - if is_negated is None: - is_negated = False - - expr = PolyMatrixImpl( - name=name, - degrees=degrees, - subs=subs, - re_index=re_index, - is_constant=is_constant, - shape=shape, - is_negated=is_negated, - ) - - return init_expr_container(expr=expr) diff --git a/polymatrix/mixins/addexprmixin.py b/polymatrix/mixins/addexprmixin.py deleted file mode 100644 index 4156829..0000000 --- a/polymatrix/mixins/addexprmixin.py +++ /dev/null @@ -1,16 +0,0 @@ -import abc - -from polymatrix.mixins.oldpolymatrixexprmixin import OldPolyMatrixExprMixin - - -class AddExprMixin(OldPolyMatrixExprMixin, abc.ABC): - @property - @abc.abstractmethod - def left(self) -> OldPolyMatrixExprMixin: - ... - - @property - @abc.abstractmethod - def right(self) -> OldPolyMatrixExprMixin: - ... - diff --git a/polymatrix/mixins/exprcontainermixin.py b/polymatrix/mixins/exprcontainermixin.py deleted file mode 100644 index b5dff91..0000000 --- a/polymatrix/mixins/exprcontainermixin.py +++ /dev/null @@ -1,82 +0,0 @@ -import abc -import dataclasses -import typing -from polymatrix.addexpr import AddExpr -from polymatrix.init.initaddexpr import init_add_expr -from polymatrix.init.initmultexpr import init_mult_expr -from polymatrix.init.initscalarmultexpr import init_scalar_mult_expr -from polymatrix.mixins.oldpolymatrixexprmixin import OldPolyMatrixExprMixin - -from polymatrix.multexpr import MultExpr -from polymatrix.scalarmultexpr import ScalarMultExpr - - -ExprType = typing.TypeVar('ExprType', bound=OldPolyMatrixExprMixin) - -class ExprContainerMixin( - typing.Generic[ExprType], - abc.ABC, -): - @property - @abc.abstractmethod - def expr(self) -> ExprType: - ... - - def __add__(self, other: 'ExprContainerMixin[ExprType]'): - return dataclasses.replace( - self, - expr=init_add_expr( - left=self.expr, - right=other.expr, - ), - ) - - def __sub__(self, other: 'ExprContainerMixin[ExprType]'): - return self + (-1.0) * other - - def __mul__(self, other: 'ExprContainerMixin[ExprType]'): - return dataclasses.replace( - self, - expr=init_mult_expr( - left=self.expr, - right=other.expr, - ), - ) - - def __rmul__(self, val: float): - return dataclasses.replace( - self, - expr=init_scalar_mult_expr( - left=val, - right=self.expr, - ), - ) - - def to_list(self) -> tuple[tuple[float, tuple[OldPolyMatrixExprMixin, ...]]]: - def get_mult_expr_list(expr: MultExpr) -> tuple[OldPolyMatrixExprMixin, ...]: - match expr: - case MultExpr(): - left = get_mult_expr_list(expr.left) - right = get_mult_expr_list(expr.right) - return left + right - case _: - return (expr,) - - def get_add_expr_list(expr: AddExpr) -> tuple[tuple[float, tuple[OldPolyMatrixExprMixin, ...]]]: - match expr: - case AddExpr(): - left = get_add_expr_list(expr.left) - right = get_add_expr_list(expr.right) - return left + right - case ScalarMultExpr(): - right = get_mult_expr_list(expr.right) - return ((expr.left, right),) - case MultExpr(): - left = get_mult_expr_list(expr.left) - right = get_mult_expr_list(expr.right) - return ((1.0, left + right),) - case _: - return ((1.0, (expr,)),) - - return get_add_expr_list(self.expr) -
\ No newline at end of file diff --git a/polymatrix/mixins/multexprmixin.py b/polymatrix/mixins/multexprmixin.py deleted file mode 100644 index c9c45c2..0000000 --- a/polymatrix/mixins/multexprmixin.py +++ /dev/null @@ -1,15 +0,0 @@ -import abc - -from polymatrix.mixins.oldpolymatrixexprmixin import OldPolyMatrixExprMixin - - -class MultExprMixin(OldPolyMatrixExprMixin, abc.ABC): - @property - @abc.abstractmethod - def left(self) -> OldPolyMatrixExprMixin: - ... - - @property - @abc.abstractmethod - def right(self) -> OldPolyMatrixExprMixin: - ... diff --git a/polymatrix/mixins/oldpolymatrixexprmixin.py b/polymatrix/mixins/oldpolymatrixexprmixin.py deleted file mode 100644 index 6d5a029..0000000 --- a/polymatrix/mixins/oldpolymatrixexprmixin.py +++ /dev/null @@ -1,2 +0,0 @@ -class OldPolyMatrixExprMixin: - pass diff --git a/polymatrix/mixins/oldpolymatrixexprstatemixin.py b/polymatrix/mixins/oldpolymatrixexprstatemixin.py deleted file mode 100644 index 0b95092..0000000 --- a/polymatrix/mixins/oldpolymatrixexprstatemixin.py +++ /dev/null @@ -1,78 +0,0 @@ -import abc -import itertools -import dataclasses - -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin -from polymatrix.utils import monomial_to_index - - -class OldPolyMatrixExprStateMixin(abc.ABC): - @property - @abc.abstractmethod - def n_var(self) -> int: - """ - dimension of x - """ - - ... - - @property - @abc.abstractmethod - def n_param(self) -> int: - """ - current number of parameters used in polynomial matrix expressions - """ - - ... - - @property - @abc.abstractmethod - def offset_dict(self) -> dict[tuple[OldPolyMatrixMixin, int], tuple[int, int]]: - ... - - # @property - # @abc.abstractmethod - # def local_index_dict(self) -> dict[tuple[PolyMatrixMixin, int], dict[int, int]]: - # ... - - def get_polymat(self, p_index: int) -> tuple[OldPolyMatrixMixin, int, int]: - for (polymat, degree), (start, end) in self.offset_dict.items(): - if start <= p_index < end: - return polymat, degree, p_index - start - - raise Exception(f'index {p_index} not found in offset dictionary') - - # def get_num_param_for(self, polymat: PolyMatrixMixin, degree: int) -> int: - # start_idx, end_idx = self.offset_dict[(polymat, degree)] - # return end_idx - start_idx - - def update_offsets(self, polymats: tuple[OldPolyMatrixMixin]) -> 'OldPolyMatrixExprStateMixin': - registered_polymats = set(polymat for polymat, _ in self.offset_dict.keys()) - parametric_polymats = set(p for p in polymats if not p.is_constant and p not in registered_polymats) - - if len(parametric_polymats) == 0: - return self - - else: - - def gen_n_param_per_struct(): - for polymat in parametric_polymats: - for degree in polymat.degrees: - - # number of terms is given by the maximum index + 1 - number_of_terms = int(polymat.shape[0] * polymat.shape[1] * (monomial_to_index(self.n_var, degree*(self.n_var-1,)) + 1)) - - # print(f'{polymat=}, {number_of_terms=}, {polymat.shape[0] * polymat.shape[1]=}, {monomial_to_index(self.n_var, degree*(self.n_var-1,)) + 1=}') - - yield (polymat, degree), number_of_terms - - param_key, num_of_terms = tuple(zip(*gen_n_param_per_struct())) - cum_sum = tuple(itertools.accumulate((self.n_param,) + num_of_terms)) - offset_dict = dict(zip(param_key, itertools.pairwise(cum_sum))) - - return dataclasses.replace(self, offset_dict=self.offset_dict | offset_dict, n_param=cum_sum[-1]) - - def update_param(self, n): - return dataclasses.replace(self, n_param=self.n_param+n) - - diff --git a/polymatrix/mixins/oldpolymatrixmixin.py b/polymatrix/mixins/oldpolymatrixmixin.py deleted file mode 100644 index dc33449..0000000 --- a/polymatrix/mixins/oldpolymatrixmixin.py +++ /dev/null @@ -1,54 +0,0 @@ -import abc -import dataclasses -import typing - -from matplotlib.pyplot import streamplot - - -class OldPolyMatrixMixin(abc.ABC): - @property - @abc.abstractmethod - def name(self) -> str: - ... - - @property - @abc.abstractmethod - def degrees(self) -> list[int]: - ... - - @property - @abc.abstractmethod - def subs(self) -> dict[int, dict[tuple[int, int], float]]: - ... - - @property - @abc.abstractmethod - def re_index(self) -> typing.Callable[[int, int, int, tuple[int, ...]], tuple[int, int, int, tuple[int, ...], float]]: - ... - - @property - @abc.abstractmethod - def is_constant(self) -> int: - ... - - @property - @abc.abstractmethod - def shape(self) -> tuple[int, int]: - ... - - @property - @abc.abstractmethod - def is_negated(self) -> bool: - ... - - def __neg__(self): - return dataclasses.replace(self, is_negated=not self.is_negated) - - def get_parameter_name(self, degree: int, index: int) -> str: - n_square = self.shape[0] * self.shape[1] - param_index = int(index / n_square) - n_index = index - param_index * n_square - col_index = int(n_index / self.shape[0]) - row_index = n_index - col_index * self.shape[0] - - return f'{self.name}_({degree},{row_index},{col_index},{param_index})' diff --git a/polymatrix/mixins/optimizationmixin.py b/polymatrix/mixins/optimizationmixin.py deleted file mode 100644 index 1a0febe..0000000 --- a/polymatrix/mixins/optimizationmixin.py +++ /dev/null @@ -1,605 +0,0 @@ -import abc -import collections -import dataclasses -import functools -import itertools - -import more_itertools -from sympy import Equivalent -from polymatrix.exprcontainer import ExprContainer -from polymatrix.mixins.oldpolymatrixexprstatemixin import OldPolyMatrixExprStateMixin - -from polymatrix.oldpolymatrix import OldPolyMatrix -from polymatrix.utils import monomial_to_index - - -class OptimizationMixin(abc.ABC): - @property - @abc.abstractmethod - def state(self) -> OldPolyMatrixExprStateMixin: - ... - - @property - @abc.abstractmethod - def equality_constraints(self) -> dict[int, dict[tuple[int, tuple[int, ...], tuple[int, ...]], float]]: - ... - - @property - @abc.abstractmethod - def inequality_constraints(self) -> dict[int, dict[tuple[int, int], float]]: - ... - - @property - @abc.abstractmethod - def auxillary_equality(self) -> dict[int, dict[tuple[int, int], float]]: - ... - - @property - def n_var(self) -> int: - return self.state.n_var - - # def pipe(self, *operators: OptimizationPipeOpMixin): - # return functools.reduce(lambda obs, op: op(obs), operators, self) - - @staticmethod - def _get_equality_terms_from_expression( - expr: ExprContainer[OldPolyMatrix], - state: OldPolyMatrixExprStateMixin, - ) -> dict[tuple[int, tuple[int, ...], tuple[int, ...]], float]: - expr_list = expr.to_list() - - for _, term in expr_list: - for left, right in itertools.pairwise(term): - assert left.shape[1] == right.shape[0], f'{left} and {right} do not match' - - all_polymats = tuple(polymat for _, term in expr_list for polymat in term) - - # update offsets with unseen polymats - state = state.update_offsets(all_polymats) - - equality_constraint = collections.defaultdict(float) - - for _, term in expr_list: - - for degrees in itertools.product(*(polymat.degrees for polymat in term)): - - total_degree = sum(degrees) - - # precompute substitution/offsets for all polynomial matrices in term - d_subs = [polymat.subs[degree] if polymat.subs is not None and degree in polymat.subs else None for degree, polymat in zip(degrees, term)] - d_offsets = [state.offset_dict.get((polymat, degree), (0, 0)) for polymat, degree in zip(term, degrees)] - - # n_var = 2, degree = 3 - # cominations = [(0, 0, 0), (0, 0, 1), (0, 1, 1), (1, 1, 1)] - for combination in itertools.combinations_with_replacement(range(state.n_var), total_degree): - for x_monomial in more_itertools.distinct_permutations(combination): - - def acc_func(acc, v): - last, _ = acc - new = last + v - - return new, x_monomial[last:new] - - # monom=(1, 0, 1) -> monom1=(x2, x1), monom2=(x2) - # the monomial either refers to a free parameter, or to a value in the substitution table - rel_x_monomials = list(monom for _, monom in itertools.accumulate(degrees, acc_func, initial=(0, None)))[1:] - - def non_increasing(seq): - return all(y <= x for x, y in zip(seq, seq[1:])) - - # (1,0) -> x2*x1 instead of (0,1)->x1*x2 - if all(non_increasing(rel_x_monomial) for rel_x_monomial in rel_x_monomials): - - col_defaults = [monomial_to_index(state.n_var, monom) for monom in rel_x_monomials] - - n_rows_col = itertools.chain((m.shape[0] for m in term), (term[-1].shape[1],)) - - for curr_rows_col in itertools.product(*[range(e) for e in n_rows_col]): - - def gen_re_indexing(): #term=term, degrees=degrees, curr_rows_col=curr_rows_col, p_monoms=p_monoms, col_defaults=col_defaults, d_subs=d_subs, d_offsets=d_offsets): - for polymat, degree, (poly_row, poly_col), rel_x_monomial, col_default, subs, offset in zip( - term, degrees, itertools.pairwise(curr_rows_col), rel_x_monomials, col_defaults, d_subs, d_offsets): - - if polymat.re_index is not None: - re_index = polymat.re_index(degree, poly_row, poly_col, rel_x_monomial) - else: - re_index = None - - if re_index is None: - col = col_default - new_poly_row = poly_row - new_poly_col = poly_col - factor = 1.0 - - else: - new_poly_row, new_poly_col, new_monom, factor = re_index - col = monomial_to_index(state.n_var, new_monom) - - if subs is not None: - try: - subs_val = subs[(new_poly_row, new_poly_col, col)] - except KeyError: - subs_val = None - else: - subs_val = None - - if subs_val is None: - # the coefficient is selected after the reindexing - row = new_poly_row + new_poly_col * polymat.shape[0] - - # linearize parameter matrix - p_index = int(offset[0] + row + col * polymat.shape[0] * polymat.shape[1]) - - assert p_index < offset[1], f'{p_index} is bigger than {offset[1]}' - - else: - p_index = None - - yield poly_row, poly_col, p_index, factor, subs_val - - data = tuple(gen_re_indexing()) - - total_factor = functools.reduce(lambda x, y: x*y, (d[3] for d in data)) - - if total_factor == 0: - continue - - value = functools.reduce(lambda x, y: x*y, (d[4] for d in data if d[4] is not None), 1) * total_factor - - if value == 0: - continue - - poly_row = data[0][0] - poly_col = data[-1][1] - p_monomial = tuple(sorted(d[2] for d in data if d[4] is None)) - x_monomial_sorted = tuple(sorted(x_monomial)) - - equality_constraint[poly_row, poly_col, x_monomial_sorted, p_monomial] += value - - return equality_constraint, state - - @staticmethod - def _get_equality_equations_from_terms( - equality_terms: dict[tuple[int, int, tuple[int, ...], tuple[int, ...]], float], - ) -> dict[int, list[tuple[list[int], float]]]: - offset = 0 - equality_constraints = collections.defaultdict(list) - - for _, data in equality_terms.items(): - equation_set = set((poly_row, poly_col, x_monomial) for (poly_row, poly_col, x_monomial, _) in data.keys()) - to_eq_index_dict = {eq: offset + idx for idx, eq in enumerate(equation_set)} - offset += len(equation_set) - - for (poly_row, poly_col, x_monomial, p_monomial), val in data.items(): - eq_index = to_eq_index_dict[(poly_row, poly_col, x_monomial)] - equality_constraints[eq_index] += ((p_monomial, val),) - - return equality_constraints - - def add_equality_constraints(self, expr: ExprContainer[OldPolyMatrix]): - - equality_constraint, state = self._get_equality_terms_from_expression( - expr=expr, - state=self.state - ) - - if len(self.equality_constraints) == 0: - constraint_index = 0 - else: - constraint_index = max(self.equality_constraints.keys()) + 1 - - return dataclasses.replace( - self, - state=state, - equality_constraints=self.equality_constraints | {constraint_index: equality_constraint}, - ) - - def add_inequality_constraints(self, expr: ExprContainer[OldPolyMatrix]): - - # all_polymats = tuple(polymat for term in expr for polymat in term) - polymat = expr.expr - all_polymats = (polymat,) - - # update offsets with unseen polymats - state = self.state.update_offsets(all_polymats) - - aux_eq_buffer = collections.defaultdict(list) - ineq_constr_buffer = collections.defaultdict(list) - - ineq_idx = len(self.inequality_constraints) - aux_eq_idx = len(self.auxillary_equality) - - for degree in polymat.degrees: - # assert degree != 0 - assert degree % 2 == 0 - - vec = tuple(itertools.combinations_with_replacement(range(state.n_var), int(degree/2))) - n_square_mat = len(vec) - - # introduce auxillary parameters - offset_x = state.n_param - - n_param_added = int(n_square_mat*(n_square_mat-1)/2) - state = state.update_param(n_param_added) - - offset_a, end_idx = state.offset_dict[(polymat, degree)] - - def get_param_index_of_poly_mat(row, col): - x_monomial = tuple(sorted(vec[row] + vec[col], reverse=True)) - - if polymat.re_index is None: - reindex = None - else: - reindex = polymat.re_index(degree, 1, 1, x_monomial) - - if reindex is None: - n_monom, factor = x_monomial, 1.0 - - else: - _, _, n_monom, factor = polymat.re_index(degree, 1, 1, vec[row] + vec[col]) - - # assert all(e < end_idx for e in n_monom) - - new_poly_row = 0 - new_poly_col = 0 - col = monomial_to_index(self.n_var, n_monom) - row = new_poly_row + new_poly_col * polymat.shape[0] - - # linearize parameter matrix - p_index = int(offset_a + row + col * polymat.shape[0] * polymat.shape[1]) - - return (p_index,), factor - - for k in range(n_square_mat): - # f in f-v^T@x-r^2 - p_monomial, factor = get_param_index_of_poly_mat(k, k) - ineq_constr_buffer[ineq_idx] += ((p_monomial, factor),) - - for row in range(k): - - # v^T@x in f-v^T@x-r^2 - p_monomial, factor = get_param_index_of_poly_mat(k, row) - ineq_constr_buffer[ineq_idx] += ((p_monomial + (offset_x + row,), -factor),) - - for col in range(k): - # P@x in P@x-v - p_monomial, factor = get_param_index_of_poly_mat(row, col) - aux_eq_buffer[aux_eq_idx] += ((p_monomial + (offset_x + col,), factor),) - - # -v in P@x-v - p_monomial, factor = get_param_index_of_poly_mat(row, k) - aux_eq_buffer[aux_eq_idx] += ((p_monomial, -factor),) - - aux_eq_idx += 1 - - ineq_idx += 1 - offset_x += k - - return dataclasses.replace(self, inequality_constraints=ineq_constr_buffer, auxillary_equality=aux_eq_buffer, state=state) - - def minimize( - self, - cost: tuple[ExprContainer[OldPolyMatrix], ...], - t: float, - ): - """ - - assume sum of squares cost function on variables x - - introduce nu/lambda for each equality/inequality - - differentiate equality and inequality constraints for each variable - - > equality constraints - > inequality constraints - r1 - > x + nu * equality constraint + lambda * inequality constraint - > lambda - r2 - > r1 * r2 - > P @ x = v - """ - - state = [self.state] - - def gen_terms(): - for idx, expr in enumerate(cost): - sum_of_squares, state[0] = self._get_equality_terms_from_expression( - expr=expr, - state=state[0] - ) - yield idx, sum_of_squares - # print(dict(gen_terms())) - sum_of_squares = self._get_equality_equations_from_terms(dict(gen_terms())) - - equality_constraints = self._get_equality_equations_from_terms(self.equality_constraints) - - n_equality_constraints = len(equality_constraints) - - # if len(self.inequality_constraints) == 0: - # n_inequality_constraints = 0 - # else: - # n_inequality_constraints = max(eq for eq, _ in self.inequality_constraints.items()) + 1 - n_inequality_constraints = len(self.inequality_constraints) - - # remove unused parameters - # ------------------------ - - # find used parameters - used_param_indices = set(param_idx for monomial_value_pairs in itertools.chain(equality_constraints.values(), self.inequality_constraints.values(), self.auxillary_equality.values()) for monomial, _ in monomial_value_pairs for param_idx in monomial) - param_update = {monom: idx for idx, monom in enumerate(sorted(used_param_indices))} - - assert max(used_param_indices) == state[0].n_param - 1, f'{max(used_param_indices)=} is not {state[0].n_param - 1=}' - - # param_indices = tuple(start + idx for start, end in state.offset_dict.values() for idx in range(end - start)) - # print(tuple(m for m in monom_update_reverse if m not in param_indices)) - # print(param_indices) - - # the variables x are assumed to be the parameters of all registered polynomial matrices - # todo: this would later come from a specific polynomial matrices - param_indices = tuple(param_update[start + idx] for start, end in state[0].offset_dict.values() for idx in range(end - start) if start + idx in used_param_indices) - - equality_constraints = tuple((eq_index, tuple((tuple(param_update[m] for m in monom), val) for monom, val in monomial_value_pair)) for eq_index, monomial_value_pair in equality_constraints.items()) - inequality_constraints = tuple((key, tuple((tuple(param_update[m] for m in monom), val) for monom, val in monom_val)) for key, monom_val in self.inequality_constraints.items()) - auxillary_equality = tuple((key, tuple((tuple(param_update[m] for m in monom), val) for monom, val in monom_val)) for key, monom_val in self.auxillary_equality.items()) - sum_of_squares = tuple((key, tuple((tuple(param_update[m] for m in monom), val) for monom, val in monom_val)) for key, monom_val in sum_of_squares.items()) - - # introduce variable nu for each equality/inequality - idx_nu = len(used_param_indices) - idx_lambda = idx_nu + n_equality_constraints - idx_r1 = idx_lambda + n_inequality_constraints - idx_r2 = idx_r1 + n_inequality_constraints - idx_dx = [idx_r2 + n_inequality_constraints] - - eq_buffer = collections.defaultdict(list) - - current_eq_offset = 0 - - # > equality constraints - for eq, monom_val_list in equality_constraints: - eq_buffer[current_eq_offset + eq] += monom_val_list - - current_eq_offset += n_equality_constraints - assert max(eq_buffer.keys()) == current_eq_offset - 1 - - # > inequality constraints - r1^2 - for eq, monom_val_list in inequality_constraints: - eq_buffer[current_eq_offset + eq] += monom_val_list + ((2 * (idx_r1 + eq,), -1),) # (tuple(), -0.1)) - - current_eq_offset += n_inequality_constraints - # assert max(eq_buffer.keys()) < current_eq_offset - assert max(eq_buffer.keys()) == current_eq_offset - 1 - - # assert max(eq_buffer.keys()) < current_eq_offset + n_inequality_constraints, f'{max(eq_buffer.keys())=} is bigger than {current_eq_offset + n_inequality_constraints=}' - - dx_dict = collections.defaultdict(dict) - - for _, aux_eq_data in auxillary_equality: - - aux_variables = set(var for p_monomial, _ in aux_eq_data for var in p_monomial if var not in param_indices) - variables = set(var for p_monomial, _ in aux_eq_data for var in p_monomial if var in param_indices) - - # eq1 = (x, a d_a x) - # eq2 = (-1, a d_b x) - for var in variables: - - # a x - b - # ((a, x), (-b,)) - for p_monomial, val in aux_eq_data: - - p_monomial_grp = dict(collections.Counter(p_monomial)) - - if var in p_monomial_grp: - def generate_monom(): - # for each variable in the monomial - for i_var, i_counter in p_monomial_grp.items(): - if var is i_var: - sel_counter = i_counter - 1 - else: - sel_counter = i_counter - - # print(f'{(i_var, sel_counter)=}') - - for _ in range(sel_counter): - yield i_var - - der_monomial = tuple(generate_monom()) - # if len(der_monomial): - # print(f'{der_monomial=}') - eq_buffer[current_eq_offset] += ((der_monomial, val * p_monomial_grp[var]),) - - for aux_var in aux_variables: - if aux_var in p_monomial_grp: - # create d_a x, d_b x, via the product rule - def generate_monom_2(): - # for each variable in the monomial - for i_var, i_counter in p_monomial_grp.items(): - if i_var is aux_var: - sel_counter = i_counter - 1 - else: - sel_counter = i_counter - - for _ in range(sel_counter): - yield i_var - - if var not in dx_dict: - dx_dict[var][aux_var] = idx_dx[0] - idx_dx[0] += 1 - - # product rule - yield dx_dict[var][aux_var] - - der_monomial_2 = tuple(generate_monom_2()) - if len(der_monomial_2): - eq_buffer[current_eq_offset] += ((der_monomial_2, val * p_monomial_grp[aux_var]),) - - current_eq_offset += 1 - - # print(dx_dict) - - def differentiate(data, var, dx_dict): - eq_buffer = tuple() - - for p_monomial, val in data: - - # count powers for each variable - p_monomial_grp = dict(collections.Counter(p_monomial)) - - # for var, counter in p_monomial_grp.items(): - if var in p_monomial_grp: - def generate_monom(): - for i_var, i_counter in p_monomial_grp.items(): - if var is i_var: - sel_counter = i_counter - 1 - else: - sel_counter = i_counter - - for _ in range(sel_counter): - yield i_var - - eq_buffer += (( - tuple(generate_monom()), - val * p_monomial_grp[var], - ),) - - if var in dx_dict: - for aux_var in dx_dict[var]: - if aux_var in p_monomial_grp: - def generate_monom(): - for i_var, i_counter in p_monomial_grp.items(): - if i_var is aux_var: - sel_counter = i_counter - 1 - else: - sel_counter = i_counter - - for _ in range(sel_counter): - yield i_var - - yield dx_dict[var][aux_var] - - eq_buffer += (( - tuple(generate_monom()), - val * p_monomial_grp[aux_var], - ),) - - return eq_buffer - - def reduce_to_quadratic_equation(eq, index): - eq_buffer = tuple() - aux_eq_buffer = tuple() - index_buffer = [index] - - for monomial, val in eq: - n_aux = len(monomial) - 2 - - if 0 < n_aux: - eq_buffer += ((monomial[0:1] + (index_buffer[0],), val),) - - for i_aux in range(n_aux): - if i_aux == n_aux - 1: - aux_monom = monomial[i_aux+1:i_aux+3] - else: - aux_monom = monomial[i_aux+1:i_aux+2] + (index_buffer[0] + 1,) - - aux_eq_buffer += (( - (aux_monom, 1.0), ((index_buffer[0],), -1.0) - ),) - index_buffer[0] += 1 - else: - eq_buffer += ((monomial, val),) - - - return eq_buffer, aux_eq_buffer, index_buffer[0] - - - # > x + nu * equality constraint + lambda * inequality constraint - # --------------------------------------------------------------- - - # print(sum_of_squares) - - for index, param in enumerate(param_indices): - # print(param) - def gen_cost_terms(): - for _, data in sum_of_squares: - equation = differentiate(data, param, {}) - for d_monomials, d_val in equation: - for monomials, val in data: - yield monomials + d_monomials, val * d_val - - eq_buffer[current_eq_offset + index] += tuple(gen_cost_terms()) - - # eq_buffer[current_eq_offset + index] += (((param,), 1),) - - # differentiate equality constraints for each variable x - for eq_constr_idx, eq_data in equality_constraints: - for var in param_indices: - - equation = differentiate(eq_data, var, dx_dict) - equation = tuple((monomials + (idx_nu + eq_constr_idx,), val) for monomials, val in equation) - equation, aux_equations, idx_dx[0] = reduce_to_quadratic_equation(equation, idx_dx[0]) - - if len(equation): - eq_idx = param_indices.index(var) - eq_buffer[current_eq_offset + eq_idx] += equation - - next_current_eq_offset = current_eq_offset + len(param_indices) - - # differentiate inequality constraints for each variable x - for ineq_constr_idx, ineq_data in inequality_constraints: - - for var in param_indices: - - equation = differentiate(ineq_data, var, dx_dict) - equation = tuple((monomials + (idx_lambda + ineq_constr_idx,), -val) for monomials, val in equation) - equation, aux_equations, idx_dx[0] = reduce_to_quadratic_equation(equation, idx_dx[0]) - - if len(equation): - eq_idx = param_indices.index(var) - # print(f'{equations=}') - eq_buffer[current_eq_offset + eq_idx] += equation - - # print(f'{aux_equations=}') - for aux_eq in aux_equations: - eq_buffer[next_current_eq_offset] = aux_eq - next_current_eq_offset += 1 - - # current_eq_offset += len(param_indices) - current_eq_offset = next_current_eq_offset - - # assert max(eq_buffer.keys()) < current_eq_offset, f'{max(eq_buffer.keys())=} is bigger than {current_eq_offset=}' - assert max(eq_buffer.keys()) == current_eq_offset - 1 - - # > lambda - r2^2 - for idx in range(n_inequality_constraints): - eq_buffer[current_eq_offset + idx] += (((idx_lambda + idx,), 1), (2 * (idx_r2 + idx,), -1)) - - current_eq_offset += n_inequality_constraints - # assert max(eq_buffer.keys()) < current_eq_offset - assert max(eq_buffer.keys()) == current_eq_offset - 1 - - next_current_eq_offset = current_eq_offset + n_inequality_constraints - # # > r1 * r2 - # for idx in range(n_inequality_constraints): - # eq_buffer[current_eq_offset + idx] += (((idx_r1 + idx, idx_r2 + idx), 1),) - - # > r1^2 * lambda - for idx in range(n_inequality_constraints): - eq_buffer[next_current_eq_offset] += (((idx_dx[0],), 1), (2 * (idx_r1 + idx,), -1)) - next_current_eq_offset += 1 - - eq_buffer[current_eq_offset + idx] += (((idx_lambda + idx, idx_dx[0]), 1), (tuple(), -1/t)) - idx_dx[0] += 1 - - # current_eq_offset += n_inequality_constraints - current_eq_offset = next_current_eq_offset - - # assert max(eq_buffer.keys()) < current_eq_offset - assert max(eq_buffer.keys()) == current_eq_offset - 1 - - # > P @ x - v - for eq, data in auxillary_equality: - eq_buffer[current_eq_offset + eq] += data - - current_eq_offset += len(used_param_indices) - len(param_indices) - # print(f'{current_eq_offset=}') - # print(f'{len(monom_update_reverse) - len(param_indices)=}') - assert max(eq_buffer.keys()) == current_eq_offset - 1, f'{max(eq_buffer.keys())} is not {current_eq_offset - 1}' - - # print(f'{current_eq_offset=}') - - return eq_buffer, param_update, current_eq_offset diff --git a/polymatrix/mixins/scalarmultexprmixin.py b/polymatrix/mixins/scalarmultexprmixin.py deleted file mode 100644 index ec98722..0000000 --- a/polymatrix/mixins/scalarmultexprmixin.py +++ /dev/null @@ -1,14 +0,0 @@ -import abc -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin - - -class ScalarMultExprMixin(ExpressionBaseMixin, abc.ABC): - @property - @abc.abstractmethod - def left(self) -> float: - ... - - @property - @abc.abstractmethod - def right(self) -> ExpressionBaseMixin: - ... diff --git a/polymatrix/mixins/statemonadmixin.py b/polymatrix/mixins/statemonadmixin.py deleted file mode 100644 index c90ebc8..0000000 --- a/polymatrix/mixins/statemonadmixin.py +++ /dev/null @@ -1,18 +0,0 @@ -import abc -import typing - - -StateType = typing.TypeVar('StateType') -ValueType = typing.TypeVar('ValueType') - -class StateMonadMixin( - # typing.Generic[StateType, ValueType], - abc.ABC, -): - @abc.abstractmethod - def apply(self, state: StateType) -> tuple[StateType, ValueType]: - ... - - # def __apply__(self, state: StateType) -> ValueType: - # _, value = self.state_func(state) - # return value diff --git a/polymatrix/multexpr.py b/polymatrix/multexpr.py deleted file mode 100644 index 568bbf8..0000000 --- a/polymatrix/multexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.mixins.multexprmixin import MultExprMixin - -class MultExpr(MultExprMixin): - pass diff --git a/polymatrix/oldpolymatrix.py b/polymatrix/oldpolymatrix.py deleted file mode 100644 index 409680d..0000000 --- a/polymatrix/oldpolymatrix.py +++ /dev/null @@ -1,6 +0,0 @@ -from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin -from polymatrix.mixins.oldpolymatrixmixin import OldPolyMatrixMixin - - -class OldPolyMatrix(OldPolyMatrixMixin): - pass diff --git a/polymatrix/oldpolymatrixexprstate.py b/polymatrix/oldpolymatrixexprstate.py deleted file mode 100644 index 1397f9e..0000000 --- a/polymatrix/oldpolymatrixexprstate.py +++ /dev/null @@ -1,5 +0,0 @@ -from polymatrix.mixins.oldpolymatrixexprstatemixin import OldPolyMatrixExprStateMixin - - -class OldPolyMatrixExprState(OldPolyMatrixExprStateMixin): - pass diff --git a/polymatrix/optimization.py b/polymatrix/optimization.py deleted file mode 100644 index c5034a1..0000000 --- a/polymatrix/optimization.py +++ /dev/null @@ -1,5 +0,0 @@ -from polymatrix.mixins.optimizationmixin import OptimizationMixin - - -class Optimization(OptimizationMixin): - pass diff --git a/polymatrix/polymatrixaddexpr.py b/polymatrix/polymatrixaddexpr.py deleted file mode 100644 index 5d4f567..0000000 --- a/polymatrix/polymatrixaddexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.additionexprmixin import AddExprMixin - -class PolyMatrixAddExpr(AddExprMixin): - pass diff --git a/polymatrix/polymatrixarrayexpr.py b/polymatrix/polymatrixarrayexpr.py deleted file mode 100644 index c9b572d..0000000 --- a/polymatrix/polymatrixarrayexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.fromarrayexprmixin import FromArrayExprMixin - -class PolyMatrixArrayExpr(FromArrayExprMixin): - pass diff --git a/polymatrix/polymatrixexpr.py b/polymatrix/polymatrixexpr.py deleted file mode 100644 index 74e8757..0000000 --- a/polymatrix/polymatrixexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.expressionmixin import ExpressionMixin - -class PolyMatrixExpr(ExpressionMixin): - pass diff --git a/polymatrix/polymatrixexprstate.py b/polymatrix/polymatrixexprstate.py deleted file mode 100644 index 11419cb..0000000 --- a/polymatrix/polymatrixexprstate.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.expressionstatemixin import ExpressionStateMixin - -class PolyMatrixExprState(ExpressionStateMixin): - pass diff --git a/polymatrix/polymatrixparamexpr.py b/polymatrix/polymatrixparamexpr.py deleted file mode 100644 index f396bc9..0000000 --- a/polymatrix/polymatrixparamexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.expression.mixins.parametrizetermsexprmixin import ParametrizeTermsExprMixin - -class PolyMatrixParamExpr(ParametrizeTermsExprMixin): - pass diff --git a/polymatrix/polysolver.py b/polymatrix/polysolver.py index 2cb9b7b..da6fd5d 100644 --- a/polymatrix/polysolver.py +++ b/polymatrix/polysolver.py @@ -9,6 +9,48 @@ import more_itertools import time +def convert_to_complex_equations(raw_data): + n = raw_data[0].shape[0] + + data = {} + data[0] = np.vstack((raw_data[0], np.zeros((n, 1)))) + data[1] = np.kron(np.eye(2), raw_data[1]) + + def get_tuple(n, idx): + row = idx % n + col = int((idx-row) / n) + return (row, col) + + def get_index(n, tple): + return sum(idx*(n**level) for level, idx in enumerate(tple)) + + csr_data = raw_data[2].tocsr() + + def gen_complex_data(): + for row in range(csr_data.shape[0]): + + pt = slice(csr_data.indptr[row], csr_data.indptr[row+1]) + + for col, value in zip(csr_data.indices[pt], csr_data.data[pt]): + + left, right = get_tuple(n, col) + + yield row, get_index(2*n, (left, right)), value + yield row, get_index(2*n, (left+n, right+n)), -value + + yield row+n, get_index(2*n, (left, right+n)), value + yield row+n, get_index(2*n, (left+n, right)), value + + row, col, idata = list(zip(*gen_complex_data())) + + data[2] = scipy.sparse.coo_array( + (idata, (row, col)), + shape=(2*n, (2*n)**2), + ) + + return data + + def substitude_x_add_a(data, a): def kron_eye_and_a(n, data, a, indices): @@ -264,15 +306,18 @@ def inner_smart_solve(data, irange=None, idegree=None, a_init=None): return a, err -def outer_smart_solve(data, n_iter=10, a_max=1.0, irange=None, idegree=None): - def gen_a_err(): - n = data[1].shape[0] +def outer_smart_solve(data, a_init=None, n_iter=10, a_max=1.0, irange=None, idegree=None): + n = data[1].shape[0] + if a_init is None: + a_init = np.zeros(n) + + def gen_a_err(): for _ in range(n_iter): - a_init = a_max * (2 * np.random.rand(n) - 1) + a_subs = a_init + a_max * (2 * np.random.rand(n) - 1) try: - a, err = inner_smart_solve(data, irange=irange, idegree=idegree, a_init=a_init) + a, err = inner_smart_solve(data, irange=irange, idegree=idegree, a_init=a_subs) subs_data = substitude_x_add_a(data, a[-1]) sol = solve_poly_system(subs_data, 6) diff --git a/polymatrix/scalarmultexpr.py b/polymatrix/scalarmultexpr.py deleted file mode 100644 index 3d9a075..0000000 --- a/polymatrix/scalarmultexpr.py +++ /dev/null @@ -1,4 +0,0 @@ -from polymatrix.mixins.scalarmultexprmixin import ScalarMultExprMixin - -class ScalarMultExpr(ScalarMultExprMixin): - pass diff --git a/polymatrix/statemonad.py b/polymatrix/statemonad.py deleted file mode 100644 index 3054160..0000000 --- a/polymatrix/statemonad.py +++ /dev/null @@ -1,33 +0,0 @@ -from typing import Callable, Tuple, Any, TypeVar, Generic - -State = TypeVar('State') -U = TypeVar('U') -V = TypeVar('V') - - -class StateMonad(Generic[U, State]): - def __init__(self, fn: Callable[[State], Tuple[U, State]]) -> None: - self._fn = fn - - @classmethod - def unit(cls, value: Any) -> 'StateMonad[U, State]': - return cls(lambda state: (value, state)) - - def map(self, fn: Callable[[U], V]) -> 'StateMonad[V, State]': - - def internal_map(state: State) -> Tuple[U, State]: - val, n_state = self._fn(state) - return fn(val), n_state - - return StateMonad(internal_map) - - def flat_map(self, fn: Callable[[U], 'StateMonad']) -> 'StateMonad[V, State]': - - def internal_map(state: State) -> Tuple[V, State]: - val, n_state = self._fn(state) - return fn(val).run(n_state) - - return StateMonad(internal_map) - - def run(self, state: Any) -> Tuple[Any, Any]: - return self._fn(state) |