move init expr functions into a single file

103 files changed, 2491 insertions, 1649 deletions

diff --git a/polymatrix/__init__.py b/polymatrix/__init__.py
index 8cffe4b..6c927d1 100644
--- a/polymatrix/__init__.py
+++ b/polymatrix/__init__.py
@@ -1,686 +1,539 @@
-import collections
-import dataclasses
-import itertools
-import math
-import typing
-import numpy as np
-import scipy.sparse
-import sympy
-
-from polymatrix.expression.expression import Expression
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-from polymatrix.expression.mixins.parametrizeexprmixin import ParametrizeExprMixin
-from polymatrix.expression.mixins.parametrizematrixexprmixin import ParametrizeMatrixExprMixin
-from polymatrix.expressionstate.expressionstate import ExpressionState
-from polymatrix.expression.init.initblockdiagexpr import init_block_diag_expr
-from polymatrix.expression.init.initexpression import init_expression
-from polymatrix.expression.init.initeyeexpr import init_eye_expr
-from polymatrix.expression.init.initfromsympyexpr import init_from_sympy_expr
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initvstackexpr import init_v_stack_expr
-from polymatrix.polymatrix.polymatrix import PolyMatrix
-from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
-from polymatrix.statemonad.init.initstatemonad import init_state_monad
-from polymatrix.statemonad.mixins.statemonadmixin import StateMonadMixin
-from polymatrix.expression.utils.monomialtoindex import monomial_to_index
-from polymatrix.expressionstate.init.initexpressionstate import init_expression_state as original_init_expression_state
-from polymatrix.statemonad.statemonad import StateMonad
-
-
-def init_expression_state():
-    return original_init_expression_state()
-
-def from_sympy(
-    data: tuple[tuple[float]],
-):
-    return init_expression(
-        init_from_sympy_expr(data)
-    )
-
-def from_state_monad(
-    data: StateMonad,
-):
-    return init_expression(
-        data.flat_map(lambda inner_data: init_from_sympy_expr(inner_data)),
-    )
-
-def from_polymatrix(
-    polymatrix: PolyMatrix,
-):
-    return init_expression(
-        init_from_terms_expr(polymatrix)
-    )
-
-def from_(
-    data: str | tuple[tuple[float]],
-):
-    if isinstance(data, str):
-        return from_(1).parametrize(data)
-
-    return from_sympy(data)
-
-def v_stack(
-    expressions: tuple[Expression],
-):
-    return init_expression(
-        init_v_stack_expr(expressions)
-    )
-
-def h_stack(
-    expressions: tuple[Expression],
-):
-    return init_expression(
-        init_v_stack_expr(tuple(expr.T for expr in expressions))
-    ).T
-
-def block_diag(
-    expressions: tuple[Expression],
-):
-    return init_expression(
-        init_block_diag_expr(expressions)
-    )
-
-def eye(
-    variable: tuple[Expression],
-):
-    return init_expression(
-        init_eye_expr(variable=variable)
-    )
-
-def kkt_equality(
-    variables: Expression,
-    equality: Expression = None,
-):
-
-    self_variables = variables
-    self_equality = equality
-
-    def func(state: ExpressionState):
-
-        state, equality = self_equality.apply(state=state)
-
-        state, equality_der = self_equality.diff(
-            self_variables, 
-            introduce_derivatives=True,
-        ).apply(state)
-
-        def acc_nu_variables(acc, v):
-            state, nu_variables = acc
-
-            nu_variable = state.n_param
-            state = state.register(n_param=1)
-            
-            return state, nu_variables + [nu_variable]
-
-        *_, (state, nu_variables) = tuple(itertools.accumulate(
-            range(equality.shape[0]), 
-            acc_nu_variables, 
-            initial=(state, []),
-        ))
-
-        terms = {}
-
-        for row in range(equality_der.shape[1]):
-
-            monomial_terms = collections.defaultdict(float)
-
-            for eq_idx, nu_variable in enumerate(nu_variables):
-
-                underlying_terms = equality_der.get_poly(eq_idx, row)
-                if underlying_terms is None:
-                    continue
+from polymatrix.expressionstate.expressionstate import ExpressionState as internal_ExpressionState
+from polymatrix.expressionstate.init.initexpressionstate import init_expression_state as internal_init_expression_state
+from polymatrix.expression.expression import Expression as internal_Expression
+from polymatrix.expression.from_ import from_ as internal_from
+# from polymatrix.expression.from_ import from_sympy as internal_from_sympy
+from polymatrix.expression import v_stack as internal_v_stack
+from polymatrix.expression import h_stack as internal_h_stack
+from polymatrix.expression import product as internal_product
+from polymatrix.denserepr.from_ import from_polymatrix_expr
+from polymatrix.expression.to import shape as internal_shape
+from polymatrix.expression.to import to_constant_repr as internal_to_constant_repr
+from polymatrix.expression.to import to_degrees as internal_to_degrees
+from polymatrix.expression.to import to_sympy_repr as internal_to_sympy_repr
+
+Expression = internal_Expression
+ExpressionState = internal_ExpressionState
+
+init_expression_state = internal_init_expression_state
+from_ = internal_from
+# from_sympy = internal_from_sympy
+v_stack = internal_v_stack
+h_stack = internal_h_stack
+product = internal_product
+to_shape = internal_shape
+to_constant_repr = internal_to_constant_repr
+to_degrees = internal_to_degrees
+to_sympy_repr = internal_to_sympy_repr
+to_matrix_repr = from_polymatrix_expr
+
+# def from_sympy(
+#     data: tuple[tuple[float]],
+# ):
+#     return init_expression(
+#         init_from_sympy_expr(data)
+#     )
+
+# def from_state_monad(
+#     data: StateMonad,
+# ):
+#     return init_expression(
+#         data.flat_map(lambda inner_data: init_from_sympy_expr(inner_data)),
+#     )
+
+# def from_polymatrix(
+#     polymatrix: PolyMatrix,
+# ):
+#     return init_expression(
+#         init_from_terms_expr(polymatrix)
+#     )
+
+# def from_(
+#     data: str | tuple[tuple[float]],
+# ):
+#     if isinstance(data, str):
+#         return from_(1).parametrize(data)
+
+#     return from_sympy(data)
+
+# def v_stack(
+#     expressions: tuple[Expression],
+# ):
+#     return init_expression(
+#         init_v_stack_expr(expressions)
+#     )
+
+# def h_stack(
+#     expressions: tuple[Expression],
+# ):
+#     return init_expression(
+#         init_v_stack_expr(tuple(expr.T for expr in expressions))
+#     ).T
+
+# def block_diag(
+#     expressions: tuple[Expression],
+# ):
+#     return init_expression(
+#         init_block_diag_expr(expressions)
+#     )
+
+# def eye(
+#     variable: tuple[Expression],
+# ):
+#     return init_expression(
+#         init_eye_expr(variable=variable)
+#     )
+
+# def kkt_equality(
+#     equality: Expression,
+#     variables: Expression,
+#     name: str,
+# ):
+#     equality_der = equality.diff(
+#         variables, 
+#         introduce_derivatives=True,
+#     )
+
+#     nu = equality_der[0, :].T.parametrize(name)
+
+#     return nu, equality_der @ nu
+
+# def kkt_inequality(
+#     inequality: Expression,
+#     variables: Expression,
+#     name: str,
+# ):
+    
+#     inequality_der = inequality.diff(
+#         variables, 
+#         introduce_derivatives=True,
+#     )
+
+#     nu = inequality_der[0, :].T.parametrize(f'lambda_{name}')
+#     r_nu = inequality_der[0, :].T.parametrize(f'r_lambda_{name}')
+#     r_ineq = inequality_der[0, :].T.parametrize(f'r_ineq_{name}')
+
+#     prim_expr = inequality - r_ineq * r_ineq
+
+#     sec_expr = nu - r_nu * r_nu
+
+#     compl_expr = nu * inequality
+
+#     cost_expr = inequality_der @ nu
+
+#     return h_stack((nu, r_nu, r_ineq)), (cost_expr, prim_expr, sec_expr, compl_expr)
+
+
+# def rows(
+#     expr: Expression,
+# ) -> StateMonadMixin[ExpressionState, np.ndarray]:
+
+#     def func(state: ExpressionState):
+#         state, underlying = expr.apply(state)
+
+#         def gen_row_terms():
+#             for row in range(underlying.shape[0]):
 
-                for monomial, value in underlying_terms.items():
-                    new_monomial = monomial + (nu_variable,)
+#                 terms = {}
 
-                    monomial_terms[new_monomial] += value
+#                 for col in range(underlying.shape[1]):
+#                     polynomial = underlying.get_poly(row, col)
+#                     if polynomial == None:
+#                         continue
 
-            # terms[row, 0] = dict(monomial_terms)
-            terms[row, 0] = monomial_terms
+#                     terms[0, col] = polynomial
 
-        cost_expr = init_expression(init_from_terms_expr(
-            terms=terms,
-            shape=(equality_der.shape[1], 1),
-        ))
+#                 yield init_expression(underlying=init_from_terms_expr(
+#                     terms=terms,
+#                     shape=(1, underlying.shape[1])
+#                 ))
 
-        nu_terms = {}
+#         row_terms = tuple(gen_row_terms())
 
-        for eq_idx, nu_variable in enumerate(nu_variables):
-            nu_terms[eq_idx, 0] = {(nu_variable,): 1}
+#         return state, row_terms
 
-        nu_expr = init_expression(init_from_terms_expr(
-            terms=nu_terms,
-            shape=(len(nu_variables), 1),
-        ))
+#     return init_state_monad(func)
 
-        return state, (nu_expr, cost_expr)   
 
-    return init_state_monad(func)
+# @dataclasses.dataclass
+# class MatrixBuffer:
+#     data: dict[int, np.ndarray]
+#     n_row: int
+#     n_param: int
 
-def kkt_inequality(
-    variables: Expression,
-    inequality: Expression = None,
-):
+#     def get_max_degree(self):
+#         return max(degree for degree in self.data.keys())
 
-    self_variables = variables
-    self_inequality = inequality
+#     def add_buffer(self, index: int):
+#         if index <= 1:
+#             buffer = np.zeros((self.n_row, self.n_param**index), dtype=np.double)
 
-    def func(state: ExpressionState):
-
-        state, inequality = self_inequality.apply(state=state)
-
-        state, inequality_der = self_inequality.diff(
-            self_variables, 
-            introduce_derivatives=True,
-        ).apply(state)
-
-        def acc_lambda_variables(acc, v):
-            state, lambda_variables = acc
-
-            lambda_variable = state.n_param
-            state = state.register(n_param=3)
+#         else:
+#             buffer = scipy.sparse.dok_array((self.n_row, self.n_param**index), dtype=np.double)
             
-            return state, lambda_variables + [lambda_variable]
-
-        *_, (state, lambda_variables) = tuple(itertools.accumulate(
-            range(inequality.shape[0]), 
-            acc_lambda_variables, 
-            initial=(state, []),
-        ))
-
-        terms = {}
-
-        for row in range(inequality_der.shape[1]):
-
-            monomial_terms = collections.defaultdict(float)
-
-            for inequality_idx, lambda_variable in enumerate(lambda_variables):
-
-                polynomial = inequality_der.get_poly(inequality_idx, row)
-                if polynomial is None:
-                    continue
-
-                for monomial, value in polynomial.items():
-                    new_monomial = monomial + (lambda_variable,)
-
-                    monomial_terms[new_monomial] += value
-
-            # terms[row, 0] = dict(monomial_terms)
-            terms[row, 0] = monomial_terms
-
-        cost_expr = init_expression(init_from_terms_expr(
-            terms=terms,
-            shape=(inequality_der.shape[1], 1),
-        ))
-
-        # inequality, dual feasibility, complementary slackness
-        # -----------------------------------------------------
-
-        inequality_terms = {}
-        feasibility_terms = {}
-        complementary_terms = {}
-
-        for inequality_idx, lambda_variable in enumerate(lambda_variables):
-            r_lambda = lambda_variable + 1
-            r_inequality = lambda_variable + 2
-
-            polynomial = inequality.get_poly(inequality_idx, 0)
-            if polynomial is None:
-                continue
-
-            # f(x) <= -0.01
-            inequality_terms[inequality_idx, 0] = polynomial | {(r_inequality, r_inequality): 1}
-
-            # dual feasibility, lambda >= 0
-            feasibility_terms[inequality_idx, 0] = {(lambda_variable,): 1, (r_lambda, r_lambda): -1}
+#         self.data[index] = buffer
 
-            # complementary slackness
-            complementary_terms[inequality_idx, 0] = {(r_lambda, r_inequality): 1}
+#     def add(self, row: int, col: int, index: int, value: float):
+#         if index not in self.data:
+#             self.add_buffer(index)
 
-        inequality_expr = init_expression(init_from_terms_expr(
-            terms=inequality_terms,
-            shape=(len(lambda_variables), 1),
-        ))
+#         self.data[index][row, col] = value
 
-        feasibility_expr = init_expression(init_from_terms_expr(
-            terms=feasibility_terms,
-            shape=(len(lambda_variables), 1),
-        ))
+#     def __getitem__(self, key):
+#         if key not in self.data:
+#             self.add_buffer(key)
 
-        complementary_expr = init_expression(init_from_terms_expr(
-            terms=complementary_terms,
-            shape=(len(lambda_variables), 1),
-        ))
+#         return self.data[key]
 
-        # lambda expression
-        # -----------------
 
-        terms = {}
-        for inequality_idx, lambda_variable in enumerate(lambda_variables):
-            terms[inequality_idx, 0] = {(lambda_variable,): 1}
+# @dataclasses.dataclass
+# class MatrixRepresentations:
+#     data: tuple[MatrixBuffer, ...]
+#     aux_data: typing.Optional[MatrixBuffer]
+#     variable_mapping: tuple[int, ...]
+#     state: ExpressionState
 
-        lambda_expr = init_expression(init_from_terms_expr(
-            terms=terms,
-            shape=(len(lambda_variables), 1),
-        ))
+#     def merge_matrix_equations(self):
+#         def gen_matrices(index: int):
+#             for equations in self.data:
+#                 if index < len(equations):
+#                     yield equations[index]
 
-        return state, (lambda_expr, cost_expr, inequality_expr, feasibility_expr, complementary_expr)   
+#             if index < len(self.aux_data):
+#                 yield self.aux_data[index]
 
-    return init_state_monad(func)
+#         indices = set(key for equations in self.data + (self.aux_data,) for key in equations.keys())
 
-def rows(
-    expr: Expression,
-) -> StateMonadMixin[ExpressionState, np.ndarray]:
+#         def gen_matrices():
+#             for index in indices:
+#                 if index <= 1:
+#                     yield index, np.vstack(tuple(gen_matrices(index)))
+#                 else:
+#                     yield index, scipy.sparse.vstack(tuple(gen_matrices(index)))
 
-    def func(state: ExpressionState):
-        state, underlying = expr.apply(state)
+#         return dict(gen_matrices())
 
-        def gen_row_terms():
-            for row in range(underlying.shape[0]):
+#     def get_value(self, variable, value):
+#         variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
 
-                terms = {}
+#         def gen_value_index():
+#             for variable_index in variable_indices:
+#                 try:
+#                     yield self.variable_mapping.index(variable_index)
+#                 except ValueError:
+#                     raise ValueError(f'{variable_index} not found in {self.variable_mapping}')
 
-                for col in range(underlying.shape[1]):
-                    polynomial = underlying.get_poly(row, col)
-                    if polynomial == None:
-                        continue
+#         value_index = list(gen_value_index())
 
-                    terms[0, col] = polynomial
+#         return value[value_index]
 
-                yield init_expression(underlying=init_from_terms_expr(
-                    terms=terms,
-                    shape=(1, underlying.shape[1])
-                ))
+#     def set_value(self, variable, value):
+#         variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
+#         value_index = list(self.variable_mapping.index(variable_index) for variable_index in variable_indices)
+#         vec = np.zeros(len(self.variable_mapping))
+#         vec[value_index] = value
+#         return vec
 
-        row_terms = tuple(gen_row_terms())
+#     # def get_matrix(self, eq_idx: int):
+#     #     equations = self.data[eq_idx].data
 
-        return state, row_terms
+#     def get_func(self, eq_idx: int):
+#         equations = self.data[eq_idx].data
+#         max_idx = max(equations.keys())
 
-    return init_state_monad(func)
+#         if 2 <= max_idx:
+#             def func(x: np.ndarray) -> np.ndarray:
+#                 if isinstance(x, tuple) or isinstance(x, list):
+#                     x = np.array(x).reshape(-1, 1)
 
-def shape(
-    expr: Expression,
-) -> StateMonadMixin[ExpressionState, tuple[int, ...]]:
-    def func(state: ExpressionState):
-        state, polymatrix = expr.apply(state)
+#                 elif x.shape[0] == 1:
+#                     x = x.reshape(-1, 1)
 
-        return state, polymatrix.shape
+#                 def acc_x_powers(acc, _):
+#                     next = (acc @ x.T).reshape(-1, 1)
+#                     return next
 
-    return init_state_monad(func)
+#                 x_powers = tuple(itertools.accumulate(
+#                     range(max_idx - 1),
+#                     acc_x_powers,
+#                     initial=x,
+#                 ))[1:]
 
+#                 def gen_value():
+#                     for idx, equation in equations.items():
+#                         if idx == 0:
+#                             yield equation
 
-@dataclasses.dataclass
-class MatrixBuffer:
-    data: dict[int, np.ndarray]
-    n_row: int
-    n_param: int
+#                         elif idx == 1:
+#                             yield equation @ x
 
-    def get_max_degree(self):
-        return max(degree for degree in self.data.keys())
+#                         else:
+#                             yield equation @ x_powers[idx-2]
 
-    def add_buffer(self, index: int):
-        if index <= 1:
-            buffer = np.zeros((self.n_row, self.n_param**index), dtype=np.double)
+#                 return sum(gen_value())
 
-        else:
-            buffer = scipy.sparse.dok_array((self.n_row, self.n_param**index), dtype=np.double)
+#         else:
+#             def func(x: np.ndarray) -> np.ndarray:
+#                 if isinstance(x, tuple) or isinstance(x, list):
+#                     x = np.array(x).reshape(-1, 1)
             
-        self.data[index] = buffer
-
-    def add(self, row: int, col: int, index: int, value: float):
-        if index not in self.data:
-            self.add_buffer(index)
-
-        self.data[index][row, col] = value
-
-    def __getitem__(self, key):
-        if key not in self.data:
-            self.add_buffer(key)
-
-        return self.data[key]
-
-
-@dataclasses.dataclass
-class MatrixRepresentations:
-    data: tuple[MatrixBuffer, ...]
-    aux_data: typing.Optional[MatrixBuffer]
-    variable_mapping: tuple[int, ...]
-    state: ExpressionState
-
-    def merge_matrix_equations(self):
-        def gen_matrices(index: int):
-            for equations in self.data:
-                if index < len(equations):
-                    yield equations[index]
-
-            if index < len(self.aux_data):
-                yield self.aux_data[index]
-
-        indices = set(key for equations in self.data + (self.aux_data,) for key in equations.keys())
-
-        def gen_matrices():
-            for index in indices:
-                if index <= 1:
-                    yield index, np.vstack(tuple(gen_matrices(index)))
-                else:
-                    yield index, scipy.sparse.vstack(tuple(gen_matrices(index)))
-
-        return dict(gen_matrices())
-
-    def get_value(self, variable, value):
-        variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
-
-        def gen_value_index():
-            for variable_index in variable_indices:
-                try:
-                    yield self.variable_mapping.index(variable_index)
-                except ValueError:
-                    raise ValueError(f'{variable_index} not found in {self.variable_mapping}')
-
-        value_index = list(gen_value_index())
-
-        return value[value_index]
-
-    def set_value(self, variable, value):
-        variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
-        value_index = list(self.variable_mapping.index(variable_index) for variable_index in variable_indices)
-        vec = np.zeros(len(self.variable_mapping))
-        vec[value_index] = value
-        return vec
-
-    def get_func(self, eq_idx: int):
-        equations = self.data[eq_idx].data
-        max_idx = max(equations.keys())
-
-        if 2 <= max_idx:
-            def func(x: np.ndarray) -> np.ndarray:
-                if isinstance(x, tuple) or isinstance(x, list):
-                    x = np.array(x).reshape(-1, 1)
-
-                elif x.shape[0] == 1:
-                    x = x.reshape(-1, 1)
-
-                def acc_x_powers(acc, _):
-                    next = (acc @ x.T).reshape(-1, 1)
-                    return next
-
-                x_powers = tuple(itertools.accumulate(
-                    range(max_idx - 1),
-                    acc_x_powers,
-                    initial=x,
-                ))[1:]
-
-                def gen_value():
-                    for idx, equation in equations.items():
-                        if idx == 0:
-                            yield equation
-
-                        elif idx == 1:
-                            yield equation @ x
-
-                        else:
-                            yield equation @ x_powers[idx-2]
-
-                return sum(gen_value())
-
-        else:
-            def func(x: np.ndarray) -> np.ndarray:
-                def gen_value():
-                    for idx, equation in equations.items():
-                        if idx == 0:
-                            yield equation
+#                 def gen_value():
+#                     for idx, equation in equations.items():
+#                         if idx == 0:
+#                             yield equation
 
-                        else:
-                            yield equation @ x
+#                         else:
+#                             yield equation @ x
 
-                return sum(gen_value())
+#                 return sum(gen_value())
 
-        return func
+#         return func
 
-def to_matrix_repr(
-    expressions: Expression | tuple[Expression],
-    variables: Expression = None,
-    sorted: bool = None,
-) -> StateMonadMixin[ExpressionState, tuple[tuple[tuple[np.ndarray, ...], ...], tuple[int, ...]]]:
+# def to_matrix_repr(
+#     expressions: Expression | tuple[Expression],
+#     variables: Expression = None,
+#     sorted: bool = None,
+# ) -> StateMonadMixin[ExpressionState, tuple[tuple[tuple[np.ndarray, ...], ...], tuple[int, ...]]]:
 
-    if isinstance(expressions, Expression):
-        expressions = (expressions,)
+#     if isinstance(expressions, Expression):
+#         expressions = (expressions,)
 
-    assert isinstance(variables, ExpressionBaseMixin) or variables is None, f'{variables=}'
+#     assert isinstance(variables, ExpressionBaseMixin) or variables is None, f'{variables=}'
 
-    def func(state: ExpressionState):
+#     def func(state: ExpressionState):
 
-        def acc_underlying_application(acc, v):
-            state, underlying_list = acc
+#         def acc_underlying_application(acc, v):
+#             state, underlying_list = acc
 
