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)