diff options
Diffstat (limited to '')
| -rw-r--r-- | .gitignore | 110 | ||||
| -rw-r--r-- | LICENSE | 191 | ||||
| -rw-r--r-- | README.md | 1 | ||||
| -rw-r--r-- | main.py | 101 | ||||
| -rw-r--r-- | polymatrix/__init__.py | 0 | ||||
| -rw-r--r-- | polymatrix/polysolver.py | 41 | ||||
| -rw-r--r-- | polymatrix/polystruct.py | 358 | ||||
| -rw-r--r-- | polymatrix/sympyutils.py | 82 | ||||
| -rw-r--r-- | polymatrix/utils.py | 34 | ||||
| -rw-r--r-- | setup.py | 5 | ||||
| -rw-r--r-- | test_polymatrix/test_polymatrix.py | 176 |
11 files changed, 1099 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d7f42da --- /dev/null +++ b/.gitignore @@ -0,0 +1,110 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +.hypothesis/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +.static_storage/ +.media/ +local_settings.py + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# celery beat schedule file +celerybeat-schedule + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ + +.idea/ + +.vscode +.pytest_cache + @@ -0,0 +1,191 @@ +Apache License +Version 2.0, January 2004 +http://www.apache.org/licenses/ + +TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + +1. Definitions. + +"License" shall mean the terms and conditions for use, reproduction, and +distribution as defined by Sections 1 through 9 of this document. + +"Licensor" shall mean the copyright owner or entity authorized by the copyright +owner that is granting the License. + +"Legal Entity" shall mean the union of the acting entity and all other entities +that control, are controlled by, or are under common control with that entity. +For the purposes of this definition, "control" means (i) the power, direct or +indirect, to cause the direction or management of such entity, whether by +contract or otherwise, or (ii) ownership of fifty percent (50%) or more of the +outstanding shares, or (iii) beneficial ownership of such entity. + +"You" (or "Your") shall mean an individual or Legal Entity exercising +permissions granted by this License. + +"Source" form shall mean the preferred form for making modifications, including +but not limited to software source code, documentation source, and configuration +files. + +"Object" form shall mean any form resulting from mechanical transformation or +translation of a Source form, including but not limited to compiled object code, +generated documentation, and conversions to other media types. + +"Work" shall mean the work of authorship, whether in Source or Object form, made +available under the License, as indicated by a copyright notice that is included +in or attached to the work (an example is provided in the Appendix below). + +"Derivative Works" shall mean any work, whether in Source or Object form, that +is based on (or derived from) the Work and for which the editorial revisions, +annotations, elaborations, or other modifications represent, as a whole, an +original work of authorship. For the purposes of this License, Derivative Works +shall not include works that remain separable from, or merely link (or bind by +name) to the interfaces of, the Work and Derivative Works thereof. + +"Contribution" shall mean any work of authorship, including the original version +of the Work and any modifications or additions to that Work or Derivative Works +thereof, that is intentionally submitted to Licensor for inclusion in the Work +by the copyright owner or by an individual or Legal Entity authorized to submit +on behalf of the copyright owner. For the purposes of this definition, +"submitted" means any form of electronic, verbal, or written communication sent +to the Licensor or its representatives, including but not limited to +communication on electronic mailing lists, source code control systems, and +issue tracking systems that are managed by, or on behalf of, the Licensor for +the purpose of discussing and improving the Work, but excluding communication +that is conspicuously marked or otherwise designated in writing by the copyright +owner as "Not a Contribution." + +"Contributor" shall mean Licensor and any individual or Legal Entity on behalf +of whom a Contribution has been received by Licensor and subsequently +incorporated within the Work. + +2. Grant of Copyright License. + +Subject to the terms and conditions of this License, each Contributor hereby +grants to You a perpetual, worldwide, non-exclusive, no-charge, royalty-free, +irrevocable copyright license to reproduce, prepare Derivative Works of, +publicly display, publicly perform, sublicense, and distribute the Work and such +Derivative Works in Source or Object form. + +3. Grant of Patent License. + +Subject to the terms and conditions of this License, each Contributor hereby +grants to You a perpetual, worldwide, non-exclusive, no-charge, royalty-free, +irrevocable (except as stated in this section) patent license to make, have +made, use, offer to sell, sell, import, and otherwise transfer the Work, where +such license applies only to those patent claims licensable by such Contributor +that are necessarily infringed by their Contribution(s) alone or by combination +of their Contribution(s) with the Work to which such Contribution(s) was +submitted. If You institute patent litigation against any entity (including a +cross-claim or counterclaim in a lawsuit) alleging that the Work or a +Contribution incorporated within the Work constitutes direct or contributory +patent infringement, then any patent licenses granted to You under this License +for that Work shall terminate as of the date such litigation is filed. + +4. Redistribution. + +You may reproduce and distribute copies of the Work or Derivative Works thereof +in any medium, with or without modifications, and in Source or Object form, +provided that You meet the following conditions: + +You must give any other recipients of the Work or Derivative Works a copy of +this License; and +You must cause any modified files to carry prominent notices stating that You +changed the files; and +You must retain, in the Source form of any Derivative Works that You distribute, +all copyright, patent, trademark, and attribution notices from the Source form +of the Work, excluding those notices that do not pertain to any part of the +Derivative Works; and +If the Work includes a "NOTICE" text file as part of its distribution, then any +Derivative Works that You distribute must include a readable copy of the +attribution notices contained within such NOTICE file, excluding those notices +that do not pertain to any part of the Derivative Works, in at least one of the +following places: within a NOTICE text file distributed as part of the +Derivative Works; within the Source form or documentation, if provided along +with the Derivative Works; or, within a display generated by the Derivative +Works, if and wherever such third-party notices normally appear. The contents of +the NOTICE file are for informational purposes only and do not modify the +License. You may add Your own attribution notices within Derivative Works that +You distribute, alongside or as an addendum to the NOTICE text from the Work, +provided that such additional attribution notices cannot be construed as +modifying the License. +You may add Your own copyright statement to Your modifications and may provide +additional or different license terms and conditions for use, reproduction, or +distribution of Your modifications, or for any such Derivative Works as a whole, +provided Your use, reproduction, and distribution of the Work otherwise complies +with the conditions stated in this License. + +5. Submission of Contributions. + +Unless You explicitly state otherwise, any Contribution intentionally submitted +for inclusion in the Work by You to the Licensor shall be under the terms and +conditions of this License, without any additional terms or conditions. +Notwithstanding the above, nothing herein shall supersede or modify the terms of +any separate license agreement you may have executed with Licensor regarding +such Contributions. + +6. Trademarks. + +This License does not grant permission to use the trade names, trademarks, +service marks, or product names of the Licensor, except as required for +reasonable and customary use in describing the origin of the Work and +reproducing the content of the NOTICE file. + +7. Disclaimer of Warranty. + +Unless required by applicable law or agreed to in writing, Licensor provides the +Work (and each Contributor provides its Contributions) on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied, +including, without limitation, any warranties or conditions of TITLE, +NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A PARTICULAR PURPOSE. You are +solely responsible for determining the appropriateness of using or +redistributing the Work and assume any risks associated with Your exercise of +permissions under this License. + +8. Limitation of Liability. + +In no event and under no legal theory, whether in tort (including negligence), +contract, or otherwise, unless required by applicable law (such as deliberate +and grossly negligent acts) or agreed to in writing, shall any Contributor be +liable to You for damages, including any direct, indirect, special, incidental, +or consequential damages of any character arising as a result of this License or +out of the use or inability to use the Work (including but not limited to +damages for loss of goodwill, work stoppage, computer failure or malfunction, or +any and all other commercial damages or losses), even if such Contributor has +been advised of the possibility of such damages. + +9. Accepting Warranty or Additional Liability. + +While redistributing the Work or Derivative Works thereof, You may choose to +offer, and charge a fee for, acceptance of support, warranty, indemnity, or +other liability obligations and/or rights consistent with this License. However, +in accepting such obligations, You may act only on Your own behalf and on Your +sole responsibility, not on behalf of any other Contributor, and only if You +agree to indemnify, defend, and hold each Contributor harmless for any liability +incurred by, or claims asserted against, such Contributor by reason of your +accepting any such warranty or additional liability. + +END OF TERMS AND CONDITIONS + +APPENDIX: How to apply the Apache License to your work + +To apply the Apache License to your work, attach the following boilerplate +notice, with the fields enclosed by brackets "[]" replaced with your own +identifying information. (Don't include the brackets!) The text should be +enclosed in the appropriate comment syntax for the file format. We also +recommend that a file or class name and description of purpose be included on +the same "printed page" as the copyright notice for easier identification within +third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/README.md b/README.md new file mode 100644 index 0000000..3ba3ae5 --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +# Polynomial matrix library @@ -0,0 +1,101 @@ +import sympy +import scipy.special +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 + +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) + +def skew_symmetric(idx1, idx2): + if idx2 < idx1: + return idx2, idx1, -1 + +def zero_diagonal_for_uneven_degree(degree, idx1, idx2): + if degree % 2 == 1 and idx1 == idx2: + return 0 + + + +nabla_ha = init_poly_vector(degrees=(1,2)) +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)) + +# 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(n_var*binom(n_var+1-1, 1)) +# print(n_var*binom(n_var+2-1, 2)) + +# print(binom(total+2-1, 2)) + +# eq = init_equation( +# # terms = [(JR, nabla_ha), (JRa, nabla_ha), (JRa, nabla_h), (mg, u)], +# # n_var = 3, +# terms = [(JR, nabla_ha)], +# n_var = 2, +# ) + +# print(eq.create()) + +mat = init_poly_matrix(degrees=(0,), re_index_func=skew_symmetric) +# vec = init_poly_vector(degrees=(0,1)) +vec = init_poly_vector(subs={0: {(0, 0): 1, (1, 0): 1}}) + +eq = init_equation( + terms = [(mat, vec)], + n_var = 2, +) +eq_tuples, offset_dict = eq.create() + +# print(eq_tuples) +# print(eq_tuples[1][(1, tuple())]) +# print(offset_dict[(mat, 0)]) +# print(offset_dict[(vec, 0)]) +print(offset_dict) +# # print(variable_to_index(10, (offset_dict[(mat, 0)], offset_dict[(vec, 1)]))) diff --git a/polymatrix/__init__.py b/polymatrix/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/polymatrix/__init__.py diff --git a/polymatrix/polysolver.py b/polymatrix/polysolver.py new file mode 100644 index 0000000..f611260 --- /dev/null +++ b/polymatrix/polysolver.py @@ -0,0 +1,41 @@ +import numpy as np +import itertools + +def quadratic_poly_solver(b, c, m): + n = b.shape[0] + nn = n*n + + b_inv = np.linalg.inv(b) + + p = b_inv @ np.ones((n, 1)) + + def square_matrix(p): + return (p @ p.T).reshape(nn, 1) + + + def merge_matrix(p1, p2): + return (np.kron(np.eye(n), p1) + np.kron(p1, np.eye(n))) @ p2 + + + def func(acc, _): + p = acc + + k = p.shape[1] + + uneven = k % 2 + half = int((k - uneven) / 2) + + def gen_pvector(): + if uneven: + yield square_matrix(p[:, half:half+1]) + + for i in range(half): + yield merge_matrix(p[:, i:i+1], p[:, k-i-1:k-i]) + + p_vector = np.add.reduce(list(gen_pvector())) + p_new = b_inv @ c @ p_vector + return np.concatenate((p, p_new), axis=1) + + *_, sol_subs = itertools.accumulate(range(m-1), func, initial=p) + + return np.sum(-sol_subs, axis=1, keepdims=True) diff --git a/polymatrix/polystruct.py b/polymatrix/polystruct.py new file mode 100644 index 0000000..69a05d6 --- /dev/null +++ b/polymatrix/polystruct.py @@ -0,0 +1,358 @@ +import abc +import collections +import typing +import numpy as np +import dataclasses +import dataclass_abc +import scipy.special +import itertools + +from polymatrix.utils import variable_to_index + +DegreeType = int +# CoordType = tuple[int, int] +# SubsDictType = dict[DegreeType, dict[CoordType, float]] + +######################################## +# Mixins +######################################## + +class PolyMatrixMixin(abc.ABC): + @property + @abc.abstractmethod + def degrees(self) -> list[int]: + ... + + @property + @abc.abstractmethod + def subs(self) -> dict[DegreeType, dict[int, int, float]]: + ... + + @property + @abc.abstractmethod + def subs_func(self) -> typing.Callable[[DegreeType, int, int], tuple[float]]: + ... + + @property + @abc.abstractmethod + def re_index(self) -> dict[DegreeType, dict[int, int, tuple[int, int, float]]]: + ... + + @property + @abc.abstractmethod + def re_index_func(self) -> typing.Callable[[int, int], tuple[int, int, float]]: + ... + + @property + @abc.abstractmethod + def re_index_func_2(self) -> typing.Callable[[int, int], tuple[int, int, float]]: + ... + + @property + @abc.abstractmethod + def is_vector(self) -> bool: + ... + + +class EquationMixin(abc.ABC): + @property + @abc.abstractmethod + def terms(self) -> list[tuple[PolyMatrixMixin, PolyMatrixMixin]]: + ... + + @property + @abc.abstractmethod + def n_var(self) -> int: + ... + + def create( + self, + subs: dict[PolyMatrixMixin, dict[DegreeType, dict[int, int, float]]] = None, + ): + if subs is None: + subs = {} + + binom = scipy.special.binom + + # create parameter offset + all_structs = set(indexed_poly_mat for term in self.terms for indexed_poly_mat in term) + + def gen_n_param_per_struct(): + acc = 0 + + for struct in all_structs: + added_subs = subs.get(struct, {}) + + for degree in struct.degrees: + + if struct.subs is not None and degree in struct.subs: + struct_subs = struct.subs[degree] + + elif struct.subs_func is not None: + def gen_coords(): + for eq in range(self.n_var): + + if struct.is_vector: + col_indices = (0,) + else: + col_indices = range(self.n_var) + + for col in col_indices: + coord = (eq, col) + + subs_val = struct.subs_func(*coord) + + if subs_val is not None: + yield coord, subs_val + + struct_subs = dict(gen_coords()) + + else: + struct_subs = {} + + if struct.re_index is not None and degree in struct.re_index: + struct_re_index = struct.re_index[degree] + + elif struct.re_index_func is not None: + def gen_coords(): + for eq in range(self.n_var): + + if struct.is_vector: + col_indices = (0,) + else: + col_indices = range(self.n_var) + + for col in col_indices: + coord = (eq, col) + + re_index_val = struct.re_index_func(*coord) + + match re_index_val: + case (row, col, value): + yield coord, ((row, col), value) + struct_re_index = dict(gen_coords()) + + else: + struct_re_index = {} + + all_subs = struct_subs | added_subs.get(degree, {}) + + yield (struct, degree), acc, (all_subs, struct_re_index) + + # def gen_predefined_coord(): + # re_index = re_index_dict[struct].get(degree, {}) + + # for coord_dict in (all_subs, re_index): + # for coord in coord_dict.keys(): + # # in case the n_var is smaller than expected + # if max(coord) <= self.n_var: + # yield coord + + # n_predefined = len(set(gen_predefined_coord())) + + if struct.is_vector: + number_of_terms = int(self.n_var * binom(self.n_var+degree-1, degree)) + else: + number_of_terms = int(self.n_var**2 * binom(self.n_var+degree-1, degree)) + + acc += number_of_terms + + yield None, acc, None + + param_list = list(gen_n_param_per_struct()) + n_param = param_list[-1][1] + offset_dict = dict(((e[0], e[1]) for e in param_list[:-1])) + subs_dict = dict(((e[0], e[2][0]) for e in param_list[:-1])) + re_index_dict = dict(((e[0], e[2][1]) for e in param_list[:-1])) + + terms = collections.defaultdict(lambda: collections.defaultdict(lambda: collections.defaultdict(int))) + + for left, right in self.terms: + for d1 in left.degrees: + + n_param_1 = binom(self.n_var+d1-1, d1) + offset_1 = offset_dict[(left, d1)] + subs_1 = subs_dict[(left, d1)] + re_index_1 = re_index_dict[(left, d1)] + + for d2 in right.degrees: + + n_param_2 = binom(self.n_var+d2-1, d2) + offset_2 = offset_dict[(right, d2)] + subs_2 = subs_dict[(right, d2)] + re_index_2 = re_index_dict[(right, d2)] + + total_degree = d1 + d2 + + for combination in itertools.combinations_with_replacement(range(self.n_var), total_degree): + # n_var = 2, degree = 3 + # cominations = [(0, 0, 0), (0, 0, 1), (0, 1, 1), (1, 1, 1)] + + # def gen_coord_for_single_equation(): + for perm in set(itertools.permutations(combination)): + + # grp1 = (1, 1, 0) -> x1 x2^2 + grp1 = perm[0:d1] + grp2 = perm[d1:] + + 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 non_increasing(grp1) and non_increasing(grp2): + + left_idx = variable_to_index(self.n_var, grp1) + right_idx = variable_to_index(self.n_var, grp2) + + # for each column of the poly matrix, and row of the poly vector + for idx_col in range(self.n_var): + + coord_2 = (idx_col, 0) + n_coord_2, factor_2 = re_index_2.get(coord_2, (coord_2, 1)) + + if factor_2 == 0: + continue + + if n_coord_2 in subs_2: + subs_val_2 = factor_2 * subs_2[n_coord_2] + if subs_val_2 == 0: + continue + else: + subs_val_2 = None + + # for each polynomial equation + for idx_eq in range(self.n_var): + + coord_1 = (idx_eq, idx_col) + n_coord_1, factor_1 = re_index_1.get(coord_1, (coord_1, 1)) + + if factor_1 == 0: + continue + + if n_coord_1 in subs_1: + subs_val_1 = factor_1 * subs_1[n_coord_1] + if subs_val_1 == 0: + continue + else: + subs_val_1 = None + + left_param_idx = int(offset_1 + left_idx + (self.n_var * n_coord_1[0] + n_coord_1[1]) * n_param_1) + right_param_idx = int(offset_2 + right_idx + n_coord_2[0] * n_param_2) + + match (subs_val_1, subs_val_2): + case (None, None): + row_idx = variable_to_index(n_param, (left_param_idx, right_param_idx)) + degree = 2 + value = factor_1*factor_2 + + case (subs_val, None): + # row_idx = variable_to_index(n_param, (right_param_idx,)) + row_idx = right_param_idx + degree = 1 + value = subs_val_1*factor_2 + + case (None, subs_val): + # row_idx = variable_to_index(n_param, (left_param_idx,)) + row_idx = left_param_idx + degree = 1 + value = subs_val_2*factor_1 + + case _: + degree, row_idx, value = 0, 0, subs_val_1*subs_val_2 + + terms[degree][(idx_eq, perm)][row_idx] += value + + return terms, offset_dict + +######################################## +# Classes +######################################## + +class PolyMatrix(PolyMatrixMixin): + pass + + +class Equation(EquationMixin): + pass + +######################################## +# Implementations +######################################## + +@dataclass_abc.dataclass_abc(frozen=True, eq=False) +class PolyMatrixImpl(PolyMatrix): + degrees: list[int] + subs: dict[int, dict[tuple[int, int], float]] + subs_func: typing.Callable[[int, int, int], tuple[float]] + re_index: dict[int, dict[tuple[int, int], tuple[int, int, float]]] + re_index_func: typing.Callable[[int, int], tuple[int, int, float]] + is_vector: bool + + +@dataclass_abc.dataclass_abc(frozen=True) +class EquationImpl(Equation): + terms: list[tuple[PolyMatrix, PolyMatrix]] + n_var: int + +######################################## +# init functions +######################################## + +def init_poly_vector( + degrees: list[int] = None, + subs: dict[int, dict[tuple[int, int], float]] = None, + subs_func: typing.Callable[[int, int, int], tuple[float]] = None, + re_index: dict[int, dict[tuple[int, int], tuple[int, int, float]]] = None, + re_index_func: typing.Callable[[int, int], tuple[int, int, float]] = None, +): + return init_poly_matrix( + degrees=degrees, + subs=subs, + subs_func=subs_func, + re_index=re_index, + re_index_func=re_index_func, + is_vector=True, + ) + + +def init_poly_matrix( + degrees: list[int] = None, + subs: dict[int, dict[tuple[int, int], float]] = None, + subs_func: typing.Callable[[int, int, int], tuple[float]] = None, + re_index: dict[int, dict[tuple[int, int], tuple[int, int, float]]] = None, + re_index_func: typing.Callable[[int, int], tuple[int, int, float]] = None, + is_vector: bool = None, +): + if degrees is None: + assert isinstance(subs, dict) + degrees = list(subs.keys()) + + if subs is None and subs_func is None: + subs = {} + + if re_index is None and re_index_func is None: + re_index = {} + + if is_vector is None: + is_vector = False + + return PolyMatrixImpl( + degrees=degrees, + subs=subs, + subs_func=subs_func, + re_index = re_index, + re_index_func = re_index_func, + is_vector = is_vector, + ) + + +def init_equation( + n_var: int, + terms: list[tuple[PolyMatrix, PolyMatrix]], +): + assert all(not left.is_vector and right.is_vector for left, right in terms) + + return EquationImpl( + n_var=n_var, + terms=terms, + ) diff --git a/polymatrix/sympyutils.py b/polymatrix/sympyutils.py new file mode 100644 index 0000000..e2c37d8 --- /dev/null +++ b/polymatrix/sympyutils.py @@ -0,0 +1,82 @@ +import itertools +import numpy as np +import scipy.sparse + +from polymatrix.utils import variable_powers_to_index + + +def poly_to_data_coord(poly_list, power = None): + """ + poly_list = [ + poly(x1*x3**2, x) + ] + power: up to which power + """ + + if power is None: + power = max(degree for poly in poly_list for degree in poly.degree_list()) + + def all_equal(iterable): + g = itertools.groupby(iterable) + return next(g, True) and not next(g, False) + + assert all_equal((p.gens for p in poly_list)), 'all polynomials need to have identical generators' + + def gen_power_mat(): + + # for all powers generate a matrix + for current_power in range(power + 1): + + def gen_value_index(): + + # each polynomial defines a row in the matrix + for row, p in enumerate(poly_list): + + # a5 x1 x3**2 -> c=a5, m=(1, 0, 2) + for c, m in zip(p.coeffs(), p.monoms()): + + if sum(m) == current_power: + + index = variable_powers_to_index(m) + yield (row, index), c + + # yield list(zip(*gen_value_index())) + data = dict(gen_value_index()) + + if len(data) > 0: + yield current_power, data + + return dict(gen_power_mat()) + + +def poly_to_matrix(poly_list, power = None): + """ + + """ + + data_coord_dict = poly_to_data_coord(poly_list, power) + + n_free_symbols = len(poly_list[0].gens) + + def gen_power_mat(): + + # for all powers generate a matrix + for current_power, data_coord in data_coord_dict: + + # empty matrix + shape = (len(poly_list), n_free_symbols**current_power) + # m = np.zeros((len(poly_list), n_free_symbols**current_power)) + + # fill matrix + if len(data_coord) == 0: + yield np.zeros((len(poly_list), n_free_symbols**current_power)) + + else: + rows, cols, data = list(zip(*data_coord)) + + yield scipy.sparse.coo_matrix((data, (rows, cols)), dtype=np.double, shape=shape) + + # m[rows, cols] = data + # yield m + + return list(gen_power_mat()) diff --git a/polymatrix/utils.py b/polymatrix/utils.py new file mode 100644 index 0000000..e53db35 --- /dev/null +++ b/polymatrix/utils.py @@ -0,0 +1,34 @@ +import scipy.special + + +def variable_to_index(n_var, combination): + """ + example: + expr = x1 * x3**2 + n_var = 3 + combination = (2, 2, 0) + + returns 5 + """ + + def gen_index_sum(): + for k, x in enumerate(combination): + yield scipy.special.binom(n_var+k-1, k+1) - scipy.special.binom(n_var+k-1-x, k+1) + + return int(sum(gen_index_sum())) + + +def variable_powers_to_index(powers): + """ + example: + expr = x1 * x3**2 + powers = (1, 0, 2) + + returns 5 + """ + + n_var = len(powers) + + indices = sorted((idx for idx, p in enumerate(powers) for _ in range(p)), reverse=True) + + return variable_to_index(n_var, indices)
\ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..e2da6c2 --- /dev/null +++ b/setup.py @@ -0,0 +1,5 @@ +from setuptools import setup + +setup( + name='polymatrix' +)
\ No newline at end of file diff --git a/test_polymatrix/test_polymatrix.py b/test_polymatrix/test_polymatrix.py new file mode 100644 index 0000000..1f0bb73 --- /dev/null +++ b/test_polymatrix/test_polymatrix.py @@ -0,0 +1,176 @@ +import unittest + +from polymatrix.polystruct import init_equation, init_poly_matrix, init_poly_vector +from polymatrix.utils import variable_to_index + +class TestPolyMatrix(unittest.TestCase): + @staticmethod + def assert_term_in_eq(terms, degree, eq_idx, row_idx, value, monoms=None): + if monoms is None: + monoms = tuple() + + assert degree in terms, f'could not find {degree} in {terms}' + degree_terms = terms[degree] + + key = eq_idx, monoms + assert key in degree_terms, f'could not find {key} in {degree_terms}' + eq_terms = degree_terms[key] + + assert row_idx in eq_terms, f'could not find {row_idx} in {eq_terms}' + assert eq_terms[row_idx] == value, f'value {eq_terms[row_idx]} and {value} do not match' + + def test_param_matrix_param_vector_degree_0(self): + """ + param = [a11 a21 a31 a41 v11 v21] + """ + + mat = init_poly_matrix(degrees=(0,)) + vec = init_poly_vector(degrees=(0,)) + n_param = 6 + + eq = init_equation( + terms = [(mat, vec)], + n_var = 2, + ) + + terms, offset_dict = list(eq.create()) + + # a11 v11 + self.assert_term_in_eq( + terms = terms, + degree = 2, + eq_idx = 0, + row_idx = variable_to_index(n_param, (offset_dict[(mat, 0)], offset_dict[(vec, 0)])), + value = 1, + ) + + # a21 v11 + self.assert_term_in_eq( + terms = terms, + degree = 2, + eq_idx = 1, + row_idx = variable_to_index(n_param, (offset_dict[(mat, 0)]+2, offset_dict[(vec, 0)])), + value = 1, + ) + + def test_param_matrix_param_vector_degree_01(self): + """ + param = [a11 a21 a31 a41 v011 v021 v111 v112 v121 v122] + """ + + mat = init_poly_matrix(degrees=(0,)) + vec = init_poly_vector(degrees=(0,1)) + n_param = 10 + + eq = init_equation( + terms = [(mat, vec)], + n_var = 2, + ) + + terms, offset_dict = list(eq.create()) + + # a11 v011 + self.assert_term_in_eq( + terms = terms, + degree = 2, + eq_idx = 0, + row_idx = variable_to_index(n_param, (offset_dict[(mat, 0)], offset_dict[(vec, 0)])), + value = 1, + ) + + # a11 v011 + self.assert_term_in_eq( + terms = terms, + degree = 2, + eq_idx = 0, + monoms=(0,), + row_idx = variable_to_index(n_param, (offset_dict[(mat, 0)], offset_dict[(vec, 1)])), + value = 1, + ) + + def test_param_matrix_const_vector_degree_0(self): + """ + param = [a11 a21 a31 a41] + """ + + mat = init_poly_matrix(degrees=(0,)) + vec = init_poly_vector(subs={0: {(0, 0): 1, (1, 0): 1}}) + + eq = init_equation( + terms = [(mat, vec)], + n_var = 2, + ) + + terms, offset_dict = list(eq.create()) + + # a11 + self.assert_term_in_eq( + terms = terms, + degree = 1, + eq_idx = 0, + row_idx = offset_dict[(mat, 0)], + value = 1, + ) + + def test_const_matrix_const_vector_degree_0(self): + """ + param = [a11 a21 a31 a41] + """ + + mat = init_poly_matrix(subs={0: {(0, 0): 1, (1, 0): 1, (0, 1): 1, (1, 1): 1}}) + vec = init_poly_vector(subs={0: {(0, 0): 1, (1, 0): 1}}) + + eq = init_equation( + terms = [(mat, vec)], + n_var = 2, + ) + + terms, offset_dict = list(eq.create()) + + self.assert_term_in_eq( + terms = terms, + degree = 0, + eq_idx = 0, + row_idx = 0, + value = 2, + ) + + def test_param_matrix_const_vector_skew_symmetric(self): + """ + param = [a11 a21 a31 a41] + """ + + def skew_symmetric(idx1, idx2): + if idx1 == idx2: + return idx1, idx2, 0 + elif idx2 < idx1: + return idx2, idx1, -1 + + mat = init_poly_matrix( + degrees=(0,), + re_index_func=skew_symmetric, + ) + vec = init_poly_vector(subs={0: {(0, 0): 1, (1, 0): 1}}) + + eq = init_equation( + terms = [(mat, vec)], + n_var = 2, + ) + + terms, offset_dict = list(eq.create()) + + self.assert_term_in_eq( + terms = terms, + degree = 1, + eq_idx = 0, + row_idx = offset_dict[(mat, 0)] + 1, + value = 1, + ) + + self.assert_term_in_eq( + terms = terms, + degree = 1, + eq_idx = 1, + row_idx = offset_dict[(mat, 0)] + 1, + value = -1, + )
\ No newline at end of file |