-            state, underlying = v.apply(state)
+#             state, underlying = v.apply(state)
 
-            assert underlying.shape[1] == 1, f'{underlying.shape[1]=} is not 1'
+#             assert underlying.shape[1] == 1, f'{underlying.shape[1]=} is not 1'
             
-            return state, underlying_list + (underlying,)
+#             return state, underlying_list + (underlying,)
 
-        *_, (state, polymatrix_list) = tuple(itertools.accumulate(
-            expressions, 
-            acc_underlying_application, 
-            initial=(state, tuple()),
-        ))
+#         *_, (state, polymatrix_list) = tuple(itertools.accumulate(
+#             expressions, 
+#             acc_underlying_application, 
+#             initial=(state, tuple()),
+#         ))
 
-        if variables is None:
-            sorted_variable_index = tuple()
+#         if variables is None:
+#             sorted_variable_index = tuple()
             
-        else:
-            state, variable_index = get_variable_indices_from_variable(state, variables)
+#         else:
+#             state, variable_index = get_variable_indices_from_variable(state, variables)
 
-            if sorted:
-                tagged_variable_index = tuple((offset, state.get_name_from_offset(offset)) for offset in variable_index)
+#             if sorted:
+#                 tagged_variable_index = tuple((offset, state.get_name_from_offset(offset)) for offset in variable_index)
 
-                sorted_variable_index = tuple(v[0] for v in sorted(tagged_variable_index, key=lambda v: v[1]))
+#                 sorted_variable_index = tuple(v[0] for v in sorted(tagged_variable_index, key=lambda v: v[1]))
                 
-            else:
-                sorted_variable_index = variable_index
+#             else:
+#                 sorted_variable_index = variable_index
 
-            sorted_variable_index_set = set(sorted_variable_index)
-            if len(sorted_variable_index) != len(sorted_variable_index_set):
-                duplicates = tuple(state.get_name_from_offset(var) for var in sorted_variable_index_set if 1 < sorted_variable_index.count(var))
-                raise Exception(f'{duplicates=}. Make sure you give a unique name for each variables.')
+#             sorted_variable_index_set = set(sorted_variable_index)
+#             if len(sorted_variable_index) != len(sorted_variable_index_set):
+#                 duplicates = tuple(state.get_name_from_offset(var) for var in sorted_variable_index_set if 1 < sorted_variable_index.count(var))
+#                 raise Exception(f'{duplicates=}. Make sure you give a unique name for each variables.')
         
-        variable_index_map = {old: new for new, old in enumerate(sorted_variable_index)}
+#         variable_index_map = {old: new for new, old in enumerate(sorted_variable_index)}
 
-        n_param = len(sorted_variable_index)
+#         n_param = len(sorted_variable_index)
 
-        def gen_numpy_matrices():
-            for polymatrix in polymatrix_list:
+#         def gen_numpy_matrices():
+#             for polymatrix in polymatrix_list:
 
-                n_row = polymatrix.shape[0]
+#                 n_row = polymatrix.shape[0]
 
-                buffer = MatrixBuffer(
-                    data={},
-                    n_row=n_row,
-                    n_param=n_param,
-                )
+#                 buffer = MatrixBuffer(
+#                     data={},
+#                     n_row=n_row,
+#                     n_param=n_param,
+#                 )
 
-                for row in range(n_row):
-                    polymatrix_terms = polymatrix.get_poly(row, 0)
+#                 for row in range(n_row):
+#                     polymatrix_terms = polymatrix.get_poly(row, 0)
                     
-                    if polymatrix_terms is None:
-                        continue
+#                     if polymatrix_terms is None:
+#                         continue
 
-                    if len(polymatrix_terms) == 0:
-                        buffer.add(row, 0, 0, 0)
+#                     if len(polymatrix_terms) == 0:
+#                         buffer.add(row, 0, 0, 0)
                     
-                    else:
-                        for monomial, value in polymatrix_terms.items():
+#                     else:
+#                         for monomial, value in polymatrix_terms.items():
 
-                            def gen_new_monomial():
-                                for var, count in monomial:
-                                    try:
-                                        index = variable_index_map[var]
-                                    except KeyError:
-                                        raise KeyError(f'{var=} ({state.get_key_from_offset(var)}) is incompatible with {variable_index_map=}')
+#                             def gen_new_monomial():
+#                                 for var, count in monomial:
+#                                     try:
+#                                         index = variable_index_map[var]
+#                                     except KeyError:
+#                                         raise KeyError(f'{var=} ({state.get_key_from_offset(var)}) is incompatible with {variable_index_map=}')
 
-                                    for _ in range(count):
-                                        yield index
+#                                     for _ in range(count):
+#                                         yield index
 
-                            new_monomial = tuple(gen_new_monomial())
+#                             new_monomial = tuple(gen_new_monomial())
 
-                            cols = monomial_to_index(n_param, new_monomial)
+#                             cols = monomial_to_index(n_param, new_monomial)
 
-                            col_value = value / len(cols)
+#                             col_value = value / len(cols)
 
-                            for col in cols:
-                                buffer.add(row, col, sum(count for _, count in monomial), col_value)
+#                             for col in cols:
+#                                 degree = sum(count for _, count in monomial)
+#                                 buffer.add(row, col, degree, col_value)
                 
-                yield buffer
+#                 yield buffer
 
-        underlying_matrices = tuple(gen_numpy_matrices())
+#         underlying_matrices = tuple(gen_numpy_matrices())
 
-        def gen_auxillary_equations():
-            for key, monomial_terms in state.auxillary_equations.items():
-                if key in sorted_variable_index:
-                    yield key, monomial_terms
+#         def gen_auxillary_equations():
+#             for key, monomial_terms in state.auxillary_equations.items():
+#                 if key in sorted_variable_index:
+#                     yield key, monomial_terms
 
-        auxillary_equations = tuple(gen_auxillary_equations())
+#         auxillary_equations = tuple(gen_auxillary_equations())
 
-        n_row = len(auxillary_equations)
+#         n_row = len(auxillary_equations)
 
-        if n_row == 0:
-            auxillary_matrix_equations = None
+#         if n_row == 0:
+#             auxillary_matrix_equations = None
 
-        else:
-            buffer = MatrixBuffer(
-                data={},
-                n_row=n_row,
-                n_param=n_param,
-            )
+#         else:
+#             buffer = MatrixBuffer(
+#                 data={},
+#                 n_row=n_row,
+#                 n_param=n_param,
+#             )
 
-            for row, (key, monomial_terms) in enumerate(auxillary_equations):
-                for monomial, value in monomial_terms.items():
-                    new_monomial = tuple(variable_index_map[var] for var, count in monomial for _ in range(count))
+#             for row, (key, monomial_terms) in enumerate(auxillary_equations):
+#                 for monomial, value in monomial_terms.items():
+#                     new_monomial = tuple(variable_index_map[var] for var, count in monomial for _ in range(count))
 
-                    cols = monomial_to_index(n_param, new_monomial)
+#                     cols = monomial_to_index(n_param, new_monomial)
 
-                    col_value = value / len(cols)
+#                     col_value = value / len(cols)
 
-                    for col in cols:
-                        buffer.add(row, col, sum(count for _, count in monomial), col_value)
+#                     for col in cols:
+#                         buffer.add(row, col, sum(count for _, count in monomial), col_value)
 
-            auxillary_matrix_equations = buffer
+#             auxillary_matrix_equations = buffer
 
-        result = MatrixRepresentations(
-            data=underlying_matrices,
-            aux_data=auxillary_matrix_equations,
-            variable_mapping=sorted_variable_index,
-            state=state,
-        )
+#         result = MatrixRepresentations(
+#             data=underlying_matrices,
+#             aux_data=auxillary_matrix_equations,
+#             variable_mapping=sorted_variable_index,
+#             state=state,
+#         )
 
-        return state, result
+#         return state, result
 
-    return init_state_monad(func)
+#     return init_state_monad(func)
 
 
-def to_constant_repr(
-    expr: Expression,
-    assert_constant: bool = True,
-) -> StateMonadMixin[ExpressionState, np.ndarray]:
+# def to_constant_repr(
+#     expr: Expression,
+#     assert_constant: bool = True,
+# ) -> StateMonadMixin[ExpressionState, np.ndarray]:
 
-    def func(state: ExpressionState):
-        state, underlying = expr.apply(state)
+#     def func(state: ExpressionState):
+#         state, underlying = expr.apply(state)
 
-        A = np.zeros(underlying.shape, dtype=np.double)
+#         A = np.zeros(underlying.shape, dtype=np.double)
 
-        for (row, col), polynomial in underlying.gen_terms():
-            for monomial, value in polynomial.items():
-                if len(monomial) == 0:
-                    A[row, col] = value
+#         for (row, col), polynomial in underlying.gen_terms():
+#             for monomial, value in polynomial.items():
+#                 if len(monomial) == 0:
+#                     A[row, col] = value
                 
-                elif assert_constant:
-                    raise Exception(f'non-constant term {monomial=}')
+#                 elif assert_constant:
+#                     raise Exception(f'non-constant term {monomial=}')
 
-        return state, A
+#         return state, A
 
-    return init_state_monad(func)
+#     return init_state_monad(func)
 
 
-def degrees(
-    expr: Expression,
-    variables: Expression,
-) -> StateMonadMixin[ExpressionState, np.ndarray]:
+# def degrees(
+#     expr: Expression,
+#     variables: Expression,
+# ) -> StateMonadMixin[ExpressionState, np.ndarray]:
 
-    def func(state: ExpressionState):
-        state, underlying = expr.apply(state)
-        state, variable_indices = get_variable_indices_from_variable(state, variables)
+#     def func(state: ExpressionState):
+#         state, underlying = expr.apply(state)
+#         state, variable_indices = get_variable_indices_from_variable(state, variables)
 
-        def gen_rows():
-            for row in range(underlying.shape[0]):
-                def gen_cols():
-                    for col in range(underlying.shape[1]):
+#         def gen_rows():
+#             for row in range(underlying.shape[0]):
+#                 def gen_cols():
+#                     for col in range(underlying.shape[1]):
                         
-                        def gen_degrees():
-                            polynomial = underlying.get_poly(row, col)
+#                         def gen_degrees():
+#                             polynomial = underlying.get_poly(row, col)
 
-                            if polynomial is None:
-                                yield 0
+#                             if polynomial is None:
+#                                 yield 0
 
-                            else:
-                                for monomial, _ in polynomial.items():
-                                    yield sum(count for var, count in monomial if var in variable_indices)
+#                             else:
+#                                 for monomial, _ in polynomial.items():
+#                                     yield sum(count for var, count in monomial if var in variable_indices)
 
-                        yield tuple(set(gen_degrees()))
+#                         yield tuple(set(gen_degrees()))
 
-                yield tuple(gen_cols())
+#                 yield tuple(gen_cols())
 
-        return state, tuple(gen_rows())
+#         return state, tuple(gen_rows())
 
-    return init_state_monad(func)
+#     return init_state_monad(func)
 
 
-def to_sympy_repr(
-    expr: Expression,
-) -> StateMonadMixin[ExpressionState, sympy.Expr]:
+# def to_sympy_repr(
+#     expr: Expression,
+# ) -> StateMonadMixin[ExpressionState, sympy.Expr]:
 
-    def func(state: ExpressionState):
-        state, underlying = expr.apply(state)
+#     def func(state: ExpressionState):
+#         state, underlying = expr.apply(state)
 
-        A = np.zeros(underlying.shape, dtype=object)
+#         A = np.zeros(underlying.shape, dtype=object)
 
-        for (row, col), polynomial in underlying.gen_terms():
+#         for (row, col), polynomial in underlying.gen_terms():
 
-            sympy_polynomial = 0
+#             sympy_polynomial = 0
 
-            for monomial, value in polynomial.items():
-                sympy_monomial = 1
+#             for monomial, value in polynomial.items():
+#                 sympy_monomial = 1
 
-                for offset, count in monomial:
+#                 for offset, count in monomial:
 
-                    variable = state.get_key_from_offset(offset)
-                    # def get_variable_from_offset(offset: int):
-                    #     for variable, (start, end) in state.offset_dict.items():
-                    #         if start <= offset < end:
-                                # assert end - start == 1, f'{start=}, {end=}, {variable=}'
+#                     variable = state.get_key_from_offset(offset)
+#                     # def get_variable_from_offset(offset: int):
+#                     #     for variable, (start, end) in state.offset_dict.items():
+#                     #         if start <= offset < end:
+#                                 # assert end - start == 1, f'{start=}, {end=}, {variable=}'
 
-                    if isinstance(variable, sympy.core.symbol.Symbol):
-                        variable_name = variable.name
-                    elif isinstance(variable, (ParametrizeExprMixin, ParametrizeMatrixExprMixin)):
-                        variable_name = variable.name
-                    elif isinstance(variable, str):
-                        variable_name = variable
-                    else:
-                        raise Exception(f'{variable=}')
+#                     if isinstance(variable, sympy.core.symbol.Symbol):
+#                         variable_name = variable.name
+#                     elif isinstance(variable, (ParametrizeExprMixin, ParametrizeMatrixExprMixin)):
+#                         variable_name = variable.name
+#                     elif isinstance(variable, str):
+#                         variable_name = variable
+#                     else:
+#                         raise Exception(f'{variable=}')
 
-                    start, end = state.offset_dict[variable]
+#                     start, end = state.offset_dict[variable]
 
-                    if end - start == 1:
-                        sympy_var = sympy.Symbol(variable_name)
-                    else:
-                        sympy_var = sympy.Symbol(f'{variable_name}_{offset - start + 1}')
+#                     if end - start == 1:
+#                         sympy_var = sympy.Symbol(variable_name)
+#                     else:
+#                         sympy_var = sympy.Symbol(f'{variable_name}_{offset - start + 1}')
 
-                    # var = get_variable_from_offset(offset)
-                    sympy_monomial *= sympy_var**count
+#                     # var = get_variable_from_offset(offset)
+#                     sympy_monomial *= sympy_var**count
 
-                sympy_polynomial += value * sympy_monomial
+#                 sympy_polynomial += value * sympy_monomial
 
-            A[row, col] = sympy_polynomial
+#             A[row, col] = sympy_polynomial
 
-        return state, A
+#         return state, A
 
-    return init_state_monad(func)
+#     return init_state_monad(func)
diff --git a/polymatrix/expression/init/__init__.py b/polymatrix/denserepr/__init__.py
index e69de29..e69de29 100644
--- a/polymatrix/expression/init/__init__.py
+++ b/polymatrix/denserepr/__init__.py
diff --git a/polymatrix/denserepr/from_.py b/polymatrix/denserepr/from_.py
new file mode 100644
index 0000000..52fda35
--- /dev/null
+++ b/polymatrix/denserepr/from_.py
@@ -0,0 +1,156 @@
+import dataclasses
+import itertools
+import typing
+import numpy as np
+import scipy.sparse
+from polymatrix.denserepr.impl import DenseReprBufferImpl, DenseReprImpl
+
+from polymatrix.expression.expression import Expression
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
+from polymatrix.statemonad.init.initstatemonad import init_state_monad
+from polymatrix.statemonad.mixins.statemonadmixin import StateMonadMixin
+from polymatrix.expression.utils.monomialtoindex import monomial_to_index
+
+
+def from_polymatrix_expr(
+    expressions: Expression | tuple[Expression],
+    variables: Expression = None,
+    sorted: bool = None,
+) -> StateMonadMixin[ExpressionState, tuple[tuple[tuple[np.ndarray, ...], ...], tuple[int, ...]]]:
+
+    if isinstance(expressions, Expression):
+        expressions = (expressions,)
+
+    assert isinstance(variables, ExpressionBaseMixin) or variables is None, f'{variables=}'
+
+    def func(state: ExpressionState):
+
+        def acc_underlying_application(acc, v):
+            state, underlying_list = acc
+
+            state, underlying = v.apply(state)
+
+            assert underlying.shape[1] == 1, f'{underlying.shape[1]=} is not 1'
+            
+            return state, underlying_list + (underlying,)
+
+        *_, (state, polymatrix_list) = tuple(itertools.accumulate(
+            expressions, 
+            acc_underlying_application, 
+            initial=(state, tuple()),
+        ))
+
+        if variables is None:
+            sorted_variable_index = tuple()
+            
+        else:
+            state, variable_index = get_variable_indices_from_variable(state, variables)
+
+            if sorted:
+                tagged_variable_index = tuple((offset, state.get_name_from_offset(offset)) for offset in variable_index)
+
+                sorted_variable_index = tuple(v[0] for v in sorted(tagged_variable_index, key=lambda v: v[1]))
+                
+            else:
+                sorted_variable_index = variable_index
+
+            sorted_variable_index_set = set(sorted_variable_index)
+            if len(sorted_variable_index) != len(sorted_variable_index_set):
+                duplicates = tuple(state.get_name_from_offset(var) for var in sorted_variable_index_set if 1 < sorted_variable_index.count(var))
+                raise Exception(f'{duplicates=}. Make sure you give a unique name for each variables.')
+        
+        variable_index_map = {old: new for new, old in enumerate(sorted_variable_index)}
+
+        n_param = len(sorted_variable_index)
+
+        def gen_numpy_matrices():
+            for polymatrix in polymatrix_list:
+
+                n_row = polymatrix.shape[0]
+
+                buffer = DenseReprBufferImpl(
+                    data={},
+                    n_row=n_row,
+                    n_param=n_param,
+                )
+
+                for row in range(n_row):
+                    polymatrix_terms = polymatrix.get_poly(row, 0)
+                    
+                    if polymatrix_terms is None:
+                        continue
+
+                    if len(polymatrix_terms) == 0:
+                        buffer.add(row, 0, 0, 0)
+                    
+                    else:
+                        for monomial, value in polymatrix_terms.items():
+
+                            def gen_new_monomial():
+                                for var, count in monomial:
+                                    try:
+                                        index = variable_index_map[var]
+                                    except KeyError:
+                                        raise KeyError(f'{var=} ({state.get_key_from_offset(var)}) is incompatible with {variable_index_map=}')
+
+                                    for _ in range(count):
+                                        yield index
+
+                            new_monomial = tuple(gen_new_monomial())
+
+                            cols = monomial_to_index(n_param, new_monomial)
+
+                            col_value = value / len(cols)
+
+                            for col in cols:
+                                degree = sum(count for _, count in monomial)
+                                buffer.add(row, col, degree, col_value)
+                
+                yield buffer
+
+        underlying_matrices = tuple(gen_numpy_matrices())
+
+        def gen_auxillary_equations():
+            for key, monomial_terms in state.auxillary_equations.items():
+                if key in sorted_variable_index:
+                    yield key, monomial_terms
+
+        auxillary_equations = tuple(gen_auxillary_equations())
+
+        n_row = len(auxillary_equations)
+
+        if n_row == 0:
+            auxillary_matrix_equations = None
+
+        else:
+            buffer = DenseReprBufferImpl(
+                data={},
+                n_row=n_row,
+                n_param=n_param,
+            )
+
+            for row, (key, monomial_terms) in enumerate(auxillary_equations):
+                for monomial, value in monomial_terms.items():
+                    new_monomial = tuple(variable_index_map[var] for var, count in monomial for _ in range(count))
+
+                    cols = monomial_to_index(n_param, new_monomial)
+
+                    col_value = value / len(cols)
+
+                    for col in cols:
+                        buffer.add(row, col, sum(count for _, count in monomial), col_value)
+
+            auxillary_matrix_equations = buffer
+
+        result = DenseReprImpl(
+            data=underlying_matrices,
+            aux_data=auxillary_matrix_equations,
+            variable_mapping=sorted_variable_index,
+            state=state,
+        )
+
+        return state, result
+
+    return init_state_monad(func)
\ No newline at end of file
diff --git a/polymatrix/denserepr/impl.py b/polymatrix/denserepr/impl.py
new file mode 100644
index 0000000..2b57e2a
--- /dev/null
+++ b/polymatrix/denserepr/impl.py
@@ -0,0 +1,143 @@
+import dataclasses
+import itertools
+import typing
+import numpy as np
+import scipy.sparse
+
+from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
+
+
+@dataclasses.dataclass
+class DenseReprBufferImpl:
+    data: dict[int, np.ndarray]
+    n_row: int
+    n_param: int
+
+    def get_max_degree(self):
+        return max(degree for degree in self.data.keys())
+
+    def add_buffer(self, index: int):
+        if index <= 1:
+            buffer = np.zeros((self.n_row, self.n_param**index), dtype=np.double)
+
+        else:
+            buffer = scipy.sparse.dok_array((self.n_row, self.n_param**index), dtype=np.double)
+            
+        self.data[index] = buffer
+
+    def add(self, row: int, col: int, index: int, value: float):
+        if index not in self.data:
+            self.add_buffer(index)
+
+        self.data[index][row, col] = value
+
+    def __getitem__(self, key):
+        if key not in self.data:
+            self.add_buffer(key)
+
+        return self.data[key]
+
+
+@dataclasses.dataclass
+class DenseReprImpl:
+    data: tuple[DenseReprBufferImpl, ...]
+    aux_data: typing.Optional[DenseReprBufferImpl]
+    variable_mapping: tuple[int, ...]
+    state: ExpressionState
+
+    def merge_matrix_equations(self):
+        def gen_matrices(index: int):
+            for equations in self.data:
+                if index < len(equations):
+                    yield equations[index]
+
+            if index < len(self.aux_data):
+                yield self.aux_data[index]
+
+        indices = set(key for equations in self.data + (self.aux_data,) for key in equations.keys())
+
+        def gen_matrices():
+            for index in indices:
+                if index <= 1:
+                    yield index, np.vstack(tuple(gen_matrices(index)))
+                else:
+                    yield index, scipy.sparse.vstack(tuple(gen_matrices(index)))
+
+        return dict(gen_matrices())
+
+    def get_value(self, variable, value):
+        variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
+
+        def gen_value_index():
+            for variable_index in variable_indices:
+                try:
+                    yield self.variable_mapping.index(variable_index)
+                except ValueError:
+                    raise ValueError(f'{variable_index} not found in {self.variable_mapping}')
+
+        value_index = list(gen_value_index())
+
+        return value[value_index]
+
+    def set_value(self, variable, value):
+        variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
+        value_index = list(self.variable_mapping.index(variable_index) for variable_index in variable_indices)
+        vec = np.zeros(len(self.variable_mapping))
+        vec[value_index] = value
+        return vec
+
+    # def get_matrix(self, eq_idx: int):
+    #     equations = self.data[eq_idx].data
+
+    def get_func(self, eq_idx: int):
+        equations = self.data[eq_idx].data
+        max_idx = max(equations.keys())
+
+        if 2 <= max_idx:
+            def func(x: np.ndarray) -> np.ndarray:
+                if isinstance(x, tuple) or isinstance(x, list):
+                    x = np.array(x).reshape(-1, 1)
+
+                elif x.shape[0] == 1:
+                    x = x.reshape(-1, 1)
+
+                def acc_x_powers(acc, _):
+                    next = (acc @ x.T).reshape(-1, 1)
+                    return next
+
+                x_powers = tuple(itertools.accumulate(
+                    range(max_idx - 1),
+                    acc_x_powers,
+                    initial=x,
+                ))[1:]
+
+                def gen_value():
+                    for idx, equation in equations.items():
+                        if idx == 0:
+                            yield equation
+
+                        elif idx == 1:
+                            yield equation @ x
+
+                        else:
+                            yield equation @ x_powers[idx-2]
+
+                return sum(gen_value())
+
+        else:
+            def func(x: np.ndarray) -> np.ndarray:
+                if isinstance(x, tuple) or isinstance(x, list):
+                    x = np.array(x).reshape(-1, 1)
+            
+                def gen_value():
+                    for idx, equation in equations.items():
+                        if idx == 0:
+                            yield equation
+
+                        else:
+                            yield equation @ x
+
+                return sum(gen_value())
+
+        return func
\ No newline at end of file
diff --git a/polymatrix/expression/__init__.py b/polymatrix/expression/__init__.py
index e69de29..75704ff 100644
--- a/polymatrix/expression/__init__.py
+++ b/polymatrix/expression/__init__.py
@@ -0,0 +1,52 @@
+import polymatrix.expression.from_
+import polymatrix.expression.implexpressionbase
+
+from collections.abc import Iterable
+
+from polymatrix.utils.getstacklines import get_stack_lines
+from polymatrix.expression.expression import init_expression, Expression
+
+
+def v_stack(
+    expressions: Iterable[Expression],
+) -> Expression:
+	
+	def gen_underlying():
+		for expr in expressions:
+			if isinstance(expr, Expression):
+				yield expr
+			else:
+				yield polymatrix.expression.from_.from_(expr)
+
+	return init_expression(
+		underlying=polymatrix.expression.implexpressionbase.VStackExprImpl(
+            underlying=tuple(gen_underlying()),
+        ),
+    )
+
+def h_stack(
+    expressions: Iterable[Expression],
+) -> Expression:
+    return v_stack((expr.T for expr in expressions)).T
+
+def block_diag(
+    expressions: tuple[Expression],
+) -> Expression:
+    return init_expression(
+        underlying=polymatrix.expression.implexpressionbase.BlockDiagExprImpl(
+            underlying=expressions,
+        )
+    )
+
+def product(
+	expressions: Iterable[Expression],
+	degrees: tuple[int, ...] = None,
+):
+    return init_expression(
+        underlying=polymatrix.expression.implexpressionbase.ProductExprImpl(
+            underlying=tuple(expressions),
+			degrees=degrees,
+            stack=get_stack_lines(),
+        )
+    )
+
diff --git a/polymatrix/expression/expression.py b/polymatrix/expression/expression.py
index 25a8c9b..ca7d474 100644
--- a/polymatrix/expression/expression.py
+++ b/polymatrix/expression/expression.py
@@ -4,48 +4,12 @@ import dataclassabc
 import typing
 import numpy as np
 
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-from polymatrix.expression.init.initadditionexpr import init_addition_expr
-from polymatrix.expression.init.initcacheexpr import init_cache_expr
-from polymatrix.expression.init.initcombinationsexpr import init_combinations_expr
-from polymatrix.expression.init.initderivativeexpr import init_derivative_expr
-from polymatrix.expression.init.initdeterminantexpr import init_determinant_expr
-from polymatrix.expression.init.initdivergenceexpr import init_divergence_expr
-from polymatrix.expression.init.initdivisionexpr import init_division_expr
-from polymatrix.expression.init.initelemmultexpr import init_elem_mult_expr
-from polymatrix.expression.init.initevalexpr import init_eval_expr
-from polymatrix.expression.init.initfilterexpr import init_filter_expr
-from polymatrix.expression.init.initfromsymmetricmatrixexpr import init_from_symmetric_matrix_expr
-from polymatrix.expression.init.inithalfnewtonpolytopeexpr import init_half_newton_polytope_expr
-from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr
-from polymatrix.expression.init.initfromsympyexpr import init_from_expr_or_none, init_from_sympy_expr
-from polymatrix.expression.init.initgetitemexpr import init_get_item_expr
-from polymatrix.expression.init.initlinearmatrixinexpr import init_linear_matrix_in_expr
-from polymatrix.expression.init.initlinearmonomialsexpr import init_linear_monomials_expr
-from polymatrix.expression.init.initmatrixmultexpr import init_matrix_mult_expr
-from polymatrix.expression.init.initmaxdegreeexpr import init_max_degree_expr
-from polymatrix.expression.init.initmaxexpr import init_max_expr
-from polymatrix.expression.init.initparametrizeexpr import init_parametrize_expr
-from polymatrix.expression.init.initparametrizematrixexpr import init_parametrize_matrix_expr
-from polymatrix.expression.init.initquadraticinexpr import init_quadratic_in_expr
-from polymatrix.expression.init.initrepmatexpr import init_rep_mat_expr
-from polymatrix.expression.init.initreshapeexpr import init_reshape_expr
-from polymatrix.expression.init.initsetelementatexpr import init_set_element_at_expr
-from polymatrix.expression.init.initquadraticmonomialsexpr import init_quadratic_monomials_expr
-from polymatrix.expression.init.initsqueezeexpr import init_squeeze_expr
-from polymatrix.expression.init.initsubstituteexpr import init_substitute_expr
-from polymatrix.expression.init.initsubtractmonomialsexpr import init_subtract_monomials_expr
-from polymatrix.expression.init.initsumexpr import init_sum_expr
-from polymatrix.expression.init.initsymmetricexpr import init_symmetric_expr
-from polymatrix.expression.init.inittoconstantexpr import init_to_constant_expr
-from polymatrix.expression.init.inittoquadraticexpr import init_to_quadratic_expr
-from polymatrix.expression.init.initdiagexpr import init_diag_expr
-from polymatrix.expression.init.inittosortedvariables import init_to_sorted_variables
-from polymatrix.expression.init.inittosymmetricmatrixexpr import init_to_symmetric_matrix_expr
-from polymatrix.expression.init.inittransposeexpr import init_transpose_expr
-from polymatrix.expression.init.inittruncateexpr import init_truncate_expr
+import polymatrix.expression.initexpressionbase
 
