Initial commit

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
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..37ec93a
--- /dev/null
+++ b/LICENSE
@@ -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
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..806d1a1
--- /dev/null
+++ b/main.py
@@ -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