+from polymatrix.utils.getstacklines import get_stack_lines
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+from polymatrix.expression.op import diff, linear_in, linear_monomials, legendre
 from polymatrix.polymatrix.polymatrix import PolyMatrix
 from polymatrix.expressionstate.expressionstate import ExpressionState
 
@@ -65,7 +29,7 @@ class Expression(
         return self.underlying.apply(state)
 
     def __add__(self, other: ExpressionBaseMixin) -> 'Expression':
-        return self._binary(init_addition_expr, self, other)
+        return self._binary(polymatrix.expression.initexpressionbase.init_addition_expr, self, other)
 
     def __getattr__(self, name):
         attr = getattr(self.underlying, name)
@@ -82,17 +46,17 @@ class Expression(
     def __getitem__(self, key: tuple[int, int]):
         return dataclasses.replace(
             self, 
-            underlying=init_get_item_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_get_item_expr(
                 underlying=self.underlying,
                 index=key,
             ),
         )     
 
     def __matmul__(self, other: typing.Union[ExpressionBaseMixin, np.ndarray]) -> 'Expression':
-        return self._binary(init_matrix_mult_expr, self, other)
+        return self._binary(polymatrix.expression.initexpressionbase.init_matrix_mult_expr, self, other)
 
     def __mul__(self, other) -> 'Expression':
-        return self._binary(init_elem_mult_expr, self, other)
+        return self._binary(polymatrix.expression.initexpressionbase.init_elem_mult_expr, self, other)
 
     def __pow__(self, num):
         curr = 1
@@ -106,10 +70,10 @@ class Expression(
         return self * (-1)
 
     def __radd__(self, other):
-        return self._binary(init_addition_expr, other, self)
+        return self._binary(polymatrix.expression.initexpressionbase.init_addition_expr, other, self)
 
     def __rmatmul__(self, other):
-        return self._binary(init_matrix_mult_expr, other, self)
+        return self._binary(polymatrix.expression.initexpressionbase.init_matrix_mult_expr, other, self)
 
     def __rmul__(self, other):
         return self * other
@@ -121,30 +85,37 @@ class Expression(
         return self + (-other)
 
     def __truediv__(self, other: ExpressionBaseMixin):
-        return self._binary(init_division_expr, self, other)
+        return self._binary(polymatrix.expression.initexpressionbase.init_division_expr, self, other)
+
+    @abc.abstractmethod
+    def copy(self, underlying: ExpressionBaseMixin) -> 'Expression':
+        ...
 
     @staticmethod
     def _binary(op, left, right):
-        if not isinstance(left, Expression):
-            left = init_from_expr_or_none(left)
 
-            if left is None:
+        stack = get_stack_lines()
+
+        if isinstance(left, Expression):
+            right = polymatrix.expression.initexpressionbase.init_from_expr_or_none(right)
+
+            # delegate to upper level
+            if right is None:
                 return NotImplemented
 
-            return dataclasses.replace(
-                right, 
-                underlying=op(left, right),
+            return left.copy(
+                underlying=op(left, right, stack),
             )
 
         else:
-            right = init_from_expr_or_none(right)
-            
-            if right is None:
+            left = polymatrix.expression.initexpressionbase.init_from_expr_or_none(left)
+
+            # delegate to upper level
+            if left is None:
                 return NotImplemented
 
-            return dataclasses.replace(
-                left, 
-                underlying=op(left, right),
+            return right.copy(
+                underlying=op(left, right, stack),
             )
 
     # def _convert_to_expression(self, other):       
@@ -156,27 +127,27 @@ class Expression(
     #     return result
 
     def cache(self) -> 'Expression':
-        return dataclasses.replace(
-            self, 
-            underlying=init_cache_expr(
+        return self.copy(
+            underlying=polymatrix.expression.initexpressionbase.init_cache_expr(
                 underlying=self.underlying,
             ),
         )
 
-    # only applies to monomials and variables
-    def combinations(self, number: int):
-        return dataclasses.replace(
+    def combinations(
             self, 
-            underlying=init_combinations_expr(
-                monomials=self.underlying,
-                number=number,
+            degrees: tuple[int, ...] | int,
+        ):
+        return self.copy(
+            underlying=polymatrix.expression.initexpressionbase.init_combinations_expr(
+                expression=self.underlying,
+                degrees=degrees,
             ),
         )
 
     def determinant(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_determinant_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_determinant_expr(
                 underlying=self.underlying,
             ),
         )
@@ -184,7 +155,7 @@ class Expression(
     def diag(self):
         return dataclasses.replace(
             self, 
-            underlying=init_diag_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_diag_expr(
                 underlying=self.underlying,
             ),
         )
@@ -196,8 +167,8 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_derivative_expr(
-                underlying=self.underlying,
+            underlying=diff(
+                expression=self,
                 variables=variables,
                 introduce_derivatives=introduce_derivatives,
             ),
@@ -209,7 +180,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_divergence_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_divergence_expr(
                 underlying=self.underlying,
                 variables=variables,
             ),
@@ -222,7 +193,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_eval_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_eval_expr(
                 underlying=self.underlying,
                 variables=variable,
                 values=value,
@@ -237,7 +208,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_filter_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_filter_expr(
                 underlying=self.underlying,
                 predicator=predicator,
                 inverse=inverse,
@@ -248,7 +219,7 @@ class Expression(
     def from_symmetric_matrix(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_from_symmetric_matrix_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_symmetric_matrix_expr(
                 underlying=self.underlying,
             ),
         )
@@ -261,7 +232,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_half_newton_polytope_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_half_newton_polytope_expr(
                 monomials=self.underlying,
                 variables=variables,
                 filter=filter,
@@ -271,17 +242,21 @@ class Expression(
     def linear_matrix_in(self, variable: 'Expression') -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_linear_matrix_in_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_linear_matrix_in_expr(
                 underlying=self.underlying,
                 variable=variable,
             ),
         )
 
-    def linear_monomials(self, variables: 'Expression') -> 'Expression':
+    def linear_monomials(
+            self, 
+            variables: 'Expression',
+        ) -> 'Expression':
+        
         return dataclasses.replace(
             self, 
-            underlying=init_linear_monomials_expr(
-                underlying=self.underlying,
+            underlying=linear_monomials(
+                expression=self.underlying,
                 variables=variables,
             ),
         )
@@ -292,23 +267,33 @@ class Expression(
         monomials: 'Expression' = None,
         ignore_unmatched: bool = None,    
     ) -> 'Expression':
-        if monomials is None:
-            monomials = self.linear_monomials(variables)
 
         return dataclasses.replace(
             self, 
-            underlying=init_linear_in_expr(
-                underlying=self.underlying,
+            underlying=linear_in(
+                expression=self.underlying,
                 monomials=monomials,
                 variables=variables,
                 ignore_unmatched=ignore_unmatched,
             ),
         )
+    
+    def legendre(
+        self,
+        degrees: tuple[int, ...] = None,
+    ) -> 'Expression':
+        return dataclasses.replace(
+            self, 
+            underlying=legendre(
+                expression=self.underlying,
+                degrees=degrees,
+            ),
+        )
 
     def max(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_max_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_max_expr(
                 underlying=self.underlying,
             ),
         )
@@ -316,7 +301,7 @@ class Expression(
     def max_degree(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_max_degree_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_max_degree_expr(
                 underlying=self.underlying,
             ),
         )
@@ -324,7 +309,7 @@ class Expression(
     def parametrize(self, name: str = None) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_parametrize_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_parametrize_expr(
                 underlying=self.underlying,
                 name=name,
             ),
@@ -343,12 +328,15 @@ class Expression(
         if monomials is None:
             monomials = self.quadratic_monomials(variables)
 
+        stack = get_stack_lines()
+
         return dataclasses.replace(
             self, 
-            underlying=init_quadratic_in_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_quadratic_in_expr(
                 underlying=self.underlying,
                 monomials=monomials,
                 variables=variables,
+                stack=stack,
             ),
         )
 
@@ -358,7 +346,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_quadratic_monomials_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_quadratic_monomials_expr(
                 underlying=self.underlying,
                 variables=variables,
             ),
@@ -367,7 +355,7 @@ class Expression(
     def reshape(self, n: int, m: int) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_reshape_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_reshape_expr(
                 underlying=self.underlying,
                 new_shape=(n, m),
             ),
@@ -376,7 +364,7 @@ class Expression(
     def rep_mat(self, n: int, m: int) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_rep_mat_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_rep_mat_expr(
                 underlying=self.underlying,
                 repetition=(n, m),
             ),
@@ -391,11 +379,11 @@ class Expression(
         if isinstance(value, Expression):
             value = value.underlying
         else:
-            value = init_from_sympy_expr(value)
+            value = polymatrix.expression.initexpressionbase.init_from_expr(value)
 
         return dataclasses.replace(
             self, 
-            underlying=init_set_element_at_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_set_element_at_expr(
                 underlying=self.underlying,
                 index=(row, col),
                 value=value,
@@ -408,7 +396,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_squeeze_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_squeeze_expr(
                 underlying=self.underlying,
             ),
         )
@@ -420,7 +408,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_subtract_monomials_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_subtract_monomials_expr(
                 underlying=self.underlying,
                 monomials=monomials,
             ),
@@ -433,7 +421,7 @@ class Expression(
     ) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_substitute_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_substitute_expr(
                 underlying=self.underlying,
                 variables=variable,
                 values=values,
@@ -453,7 +441,7 @@ class Expression(
     def sum(self):
         return dataclasses.replace(
             self, 
-            underlying=init_sum_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_sum_expr(
                 underlying=self.underlying,
             ),
         )
@@ -461,7 +449,7 @@ class Expression(
     def symmetric(self):
         return dataclasses.replace(
             self, 
-            underlying=init_symmetric_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_symmetric_expr(
                 underlying=self.underlying,
             ),
         )
@@ -469,7 +457,7 @@ class Expression(
     def transpose(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_transpose_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_transpose_expr(
                 underlying=self.underlying,
             ),
         )
@@ -481,7 +469,7 @@ class Expression(
     def to_constant(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_to_constant_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_to_constant_expr(
                 underlying=self.underlying,
             ),
         )
@@ -497,7 +485,7 @@ class Expression(
     def to_symmetric_matrix(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_to_symmetric_matrix_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_to_symmetric_matrix_expr(
                 underlying=self.underlying,
             ),
         )
@@ -506,7 +494,7 @@ class Expression(
     def to_sorted_variables(self) -> 'Expression':
         return dataclasses.replace(
             self, 
-            underlying=init_to_sorted_variables(
+            underlying=polymatrix.expression.initexpressionbase.init_to_sorted_variables(
                 underlying=self.underlying,
             ),
         )
@@ -520,7 +508,7 @@ class Expression(
     ):
         return dataclasses.replace(
             self, 
-            underlying=init_truncate_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_truncate_expr(
                 underlying=self.underlying,
                 variables=variables,
                 degrees=degrees,
@@ -531,7 +519,21 @@ class Expression(
 
 @dataclassabc.dataclassabc(frozen=True, repr=False)
 class ExpressionImpl(Expression):
-	underlying: ExpressionBaseMixin
+    underlying: ExpressionBaseMixin
+
+    def __repr__(self) -> str:
+        return self.underlying.__repr__()
+    
+    def copy(self, underlying: ExpressionBaseMixin) -> 'Expression':
+        return dataclasses.replace(
+            self, 
+            underlying=underlying,
+        )
 
-	def __repr__(self) -> str:
-		return self.underlying.__repr__()
+
+def init_expression(
+    underlying: ExpressionBaseMixin,
+):
+    return ExpressionImpl(
+        underlying=underlying,
+    )
diff --git a/polymatrix/expression/from_.py b/polymatrix/expression/from_.py
new file mode 100644
index 0000000..637c9d1
--- /dev/null
+++ b/polymatrix/expression/from_.py
@@ -0,0 +1,56 @@
+import numpy as np
+import sympy
+
+import polymatrix.expression.initexpressionbase
+
+from polymatrix.expression.expression import init_expression, Expression
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+
+
+DATA_TYPE = str | np.ndarray | sympy.Matrix | sympy.Expr | tuple | ExpressionBaseMixin
+
+def from_expr_or_none(
+	data: DATA_TYPE,
+) -> Expression | None:
+	
+	return init_expression(
+        underlying=polymatrix.expression.initexpressionbase.init_from_expr_or_none(
+            data=data,
+        ),
+    )
+
+def from_(
+    data: DATA_TYPE,
+) -> Expression:
+	
+    return init_expression(
+        underlying=polymatrix.expression.initexpressionbase.init_from_expr(
+            data=data,
+        ),
+    )
+
+# def from_expr(
+#     data: DATA_TYPE,
+# ) -> Expression:
+# 	return from_(data=data)
+
+# def from_sympy(
+#     data: tuple[tuple[float]],
+# ):
+#     return init_expression(
+#         polymatrix.expression.init.init_from_sympy_expr(data)
+#     )
+
+# def from_state_monad(
+#     data: StateMonad,
+# ):
+#     return init_expression(
+#         data.flat_map(lambda inner_data: polymatrix.expression.init.init_from_sympy_expr(inner_data)),
+#     )
+
+# def from_polymatrix(
+#     polymatrix: PolyMatrix,
+# ):
+#     return init_expression(
+#         polymatrix.expression.init.init_from_terms_expr(polymatrix)
+#     )
diff --git a/polymatrix/expression/impl/reshapeexprmixin.py b/polymatrix/expression/impl/reshapeexprmixin.py
deleted file mode 100644
index 00a93b8..0000000
--- a/polymatrix/expression/impl/reshapeexprmixin.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import dataclassabc
-
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-from polymatrix.expression.mixins.reshapeexprmixin import ReshapeExprMixin
-
-
-@dataclassabc.dataclassabc(frozen=True)
-class ReshapeExprImpl(ReshapeExprMixin):
-	underlying: ExpressionBaseMixin
-	new_shape: tuple
diff --git a/polymatrix/expression/impl/impl.py b/polymatrix/expression/implexpressionbase.py
index 7a919e5..79761f2 100644
--- a/polymatrix/expression/impl/impl.py
+++ b/polymatrix/expression/implexpressionbase.py
@@ -1,3 +1,6 @@
+from polymatrix.expression.mixins.legendreseriesmixin import LegendreSeriesMixin
+from polymatrix.expression.mixins.productexprmixin import ProductExprMixin
+from polymatrix.utils.getstacklines import FrameSummary
 import dataclassabc
 
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
@@ -20,7 +23,7 @@ from polymatrix.expression.mixins.eyeexprmixin import EyeExprMixin
 from polymatrix.expression.mixins.filterexprmixin import FilterExprMixin
 from polymatrix.expression.mixins.fromsymmetricmatrixexprmixin import \
     FromSymmetricMatrixExprMixin
-from polymatrix.expression.mixins.fromsympyexprmixin import FromSympyExprMixin
+from polymatrix.expression.mixins.fromtupleexprmixin import FromTupleExprMixin
 from polymatrix.expression.mixins.fromtermsexprmixin import FromTermsExprMixin
 from polymatrix.expression.mixins.getitemexprmixin import GetItemExprMixin
 from polymatrix.expression.mixins.halfnewtonpolytopeexprmixin import \
@@ -69,14 +72,16 @@ from polymatrix.expression.mixins.tosortedvariablesmixin import ToSortedVariable
 class AdditionExprImpl(AdditionExprMixin):
 	left: ExpressionBaseMixin
 	right: ExpressionBaseMixin
+	stack: tuple[FrameSummary]
 
+	def __repr__(self):
+		return f'{self.__class__.__name__}(left={self.left}, right={self.right})'
 
 @dataclassabc.dataclassabc(frozen=True)
 class BlockDiagExprImpl(BlockDiagExprMixin):
-	underlying: tuple
+	underlying: tuple[ExpressionBaseMixin]
 
 
-# @dataclassabc.dataclassabc(frozen=True, eq=False)
 @dataclassabc.dataclassabc(frozen=True)
 class CacheExprImpl(CacheExprMixin):
 	underlying: ExpressionBaseMixin
@@ -84,8 +89,8 @@ class CacheExprImpl(CacheExprMixin):
 
 @dataclassabc.dataclassabc(frozen=True)
 class CombinationsExprImpl(CombinationsExprMixin):
-	monomials: ExpressionBaseMixin
-	number: int
+	expression: ExpressionBaseMixin
+	degrees: tuple[int, ...]
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -93,6 +98,10 @@ class DerivativeExprImpl(DerivativeExprMixin):
 	underlying: ExpressionBaseMixin
 	variables: tuple
 	introduce_derivatives: bool
+	stack: tuple[FrameSummary]
+
+	def __repr__(self):
+		return f'{self.__class__.__name__}(variables={self.variables}, underlying={self.underlying})'
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -110,11 +119,11 @@ class DivergenceExprImpl(DivergenceExprMixin):
 	underlying: ExpressionBaseMixin
 	variables: tuple
 
-
-@dataclassabc.dataclassabc(frozen=True, eq=False)
+@dataclassabc.dataclassabc(frozen=True)
 class DivisionExprImpl(DivisionExprMixin):
 	left: ExpressionBaseMixin
 	right: ExpressionBaseMixin
+	stack: tuple[FrameSummary]
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -147,9 +156,12 @@ class FromSymmetricMatrixExprImpl(FromSymmetricMatrixExprMixin):
 
 
 @dataclassabc.dataclassabc(frozen=True)
-class FromSympyExprImpl(FromSympyExprMixin):
+class FromTupleExprImpl(FromTupleExprMixin):
 	data: tuple[tuple[float]]
+	stack: tuple[FrameSummary]
 
+	def __repr__(self):
+		return f'{self.__class__.__name__}(data={self.data})'
 
 @dataclassabc.dataclassabc(frozen=True)
 class FromTermsExprImpl(FromTermsExprMixin):
@@ -191,9 +203,20 @@ class LinearMonomialsExprImpl(LinearMonomialsExprMixin):
 
 
 @dataclassabc.dataclassabc(frozen=True)
+class LegendreSeriesImpl(LegendreSeriesMixin):
+	underlying: ExpressionBaseMixin
+	degrees: tuple[int, ...] | None
+	stack: tuple[FrameSummary]
+
+
+@dataclassabc.dataclassabc(frozen=True)
 class MatrixMultExprImpl(MatrixMultExprMixin):
 	left: ExpressionBaseMixin
 	right: ExpressionBaseMixin
+	stack: tuple[FrameSummary]
+
+	def __repr__(self):
+		return f'{self.__class__.__name__}(left={self.left}, right={self.right})'
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -206,14 +229,13 @@ class MaxExprImpl(MaxExprMixin):
 	underlying: ExpressionBaseMixin
 
 
-# @dataclassabc.dataclassabc(frozen=True, eq=False, repr=False)
 @dataclassabc.dataclassabc(frozen=True, repr=False)
 class ParametrizeExprImpl(ParametrizeExprMixin):
 	underlying: ExpressionBaseMixin
 	name: str
 
 	def __repr__(self) -> str:
-		return f'{self.__class__.__name__}(name={repr(self.name)})'
+		return f'{self.__class__.__name__}(name={self.name})'
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -223,10 +245,21 @@ class ParametrizeMatrixExprImpl(ParametrizeMatrixExprMixin):
 
 
 @dataclassabc.dataclassabc(frozen=True)
+class ProductExprImpl(ProductExprMixin):
+	underlying: tuple[ExpressionBaseMixin]
+	degrees: tuple[int, ...] | None
+	stack: tuple[FrameSummary]
+
+
+@dataclassabc.dataclassabc(frozen=True)
 class QuadraticInExprImpl(QuadraticInExprMixin):
 	underlying: ExpressionBaseMixin
 	monomials: ExpressionBaseMixin
 	variables: tuple
+	stack: tuple[FrameSummary]
+
+	def __repr__(self):
+		return f'{self.__class__.__name__}(variables={self.variables}, monomials={self.monomials}, underlying={self.underlying})'
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -241,10 +274,10 @@ class RepMatExprImpl(RepMatExprMixin):
 	repetition: tuple
 
 
-# @dataclassabc.dataclassabc(frozen=True)
-# class ReshapeExprImpl(ReshapeExprMixin):
-# 	underlying: ExpressionBaseMixin
-# 	new_shape: tuple
+@dataclassabc.dataclassabc(frozen=True)
+class ReshapeExprImpl(ReshapeExprMixin):
+	underlying: ExpressionBaseMixin
+	new_shape: tuple
 
 
 @dataclassabc.dataclassabc(frozen=True)
@@ -262,7 +295,6 @@ class SqueezeExprImpl(SqueezeExprMixin):
 @dataclassabc.dataclassabc(frozen=True)
 class SubstituteExprImpl(SubstituteExprMixin):
 	underlying: ExpressionBaseMixin
-	# variables: tuple
 	substitutions: tuple
 
 
diff --git a/polymatrix/expression/init/initadditionexpr.py b/polymatrix/expression/init/initadditionexpr.py
deleted file mode 100644
index cdc97b1..0000000
--- a/polymatrix/expression/init/initadditionexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_addition_expr(
-		left: ExpressionBaseMixin,
-		right: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.AdditionExprImpl(
-		left=left,
-		right=right,
-)
diff --git a/polymatrix/expression/init/initblockdiagexpr.py b/polymatrix/expression/init/initblockdiagexpr.py
deleted file mode 100644
index a9ad8a6..0000000
--- a/polymatrix/expression/init/initblockdiagexpr.py
+++ /dev/null
@@ -1,7 +0,0 @@
-import polymatrix.expression.impl.impl
-def init_block_diag_expr(
-		underlying: tuple,
-):
-	return polymatrix.expression.impl.impl.BlockDiagExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initcacheexpr.py b/polymatrix/expression/init/initcacheexpr.py
deleted file mode 100644
index 861ffa1..0000000
--- a/polymatrix/expression/init/initcacheexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_cache_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.CacheExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initcombinationsexpr.py b/polymatrix/expression/init/initcombinationsexpr.py
deleted file mode 100644
index 10931e6..0000000
--- a/polymatrix/expression/init/initcombinationsexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_combinations_expr(
-		monomials: ExpressionBaseMixin,
-		number: int,
-):
-	return polymatrix.expression.impl.impl.CombinationsExprImpl(
-		monomials=monomials,
-		number=number,
-)
diff --git a/polymatrix/expression/init/initderivativeexpr.py b/polymatrix/expression/init/initderivativeexpr.py
deleted file mode 100644
index b1c9405..0000000
--- a/polymatrix/expression/init/initderivativeexpr.py
+++ /dev/null
@@ -1,20 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_derivative_expr(
-		underlying: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-		introduce_derivatives: bool = None,
-):
-
-	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
-
-	if introduce_derivatives is None:
-		introduce_derivatives = False
-
-	return polymatrix.expression.impl.impl.DerivativeExprImpl(
-		underlying=underlying,
-		variables=variables,
-		introduce_derivatives=introduce_derivatives,
-)
diff --git a/polymatrix/expression/init/initdeterminantexpr.py b/polymatrix/expression/init/initdeterminantexpr.py
deleted file mode 100644
index 660f3ae..0000000
--- a/polymatrix/expression/init/initdeterminantexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_determinant_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.DeterminantExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initdiagexpr.py b/polymatrix/expression/init/initdiagexpr.py
deleted file mode 100644
index 794c5bf..0000000
--- a/polymatrix/expression/init/initdiagexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_diag_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.DiagExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initdivergenceexpr.py b/polymatrix/expression/init/initdivergenceexpr.py
deleted file mode 100644
index 6f5d99d..0000000
--- a/polymatrix/expression/init/initdivergenceexpr.py
+++ /dev/null
@@ -1,13 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-from typing import Union
-import polymatrix.expression.impl.impl
-
-
-def init_divergence_expr(
-		underlying: ExpressionBaseMixin,
-		variables: tuple,
-):
-	return polymatrix.expression.impl.impl.DivergenceExprImpl(
-		underlying=underlying,
-		variables=variables,
-)
diff --git a/polymatrix/expression/init/initdivisionexpr.py b/polymatrix/expression/init/initdivisionexpr.py
deleted file mode 100644
index d7870bc..0000000
--- a/polymatrix/expression/init/initdivisionexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_division_expr(
-		left: ExpressionBaseMixin,
-		right: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.DivisionExprImpl(
-		left=left,
-		right=right,
-)
diff --git a/polymatrix/expression/init/initelemmultexpr.py b/polymatrix/expression/init/initelemmultexpr.py
deleted file mode 100644
index ef66c68..0000000
--- a/polymatrix/expression/init/initelemmultexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_elem_mult_expr(
-		left: ExpressionBaseMixin,
-		right: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.ElemMultExprImpl(
-		left=left,
-		right=right,
-)
diff --git a/polymatrix/expression/init/initevalexpr.py b/polymatrix/expression/init/initevalexpr.py
deleted file mode 100644
index 29359d6..0000000
--- a/polymatrix/expression/init/initevalexpr.py
+++ /dev/null
@@ -1,37 +0,0 @@
-import typing
-import numpy as np
-from polymatrix.expression.init.initsubstituteexpr import format_substitutions
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_eval_expr(
-		underlying: ExpressionBaseMixin,
-		variables: typing.Union[typing.Any, tuple, dict],
-		values: typing.Union[float, tuple] = None,
-):
-
-	substitutions = format_substitutions(
-		variables=variables, 
-		values=values,
-	)
-
-	def formatted_values(value):
-		if isinstance(value, np.ndarray):
-			return tuple(value.reshape(-1))
-
-		elif isinstance(value, tuple):
-			return value
-
-		elif isinstance(value, int) or isinstance(value, float):
-			return (value,)
-
-		else:
-			return (float(value),)
-
-	substitutions = tuple((variable, formatted_values(value)) for variable, value in substitutions)
-
-	return polymatrix.expression.impl.impl.EvalExprImpl(
-		underlying=underlying,
-		substitutions=substitutions,
-	)
diff --git a/polymatrix/expression/init/initexpression.py b/polymatrix/expression/init/initexpression.py
deleted file mode 100644
index 3c3258c..0000000
--- a/polymatrix/expression/init/initexpression.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.expression import ExpressionImpl
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-
-
-def init_expression(
-		underlying: ExpressionBaseMixin,
-):
-	return ExpressionImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initeyeexpr.py b/polymatrix/expression/init/initeyeexpr.py
deleted file mode 100644
index 642e4c6..0000000
--- a/polymatrix/expression/init/initeyeexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_eye_expr(
-		variable: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.EyeExprImpl(
-		variable=variable,
-)
diff --git a/polymatrix/expression/init/initfilterexpr.py b/polymatrix/expression/init/initfilterexpr.py
deleted file mode 100644
index 7164183..0000000
--- a/polymatrix/expression/init/initfilterexpr.py
+++ /dev/null
@@ -1,17 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_filter_expr(
-		underlying: ExpressionBaseMixin,
-		predicator: ExpressionBaseMixin,
-		inverse: bool == None,
-):
-	if inverse is None:
-		inverse = False
-
-	return polymatrix.expression.impl.impl.FilterExprImpl(
-		underlying=underlying,
-		predicator=predicator,
-		inverse=inverse,
-)
diff --git a/polymatrix/expression/init/initfromsymmetricmatrixexpr.py b/polymatrix/expression/init/initfromsymmetricmatrixexpr.py
deleted file mode 100644
index b3ef782..0000000
--- a/polymatrix/expression/init/initfromsymmetricmatrixexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_from_symmetric_matrix_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.FromSymmetricMatrixExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initfromsympyexpr.py b/polymatrix/expression/init/initfromsympyexpr.py
deleted file mode 100644
index 4ca6fcf..0000000
--- a/polymatrix/expression/init/initfromsympyexpr.py
+++ /dev/null
@@ -1,68 +0,0 @@
-import typing
-import numpy as np
-import sympy
-
-from polymatrix.expression.impl.impl import FromSympyExprImpl
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-
-
-def init_from_expr_or_none(
-	data: typing.Union[np.ndarray, tuple[tuple[float]]],
-) -> ExpressionBaseMixin | None:
-	match data:
-		case np.ndarray():
-			assert len(data.shape) <= 2
-
-			def gen_elements():
-				for row in data:
-					if isinstance(row, np.ndarray):
-						yield tuple(row)
-					else:
-						yield (row,)
-
-			data = tuple(gen_elements())
-
-		case sympy.Matrix():
-			data = tuple(tuple(i for i in data.row(row)) for row in range(data.rows))
-
-		case tuple():
-
-			match data[0]:
-
-				case tuple():
-					n_col = len(data[0])
-					assert all(len(col) == n_col for col in data)
-
-				case _:
-					data = tuple((e,) for e in data)
-
-		case sympy.Expr():
-			data = ((sympy.expand(data),),)
-
-		case ExpressionBaseMixin():
-			return data
-
-		case _:
-			if not isinstance(data, (float, int, np.number)):
-			# 	raise Exception(f'{data=}, {type(data)=}')
-				return None
-
-			data = ((data,),)
-
-	return FromSympyExprImpl(
-			data=data,
-	)
-
-
-def init_from_sympy_expr(
-		data: typing.Union[np.ndarray, tuple[tuple[float]]],
-):
-	# if isinstance(data, ExpressionBaseMixin):
-	# 	return data
-
-	expr = init_from_expr_or_none(data)
-
-	if expr is None:
-		raise Exception(f'{data=}')
-
-	return expr
diff --git a/polymatrix/expression/init/initfromtermsexpr.py b/polymatrix/expression/init/initfromtermsexpr.py
deleted file mode 100644
index b2c1dc3..0000000
--- a/polymatrix/expression/init/initfromtermsexpr.py
+++ /dev/null
@@ -1,33 +0,0 @@
-import typing
-from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
-import polymatrix.expression.impl.impl
-
-
-def init_from_terms_expr(
-	terms: typing.Union[tuple, PolyMatrixMixin],
-	shape: tuple[int, int] = None,
-):
-
-	if isinstance(terms, PolyMatrixMixin):
-		shape = terms.shape
-		gen_terms = terms.gen_terms()
-
-	else:
-		assert shape is not None
-
-		if isinstance(terms, tuple):
-			gen_terms = terms
-
-		elif isinstance(terms, dict):
-			gen_terms = terms.items()
-
-		else:
-			raise Exception(f'{terms=}')
-
-	# Expression needs to be hashable
-	terms_formatted = tuple((key, tuple(monomials.items())) for key, monomials in gen_terms)
-
-	return polymatrix.expression.impl.impl.FromTermsExprImpl(
-		terms=terms_formatted,
-		shape=shape,
-	)
diff --git a/polymatrix/expression/init/initgetitemexpr.py b/polymatrix/expression/init/initgetitemexpr.py
deleted file mode 100644
index f7fc0c9..0000000
--- a/polymatrix/expression/init/initgetitemexpr.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from numpy import isin
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_get_item_expr(
-		underlying: ExpressionBaseMixin,
-		index: tuple[tuple[int, ...], tuple[int, ...]],
-):
-
-	def get_hashable_slice(index):
-		if isinstance(index, slice):
-			return polymatrix.expression.impl.impl.GetItemExprImpl.Slice(start=index.start, stop=index.stop, step=index.step)
-		else:
-			return index
-
-	proper_index = (get_hashable_slice(index[0]), get_hashable_slice(index[1]))
-
-	return polymatrix.expression.impl.impl.GetItemExprImpl(
-		underlying=underlying,
-		index=proper_index,
-)
diff --git a/polymatrix/expression/init/inithalfnewtonpolytopeexpr.py b/polymatrix/expression/init/inithalfnewtonpolytopeexpr.py
deleted file mode 100644
index 794c938..0000000
--- a/polymatrix/expression/init/inithalfnewtonpolytopeexpr.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_half_newton_polytope_expr(
-		monomials: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-		filter: ExpressionBaseMixin | None = None,
-):
-	return polymatrix.expression.impl.impl.HalfNewtonPolytopeExprImpl(
-		monomials=monomials,
-		variables=variables,
-		filter=filter
-)
diff --git a/polymatrix/expression/init/initlinearinexpr.py b/polymatrix/expression/init/initlinearinexpr.py
deleted file mode 100644
index 30705e2..0000000
--- a/polymatrix/expression/init/initlinearinexpr.py
+++ /dev/null
@@ -1,18 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_linear_in_expr(
-		underlying: ExpressionBaseMixin,
-		monomials: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-		ignore_unmatched: bool = None,
-):
-	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
-
-	return polymatrix.expression.impl.impl.LinearInExprImpl(
-		underlying=underlying,
-		monomials=monomials,
-		variables=variables,
-		ignore_unmatched = ignore_unmatched,
-)
diff --git a/polymatrix/expression/init/initlinearmatrixinexpr.py b/polymatrix/expression/init/initlinearmatrixinexpr.py
deleted file mode 100644
index 83e9aea..0000000
--- a/polymatrix/expression/init/initlinearmatrixinexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_linear_matrix_in_expr(
-		underlying: ExpressionBaseMixin,
-		variable: int,
-):
-	return polymatrix.expression.impl.impl.LinearMatrixInExprImpl(
-		underlying=underlying,
-		variable=variable,
-)
diff --git a/polymatrix/expression/init/initlinearmonomialsexpr.py b/polymatrix/expression/init/initlinearmonomialsexpr.py
deleted file mode 100644
index f830542..0000000
--- a/polymatrix/expression/init/initlinearmonomialsexpr.py
+++ /dev/null
@@ -1,15 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_linear_monomials_expr(
-		underlying: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-):
-
-	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
-
-	return polymatrix.expression.impl.impl.LinearMonomialsExprImpl(
-		underlying=underlying,
-		variables=variables,
-	)
diff --git a/polymatrix/expression/init/initmatrixmultexpr.py b/polymatrix/expression/init/initmatrixmultexpr.py
deleted file mode 100644
index 20bffec..0000000
--- a/polymatrix/expression/init/initmatrixmultexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_matrix_mult_expr(
-		left: ExpressionBaseMixin,
-		right: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.MatrixMultExprImpl(
-		left=left,
-		right=right,
-)
diff --git a/polymatrix/expression/init/initmaxdegreeexpr.py b/polymatrix/expression/init/initmaxdegreeexpr.py
deleted file mode 100644
index d232163..0000000
--- a/polymatrix/expression/init/initmaxdegreeexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_max_degree_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.MaxDegreeExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initmaxexpr.py b/polymatrix/expression/init/initmaxexpr.py
deleted file mode 100644
index 41dde90..0000000
--- a/polymatrix/expression/init/initmaxexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_max_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.MaxExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initparametrizeexpr.py b/polymatrix/expression/init/initparametrizeexpr.py
deleted file mode 100644
index 178d963..0000000
--- a/polymatrix/expression/init/initparametrizeexpr.py
+++ /dev/null
@@ -1,15 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_parametrize_expr(
-		underlying: ExpressionBaseMixin,
-		name: str = None,
-):
-	if name is None:
-		name = 'undefined'
-
-	return polymatrix.expression.impl.impl.ParametrizeExprImpl(
-		underlying=underlying,
-		name=name,
-)
diff --git a/polymatrix/expression/init/initparametrizematrixexpr.py b/polymatrix/expression/init/initparametrizematrixexpr.py
deleted file mode 100644
index bf90510..0000000
--- a/polymatrix/expression/init/initparametrizematrixexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_parametrize_matrix_expr(
-		underlying: ExpressionBaseMixin,
-		name: str,
-):
-	return polymatrix.expression.impl.impl.ParametrizeMatrixExprImpl(
-		underlying=underlying,
-		name=name,
-)
diff --git a/polymatrix/expression/init/initquadraticinexpr.py b/polymatrix/expression/init/initquadraticinexpr.py
deleted file mode 100644
index 9569cd3..0000000
--- a/polymatrix/expression/init/initquadraticinexpr.py
+++ /dev/null
@@ -1,17 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_quadratic_in_expr(
-		underlying: ExpressionBaseMixin,
-		monomials: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-):
-
-	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
-
-	return polymatrix.expression.impl.impl.QuadraticInExprImpl(
-		underlying=underlying,
-		monomials=monomials,
-		variables=variables,
-	)
diff --git a/polymatrix/expression/init/initquadraticmonomialsexpr.py b/polymatrix/expression/init/initquadraticmonomialsexpr.py
deleted file mode 100644
index c965312..0000000
--- a/polymatrix/expression/init/initquadraticmonomialsexpr.py
+++ /dev/null
@@ -1,15 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_quadratic_monomials_expr(
-		underlying: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-):
-
-	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
-
-	return polymatrix.expression.impl.impl.QuadraticMonomialsExprImpl(
-		underlying=underlying,
-		variables=variables,
-	)
diff --git a/polymatrix/expression/init/initrepmatexpr.py b/polymatrix/expression/init/initrepmatexpr.py
deleted file mode 100644
index 152fff8..0000000
--- a/polymatrix/expression/init/initrepmatexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_rep_mat_expr(
-		underlying: ExpressionBaseMixin,
-		repetition: tuple,
-):
-	return polymatrix.expression.impl.impl.RepMatExprImpl(
-		underlying=underlying,
-		repetition=repetition,
-)
diff --git a/polymatrix/expression/init/initreshapeexpr.py b/polymatrix/expression/init/initreshapeexpr.py
deleted file mode 100644
index 4e6ab57..0000000
--- a/polymatrix/expression/init/initreshapeexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.impl.reshapeexprmixin import ReshapeExprImpl
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-
-
-def init_reshape_expr(
-		underlying: ExpressionBaseMixin,
-		new_shape: tuple,
-):
-	return ReshapeExprImpl(
-		underlying=underlying,
-		new_shape=new_shape,
-)
diff --git a/polymatrix/expression/init/initsetelementatexpr.py b/polymatrix/expression/init/initsetelementatexpr.py
deleted file mode 100644
index 250c35b..0000000
--- a/polymatrix/expression/init/initsetelementatexpr.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_set_element_at_expr(
-		underlying: ExpressionBaseMixin,
-		index: tuple,
-		value: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.SetElementAtExprImpl(
-		underlying=underlying,
-		index=index,
-		value=value,
-)
diff --git a/polymatrix/expression/init/initsqueezeexpr.py b/polymatrix/expression/init/initsqueezeexpr.py
deleted file mode 100644
index d539453..0000000
--- a/polymatrix/expression/init/initsqueezeexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_squeeze_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.SqueezeExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initsubtractmonomialsexpr.py b/polymatrix/expression/init/initsubtractmonomialsexpr.py
deleted file mode 100644
index de720bd..0000000
--- a/polymatrix/expression/init/initsubtractmonomialsexpr.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_subtract_monomials_expr(
-		underlying: ExpressionBaseMixin,
-		monomials: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.SubtractMonomialsExprImpl(
-		underlying=underlying,
-		monomials=monomials,
-)
diff --git a/polymatrix/expression/init/initsumexpr.py b/polymatrix/expression/init/initsumexpr.py
deleted file mode 100644
index 288ae72..0000000
--- a/polymatrix/expression/init/initsumexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_sum_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.SumExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/initsymmetricexpr.py b/polymatrix/expression/init/initsymmetricexpr.py
deleted file mode 100644
index dc321e5..0000000
--- a/polymatrix/expression/init/initsymmetricexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_symmetric_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.SymmetricExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/inittoconstantexpr.py b/polymatrix/expression/init/inittoconstantexpr.py
deleted file mode 100644
index a0084a1..0000000
--- a/polymatrix/expression/init/inittoconstantexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_to_constant_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.ToConstantExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/inittoquadraticexpr.py b/polymatrix/expression/init/inittoquadraticexpr.py
deleted file mode 100644
index 77e6943..0000000
--- a/polymatrix/expression/init/inittoquadraticexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_to_quadratic_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.ToQuadraticExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/inittosortedvariables.py b/polymatrix/expression/init/inittosortedvariables.py
deleted file mode 100644
index 3f5222a..0000000
--- a/polymatrix/expression/init/inittosortedvariables.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_to_sorted_variables(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.ToSortedVariablesImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/inittosymmetricmatrixexpr.py b/polymatrix/expression/init/inittosymmetricmatrixexpr.py
deleted file mode 100644
index 68b27fe..0000000
--- a/polymatrix/expression/init/inittosymmetricmatrixexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_to_symmetric_matrix_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.ToSymmetricMatrixExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/inittransposeexpr.py b/polymatrix/expression/init/inittransposeexpr.py
deleted file mode 100644
index 7807e4d..0000000
--- a/polymatrix/expression/init/inittransposeexpr.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_transpose_expr(
-		underlying: ExpressionBaseMixin,
-):
-	return polymatrix.expression.impl.impl.TransposeExprImpl(
-		underlying=underlying,
-)
diff --git a/polymatrix/expression/init/inittruncateexpr.py b/polymatrix/expression/init/inittruncateexpr.py
deleted file mode 100644
index d36db03..0000000
--- a/polymatrix/expression/init/inittruncateexpr.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_truncate_expr(
-		underlying: ExpressionBaseMixin,
-		variables: ExpressionBaseMixin,
-		degrees: tuple[int],
-		inverse: bool = None,
-):
-	if isinstance(degrees, int):
-		degrees = (degrees,)
-
-	if inverse is None:
-		inverse = False
-
-	return polymatrix.expression.impl.impl.TruncateExprImpl(
-		underlying=underlying,
-		variables=variables,
-		degrees=degrees,
-		inverse=inverse,
-	)
diff --git a/polymatrix/expression/init/initvstackexpr.py b/polymatrix/expression/init/initvstackexpr.py
deleted file mode 100644
index 3d47979..0000000
--- a/polymatrix/expression/init/initvstackexpr.py
+++ /dev/null
@@ -1,20 +0,0 @@
-from polymatrix.expression.init.initfromsympyexpr import init_from_sympy_expr
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
-
-def init_v_stack_expr(
-		underlying: tuple,
-):
-
-	def gen_underlying():
-	
-		for e in underlying:
-			if isinstance(e, ExpressionBaseMixin):
-				yield e
-			else:
-				yield init_from_sympy_expr(e)
-
-	return polymatrix.expression.impl.impl.VStackExprImpl(
-		underlying=tuple(gen_underlying()),
-)
diff --git a/polymatrix/expression/initexpressionbase.py b/polymatrix/expression/initexpressionbase.py
new file mode 100644
index 0000000..1a3e52c
--- /dev/null
+++ b/polymatrix/expression/initexpressionbase.py
@@ -0,0 +1,604 @@
+import typing
+import numpy as np
+import sympy
+
+import polymatrix.expression.implexpressionbase
+
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.getstacklines import get_stack_lines
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+from polymatrix.expression.utils.formatsubstitutions import format_substitutions
+from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
+from polymatrix.expression.implexpressionbase import FromTupleExprImpl, AdditionExprImpl
+
+
+def init_addition_expr(
+		left: ExpressionBaseMixin,
+		right: ExpressionBaseMixin,
+		stack: tuple[FrameSummary],
+):
+	return AdditionExprImpl(
+		left=left,
+		right=right,
+		stack=stack,
+)
+
+
+def init_block_diag_expr(
+		underlying: tuple,
+):
+	return polymatrix.expression.implexpressionbase.BlockDiagExprImpl(
+		underlying=underlying,
+	)
+
+
+def init_cache_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.CacheExprImpl(
+		underlying=underlying,
+)
+
+
+def init_combinations_expr(
+		expression: ExpressionBaseMixin,
+		degrees: tuple[int, ...] | int,
+):
+	if isinstance(degrees, int):
+		degrees = (degrees,)
+
+	return polymatrix.expression.implexpressionbase.CombinationsExprImpl(
+		expression=expression,
+		degrees=degrees,
+)
+
+
+# def init_derivative_expr(
+# 		underlying: ExpressionBaseMixin,
+# 		variables: ExpressionBaseMixin,
+# 		stack: tuple[FrameSummary],
+# 		introduce_derivatives: bool = None,
+# ):
+
+# 	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
+
+# 	if introduce_derivatives is None:
+# 		introduce_derivatives = False
+
+# 	return polymatrix.expression.implexpressionbase.DerivativeExprImpl(
+# 		underlying=underlying,
+# 		variables=variables,
+# 		introduce_derivatives=introduce_derivatives,
+# 		stack=stack,
+# 	)
+
+
+def init_determinant_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.DeterminantExprImpl(
+		underlying=underlying,
+)
+
+
+def init_diag_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.DiagExprImpl(
+		underlying=underlying,
+)
+
+
+def init_divergence_expr(
+		underlying: ExpressionBaseMixin,
+		variables: tuple,
+):
+	return polymatrix.expression.implexpressionbase.DivergenceExprImpl(
+		underlying=underlying,
+		variables=variables,
+)
+
+
+def init_division_expr(
+		left: ExpressionBaseMixin,
+		right: ExpressionBaseMixin,
+		stack: tuple[FrameSummary],
+):
+	return polymatrix.expression.implexpressionbase.DivisionExprImpl(
+		left=left,
+		right=right,
+		stack=stack,
+)
+
+
+def init_elem_mult_expr(
+		left: ExpressionBaseMixin,
+		right: ExpressionBaseMixin,
+		stack: tuple[FrameSummary],
+):
+	return polymatrix.expression.implexpressionbase.ElemMultExprImpl(
+		left=left,
+		right=right,
+)
+
+
+def init_eval_expr(
+		underlying: ExpressionBaseMixin,
+		variables: typing.Union[typing.Any, tuple, dict],
+		values: typing.Union[float, tuple] = None,
+):
+
+	substitutions = format_substitutions(
+		variables=variables, 
+		values=values,
+	)
+
+	def formatted_values(value):
+		if isinstance(value, np.ndarray):
+			return tuple(value.reshape(-1))
+
+		elif isinstance(value, tuple):
+			return value
+
+		elif isinstance(value, int) or isinstance(value, float):
+			return (value,)
+
+		else:
+			return (float(value),)
+
+	substitutions = tuple((variable, formatted_values(value)) for variable, value in substitutions)
+
+	return polymatrix.expression.implexpressionbase.EvalExprImpl(
+		underlying=underlying,
+		substitutions=substitutions,
+	)
+
+
+def init_eye_expr(
+		variable: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.EyeExprImpl(
+		variable=variable,
+)
+
+
+def init_filter_expr(
+		underlying: ExpressionBaseMixin,
+		predicator: ExpressionBaseMixin,
+		inverse: bool = None,
+):
+	if inverse is None:
+		inverse = False
+
+	return polymatrix.expression.implexpressionbase.FilterExprImpl(
+		underlying=underlying,
+		predicator=predicator,
+		inverse=inverse,
+)
+
+
+def init_from_symmetric_matrix_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.FromSymmetricMatrixExprImpl(
+		underlying=underlying,
+)
+
+
+DATA_TYPE = str | np.ndarray | sympy.Matrix | sympy.Expr | tuple | ExpressionBaseMixin
+
+def init_from_expr_or_none(
+	data: DATA_TYPE,
+) -> ExpressionBaseMixin | None:
+	if isinstance(data, str):
+		return init_parametrize_expr(
+			underlying=init_from_expr_or_none(1),
+			name=data,
+		)
+
+	elif isinstance(data, np.ndarray):
+		assert len(data.shape) <= 2
+
+		def gen_elements():
+			for row in data:
+				if isinstance(row, np.ndarray):
+					yield tuple(row)
+				else:
+					yield (row,)
+
+		data = tuple(gen_elements())
+
+	elif isinstance(data, sympy.Matrix):
+			data = tuple(tuple(i for i in data.row(row)) for row in range(data.rows))
+
+	elif isinstance(data, sympy.Expr):
+			data = ((sympy.expand(data),),)
+
+	elif isinstance(data, tuple):
+
+		if isinstance(data[0], tuple):
+			n_col = len(data[0])
+			assert all(len(col) == n_col for col in data)
+
+		else:
+			data = tuple((e,) for e in data)
+
+	elif isinstance(data, ExpressionBaseMixin):
+		return data
+
+	else:
+		if not isinstance(data, (float, int, np.number)):
+			return None
+
+		data = ((data,),)
+
+	return FromTupleExprImpl(
+			data=data,
+			stack=get_stack_lines(),
+	)
+
+
+def init_from_expr(
+	data: DATA_TYPE
+):
+	expr = init_from_expr_or_none(data)
+
+	if expr is None:
+		raise Exception(f'{data=}')
+
+	return expr
+
+
+def init_from_terms_expr(
+	terms: typing.Union[tuple, PolyMatrixMixin],
+	shape: tuple[int, int] = None,
+):
+
+	if isinstance(terms, PolyMatrixMixin):
+		shape = terms.shape
+		gen_terms = terms.gen_terms()
+
+	else:
+		assert shape is not None
+
+		if isinstance(terms, tuple):
+			gen_terms = terms
+
+		elif isinstance(terms, dict):
+			gen_terms = terms.items()
+
+		else:
+			raise Exception(f'{terms=}')
+
+	# Expression needs to be hashable
+	terms_formatted = tuple((key, tuple(monomials.items())) for key, monomials in gen_terms)
+
+	return polymatrix.expression.implexpressionbase.FromTermsExprImpl(
+		terms=terms_formatted,
+		shape=shape,
+	)
+
+
+def init_get_item_expr(
+		underlying: ExpressionBaseMixin,
+		index: tuple[tuple[int, ...], tuple[int, ...]],
+):
+
+	def get_hashable_slice(index):
+		if isinstance(index, slice):
+			return polymatrix.expression.implexpressionbase.GetItemExprImpl.Slice(start=index.start, stop=index.stop, step=index.step)
+		else:
+			return index
+
+	proper_index = (get_hashable_slice(index[0]), get_hashable_slice(index[1]))
+
+	return polymatrix.expression.implexpressionbase.GetItemExprImpl(
+		underlying=underlying,
+		index=proper_index,
+)
+
+
+def init_half_newton_polytope_expr(
+		monomials: ExpressionBaseMixin,
+		variables: ExpressionBaseMixin,
+		filter: ExpressionBaseMixin | None = None,
+):
+	return polymatrix.expression.implexpressionbase.HalfNewtonPolytopeExprImpl(
+		monomials=monomials,
+		variables=variables,
+		filter=filter
+)
+
+
+# def init_linear_in_expr(
+# 		underlying: ExpressionBaseMixin,
+# 		monomials: ExpressionBaseMixin,
+# 		variables: ExpressionBaseMixin,
+# 		ignore_unmatched: bool = None,
+# ):
+# 	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
+
+# 	return polymatrix.expression.implexpressionbase.LinearInExprImpl(
+# 		underlying=underlying,
+# 		monomials=monomials,
+# 		variables=variables,
+# 		ignore_unmatched = ignore_unmatched,
+# )
+
+
+def init_linear_matrix_in_expr(
+		underlying: ExpressionBaseMixin,
+		variable: int,
+):
+	return polymatrix.expression.implexpressionbase.LinearMatrixInExprImpl(
+		underlying=underlying,
+		variable=variable,
+)
+
+
+# def init_linear_monomials_expr(
+# 		underlying: ExpressionBaseMixin,
+# 		variables: ExpressionBaseMixin,
+# ):
+
+# 	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
+
+# 	return polymatrix.expression.implexpressionbase.LinearMonomialsExprImpl(
+# 		underlying=underlying,
+# 		variables=variables,
+# 	)
+
+
+def init_matrix_mult_expr(
+		left: ExpressionBaseMixin,
+		right: ExpressionBaseMixin,
+		stack: tuple[FrameSummary],
+):
+	return polymatrix.expression.implexpressionbase.MatrixMultExprImpl(
+		left=left,
+		right=right,
+		stack=stack,
+)
+
+
+def init_max_degree_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.MaxDegreeExprImpl(
+		underlying=underlying,
+)
+
+
+def init_max_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.MaxExprImpl(
+		underlying=underlying,
+)
+
+
+def init_parametrize_expr(
+		underlying: ExpressionBaseMixin,
+		name: str = None,
+):
+	if name is None:
+		name = 'undefined'
+
+	return polymatrix.expression.implexpressionbase.ParametrizeExprImpl(
+		underlying=underlying,
+		name=name,
+)
+
+
+def init_parametrize_matrix_expr(
+		underlying: ExpressionBaseMixin,
+		name: str,
+):
+	return polymatrix.expression.implexpressionbase.ParametrizeMatrixExprImpl(
+		underlying=underlying,
+		name=name,
+	)
+
+
+def init_quadratic_in_expr(
+		underlying: ExpressionBaseMixin,
+		monomials: ExpressionBaseMixin,
+		variables: ExpressionBaseMixin,
+		stack: tuple[FrameSummary],
+):
+
+	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
+
+	return polymatrix.expression.implexpressionbase.QuadraticInExprImpl(
+		underlying=underlying,
+		monomials=monomials,
+		variables=variables,
+		stack=stack,
+	)
+
+
+def init_quadratic_monomials_expr(
+		underlying: ExpressionBaseMixin,
+		variables: ExpressionBaseMixin,
+):
+
+	assert isinstance(variables, ExpressionBaseMixin), f'{variables=}'
+
+	return polymatrix.expression.implexpressionbase.QuadraticMonomialsExprImpl(
+		underlying=underlying,
+		variables=variables,
+	)
+
+
+def init_rep_mat_expr(
+		underlying: ExpressionBaseMixin,
+		repetition: tuple,
+):
+	return polymatrix.expression.implexpressionbase.RepMatExprImpl(
+		underlying=underlying,
+		repetition=repetition,
+)
+
+
+def init_reshape_expr(
+		underlying: ExpressionBaseMixin,
+		new_shape: tuple,
+):
+	return polymatrix.expression.implexpressionbase.ReshapeExprImpl(
+		underlying=underlying,
+		new_shape=new_shape,
+)
+
+
+def init_set_element_at_expr(
+		underlying: ExpressionBaseMixin,
+		index: tuple,
+		value: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.SetElementAtExprImpl(
+		underlying=underlying,
+		index=index,
+		value=value,
+)
+
+
+def init_squeeze_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.SqueezeExprImpl(
+		underlying=underlying,
+)
+
+
+def init_substitute_expr(
+		underlying: ExpressionBaseMixin,
+		variables: tuple,
+		values: tuple = None,
+):
+
+	substitutions = format_substitutions(
+		variables=variables, 
+		values=values,
+	)
+
+	def formatted_values(value) -> ExpressionBaseMixin:
+		if isinstance(value, ExpressionBaseMixin):
+			expr = value
+		
+		else:
+			expr = init_from_expr(value)
+
+		return polymatrix.expression.implexpressionbase.ReshapeExprImpl(
+			underlying=expr,
+			new_shape=(-1, 1),
+		)
+
+	substitutions = tuple((variable, formatted_values(value)) for variable, value in substitutions)
+
+	return polymatrix.expression.implexpressionbase.SubstituteExprImpl(
+		underlying=underlying,
+		substitutions=substitutions,
+	)
+
+
+def init_subtract_monomials_expr(
+		underlying: ExpressionBaseMixin,
+		monomials: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.SubtractMonomialsExprImpl(
+		underlying=underlying,
+		monomials=monomials,
+)
+
+
+def init_sum_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.SumExprImpl(
+		underlying=underlying,
+)
+
+
+def init_symmetric_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.SymmetricExprImpl(
+		underlying=underlying,
+)
+
+
+def init_to_constant_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.ToConstantExprImpl(
+		underlying=underlying,
+)
+
+
+def init_to_quadratic_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.ToQuadraticExprImpl(
+		underlying=underlying,
+)
+
+
+def init_to_sorted_variables(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.ToSortedVariablesImpl(
+		underlying=underlying,
+)
+
+
+def init_to_symmetric_matrix_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.ToSymmetricMatrixExprImpl(
+		underlying=underlying,
+)
+
+
+def init_transpose_expr(
+		underlying: ExpressionBaseMixin,
+):
+	return polymatrix.expression.implexpressionbase.TransposeExprImpl(
+		underlying=underlying,
+)
+
+
+def init_truncate_expr(
+		underlying: ExpressionBaseMixin,
+		variables: ExpressionBaseMixin,
+		degrees: tuple[int],
+		inverse: bool = None,
+):
+	if isinstance(degrees, int):
+		degrees = (degrees,)
+
+	if inverse is None:
+		inverse = False
+
+	return polymatrix.expression.implexpressionbase.TruncateExprImpl(
+		underlying=underlying,
+		variables=variables,
+		degrees=degrees,
+		inverse=inverse,
+	)
+
+
+# def init_v_stack_expr(
+# 		underlying: tuple,
+# ):
+
+# 	def gen_underlying():
+	
+# 		for e in underlying:
+# 			if isinstance(e, ExpressionBaseMixin):
+# 				yield e
+# 			else:
+# 				yield init_from_(e)
+
+# 	return polymatrix.expression.implexpressionbase.VStackExprImpl(
+# 		underlying=tuple(gen_underlying()),
+# )
diff --git a/polymatrix/expression/mixins/additionexprmixin.py b/polymatrix/expression/mixins/additionexprmixin.py
index e247ebd..6257d55 100644
--- a/polymatrix/expression/mixins/additionexprmixin.py
+++ b/polymatrix/expression/mixins/additionexprmixin.py
@@ -3,11 +3,13 @@ import math
 import typing
 import dataclassabc
 
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.tooperatorexception import to_operator_exception
 from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 from polymatrix.polymatrix.polymatrix import PolyMatrix
 from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 
 
 class AdditionExprMixin(ExpressionBaseMixin):
@@ -21,6 +23,11 @@ class AdditionExprMixin(ExpressionBaseMixin):
     def right(self) -> ExpressionBaseMixin:
         ...
 
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
     # overwrites abstract method of `ExpressionBaseMixin`
     def apply(
         self, 
@@ -54,7 +61,11 @@ class AdditionExprMixin(ExpressionBaseMixin):
                 all_underlying = (left, broadcasted_right)
 
         else:
-            assert left.shape == right.shape, f'{left.shape} != {right.shape}'
+            if not (left.shape == right.shape):
+                raise AssertionError(to_operator_exception(
+                    message=f'{left.shape} != {right.shape}',
+                    stack=self.stack,
+                ))
 
             all_underlying = (left, right)
 
diff --git a/polymatrix/expression/mixins/assertdegreeexprmixin.py b/polymatrix/expression/mixins/assertdegreeexprmixin.py
deleted file mode 100644
index e57de1c..0000000
--- a/polymatrix/expression/mixins/assertdegreeexprmixin.py
+++ /dev/null
@@ -1,49 +0,0 @@
-
-import abc
-import itertools
-import dataclassabc
-from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
-
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-from polymatrix.polymatrix.polymatrix import PolyMatrix
-from polymatrix.expressionstate.expressionstate import ExpressionState
-
-
-class AssertDegreeExprMixin(ExpressionBaseMixin):
-    @property
-    @abc.abstractmethod
-    def underlying(self) -> tuple[ExpressionBaseMixin, ...]:
-        ...
-
-    @property
-    @abc.abstractmethod
-    def degrees(self) -> tuple[int, ...]:
-        ...
-
-    # overwrites abstract method of `ExpressionBaseMixin`
-    def apply(
-        self, 
-        state: ExpressionState,
-    ) -> tuple[ExpressionState, PolyMatrix]:
-    
-        state, underlying = self.underlying.apply(state)
-
-        @dataclassabc.dataclassabc(frozen=True)
-        class AssertDegreePolyMatrix(PolyMatrixMixin):
-            underlying: PolyMatrixMixin
-            degrees: tuple[int, ...]
-
-            def get_poly(self, row: int, col: int) -> dict[tuple[int, ...], float]:
-                terms = self.underlying.get_poly(row, col)
-
-                for monomial in terms.keys():
-                    pass
-
-                return terms
-
-        polymatrix = AssertDegreePolyMatrix(
-            underlying=underlying,
-            degrees=self.degrees,
-        )
-
-        return state, polymatrix
diff --git a/polymatrix/expression/mixins/blockdiagexprmixin.py b/polymatrix/expression/mixins/blockdiagexprmixin.py
index b210d27..0bafd0d 100644
--- a/polymatrix/expression/mixins/blockdiagexprmixin.py
+++ b/polymatrix/expression/mixins/blockdiagexprmixin.py
@@ -26,8 +26,6 @@ class BlockDiagExprMixin(ExpressionBaseMixin):
             state, polymat = expr.apply(state=state)
             all_underlying.append(polymat)
 
-        # assert all(underlying.shape[0] == underlying.shape[1] for underlying in all_underlying)
-
         @dataclassabc.dataclassabc(frozen=True)
         class BlockDiagPolyMatrix(PolyMatrixMixin):
             all_underlying: tuple[PolyMatrixMixin]
diff --git a/polymatrix/expression/mixins/cacheexprmixin.py b/polymatrix/expression/mixins/cacheexprmixin.py
index fe6dbb7..7a92259 100644
--- a/polymatrix/expression/mixins/cacheexprmixin.py
+++ b/polymatrix/expression/mixins/cacheexprmixin.py
@@ -10,6 +10,8 @@ from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 
 
 class CacheExprMixin(ExpressionBaseMixin):
+    """ Caches the polynomial matrix using the state """
+    
     @property
     @abc.abstractclassmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/combinationsexprmixin.py b/polymatrix/expression/mixins/combinationsexprmixin.py
index 4841e8e..a9cc247 100644
--- a/polymatrix/expression/mixins/combinationsexprmixin.py
+++ b/polymatrix/expression/mixins/combinationsexprmixin.py
@@ -1,26 +1,29 @@
 
 import abc
 import itertools
-import math
 
 from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.polymatrix.polymatrix import PolyMatrix
 from polymatrix.expressionstate.expressionstate import ExpressionState
-from polymatrix.expression.utils.getmonomialindices import get_monomial_indices
-from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices
-from polymatrix.expression.utils.sortmonomials import sort_monomials
+from polymatrix.expression.utils.multiplymonomials import multiply_monomials
 
 
 class CombinationsExprMixin(ExpressionBaseMixin):
+    """
+    combination using degrees=(0, 1, 2, 3):
+
+    [[x]]  ->  [[1], [x], [x**2], [x**3]]
+    """
+     
     @property
     @abc.abstractmethod
-    def monomials(self) -> ExpressionBaseMixin:
+    def expression(self) -> ExpressionBaseMixin:
         ...
 
     @property
     @abc.abstractmethod
-    def number(self) -> int:
+    def degrees(self) -> tuple[int, ...]:
         ...
 
     # overwrites abstract method of `ExpressionBaseMixin`
@@ -28,34 +31,82 @@ class CombinationsExprMixin(ExpressionBaseMixin):
         self, 
         state: ExpressionState,
     ) -> tuple[ExpressionState, PolyMatrix]:
-        if self.number == 0:
-            terms = {(0, 0): {tuple(): 1.0}}
+        # if self.degree == 0:
+        #     terms = {(0, 0): {tuple(): 1.0}}
 
-            poly_matrix = init_poly_matrix(
-                terms=terms,
-                shape=(1, 1),
-            )
+        #     poly_matrix = init_poly_matrix(
+        #         terms=terms,
+        #         shape=(1, 1),
+        #     )
+
+        # elif self.degree == 1:
+        #     state, monomials = self.expression.apply(state=state) 
+        #     poly_matrix = monomials
+
+        # else:
+
+        state, poly_matrix = self.expression.apply(state=state)
+
+        assert poly_matrix.shape[1] == 1
+
+        def gen_indices():
+            for degree in self.degrees:
+                yield from itertools.combinations_with_replacement(range(poly_matrix.shape[0]), degree)
+
+        indices = tuple(gen_indices())
+
+        terms = {}
 
-        elif self.number == 1:
-            state, monomials = self.monomials.apply(state=state) 
-            poly_matrix = monomials
+        for row, indexing in enumerate(indices):
 
-        else:
+            # print(indexing)
 
-            state, monomials = get_monomial_indices(state, self.monomials)
+            if indexing is tuple():
+                terms[row, 0] = {tuple(): 1.0}
+                continue
 
-            combinations = tuple(itertools.combinations_with_replacement(monomials, self.number))
+            def acc_product(acc, v):
+                left_monomials = acc
+                row = v
 
-            terms = {}
+                right_monomials = poly_matrix.get_poly(row, 0).keys()
+                # print(right_monomials)
 
-            for row, combination in enumerate(combinations):
-                combination_monomial = merge_monomial_indices(combination)
+                if left_monomials is (None,):
+                    return right_monomials
 
-                terms[row, 0] = {combination_monomial: 1.0}
+                monomials = tuple(multiply_monomials(left_monomials, right_monomials))
+                return monomials
 
-            poly_matrix = init_poly_matrix(
-                terms=terms,
-                shape=(math.comb(len(monomials) + self.number - 1, self.number), 1),
+            *_, monomials = itertools.accumulate(
+                indexing,
+                acc_product,
+                initial=(None,),
             )
 
+            terms[row, 0] = {m: 1.0 for m in monomials} 
+
+        poly_matrix = init_poly_matrix(
+            terms=terms,
+            shape=(len(terms), 1),
+        )
+
+        # indices = filter(lambda v: sum(v) <= self.degree, itertools.product(*(range(self.degree) for _ in range(dim))))
+
+        # state, monomials = get_monomial_indices(state, self.monomials)
+
+        # combinations = tuple(itertools.combinations_with_replacement(monomials, self.number))
+
+        # terms = {}
+
+        # for row, combination in enumerate(combinations):
+        #     combination_monomial = merge_monomial_indices(combination)
+
+        #     terms[row, 0] = {combination_monomial: 1.0}
+
+        # poly_matrix = init_poly_matrix(
+        #     terms=terms,
+        #     shape=(math.comb(len(monomials) + self.number - 1, self.number), 1),
+        # )
+
         return state, poly_matrix   
diff --git a/polymatrix/expression/mixins/derivativeexprmixin.py b/polymatrix/expression/mixins/derivativeexprmixin.py
index 183608b..15873b6 100644
--- a/polymatrix/expression/mixins/derivativeexprmixin.py
+++ b/polymatrix/expression/mixins/derivativeexprmixin.py
@@ -8,9 +8,19 @@ from polymatrix.polymatrix.polymatrix import PolyMatrix
 from polymatrix.expressionstate.expressionstate import ExpressionState
 from polymatrix.expression.utils.getderivativemonomials import get_derivative_monomials
 from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.tooperatorexception import to_operator_exception
 
 
 class DerivativeExprMixin(ExpressionBaseMixin):
+    """
+    differentiate w.r.t. x:
+
+    [[x**2]]  ->  [[2*x]]
+
+    introduce_derivatives: not used at the moment
+    """
+
     @property
     @abc.abstractmethod
     def underlying(self) -> ExpressionBaseMixin:
@@ -26,6 +36,11 @@ class DerivativeExprMixin(ExpressionBaseMixin):
     def introduce_derivatives(self) -> bool:
         ...
 
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
     # overwrites abstract method of `ExpressionBaseMixin`
     def apply(
         self, 
@@ -35,7 +50,11 @@ class DerivativeExprMixin(ExpressionBaseMixin):
         state, underlying = self.underlying.apply(state=state)        
         state, diff_wrt_variables = get_variable_indices_from_variable(state, self.variables)
 
-        assert underlying.shape[1] == 1, f'{underlying.shape=}'
+        if not (underlying.shape[1] == 1):
+            raise AssertionError(to_operator_exception(
+                message=f'{underlying.shape[1]=} is not 1',
+                stack=self.stack,
+            ))
 
         terms = {}
 
@@ -52,7 +71,7 @@ class DerivativeExprMixin(ExpressionBaseMixin):
                     monomial_terms=underlying_terms,
                     diff_wrt_variable=diff_wrt_variable,
                     state=state,
-                    considered_variables=set(),
+                    considered_variables=set(diff_wrt_variables),
                     introduce_derivatives=self.introduce_derivatives,
                 )
 
diff --git a/polymatrix/expression/mixins/diagexprmixin.py b/polymatrix/expression/mixins/diagexprmixin.py
index 26c75d0..6651e13 100644
--- a/polymatrix/expression/mixins/diagexprmixin.py
+++ b/polymatrix/expression/mixins/diagexprmixin.py
@@ -6,6 +6,14 @@ from polymatrix.expressionstate.mixins.expressionstatemixin import ExpressionSta
 from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 
 class DiagExprMixin(ExpressionBaseMixin):
+    """
+    [[1],[2]]  ->  [[1,0],[0,2]]
+
+    or
+
+    [[1,0],[0,2]]  ->  [[1],[2]]
+    """
+    
     @property
     @abc.abstractclassmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/divisionexprmixin.py b/polymatrix/expression/mixins/divisionexprmixin.py
index 94bfec1..e236456 100644
--- a/polymatrix/expression/mixins/divisionexprmixin.py
+++ b/polymatrix/expression/mixins/divisionexprmixin.py
@@ -2,6 +2,9 @@
 import abc
 import dataclasses
 
+from polymatrix.expression.mixins.elemmultexprmixin import ElemMultExprMixin
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.tooperatorexception import to_operator_exception
 from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.polymatrix.polymatrix import PolyMatrix
@@ -19,20 +22,40 @@ class DivisionExprMixin(ExpressionBaseMixin):
     def right(self) -> ExpressionBaseMixin:
         ...
 
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
     # overwrites abstract method of `ExpressionBaseMixin`
     def apply(
         self, 
         state: ExpressionState,
     ) -> tuple[ExpressionState, PolyMatrix]:
 
-        # add an auxillary equation and, therefore, needs to be cached
-        if self in state.cache:
-            return state, state.cache[self]
-
         state, left = self.left.apply(state=state)
         state, right = self.right.apply(state=state)
 
-        assert right.shape == (1, 1)
+        # if left.shape == (1, 1):
+        #     left, right = right, left
+
+        # assert right.shape == (1, 1)
+
+        if not (right.shape == (1, 1)):
+            raise AssertionError(to_operator_exception(
+                message=f'{right.shape=} is not (1, 1)',
+                stack=self.stack,
+            ))
+
+        right_poly = right.get_poly(0, 0)
+
+        if len(right_poly) == 1 and tuple() in right_poly:
+            right_inv = {(0, 0): {tuple(): 1/right_poly[tuple()]}}
+            return ElemMultExprMixin.elem_mult(state, left, right)
+
+        # add an auxillary equation and, therefore, needs to be cached
+        if self in state.cache:
+            return state, state.cache[self]
 
         terms = {}
 
@@ -53,7 +76,7 @@ class DivisionExprMixin(ExpressionBaseMixin):
                 terms[row, col] = dict(gen_monomial_terms())
 
         def gen_auxillary_terms():
-            for monomial, value in right.get_poly(0, 0).items():
+            for monomial, value in right_poly.items():
                 yield monomial + ((division_variable, 1),), value
 
         auxillary_terms = dict(gen_auxillary_terms())
diff --git a/polymatrix/expression/mixins/elemmultexprmixin.py b/polymatrix/expression/mixins/elemmultexprmixin.py
index 2edf9e6..37c6f6f 100644
--- a/polymatrix/expression/mixins/elemmultexprmixin.py
+++ b/polymatrix/expression/mixins/elemmultexprmixin.py
@@ -28,14 +28,12 @@ class ElemMultExprMixin(ExpressionBaseMixin):
     # def shape(self) -> tuple[int, int]:
     #     return self.left.shape
 
-    # overwrites abstract method of `ExpressionBaseMixin`
-    def apply(
-        self, 
+    @staticmethod
+    def elem_mult(
         state: ExpressionState,
-    ) -> tuple[ExpressionState, PolyMatrix]:
-        state, left = self.left.apply(state=state)
-        state, right = self.right.apply(state=state)
-
+        left: PolyMatrix,
+        right: PolyMatrix,
+    ):
         if left.shape != right.shape and left.shape == (1, 1):
             left, right = right, left
 
@@ -95,4 +93,76 @@ class ElemMultExprMixin(ExpressionBaseMixin):
             shape=left.shape,
         )
 
-        return state, poly_matrix        
+        return state, poly_matrix      
+        
+
+    # overwrites abstract method of `ExpressionBaseMixin`
+    def apply(
+        self, 
+        state: ExpressionState,
+    ) -> tuple[ExpressionState, PolyMatrix]:
+        state, left = self.left.apply(state=state)
+        state, right = self.right.apply(state=state)
+
+        return self.elem_mult(state, left, right)
+
+        # if left.shape != right.shape and left.shape == (1, 1):
+        #     left, right = right, left
+
+        # if right.shape == (1, 1):
+        #     right_poly = right.get_poly(0, 0)
+
+        #     @dataclassabc.dataclassabc(frozen=True)
+        #     class BroadCastedPolyMatrix(PolyMatrixMixin):
+        #         underlying: tuple[tuple[int], float]
+        #         shape: tuple[int, int]
+
+        #         def get_poly(self, row: int, col: int) -> typing.Optional[dict[tuple[int, ...], float]]:
+        #             return self.underlying
+
+        #     right = BroadCastedPolyMatrix(
+        #         underlying=right_poly,
+        #         shape=left.shape,
+        #     )
+
+        # terms = {}
+
+        # for poly_row in range(left.shape[0]):
+        #     for poly_col in range(left.shape[1]):
+
+        #         terms_row_col = {}
+
+        #         left_terms = left.get_poly(poly_row, poly_col)
+        #         if left_terms is None:
+        #             continue
+
+        #         right_terms = right.get_poly(poly_row, poly_col)
+        #         if right_terms is None:
+        #             continue
+
+        #         for (left_monomial, left_value), (right_monomial, right_value) \
+        #                 in itertools.product(left_terms.items(), right_terms.items()):
+
+        #             value = left_value * right_value
+
+        #             # if value == 0:
+        #             #     continue
+
+        #             # monomial = tuple(sorted(left_monomial + right_monomial))
+
+        #             new_monomial = merge_monomial_indices((left_monomial, right_monomial))
+
+        #             if new_monomial not in terms_row_col:
+        #                 terms_row_col[new_monomial] = 0
+
+        #             terms_row_col[new_monomial] += value
+
+        #         if 0 < len(terms_row_col):
+        #             terms[poly_row, poly_col] = terms_row_col
+
+        # poly_matrix = init_poly_matrix(
+        #     terms=terms,
+        #     shape=left.shape,
+        # )
+
+        # return state, poly_matrix        
diff --git a/polymatrix/expression/mixins/expressionbasemixin.py b/polymatrix/expression/mixins/expressionbasemixin.py
index f2f8507..8aedc25 100644
--- a/polymatrix/expression/mixins/expressionbasemixin.py
+++ b/polymatrix/expression/mixins/expressionbasemixin.py
@@ -1,6 +1,7 @@
 import abc
-from polymatrix.expressionstate.mixins.expressionstatemixin import ExpressionStateMixin
-from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
+
+from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.polymatrix.polymatrix import PolyMatrix
 
 
 class ExpressionBaseMixin(
@@ -8,5 +9,5 @@ class ExpressionBaseMixin(
 ):
 
     @abc.abstractmethod
-    def apply(self, state: ExpressionStateMixin) -> tuple[ExpressionStateMixin, PolyMatrixMixin]:
+    def apply(self, state: ExpressionState) -> tuple[ExpressionState, PolyMatrix]:
         ...
diff --git a/polymatrix/expression/mixins/fromsymmetricmatrixexprmixin.py b/polymatrix/expression/mixins/fromsymmetricmatrixexprmixin.py
index 906adf8..395eda9 100644
--- a/polymatrix/expression/mixins/fromsymmetricmatrixexprmixin.py
+++ b/polymatrix/expression/mixins/fromsymmetricmatrixexprmixin.py
@@ -11,6 +11,8 @@ from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 
 
 class FromSymmetricMatrixExprMixin(ExpressionBaseMixin):
+    """ inverse of ToSymmetricMatrixExprMixin """
+
     @property
     @abc.abstractclassmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/fromsympyexprmixin.py b/polymatrix/expression/mixins/fromtupleexprmixin.py
index 25edde0..4157ce7 100644
--- a/polymatrix/expression/mixins/fromsympyexprmixin.py
+++ b/polymatrix/expression/mixins/fromtupleexprmixin.py
@@ -4,16 +4,25 @@ import math
 import sympy
 import numpy as np
 
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.tooperatorexception import to_operator_exception
 from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-from polymatrix.expressionstate.mixins.expressionstatemixin import ExpressionStateMixin
 from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
+from polymatrix.expressionstate.mixins.expressionstatemixin import ExpressionStateMixin
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+
 
+class FromTupleExprMixin(ExpressionBaseMixin):
+    DATA_TYPE = int | float | np.number | sympy.Expr | ExpressionBaseMixin
 
-class FromSympyExprMixin(ExpressionBaseMixin):
     @property
     @abc.abstractmethod
-    def data(self) -> tuple[tuple[float]]:
+    def data(self) -> tuple[tuple['FromTupleExprMixin.DATA_TYPE']]:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
         ...
 
     # overwrites abstract method of `ExpressionBaseMixin`
@@ -68,18 +77,24 @@ class FromSympyExprMixin(ExpressionBaseMixin):
 
                         polynomial[monomial] = value
 
-                # elif isinstance(poly_data, np.number):
-                #     terms_row_col = {tuple(): float(poly_data)}
-
                 elif isinstance(poly_data, ExpressionBaseMixin):
                     state, instance = poly_data.apply(state)
 
-                    assert instance.shape == (1, 1)
+                    # assert instance.shape == (1, 1)
+                    if not (instance.shape == (1, 1)):
+                        raise AssertionError(to_operator_exception(
+                            message=f'{instance.shape=} is not (1, 1)',
+                            stack=self.stack,
+                        ))
 
                     polynomial = instance.get_poly(0, 0)
 
                 else:
-                    raise Exception(f'{poly_data=}, {type(poly_data)=}')
+                    # raise Exception(f'{poly_data=}, {type(poly_data)=}')
+                    raise AssertionError(to_operator_exception(
+                        message=f'unknown data type {type(poly_data)=}',
+                        stack=self.stack,
+                    ))
 
                 polynomials[poly_row, poly_col] = polynomial
 
diff --git a/polymatrix/expression/mixins/legendreseriesmixin.py b/polymatrix/expression/mixins/legendreseriesmixin.py
new file mode 100644
index 0000000..3232207
--- /dev/null
+++ b/polymatrix/expression/mixins/legendreseriesmixin.py
@@ -0,0 +1,62 @@
+import abc
+
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.polymatrix.polymatrix import PolyMatrix
+from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
+from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+
+
+class LegendreSeriesMixin(ExpressionBaseMixin):
+    @property
+    @abc.abstractmethod
+    def underlying(self) -> tuple[ExpressionBaseMixin]:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def degrees(self) -> tuple[int, ...] | None:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
+    # overwrites abstract method of `ExpressionBaseMixin`
+    def apply(
+        self, 
+        state: ExpressionState,
+    ) -> tuple[ExpressionState, PolyMatrix]:
+        
+        state, underlying = self.underlying.apply(state)
+
+        if self.degrees is None:
+            degrees = range(underlying.shape[0])
+        else:
+            degrees = self.degrees
+        
+        terms = {}
+
+        for degree in degrees:
+        # for degree in self.degree:
+            poly = underlying.get_poly(degree, 0)
+
+            terms[degree, 0] = dict(poly)
+
+            if 2 <= degree:
+                poly = underlying.get_poly(degree - 2, 0)
+                factor = - (degree - 1) / (degree + 1)
+
+                for m, v in poly.items():
+                    if m in terms[degree, 0]:
+                        terms[degree, 0][m] += v*factor
+                    else:
+                        terms[degree, 0][m] = v*factor
+
+        poly_matrix = init_poly_matrix(
+            terms=terms,
+            shape=(len(terms), 1),
+        )
+
+        return state, poly_matrix
diff --git a/polymatrix/expression/mixins/matrixmultexprmixin.py b/polymatrix/expression/mixins/matrixmultexprmixin.py
index a02b4f1..98527e7 100644
--- a/polymatrix/expression/mixins/matrixmultexprmixin.py
+++ b/polymatrix/expression/mixins/matrixmultexprmixin.py
@@ -1,14 +1,13 @@
 
 import abc
-import itertools
-import math
 
+from polymatrix.utils.getstacklines import FrameSummary
 from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.polymatrix.polymatrix import PolyMatrix
 from polymatrix.expressionstate.expressionstate import ExpressionState
-from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices
 from polymatrix.expression.utils.multiplypolynomial import multiply_polynomial
+from polymatrix.utils.tooperatorexception import to_operator_exception
 
 
 class MatrixMultExprMixin(ExpressionBaseMixin):
@@ -22,6 +21,11 @@ class MatrixMultExprMixin(ExpressionBaseMixin):
     def right(self) -> ExpressionBaseMixin:
         ...
 
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
     # overwrites abstract method of `ExpressionBaseMixin`
     def apply(
         self, 
@@ -30,7 +34,13 @@ class MatrixMultExprMixin(ExpressionBaseMixin):
         state, left = self.left.apply(state=state)
         state, right = self.right.apply(state=state)
 
-        assert left.shape[1] == right.shape[0], f'{left.shape[1]} is not equal to {right.shape[0]}'
+        # assert left.shape[1] == right.shape[0], f'{left.shape[1]} is not equal to {right.shape[0]}'
+
+        if not (left.shape[1] == right.shape[0]):
+            raise AssertionError(to_operator_exception(
+                message=f'{left.shape[1]} is not equal to {right.shape[0]}',
+                stack=self.stack,
+            ))
 
         terms = {}
 
diff --git a/polymatrix/expression/mixins/parametrizematrixexprmixin.py b/polymatrix/expression/mixins/parametrizematrixexprmixin.py
index f59ddb7..de7be88 100644
--- a/polymatrix/expression/mixins/parametrizematrixexprmixin.py
+++ b/polymatrix/expression/mixins/parametrizematrixexprmixin.py
@@ -8,7 +8,6 @@ from polymatrix.expressionstate.mixins.expressionstatemixin import ExpressionSta
 from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 
 
-# rename ParametrizeSymmetricMatrixExprMixin
 class ParametrizeMatrixExprMixin(ExpressionBaseMixin):
     @property
     @abc.abstractclassmethod
diff --git a/polymatrix/expression/mixins/productexprmixin.py b/polymatrix/expression/mixins/productexprmixin.py
new file mode 100644
index 0000000..bacd29e
--- /dev/null
+++ b/polymatrix/expression/mixins/productexprmixin.py
@@ -0,0 +1,132 @@
+
+import abc
+import itertools
+
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.tooperatorexception import to_operator_exception
+from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
+from polymatrix.polymatrix.polymatrix import PolyMatrix
+from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.expression.utils.multiplymonomials import multiply_monomials
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+
+
+class ProductExprMixin(ExpressionBaseMixin):
+    @property
+    @abc.abstractmethod
+    def underlying(self) -> tuple[ExpressionBaseMixin]:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def degrees(self) -> tuple[int, ...] | None:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
+    # overwrites abstract method of `ExpressionBaseMixin`
+    def apply(
+        self, 
+        state: ExpressionState,
+    ) -> tuple[ExpressionState, PolyMatrix]:
+        
+        # if self.number == 0:
+        #     terms = {(0, 0): {tuple(): 1.0}}
+
+        #     poly_matrix = init_poly_matrix(
+        #         terms=terms,
+        #         shape=(1, 1),
+        #     )
+
+        # elif self.number == 1:
+        #     state, monomials = self.monomials.apply(state=state) 
+        #     poly_matrix = monomials
+
+        # else:
+
+        if len(self.underlying) == 0:
+            terms = {(0,0): {tuple(): 1}}
+        
+        else:
+
+            def acc_underlying(acc, v):
+                state, acc_polymatrix = acc
+                underlying = v
+
+                state, poly_matrix = underlying.apply(state=state)
+
+                return state, acc_polymatrix + (poly_matrix,)
+
+            *_, (state, underlying) = itertools.accumulate(
+                self.underlying,
+                acc_underlying,
+                initial=(state, tuple())
+            )
+
+            # highest_degrees = tuple(e.shape[0] for e in underlying)
+
+            # if self.degrees is None:
+            #     degrees = range(sum(highest_degrees))
+
+            # else:
+            #     degrees = self.degrees
+
+            # max_degree = max(degrees)
+
+            # for poly_matrix in underlying:
+            #     if not (max_degree <= poly_matrix.shape[0]):
+            #         raise AssertionError(to_operator_exception(
+            #             message=f'{poly_matrix.shape[0]} < {max_degree}',
+            #             stack=self.stack,
+            #         ))
+
+            def gen_indices():
+                product_indices = itertools.product(*(range(e.shape[0]) for e in underlying))
+
+                if self.degrees is None:
+                    yield from product_indices
+
+                else:
+                    yield from filter(lambda v: sum(v) in self.degrees, product_indices)
+                    
+            indices = tuple(gen_indices())
+
+            # indices = filter(lambda v: sum(v) <= self.degree, itertools.product(*(range(self.degree) for _ in range(dim))))
+
+            terms = {}
+
+            for row, indexing in enumerate(indices):
+
+                # print(indexing)
+
+                def acc_product(acc, v):
+                    left_monomials = acc
+                    polymatrix, row = v
+
+                    right_monomials = polymatrix.get_poly(row, 0).keys()
+
+                    # print(f'{left_monomials=}')
+                    # print(f'{right_monomials=}')
+
+                    if left_monomials is (None,):
+                        return right_monomials
+
+                    return tuple(multiply_monomials(left_monomials, right_monomials))
+
+                *_, monomials = itertools.accumulate(
+                    zip(underlying, indexing),
+                    acc_product,
+                    initial=(None,),
+                )
+
+                terms[row, 0] = {m: 1.0 for m in monomials} 
+
+        poly_matrix = init_poly_matrix(
+            terms=terms,
+            shape=(len(terms), 1),
+        )
+
+        return state, poly_matrix   
diff --git a/polymatrix/expression/mixins/quadraticinexprmixin.py b/polymatrix/expression/mixins/quadraticinexprmixin.py
index 2b23827..052124f 100644
--- a/polymatrix/expression/mixins/quadraticinexprmixin.py
+++ b/polymatrix/expression/mixins/quadraticinexprmixin.py
@@ -9,6 +9,8 @@ from polymatrix.expressionstate.expressionstate import ExpressionState
 from polymatrix.expression.utils.getmonomialindices import get_monomial_indices
 from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
 from polymatrix.expression.utils.splitmonomialindices import split_monomial_indices
+from polymatrix.utils.getstacklines import FrameSummary
+from polymatrix.utils.tooperatorexception import to_operator_exception
 
 
 class QuadraticInExprMixin(ExpressionBaseMixin):
@@ -27,6 +29,11 @@ class QuadraticInExprMixin(ExpressionBaseMixin):
     def variables(self) -> ExpressionBaseMixin:
         ...
 
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
     # overwrites abstract method of `ExpressionBaseMixin`
     def apply(
         self, 
@@ -37,14 +44,23 @@ class QuadraticInExprMixin(ExpressionBaseMixin):
         state, sos_monomials = get_monomial_indices(state, self.monomials)
         state, variable_indices = get_variable_indices_from_variable(state, self.variables)
 
-        assert underlying.shape == (1, 1), f'underlying shape is {underlying.shape}'
+        # assert underlying.shape == (1, 1), f'underlying shape is {underlying.shape}'
+        if not (underlying.shape == (1, 1)):
+            raise AssertionError(to_operator_exception(
+                message=f'underlying shape is {underlying.shape}',
+                stack=self.stack,
+            ))
 
         terms = collections.defaultdict(lambda: collections.defaultdict(float))
 
         underlying_poly = underlying.get_poly(0, 0)
 
         if underlying_poly is None:
-            raise Exception(f'{self.underlying} is empty')
+            raise AssertionError(to_operator_exception(
+                message=f'{underlying} is empty',
+                # message=f'{self.underlying} is empty',
+                stack=self.stack,
+            ))
 
         for monomial, value in underlying_poly.items():
 
@@ -56,12 +72,20 @@ class QuadraticInExprMixin(ExpressionBaseMixin):
             try:
                 col = sos_monomials.index(left)
             except ValueError:
-                raise ValueError(f'{left=} not in {sos_monomials=}')
+                raise AssertionError(to_operator_exception(
+                    message=f'{left=} not in {sos_monomials=}',
+                    stack=self.stack,
+                ))
+                # raise ValueError(f'{left=} not in {sos_monomials=}')
 
             try:
                 row = sos_monomials.index(right)
             except ValueError:
-                raise ValueError(f'{right=} not in {sos_monomials=}')
+                raise AssertionError(to_operator_exception(
+                    message=f'{right=} not in {sos_monomials=}',
+                    stack=self.stack,
+                ))
+                # raise ValueError(f'{right=} not in {sos_monomials=}')
 
             terms[row, col][p_monomial] += value
 
diff --git a/polymatrix/expression/mixins/substituteexprmixin.py b/polymatrix/expression/mixins/substituteexprmixin.py
index e693de1..432c9c6 100644
--- a/polymatrix/expression/mixins/substituteexprmixin.py
+++ b/polymatrix/expression/mixins/substituteexprmixin.py
@@ -1,10 +1,8 @@
-
 import abc
 import collections
 import itertools
 import math
 import typing
-from polymatrix.expression.impl.reshapeexprmixin import ReshapeExprImpl
 
 from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
@@ -42,13 +40,14 @@ class SubstituteExprMixin(ExpressionBaseMixin):
             if indices is None:
                 return acc
 
-            state, substitution = ReshapeExprImpl(
-                underlying=expr,
-                new_shape=(-1, 1),
-            ).apply(state)
+            state, substitution = expr.apply(state)
+
+            # state, substitution = ReshapeExprImpl(
+            #     underlying=expr,
+            #     new_shape=(-1, 1),
+            # ).apply(state)
 
             # # todo: improve this?
-            # substitution = substitution.reshape(-1, 1)
             # assert substitution.shape[1] == 1, f'The substitutions have to be a vector {substitution=}'
 
             def gen_polynomials():
diff --git a/polymatrix/expression/mixins/sumexprmixin.py b/polymatrix/expression/mixins/sumexprmixin.py
index 9957f10..9e12e07 100644
--- a/polymatrix/expression/mixins/sumexprmixin.py
+++ b/polymatrix/expression/mixins/sumexprmixin.py
@@ -10,6 +10,12 @@ from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 
 
 class SumExprMixin(ExpressionBaseMixin):
+    """ 
+    For each row of the matrix sum the colum elements. 
+    
+    [[1, 2, 3], [4, 5, 6]]  ->  [[6], [15]]
+    """
+
     @property
     @abc.abstractclassmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/symmetricexprmixin.py b/polymatrix/expression/mixins/symmetricexprmixin.py
index 2677501..2b191db 100644
--- a/polymatrix/expression/mixins/symmetricexprmixin.py
+++ b/polymatrix/expression/mixins/symmetricexprmixin.py
@@ -11,6 +11,12 @@ from polymatrix.expressionstate.expressionstate import ExpressionState
 
 
 class SymmetricExprMixin(ExpressionBaseMixin):
+    """
+    Maps a square matrix to a symmetric square matrix by taking the average of the diagonal elements.
+
+    [[1, 2], [3, 4]]  ->  [[1, 2.5], [2.5, 4]]
+    """
+    
     @property
     @abc.abstractmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/toconstantexprmixin.py b/polymatrix/expression/mixins/toconstantexprmixin.py
index ca928c6..91fbb2b 100644
--- a/polymatrix/expression/mixins/toconstantexprmixin.py
+++ b/polymatrix/expression/mixins/toconstantexprmixin.py
@@ -9,6 +9,11 @@ from polymatrix.expressionstate.expressionstate import ExpressionState
 
 
 class ToConstantExprMixin(ExpressionBaseMixin):
+    """
+    Discards non-constant polynomial coefficients
+
+    [[a + b*x + c*x^2]]  ->  [[a]]
+    """
     @property
     @abc.abstractmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/tosortedvariablesmixin.py b/polymatrix/expression/mixins/tosortedvariablesmixin.py
index eb80597..d0a5d79 100644
--- a/polymatrix/expression/mixins/tosortedvariablesmixin.py
+++ b/polymatrix/expression/mixins/tosortedvariablesmixin.py
@@ -8,6 +8,7 @@ from polymatrix.polymatrix.polymatrix import PolyMatrix
 from polymatrix.expressionstate.expressionstate import ExpressionState
 
 
+# to be deleted?
 class ToSortedVariablesExprMixin(ExpressionBaseMixin):
     @property
     @abc.abstractmethod
diff --git a/polymatrix/expression/mixins/tosymmetricmatrixexprmixin.py b/polymatrix/expression/mixins/tosymmetricmatrixexprmixin.py
index 0d8c6ed..a34c2b3 100644
--- a/polymatrix/expression/mixins/tosymmetricmatrixexprmixin.py
+++ b/polymatrix/expression/mixins/tosymmetricmatrixexprmixin.py
@@ -9,6 +9,12 @@ from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin
 
 
 class ToSymmetricMatrixExprMixin(ExpressionBaseMixin):
+    """
+    Convert a vector to a symmetric matrix
+
+    [1, 2, 3, 4, 5, 6].T  ->  [[1, 2, 3], [2, 4, 5], [3, 5, 6]]
+    """
+
     @property
     @abc.abstractclassmethod
     def underlying(self) -> ExpressionBaseMixin:
@@ -24,8 +30,6 @@ class ToSymmetricMatrixExprMixin(ExpressionBaseMixin):
 
         assert underlying.shape[1] == 1
 
-        # state, variable_indices = get_variable_indices_from_variable(state, self.underlying)
-
         def invert_binomial_coefficient(val):
             idx = 1
             sum_val = 1
@@ -38,7 +42,6 @@ class ToSymmetricMatrixExprMixin(ExpressionBaseMixin):
 
             return idx
 
-        # n_row = invert_binomial_coefficient(len(variable_indices))
         n_row = invert_binomial_coefficient(underlying.shape[0])
 
         terms = {}
@@ -46,7 +49,6 @@ class ToSymmetricMatrixExprMixin(ExpressionBaseMixin):
 
         for row in range(n_row):
             for col in range(row, n_row):
-                # terms[row, col] = {((variable_indices[var_index], 1),): 1.0}
                 terms[row, col] = underlying.get_poly(var_index, 0)
 
                 if row != col:
diff --git a/polymatrix/expression/mixins/transposeexprmixin.py b/polymatrix/expression/mixins/transposeexprmixin.py
index 5aa59e2..7c60daa 100644
--- a/polymatrix/expression/mixins/transposeexprmixin.py
+++ b/polymatrix/expression/mixins/transposeexprmixin.py
@@ -12,6 +12,11 @@ from polymatrix.expressionstate.expressionstate import ExpressionState
 
 
 class TransposeExprMixin(ExpressionBaseMixin):
+    """
+    Transpose the polynomial matrix
+
+    [[1, 2, 3]]  ->  [[1], [2], [3]]
+    """
     @property
     @abc.abstractmethod
     def underlying(self) -> ExpressionBaseMixin:
diff --git a/polymatrix/expression/mixins/vstackexprmixin.py b/polymatrix/expression/mixins/vstackexprmixin.py
index fd68ca3..e787495 100644
--- a/polymatrix/expression/mixins/vstackexprmixin.py
+++ b/polymatrix/expression/mixins/vstackexprmixin.py
@@ -10,6 +10,12 @@ from polymatrix.expressionstate.expressionstate import ExpressionState
 
 
 class VStackExprMixin(ExpressionBaseMixin):
+    """
+    Vertical stacking of the underlying polynomial matrices
+
+    [[1, 2]], [[3, 4]]  ->  [[1, 2], [3, 4]]
+    """
+    
     @property
     @abc.abstractmethod
     def underlying(self) -> tuple[ExpressionBaseMixin, ...]:
diff --git a/polymatrix/expression/op.py b/polymatrix/expression/op.py
new file mode 100644
index 0000000..0d98401
--- /dev/null
+++ b/polymatrix/expression/op.py
@@ -0,0 +1,67 @@
+import polymatrix.expression.initexpressionbase
+
+from polymatrix.utils.getstacklines import get_stack_lines
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+
+
+def diff(
+    expression: ExpressionBaseMixin, 
+    variables: ExpressionBaseMixin,
+    introduce_derivatives: bool = None,
+) -> ExpressionBaseMixin:
+            
+    if not isinstance(variables, ExpressionBaseMixin):
+        variables=polymatrix.expression.initexpressionbase.init_from_expr(variables)
+
+    if introduce_derivatives is None:
+        introduce_derivatives = False
+
+    return polymatrix.expression.implexpressionbase.DerivativeExprImpl(
+		underlying=expression,
+		variables=variables,
+		introduce_derivatives=introduce_derivatives,
+		stack=get_stack_lines(),
+	)
+
+
+def legendre(
+    expression: ExpressionBaseMixin,
+    degrees: tuple[int, ...] = None,
+) -> ExpressionBaseMixin:
+    
+    return polymatrix.expression.implexpressionbase.LegendreSeriesImpl(
+        underlying=expression,
+        degrees=degrees,
+        stack=get_stack_lines(),
+    )
+    
+
+def linear_in(
+    expression: ExpressionBaseMixin, 
+    variables: ExpressionBaseMixin, 
+    monomials: ExpressionBaseMixin = None,
+    ignore_unmatched: bool = None,    
+) -> ExpressionBaseMixin:
+    
+    if monomials is None:
+        monomials = linear_monomials(
+            expression=expression,
+            variables=variables,
+        )
+
+    return polymatrix.expression.implexpressionbase.LinearInExprImpl(
+        underlying=expression,
+        monomials=monomials,
+        variables=variables,
+        ignore_unmatched = ignore_unmatched,
+    )
+
+def linear_monomials(
+    expression: ExpressionBaseMixin, 
+    variables: ExpressionBaseMixin,
+) -> ExpressionBaseMixin:
+    
+    return polymatrix.expression.implexpressionbase.LinearMonomialsExprImpl(
+		underlying=expression,
+		variables=variables,
+	)
diff --git a/polymatrix/expression/to.py b/polymatrix/expression/to.py
new file mode 100644
index 0000000..53cd8a8
--- /dev/null
+++ b/polymatrix/expression/to.py
@@ -0,0 +1,129 @@
+import sympy
+import numpy as np
+
+from polymatrix.expression.expression import Expression
+from polymatrix.expression.mixins.parametrizeexprmixin import ParametrizeExprMixin
+from polymatrix.expression.mixins.parametrizematrixexprmixin import ParametrizeMatrixExprMixin
+from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
+from polymatrix.expressionstate.expressionstate import ExpressionState
+from polymatrix.statemonad.init.initstatemonad import init_state_monad
+from polymatrix.statemonad.mixins.statemonadmixin import StateMonadMixin
+
+
+def shape(
+    expr: Expression,
+) -> StateMonadMixin[ExpressionState, tuple[int, ...]]:
+    def func(state: ExpressionState):
+        state, polymatrix = expr.apply(state)
+
+        return state, polymatrix.shape
+
+    return init_state_monad(func)
+
+
+def to_constant_repr(
+    expr: Expression,
+    assert_constant: bool = True,
+) -> StateMonadMixin[ExpressionState, np.ndarray]:
+
+    def func(state: ExpressionState):
+        state, underlying = expr.apply(state)
+
+        A = np.zeros(underlying.shape, dtype=np.double)
+
+        for (row, col), polynomial in underlying.gen_terms():
+            for monomial, value in polynomial.items():
+                if len(monomial) == 0:
+                    A[row, col] = value
+                
+                elif assert_constant:
+                    raise Exception(f'non-constant term {monomial=}')
+
+        return state, A
+
+    return init_state_monad(func)
+
+
+def to_degrees(
+    expr: Expression,
+    variables: Expression,
+) -> StateMonadMixin[ExpressionState, np.ndarray]:
+
+    def func(state: ExpressionState):
+        state, underlying = expr.apply(state)
+        state, variable_indices = get_variable_indices_from_variable(state, variables)
+
+        def gen_rows():
+            for row in range(underlying.shape[0]):
+                def gen_cols():
+                    for col in range(underlying.shape[1]):
+                        
+                        def gen_degrees():
+                            polynomial = underlying.get_poly(row, col)
+
+                            if polynomial is None:
+                                yield 0
+
+                            else:
+                                for monomial, _ in polynomial.items():
+                                    yield sum(count for var, count in monomial if var in variable_indices)
+
+                        yield tuple(set(gen_degrees()))
+
+                yield tuple(gen_cols())
+
+        return state, tuple(gen_rows())
+
+    return init_state_monad(func)
+
+
+def to_sympy_repr(
+    expr: Expression,
+) -> StateMonadMixin[ExpressionState, sympy.Expr]:
+
+    def func(state: ExpressionState):
+        state, underlying = expr.apply(state)
+
+        A = np.zeros(underlying.shape, dtype=object)
+
+        for (row, col), polynomial in underlying.gen_terms():
+
+            sympy_polynomial = 0
+
+            for monomial, value in polynomial.items():
+                sympy_monomial = 1
+
+                for offset, count in monomial:
+
+                    variable = state.get_key_from_offset(offset)
+                    # def get_variable_from_offset(offset: int):
+                    #     for variable, (start, end) in state.offset_dict.items():
+                    #         if start <= offset < end:
+                                # assert end - start == 1, f'{start=}, {end=}, {variable=}'
+
+                    if isinstance(variable, sympy.core.symbol.Symbol):
+                        variable_name = variable.name
+                    elif isinstance(variable, (ParametrizeExprMixin, ParametrizeMatrixExprMixin)):
+                        variable_name = variable.name
+                    elif isinstance(variable, str):
+                        variable_name = variable
+                    else:
+                        raise Exception(f'{variable=}')
+
+                    start, end = state.offset_dict[variable]
+
+                    if end - start == 1:
+                        sympy_var = sympy.Symbol(variable_name)
+                    else:
+                        sympy_var = sympy.Symbol(f'{variable_name}_{offset - start + 1}')
+
+                    # var = get_variable_from_offset(offset)
+                    sympy_monomial *= sympy_var**count
+
+                sympy_polynomial += value * sympy_monomial
+
+            A[row, col] = sympy_polynomial
+
+        return state, A
+
+    return init_state_monad(func)
diff --git a/polymatrix/expression/init/initsubstituteexpr.py b/polymatrix/expression/utils/formatsubstitutions.py
index e141d9e..c9b2972 100644
--- a/polymatrix/expression/init/initsubstituteexpr.py
+++ b/polymatrix/expression/utils/formatsubstitutions.py
@@ -1,10 +1,4 @@
 import typing
-import numpy as np
-
-from polymatrix.expression.init.initfromsympyexpr import init_from_sympy_expr
-from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
-import polymatrix.expression.impl.impl
-
 
 def format_substitutions(
 	variables: typing.Union[typing.Any, tuple, dict],
@@ -48,28 +42,3 @@ def format_substitutions(
 
 	return substitutions
 
-
-def init_substitute_expr(
-		underlying: ExpressionBaseMixin,
-		variables: tuple,
-		values: tuple = None,
-):
-
-	substitutions = format_substitutions(
-		variables=variables, 
-		values=values,
-	)
-
-	def formatted_values(value) -> ExpressionBaseMixin:
-		if isinstance(value, ExpressionBaseMixin):
-			return value
-		
-		else:
-			return init_from_sympy_expr(value)
-
-	substitutions = tuple((variable, formatted_values(value)) for variable, value in substitutions)
-
-	return polymatrix.expression.impl.impl.SubstituteExprImpl(
-		underlying=underlying,
-		substitutions=substitutions,
-	)
diff --git a/polymatrix/expression/utils/getderivativemonomials.py b/polymatrix/expression/utils/getderivativemonomials.py
index 4aa7144..9bc1e1e 100644
--- a/polymatrix/expression/utils/getderivativemonomials.py
+++ b/polymatrix/expression/utils/getderivativemonomials.py
@@ -1,4 +1,6 @@
 import collections
+import dataclasses
+import itertools
 
 from polymatrix.expressionstate.expressionstate import ExpressionState
 
@@ -10,63 +12,66 @@ def get_derivative_monomials(
     considered_variables: set,
     introduce_derivatives: bool,
 ):
+    
+    @dataclasses.dataclass(frozen=True)
+    class DerivativeKey:
+        variable: int
+        with_respect_to: int
 
-    # if self.introduce_derivatives:
-
-    #     raise Exception('not implemented')
+    if introduce_derivatives:
         
-    #     def gen_new_variables():
-    #         for monomial in monomial_terms.keys():
-    #             for var in monomial:
-    #                 if var not in diff_wrt_variables and var not in considered_variables:
-    #                     yield var
+        def gen_new_variables():
+            for monomial in monomial_terms.keys():
+                for var in monomial:
+                    if var is not diff_wrt_variable and var not in considered_variables:
+                        yield var
 
-    #     new_variables = set(gen_new_variables())
+        new_variables = set(gen_new_variables())
 
-    #     new_considered_variables = considered_variables | new_variables
+        new_considered_variables = considered_variables | new_variables
 
-    #     def acc_state_candidates(acc, new_variable):
-    #         state, candidates = acc
+        def acc_state_candidates(acc, new_variable):
+            state, candidates = acc
 
-    #         @dataclasses.dataclass
-    #         class DerivativeKey:
-    #             variable: int
-    #             with_respect_to: int
+            # key representing the derivative of the variable
+            key = DerivativeKey(
+                variable=new_variable,
+                with_respect_to=diff_wrt_variable,
+            )
+            # register if it doesn't already exist
+            state = state.register(key=key, n_param=1)
 
-    #         key = DerivativeKey(
-    #             variable=new_variable,
-    #             with_respect_to=diff_wrt_variable,
-    #         )
-    #         state = state.register(key=key, n_param=1)
+            # for each new variable we expect an auxillary equation
+            state, auxillary_derivation_terms = get_derivative_monomials(
+                monomial_terms=state.auxillary_equations[new_variable],
+                diff_wrt_variable=diff_wrt_variable,
+                state=state,
+                considered_variables=new_considered_variables,
+                introduce_derivatives=True,
+            )
 
-    #         state, auxillary_derivation_terms = get_derivative_terms(
-    #             monomial_terms=state.auxillary_equations[new_variable],
-    #             diff_wrt_variable=diff_wrt_variable,
-    #             state=state,
-    #             considered_variables=new_considered_variables,
-    #         )
+            if 1 < len(auxillary_derivation_terms):
+                derivation_variable = state.offset_dict[key][0]
 
-    #         if 1 < len(auxillary_derivation_terms):
-    #             derivation_variable = state.offset_dict[key][0]
+                state = dataclasses.replace(
+                    state,
+                    auxillary_equations=state.auxillary_equations | {derivation_variable: auxillary_derivation_terms},
+                )
 
-    #             state = dataclasses.replace(
-    #                 state,
-    #                 auxillary_equations=state.auxillary_equations | {derivation_variable: auxillary_derivation_terms},
-    #             )
+                return state, candidates + (new_variable,)
 
-    #             return state, candidates + (new_variable,)
+            else:
+                return state, candidates
 
-    #         else:
-    #             return state, candidates
+        *_, (state, confirmed_variables) = itertools.accumulate(
+            new_variables,
+            acc_state_candidates,
+            initial=(state, tuple()),
+        )
 
-    #     *_, (state, confirmed_variables) = itertools.accumulate(
-    #         new_variables,
-    #         acc_state_candidates,
-    #         initial=(state, tuple()),
-    #     )
+    else:
+        confirmed_variables = tuple()
 
-    # else:
-    #     confirmed_variables = tuple()
 
     derivation_terms = collections.defaultdict(float)
 
@@ -100,20 +105,21 @@ def get_derivative_monomials(
             diff_monomial, value = differentiate_monomial(diff_wrt_variable)
             derivation_terms[diff_monomial] += value
 
-        # # only used if introduce_derivatives == True
-        # for candidate_variable in monomial_cnt.keys():
-        #     if candidate_variable in considered_variables or candidate_variable in confirmed_variables:
-        #         key = init_derivative_key(
-        #             variable=candidate_variable,
-        #             with_respect_to=diff_wrt_variable,
-        #         )
-        #         derivation_variable = state.offset_dict[key][0]
-
-        #         diff_monomial, value = differentiate_monomial(
-        #             dependent_variable=candidate_variable, 
-        #             derivation_variable=derivation_variable,
-        #         )
-        #         derivation_terms[diff_monomial] += value
+        # only used if introduce_derivatives == True
+        if introduce_derivatives:
+            for candidate_variable in monomial_cnt.keys():
+                if candidate_variable in considered_variables or candidate_variable in confirmed_variables:
+                    key = DerivativeKey(
+                        variable=candidate_variable,
+                        with_respect_to=diff_wrt_variable,
+                    )
+                    derivation_variable = state.offset_dict[key][0]
+
+                    diff_monomial, value = differentiate_monomial(
+                        dependent_variable=candidate_variable, 
+                        derivation_variable=derivation_variable,
+                    )
+                    derivation_terms[diff_monomial] += value
 
     # return state, dict(derivation_terms)
     return state, derivation_terms
\ No newline at end of file
diff --git a/polymatrix/expression/utils/getmonomialindices.py b/polymatrix/expression/utils/getmonomialindices.py
index 6d6c599..ac6d5f7 100644
--- a/polymatrix/expression/utils/getmonomialindices.py
+++ b/polymatrix/expression/utils/getmonomialindices.py
@@ -1,7 +1,8 @@
+from polymatrix.expressionstate.expressionstate import ExpressionState
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 
 
-def get_monomial_indices(state, monomials) -> tuple:
+def get_monomial_indices(state: ExpressionState, monomials: ExpressionBaseMixin) -> tuple[ExpressionState, tuple[int, ...]]:
 
     state, monomials_obj = monomials.apply(state)
 
diff --git a/polymatrix/expression/utils/multiplymonomials.py b/polymatrix/expression/utils/multiplymonomials.py
new file mode 100644
index 0000000..a4b3273
--- /dev/null
+++ b/polymatrix/expression/utils/multiplymonomials.py
@@ -0,0 +1,15 @@
+import itertools
+import math
+
+from polymatrix.expression.utils.mergemonomialindices import merge_monomial_indices
+
+
+def multiply_monomials(left_monomials, right_monomials):
+    def gen_monomials():
+        for left_monomial, right_monomial in itertools.product(left_monomials, right_monomials):
+            
+            monomial = merge_monomial_indices((left_monomial, right_monomial))
+
+            yield monomial
+
+    return gen_monomials()
diff --git a/polymatrix/utils/__init__.py b/polymatrix/utils/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/polymatrix/utils/__init__.py
diff --git a/polymatrix/utils/getstacklines.py b/polymatrix/utils/getstacklines.py
new file mode 100644
index 0000000..6fdd69f
--- /dev/null
+++ b/polymatrix/utils/getstacklines.py
@@ -0,0 +1,23 @@
+import dataclasses
+import traceback
+
+
+@dataclasses.dataclass(frozen=True)
+class FrameSummary:
+    filename: str
+    lineno: int
+    name: str
+    line: str
+
+
+def get_stack_lines(index: int = 2) -> tuple[FrameSummary]:
+    def gen_stack_lines():
+        for obj in traceback.extract_stack()[:-index]:
+            yield FrameSummary(
+                filename=obj.filename,
+                lineno=obj.lineno,
+                name=obj.name,
+                line=obj.line,
+            )
+
+    return tuple(gen_stack_lines())
diff --git a/polymatrix/utils/tooperatorexception.py b/polymatrix/utils/tooperatorexception.py
new file mode 100644
index 0000000..5f7db0f
--- /dev/null
+++ b/polymatrix/utils/tooperatorexception.py
@@ -0,0 +1,14 @@
+from polymatrix.utils.getstacklines import FrameSummary
+
+
+def to_operator_exception(
+        message: str,
+        stack: tuple[FrameSummary],
+) -> str:
+    exception_lines = [
+        message,
+        f'  Assertion traceback (most recent call last):',
+        *(f'    File "{stack_line.filename}", line {stack_line.lineno}\n      {stack_line.line}' for stack_line in stack),
+    ]
+
+    return '\n'.join(exception_lines)
diff --git a/test_polymatrix/test_expression/test_addition.py b/test_polymatrix/test_expression/test_addition.py
index 6c8a9a2..e3e7f05 100644
--- a/test_polymatrix/test_expression/test_addition.py
+++ b/test_polymatrix/test_expression/test_addition.py
@@ -1,8 +1,7 @@
 import unittest
 
-from polymatrix.expression.init.initadditionexpr import init_addition_expr
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestAddition(unittest.TestCase):
@@ -28,17 +27,17 @@ class TestAddition(unittest.TestCase):
             },
         }
 
-        left = init_from_terms_expr(
+        left = polymatrix.expression.initexpressionbase.init_from_terms_expr(
             terms=left_terms,
             shape=(2, 2),
         )
 
-        right = init_from_terms_expr(
+        right = polymatrix.expression.initexpressionbase.init_from_terms_expr(
             terms=right_terms,
             shape=(2, 2),
         )
 
-        expr = init_addition_expr(
+        expr = polymatrix.expression.initexpressionbase.init_addition_expr(
             left=left,
             right=right,
         )
diff --git a/test_polymatrix/test_expression/test_blockdiag.py b/test_polymatrix/test_expression/test_blockdiag.py
index a66935e..4a67c77 100644
--- a/test_polymatrix/test_expression/test_blockdiag.py
+++ b/test_polymatrix/test_expression/test_blockdiag.py
@@ -1,8 +1,7 @@
 import unittest
 
-from polymatrix.expression.init.initblockdiagexpr import init_block_diag_expr
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestBlockDiag(unittest.TestCase):
@@ -26,10 +25,10 @@ class TestBlockDiag(unittest.TestCase):
             },
         }
 
-        expr = init_block_diag_expr(
+        expr = polymatrix.expression.initexpressionbase.init_block_diag_expr(
             underlying=(
-                init_from_terms_expr(terms=terms1, shape=(2, 2),), 
-                init_from_terms_expr(terms=terms2, shape=(2, 2),),
+                polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms1, shape=(2, 2),), 
+                polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms2, shape=(2, 2),),
             ),
         )
 
diff --git a/test_polymatrix/test_expression/test_derivative.py b/test_polymatrix/test_expression/test_derivative.py
index 27dfbeb..e8f551b 100644
--- a/test_polymatrix/test_expression/test_derivative.py
+++ b/test_polymatrix/test_expression/test_derivative.py
@@ -1,10 +1,7 @@
 import unittest
-from polymatrix.expression.init.initderivativeexpr import init_derivative_expr
-from polymatrix.expression.init.initdivergenceexpr import init_divergence_expr
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestDerivative(unittest.TestCase):
@@ -27,9 +24,9 @@ class TestDerivative(unittest.TestCase):
             (2, 0): {((2, 1),): 1},
         }
 
-        expr = init_derivative_expr(
-            underlying=init_from_terms_expr(terms=underlying_terms, shape=(2, 1)),
-            variables=init_from_terms_expr(terms=variable_terms, shape=(3, 1),),
+        expr = polymatrix.expression.initexpressionbase.init_derivative_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=underlying_terms, shape=(2, 1)),
+            variables=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(3, 1),),
         )
 
         state = init_expression_state(n_param=3)
diff --git a/test_polymatrix/test_expression/test_divergence.py b/test_polymatrix/test_expression/test_divergence.py
index 387a30f..2d2b6ad 100644
--- a/test_polymatrix/test_expression/test_divergence.py
+++ b/test_polymatrix/test_expression/test_divergence.py
@@ -1,9 +1,7 @@
 import unittest
-from polymatrix.expression.init.initdivergenceexpr import init_divergence_expr
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestDivergence(unittest.TestCase):
@@ -30,9 +28,9 @@ class TestDivergence(unittest.TestCase):
             (2, 0): {((2, 1),): 1},
         }
 
-        expr = init_divergence_expr(
-            underlying=init_from_terms_expr(terms=underlying_terms, shape=(3, 1)),
-            variables=init_from_terms_expr(terms=variable_terms, shape=(3, 1),),
+        expr = polymatrix.expression.initexpressionbase.init_divergence_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=underlying_terms, shape=(3, 1)),
+            variables=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(3, 1),),
         )
 
         state = init_expression_state(n_param=2)
diff --git a/test_polymatrix/test_expression/test_eval.py b/test_polymatrix/test_expression/test_eval.py
index 9b9454f..5787f2a 100644
--- a/test_polymatrix/test_expression/test_eval.py
+++ b/test_polymatrix/test_expression/test_eval.py
@@ -1,8 +1,7 @@
 import unittest
 
-from polymatrix.expression.init.initevalexpr import init_eval_expr
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestEval(unittest.TestCase):
@@ -19,8 +18,8 @@ class TestEval(unittest.TestCase):
             },
         }
 
-        expr = init_eval_expr(
-            underlying=init_from_terms_expr(terms=terms, shape=(2, 1)),
+        expr = polymatrix.expression.initexpressionbase.init_eval_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms, shape=(2, 1)),
             variables=(0, 1),
             values=(2.0, 3.0),
         )
diff --git a/test_polymatrix/test_expression/test_linearin.py b/test_polymatrix/test_expression/test_linearin.py
index 572c34c..d0fb0d9 100644
--- a/test_polymatrix/test_expression/test_linearin.py
+++ b/test_polymatrix/test_expression/test_linearin.py
@@ -1,9 +1,8 @@
 import unittest
-from polymatrix.expression.init.initfromsympyexpr import init_from_sympy_expr
+from polymatrix.expression.initexpressionbase import init_from_
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestLinearIn(unittest.TestCase):
@@ -36,10 +35,10 @@ class TestLinearIn(unittest.TestCase):
             (1, 0): {((1, 1),): 1},
         }
 
-        expr = init_linear_in_expr(
-            underlying=init_from_terms_expr(terms=underlying_terms, shape=(2, 1)),
-            monomials=init_from_terms_expr(terms=monomial_terms, shape=(4, 1),),
-            variables=init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
+        expr = polymatrix.expression.initexpressionbase.init_linear_in_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=underlying_terms, shape=(2, 1)),
+            monomials=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=monomial_terms, shape=(4, 1),),
+            variables=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
         )
 
         state = init_expression_state(n_param=2)
diff --git a/test_polymatrix/test_expression/test_matrixmult.py b/test_polymatrix/test_expression/test_matrixmult.py
index 12d0584..b15f7eb 100644
--- a/test_polymatrix/test_expression/test_matrixmult.py
+++ b/test_polymatrix/test_expression/test_matrixmult.py
@@ -1,9 +1,7 @@
 import unittest
-from polymatrix.expression.init.initadditionexpr import init_addition_expr
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
 
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initmatrixmultexpr import init_matrix_mult_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestMatrixMult(unittest.TestCase):
@@ -32,17 +30,17 @@ class TestMatrixMult(unittest.TestCase):
             },
         }
 
-        left = init_from_terms_expr(
+        left = polymatrix.expression.initexpressionbase.init_from_terms_expr(
             terms=left_terms,
             shape=(2, 2),
         )
 
-        right = init_from_terms_expr(
+        right = polymatrix.expression.initexpressionbase.init_from_terms_expr(
             terms=right_terms,
             shape=(2, 1),
         )
 
-        expr = init_matrix_mult_expr(
+        expr = polymatrix.expression.initexpressionbase.init_matrix_mult_expr(
             left=left,
             right=right,
         )
diff --git a/test_polymatrix/test_expression/test_quadraticin.py b/test_polymatrix/test_expression/test_quadraticin.py
index debf779..7f3767d 100644
--- a/test_polymatrix/test_expression/test_quadraticin.py
+++ b/test_polymatrix/test_expression/test_quadraticin.py
@@ -1,8 +1,7 @@
 import unittest
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initquadraticinexpr import init_quadratic_in_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestQuadraticIn(unittest.TestCase):
@@ -35,10 +34,10 @@ class TestQuadraticIn(unittest.TestCase):
             (1, 0): {((1, 1),): 1},
         }
 
-        expr = init_quadratic_in_expr(
-            underlying=init_from_terms_expr(terms=underlying_terms, shape=(1, 1)),
-            monomials=init_from_terms_expr(terms=monomial_terms, shape=(3, 1)),
-            variables=init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
+        expr = polymatrix.expression.initexpressionbase.init_quadratic_in_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=underlying_terms, shape=(1, 1)),
+            monomials=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=monomial_terms, shape=(3, 1)),
+            variables=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
         )
 
         state = init_expression_state(n_param=2)
diff --git a/test_polymatrix/test_expression/test_substitude.py b/test_polymatrix/test_expression/test_substitude.py
index 457d35e..4a1f224 100644
--- a/test_polymatrix/test_expression/test_substitude.py
+++ b/test_polymatrix/test_expression/test_substitude.py
@@ -1,8 +1,7 @@
 import unittest
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initsubstituteexpr import init_substitute_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestEval(unittest.TestCase):
@@ -28,10 +27,10 @@ class TestEval(unittest.TestCase):
             (0, 0): {((0, 1),): 1},
         }
 
-        expr = init_substitute_expr(
-            underlying=init_from_terms_expr(terms=terms, shape=(1, 1)),
-            variables=init_from_terms_expr(terms=variable_terms, shape=(1, 1),),
-            values=(init_from_terms_expr(terms=substitution, shape=(1, 1)),),
+        expr = polymatrix.expression.initexpressionbase.init_substitute_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms, shape=(1, 1)),
+            variables=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(1, 1),),
+            values=(polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=substitution, shape=(1, 1)),),
         )
 
         state = init_expression_state(n_param=2)
diff --git a/test_polymatrix/test_expression/test_subtractmonomials.py b/test_polymatrix/test_expression/test_subtractmonomials.py
index 8cd8585..77fcc84 100644
--- a/test_polymatrix/test_expression/test_subtractmonomials.py
+++ b/test_polymatrix/test_expression/test_subtractmonomials.py
@@ -1,11 +1,7 @@
 import unittest
-from polymatrix.expression.init.initderivativeexpr import init_derivative_expr
-from polymatrix.expression.init.initdivergenceexpr import init_divergence_expr
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr
-from polymatrix.expression.init.initsubtractmonomialsexpr import init_subtract_monomials_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestDerivative(unittest.TestCase):
@@ -29,9 +25,9 @@ class TestDerivative(unittest.TestCase):
             },
         }
 
-        expr = init_subtract_monomials_expr(
-            underlying=init_from_terms_expr(terms=monomials1, shape=(2, 1)),
-            monomials=init_from_terms_expr(terms=monomials2, shape=(2, 1)),
+        expr = polymatrix.expression.initexpressionbase.init_subtract_monomials_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=monomials1, shape=(2, 1)),
+            monomials=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=monomials2, shape=(2, 1)),
         )
 
         state = init_expression_state(n_param=3)
diff --git a/test_polymatrix/test_expression/test_sum.py b/test_polymatrix/test_expression/test_sum.py
index 053e428..27c4766 100644
--- a/test_polymatrix/test_expression/test_sum.py
+++ b/test_polymatrix/test_expression/test_sum.py
@@ -1,9 +1,7 @@
 import unittest
 
-from polymatrix.expression.init.initevalexpr import init_eval_expr
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initsumexpr import init_sum_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestSum(unittest.TestCase):
@@ -20,8 +18,8 @@ class TestSum(unittest.TestCase):
             },
         }
 
-        expr = init_sum_expr(
-            underlying=init_from_terms_expr(terms=terms, shape=(1, 2)),
+        expr = polymatrix.expression.initexpressionbase.init_sum_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms, shape=(1, 2)),
         )
 
         state = init_expression_state(n_param=2)
diff --git a/test_polymatrix/test_expression/test_symmetric.py b/test_polymatrix/test_expression/test_symmetric.py
index a93d975..5313ff1 100644
--- a/test_polymatrix/test_expression/test_symmetric.py
+++ b/test_polymatrix/test_expression/test_symmetric.py
@@ -1,9 +1,7 @@
 import unittest
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initquadraticinexpr import init_quadratic_in_expr
-from polymatrix.expression.init.initsymmetricexpr import init_symmetric_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestQuadraticIn(unittest.TestCase):
@@ -22,12 +20,12 @@ class TestQuadraticIn(unittest.TestCase):
             },
         }
 
-        underlying = init_from_terms_expr(
+        underlying = polymatrix.expression.initexpressionbase.init_from_terms_expr(
             terms=terms,
             shape=(2, 2),
         )
 
-        expr = init_symmetric_expr(
+        expr = polymatrix.expression.initexpressionbase.init_symmetric_expr(
             underlying=underlying,
         )
 
diff --git a/test_polymatrix/test_expression/test_toconstant.py b/test_polymatrix/test_expression/test_toconstant.py
index f3a0643..c85c482 100644
--- a/test_polymatrix/test_expression/test_toconstant.py
+++ b/test_polymatrix/test_expression/test_toconstant.py
@@ -1,11 +1,7 @@
 import unittest
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initquadraticinexpr import init_quadratic_in_expr
-from polymatrix.expression.init.initsymmetricexpr import init_symmetric_expr
-from polymatrix.expression.init.inittoconstantexpr import init_to_constant_expr
-from polymatrix.expression.init.inittruncateexpr import init_truncate_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestToConstant(unittest.TestCase):
@@ -27,8 +23,8 @@ class TestToConstant(unittest.TestCase):
             },
         }
 
-        expr = init_to_constant_expr(
-            underlying=init_from_terms_expr(terms=terms, shape=(2, 2)),
+        expr = polymatrix.expression.initexpressionbase.init_to_constant_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms, shape=(2, 2)),
         )
 
         state = init_expression_state(n_param=2)
diff --git a/test_polymatrix/test_expression/test_truncate.py b/test_polymatrix/test_expression/test_truncate.py
index d6cafbe..c73e98d 100644
--- a/test_polymatrix/test_expression/test_truncate.py
+++ b/test_polymatrix/test_expression/test_truncate.py
@@ -1,10 +1,7 @@
 import unittest
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initquadraticinexpr import init_quadratic_in_expr
-from polymatrix.expression.init.initsymmetricexpr import init_symmetric_expr
-from polymatrix.expression.init.inittruncateexpr import init_truncate_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestTruncate(unittest.TestCase):
@@ -30,9 +27,9 @@ class TestTruncate(unittest.TestCase):
         }
 
 
-        expr = init_truncate_expr(
-            underlying=init_from_terms_expr(terms=terms, shape=(2, 2)),
-            variables=init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
+        expr = polymatrix.expression.initexpressionbase.init_truncate_expr(
+            underlying=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms, shape=(2, 2)),
+            variables=polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
             degrees=(1, 2),
         )
 
diff --git a/test_polymatrix/test_expression/test_vstack.py b/test_polymatrix/test_expression/test_vstack.py
index 8778f8a..27e52f9 100644
--- a/test_polymatrix/test_expression/test_vstack.py
+++ b/test_polymatrix/test_expression/test_vstack.py
@@ -1,8 +1,7 @@
 import unittest
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initvstackexpr import init_v_stack_expr
+import polymatrix.expression.initexpressionbase
 
 
 class TestVStack(unittest.TestCase):
@@ -26,10 +25,10 @@ class TestVStack(unittest.TestCase):
             },
         }
 
-        expr = init_v_stack_expr(
+        expr = polymatrix.expression.initexpressionbase.init_v_stack_expr(
             underlying=(
-                init_from_terms_expr(terms=terms1, shape=(2, 2),), 
-                init_from_terms_expr(terms=terms2, shape=(2, 2),),
+                polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms1, shape=(2, 2),), 
+                polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=terms2, shape=(2, 2),),
             ),
         )
 
diff --git a/test_polymatrix/test_tomatrixrepr.py b/test_polymatrix/test_tomatrixrepr.py
index 88a35f2..dc023ca 100644
--- a/test_polymatrix/test_tomatrixrepr.py
+++ b/test_polymatrix/test_tomatrixrepr.py
@@ -1,9 +1,7 @@
 import unittest
-import polymatrix
+import polymatrix.expression.initexpressionbase
 
 from polymatrix.expressionstate.init.initexpressionstate import init_expression_state
-from polymatrix.expression.init.initfromtermsexpr import init_from_terms_expr
-from polymatrix.expression.init.initlinearinexpr import init_linear_in_expr
 
 
 class TestToMatrixRepr(unittest.TestCase):
@@ -29,12 +27,12 @@ class TestToMatrixRepr(unittest.TestCase):
             (1, 0): {((1, 1),): 1},
         }
 
-        expr = init_from_terms_expr(terms=underlying_terms, shape=(3, 1))
+        expr = polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=underlying_terms, shape=(3, 1))
 
         state = init_expression_state(n_param=2)
         state, result = polymatrix.to_matrix_repr(
             (expr,), 
-            init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
+            polymatrix.expression.initexpressionbase.init_from_terms_expr(terms=variable_terms, shape=(2, 1),),
         ).apply(state)
 
         A0 = result.data[0][0]