17 files changed, 260 insertions, 77 deletions

diff --git a/polymatrix/expression/expression.py b/polymatrix/expression/expression.py
index 45ac41b..f190482 100644
--- a/polymatrix/expression/expression.py
+++ b/polymatrix/expression/expression.py
@@ -9,7 +9,7 @@ import polymatrix.expression.init
 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, filter_, degree
+from polymatrix.expression.op import diff, integrate, linear_in, linear_monomials, legendre, filter_, degree
 from polymatrix.polymatrix.abc import PolyMatrix
 from polymatrix.expressionstate.abc import ExpressionState
 
@@ -232,6 +232,21 @@ class Expression(
             ),
         )
 
+    def integrate(
+        self, 
+        variables: 'Expression',
+        from_: tuple[float, ...],
+        to: tuple[float, ...]
+    ) -> 'Expression':
+        return self.copy(
+            underlying=integrate(
+                expression=self,
+                variables=variables,
+                from_=from_,
+                to=to,
+            ),
+        )
+
     def linear_matrix_in(self, variable: 'Expression') -> 'Expression':
         return self.copy(
             underlying=polymatrix.expression.init.init_linear_matrix_in_expr(
diff --git a/polymatrix/expression/impl.py b/polymatrix/expression/impl.py
index 937676d..22c33bc 100644
--- a/polymatrix/expression/impl.py
+++ b/polymatrix/expression/impl.py
@@ -1,5 +1,6 @@
 import typing
 import dataclassabc
+from polymatrix.expression.mixins.integrateexprmixin import IntegrateExprMixin
 
 from polymatrix.expression.mixins.legendreseriesmixin import LegendreSeriesMixin
 from polymatrix.expression.mixins.productexprmixin import ProductExprMixin
@@ -25,7 +26,7 @@ from polymatrix.expression.mixins.filterexprmixin import FilterExprMixin
 from polymatrix.expression.mixins.fromsymmetricmatrixexprmixin import \
     FromSymmetricMatrixExprMixin
 from polymatrix.expression.mixins.fromtupleexprmixin import FromTupleExprMixin
-from polymatrix.expression.mixins.fromtermsexprmixin import FromTermsExprMixin
+from polymatrix.expression.mixins.fromtermsexprmixin import FromPolynomialDataExprMixin, PolynomialMatrixTupledData
 from polymatrix.expression.mixins.getitemexprmixin import GetItemExprMixin
 from polymatrix.expression.mixins.halfnewtonpolytopeexprmixin import \
     HalfNewtonPolytopeExprMixin
@@ -75,6 +76,7 @@ class AdditionExprImpl(AdditionExprMixin):
 	right: ExpressionBaseMixin
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(left={self.left}, right={self.right})'
 
@@ -101,6 +103,7 @@ class DerivativeExprImpl(DerivativeExprMixin):
 	introduce_derivatives: bool
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(variables={self.variables}, underlying={self.underlying})'
 
@@ -126,6 +129,7 @@ class DivisionExprImpl(DivisionExprMixin):
 	right: ExpressionBaseMixin
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(left={self.left}, right={self.right})'
 
@@ -164,12 +168,13 @@ class FromTupleExprImpl(FromTupleExprMixin):
 	data: tuple[tuple[float]]
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(data={self.data})'
 
 @dataclassabc.dataclassabc(frozen=True)
-class FromTermsExprImpl(FromTermsExprMixin):
-	terms: tuple
+class FromPolynomialDataExprImpl(FromPolynomialDataExprMixin):
+	data: PolynomialMatrixTupledData
 	shape: tuple[int, int]
 
 
@@ -187,6 +192,19 @@ class HalfNewtonPolytopeExprImpl(HalfNewtonPolytopeExprMixin):
 
 
 @dataclassabc.dataclassabc(frozen=True)
+class IntegrateExprImpl(IntegrateExprMixin):
+	underlying: ExpressionBaseMixin
+	variables: ExpressionBaseMixin
+	from_: tuple[float, ...]
+	to: tuple[float, ...]
+	stack: tuple[FrameSummary]
+
+	# implement custom __repr__ method that returns a representation without the stack
+	def __repr__(self):
+		return f'{self.__class__.__name__}(underlying={self.underlying}, variables={self.variables}, from_={self.from_}, to={self.to})'
+
+
+@dataclassabc.dataclassabc(frozen=True)
 class LinearInExprImpl(LinearInExprMixin):
 	underlying: ExpressionBaseMixin
 	monomials: ExpressionBaseMixin
@@ -212,6 +230,7 @@ class LegendreSeriesImpl(LegendreSeriesMixin):
 	degrees: tuple[int, ...] | None
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(underlying={self.underlying}, degrees={self.degrees})'
 
@@ -222,6 +241,7 @@ class MatrixMultExprImpl(MatrixMultExprMixin):
 	right: ExpressionBaseMixin
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(left={self.left}, right={self.right})'
 
@@ -231,6 +251,7 @@ class DegreeExprImpl(DegreeExprMixin):
 	underlying: ExpressionBaseMixin
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(underlying={self.underlying})'
 
@@ -261,6 +282,7 @@ class ProductExprImpl(ProductExprMixin):
 	degrees: tuple[int, ...] | None
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(underlying={self.underlying}, degrees={self.degrees})'
 
@@ -272,6 +294,7 @@ class QuadraticInExprImpl(QuadraticInExprMixin):
 	variables: tuple
 	stack: tuple[FrameSummary]
 
+	# implement custom __repr__ method that returns a representation without the stack
 	def __repr__(self):
 		return f'{self.__class__.__name__}(variables={self.variables}, monomials={self.monomials}, underlying={self.underlying})'
 
diff --git a/polymatrix/expression/init.py b/polymatrix/expression/init.py
index f8092d2..7465f26 100644
--- a/polymatrix/expression/init.py
+++ b/polymatrix/expression/init.py
@@ -5,6 +5,7 @@ import sympy
 import polymatrix.expression.impl
 
 from polymatrix.polymatrix.mixins import PolyMatrixMixin
+from polymatrix.polymatrix.typing import PolynomialMatrixData
 from polymatrix.statemonad.abc import StateMonad
 from polymatrix.utils.getstacklines import FrameSummary
 from polymatrix.utils.getstacklines import get_stack_lines
@@ -220,31 +221,31 @@ def init_from_expr(
 
 
 def init_from_terms_expr(
-    terms: typing.Union[tuple, PolyMatrixMixin],
+    data: PolyMatrixMixin | PolynomialMatrixData,
     shape: tuple[int, int] = None,
 ):
 
-    if isinstance(terms, PolyMatrixMixin):
-        shape = terms.shape
-        gen_terms = terms.gen_data()
+    if isinstance(data, PolyMatrixMixin):
+        shape = data.shape
+        poly_matrix_data = data.gen_data()
 
     else:
         assert shape is not None
 
-        if isinstance(terms, tuple):
-            gen_terms = terms
+        if isinstance(data, tuple):
+            poly_matrix_data = data
 
-        elif isinstance(terms, dict):
-            gen_terms = terms.items()
+        elif isinstance(data, dict):
+            poly_matrix_data = data.items()
 
         else:
-            raise Exception(f'{terms=}')
+            raise Exception(f'{data=}')
 
     # Expression needs to be hashable
-    terms_formatted = tuple((key, tuple(monomials.items())) for key, monomials in gen_terms)
+    data_as_tuple = tuple((coord, tuple(polynomial.items())) for coord, polynomial in poly_matrix_data)
 
-    return polymatrix.expression.impl.FromTermsExprImpl(
-        terms=terms_formatted,
+    return polymatrix.expression.impl.FromPolynomialDataExprImpl(
+        terms=data_as_tuple,
         shape=shape,
     )
 
diff --git a/polymatrix/expression/mixins/derivativeexprmixin.py b/polymatrix/expression/mixins/derivativeexprmixin.py
index 0b6ec6e..64f2640 100644
--- a/polymatrix/expression/mixins/derivativeexprmixin.py
+++ b/polymatrix/expression/mixins/derivativeexprmixin.py
@@ -48,7 +48,7 @@ class DerivativeExprMixin(ExpressionBaseMixin):
     ) -> tuple[ExpressionState, PolyMatrix]:
 
         state, underlying = self.underlying.apply(state=state)        
-        state, diff_wrt_variables = get_variable_indices_from_variable(state, self.variables)
+        state, variables = get_variable_indices_from_variable(state, self.variables)
 
         if not (underlying.shape[1] == 1):
             raise AssertionError(to_operator_exception(
@@ -60,27 +60,28 @@ class DerivativeExprMixin(ExpressionBaseMixin):
 
         for row in range(underlying.shape[0]):
 
-            polynomial = underlying.get_poly(row, 0)
-            if polynomial is None:
+            underlying_poly = underlying.get_poly(row, 0)
+
+            if underlying_poly is None:
                 continue
 
             # derivate each variable and map result to the corresponding column
-            for col, diff_wrt_variable in enumerate(diff_wrt_variables):
+            for col, variable in enumerate(variables):
 
-                state, derivation = differentiate_polynomial(
-                    polynomial=polynomial,
-                    diff_wrt_variable=diff_wrt_variable,
+                state, diff_polynomial = differentiate_polynomial(
+                    polynomial=underlying_poly,
+                    diff_wrt_variable=variable,
                     state=state,
-                    considered_variables=set(diff_wrt_variables),
+                    considered_variables=set(variables),
                     introduce_derivatives=self.introduce_derivatives,
                 )
 
-                if 0 < len(derivation):
-                    poly_matrix_data[row, col] = derivation
+                if 0 < len(diff_polynomial):
+                    poly_matrix_data[row, col] = diff_polynomial
 
         poly_matrix = init_poly_matrix(
             data=poly_matrix_data,
-            shape=(underlying.shape[0], len(diff_wrt_variables)),
+            shape=(underlying.shape[0], len(variables)),
         )
 
         return state, poly_matrix   
diff --git a/polymatrix/expression/mixins/divergenceexprmixin.py b/polymatrix/expression/mixins/divergenceexprmixin.py
index fc12abd..fbe6da0 100644
--- a/polymatrix/expression/mixins/divergenceexprmixin.py
+++ b/polymatrix/expression/mixins/divergenceexprmixin.py
@@ -34,7 +34,7 @@ class DivergenceExprMixin(ExpressionBaseMixin):
         assert underlying.shape[1] == 1, f'{underlying.shape=}'
         assert len(variables) == underlying.shape[0], f'{variables=}, {underlying.shape=}'
 
-        monomial_terms = collections.defaultdict(float)
+        polynomial_data = collections.defaultdict(float)
 
         for row, variable in enumerate(variables):
 
@@ -52,10 +52,10 @@ class DivergenceExprMixin(ExpressionBaseMixin):
             )
 
             for monomial, value in derivation.items():
-                monomial_terms[monomial] += value
+                polynomial_data[monomial] += value
 
         poly_matrix = init_poly_matrix(
-            data={(0, 0): monomial_terms},
+            data={(0, 0): dict(polynomial_data)},
             shape=(1, 1),
         )
 
diff --git a/polymatrix/expression/mixins/divisionexprmixin.py b/polymatrix/expression/mixins/divisionexprmixin.py
index 65e6bc3..71bb86f 100644
--- a/polymatrix/expression/mixins/divisionexprmixin.py
+++ b/polymatrix/expression/mixins/divisionexprmixin.py
@@ -72,26 +72,26 @@ class DivisionExprMixin(ExpressionBaseMixin):
         for row in range(left.shape[0]):
             for col in range(left.shape[1]):
 
-                underlying_terms = left.get_poly(row, col)
-                if underlying_terms is None:
+                underlying_poly = left.get_poly(row, col)
+                if underlying_poly is None:
                     continue
 
-                def gen_monomial_terms():
-                    for monomial, value in underlying_terms.items():
+                def gen_polynomial():
+                    for monomial, value in underlying_poly.items():
                         yield monomial + ((division_variable, 1),), value
 
-                poly_matrix_data[row, col] = dict(gen_monomial_terms())
+                poly_matrix_data[row, col] = dict(gen_polynomial())
 
-        def gen_auxillary_terms():
+        def gen_auxillary_polynomials():
             for monomial, value in right_poly.items():
                 yield monomial + ((division_variable, 1),), value
 
-        auxillary_terms = dict(gen_auxillary_terms())
+        auxillary_poly = dict(gen_auxillary_polynomials())
 
-        if tuple() not in auxillary_terms:
-            auxillary_terms[tuple()] = 0
+        if tuple() not in auxillary_poly:
+            auxillary_poly[tuple()] = 0
 
-        auxillary_terms[tuple()] -= 1
+        auxillary_poly[tuple()] -= 1
 
         poly_matrix = init_poly_matrix(
             data=poly_matrix_data,
@@ -100,7 +100,7 @@ class DivisionExprMixin(ExpressionBaseMixin):
 
         state = dataclasses.replace(
             state, 
-            auxillary_equations=state.auxillary_equations | {division_variable: auxillary_terms},
+            auxillary_equations=state.auxillary_equations | {division_variable: auxillary_poly},
             cache=state.cache | {self: poly_matrix},
         )
 
diff --git a/polymatrix/expression/mixins/fromtermsexprmixin.py b/polymatrix/expression/mixins/fromtermsexprmixin.py
index 22512d5..c0d4dd5 100644
--- a/polymatrix/expression/mixins/fromtermsexprmixin.py
+++ b/polymatrix/expression/mixins/fromtermsexprmixin.py
@@ -1,18 +1,23 @@
 
 import abc
-from numpy import poly
-import sympy
 
 from polymatrix.polymatrix.init import init_poly_matrix
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.expressionstate.mixins import ExpressionStateMixin
 from polymatrix.polymatrix.mixins import PolyMatrixMixin
+from polymatrix.polymatrix.typing import MonomialData
 
 
-class FromTermsExprMixin(ExpressionBaseMixin):
+PolynomialTupledData = tuple[tuple[MonomialData, float], ...]
+PolynomialMatrixTupledData = tuple[tuple[tuple[int, int], PolynomialTupledData], ...]
+
+
+class FromPolynomialDataExprMixin(ExpressionBaseMixin):
+
+    # an Expression needs to be hashable
     @property
     @abc.abstractmethod
-    def terms(self) -> tuple[tuple[tuple[tuple[int, int], ...], float], ...]:
+    def data(self) -> PolynomialMatrixTupledData:
         pass
 
     @property
@@ -26,7 +31,7 @@ class FromTermsExprMixin(ExpressionBaseMixin):
         state: ExpressionStateMixin,
     ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]:
 
-        data = {coord: dict(monomials) for coord, monomials in self.terms}
+        data = {coord: dict(polynomial) for coord, polynomial in self.data}
 
         poly_matrix = init_poly_matrix(
             data=data,
diff --git a/polymatrix/expression/mixins/integrateexprmixin.py b/polymatrix/expression/mixins/integrateexprmixin.py
new file mode 100644
index 0000000..4104995
--- /dev/null
+++ b/polymatrix/expression/mixins/integrateexprmixin.py
@@ -0,0 +1,116 @@
+
+import abc
+import collections
+import itertools
+import typing
+
+from polymatrix.polymatrix.init import init_poly_matrix
+from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
+from polymatrix.polymatrix.abc import PolyMatrix
+from polymatrix.expressionstate.abc import ExpressionState
+from polymatrix.expression.utils.getderivativemonomials import differentiate_polynomial
+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 IntegrateExprMixin(ExpressionBaseMixin):
+    """
+    integrate w.r.t. x on interval (a, b):
+
+    [[c*x**2]]  ->  [[(1/3)*c*(b**3 - a**3)]]
+    """
+
+    @property
+    @abc.abstractmethod
+    def underlying(self) -> ExpressionBaseMixin:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def variables(self) -> ExpressionBaseMixin:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def from_(self) -> tuple[float, ...]:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def to(self) -> tuple[float, ...]:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def stack(self) -> tuple[FrameSummary]:
+        ...
+
+    # overwrites the abstract method of `ExpressionBaseMixin`
+    def apply(
+        self, 
+        state: ExpressionState,
+    ) -> tuple[ExpressionState, PolyMatrix]:
+
+        state, underlying = self.underlying.apply(state=state)        
+        state, variables = get_variable_indices_from_variable(state, self.variables)
+
+        if not (underlying.shape[1] == 1):
+            raise AssertionError(to_operator_exception(
+                message=f'{underlying.shape[1]=} is not 1',
+                stack=self.stack,
+            ))
+        
+        if not (len(variables) == len(self.from_)):
+            raise AssertionError(to_operator_exception(
+                message=f'length of tuples do not match: {variables=}, {self.from_=}, {self.to}',
+                stack=self.stack,
+            ))
+
+        poly_matrix_data = {}
+
+        for row in range(underlying.shape[0]):
+
+            underlying_poly = underlying.get_poly(row, 0)
+
+            if underlying_poly is None:
+                continue
+
+            # integrate each variable and map result to the corresponding column
+            for col, (variable, from_, to) in enumerate(zip(variables, self.from_, self.to)):
+
+                integrated_polynomial = dict()
+
+                for monomial, value in underlying_poly.items():
+
+                    monomial_cnt = dict(monomial)
+
+                    if variable in monomial_cnt:
+                        
+                        def acc_integrated_term(acc, v):
+                            monomial, value = acc
+                            curr_var, curr_var_cnt = v
+
+                            if curr_var is variable:
+                                exponent = curr_var_cnt + 1
+                                return monomial, value * (to**exponent - from_**exponent)/exponent
+                            
+                            else:
+                                return monomial + v, value
+
+                        *_, (integrated_monomial, value) = itertools.accumulate(
+                            monomial,
+                            acc_integrated_term,
+                            initial=(tuple(), 1)
+                        )
+                        integrated_polynomial[integrated_monomial] = value
+
+                if 0 < len(integrated_polynomial):
+                    poly_matrix_data[row, col] = integrated_polynomial
+
+        poly_matrix = init_poly_matrix(
+            data=poly_matrix_data,
+            shape=(underlying.shape[0], len(variables)),
+        )
+
+        return state, poly_matrix   
diff --git a/polymatrix/expression/mixins/linearinexprmixin.py b/polymatrix/expression/mixins/linearinexprmixin.py
index f0b57d3..a900b7f 100644
--- a/polymatrix/expression/mixins/linearinexprmixin.py
+++ b/polymatrix/expression/mixins/linearinexprmixin.py
@@ -63,7 +63,7 @@ class LinearInExprMixin(ExpressionBaseMixin):
 
         assert underlying.shape[1] == 1
 
-        terms = collections.defaultdict(dict)
+        poly_matrix_data = collections.defaultdict(dict)
 
         for row in range(underlying.shape[0]):
 
@@ -84,10 +84,10 @@ class LinearInExprMixin(ExpressionBaseMixin):
                     else:
                         raise Exception(f'{x_monomial} not in {monomials}')
 
-                terms[row, col][p_monomial] = value
+                poly_matrix_data[row, col][p_monomial] = value
         
         poly_matrix = init_poly_matrix(
-            data=dict(terms),
+            data=dict(poly_matrix_data),
             shape=(underlying.shape[0], len(monomials)),
         )
 
diff --git a/polymatrix/expression/mixins/linearmatrixinexprmixin.py b/polymatrix/expression/mixins/linearmatrixinexprmixin.py
index 3bb8dfe..a7d5f34 100644
--- a/polymatrix/expression/mixins/linearmatrixinexprmixin.py
+++ b/polymatrix/expression/mixins/linearmatrixinexprmixin.py
@@ -31,16 +31,17 @@ class LinearMatrixInExprMixin(ExpressionBaseMixin):
 
         assert len(variable_index) == 1
 
-        terms = collections.defaultdict(dict)
+        poly_matrix_data = collections.defaultdict(dict)
 
         for row in range(underlying.shape[0]):
             for col in range(underlying.shape[1]):
 
-                underlying_terms = underlying.get_poly(row, col)
-                if underlying_terms is None:
+                underlying_poly = underlying.get_poly(row, col)
+
+                if underlying_poly is None:
                     continue
 
-                for monomial, value in underlying_terms.items():
+                for monomial, value in underlying_poly.items():
 
                     if len(monomial) == 1:
 
@@ -48,12 +49,12 @@ class LinearMatrixInExprMixin(ExpressionBaseMixin):
 
                         if variable == variable_index:
 
-                            terms[row, col][tuple()] = value
+                            poly_matrix_data[row, col][tuple()] = value
 
                             break
         
         poly_matrix = init_poly_matrix(
-            data=dict(terms),
+            data=dict(poly_matrix_data),
             shape=underlying.shape,
         )
 
diff --git a/polymatrix/expression/mixins/quadraticinexprmixin.py b/polymatrix/expression/mixins/quadraticinexprmixin.py
index 2cc9e12..95d3a2d 100644
--- a/polymatrix/expression/mixins/quadraticinexprmixin.py
+++ b/polymatrix/expression/mixins/quadraticinexprmixin.py
@@ -50,7 +50,7 @@ class QuadraticInExprMixin(ExpressionBaseMixin):
                 stack=self.stack,
             ))
 
-        terms = collections.defaultdict(lambda: collections.defaultdict(float))
+        poly_matrix_data = collections.defaultdict(lambda: collections.defaultdict(float))
 
         underlying_poly = underlying.get_poly(0, 0)
 
@@ -84,10 +84,10 @@ class QuadraticInExprMixin(ExpressionBaseMixin):
                     stack=self.stack,
                 ))
 
-            terms[row, col][p_monomial] += value
+            poly_matrix_data[row, col][p_monomial] += value
 
         poly_matrix = init_poly_matrix(
-            data=dict((k, dict(v)) for k, v in terms.items()),
+            data=dict((k, dict(v)) for k, v in poly_matrix_data.items()),
             shape=2*(len(sos_monomials),),
         )
 
diff --git a/polymatrix/expression/mixins/symmetricexprmixin.py b/polymatrix/expression/mixins/symmetricexprmixin.py
index 769e8c6..624f0ee 100644
--- a/polymatrix/expression/mixins/symmetricexprmixin.py
+++ b/polymatrix/expression/mixins/symmetricexprmixin.py
@@ -55,15 +55,15 @@ class SymmetricExprMixin(ExpressionBaseMixin):
                     return None
 
                 else:
-                    terms = collections.defaultdict(float)
+                    polynomial = collections.defaultdict(float)
 
                     # merge monomials
                     for monomials in all_monomials:
                         for monomial, value in monomials.items():
-                            terms[monomial] += value / 2
+                            polynomial[monomial] += value / 2
 
-                    # return dict(terms)
-                    return terms
+                    return dict(polynomial)
+                    # return terms
 
         polymatrix = SymmetricPolyMatrix(
             underlying=underlying,
diff --git a/polymatrix/expression/mixins/toconstantexprmixin.py b/polymatrix/expression/mixins/toconstantexprmixin.py
index c3dc4e2..53db4c2 100644
--- a/polymatrix/expression/mixins/toconstantexprmixin.py
+++ b/polymatrix/expression/mixins/toconstantexprmixin.py
@@ -26,7 +26,7 @@ class ToConstantExprMixin(ExpressionBaseMixin):
     ) -> tuple[ExpressionState, PolyMatrix]:
         state, underlying = self.underlying.apply(state=state)
 
-        terms = collections.defaultdict(dict)
+        poly_matrix_data = collections.defaultdict(dict)
 
         for row in range(underlying.shape[0]):
             for col in range(underlying.shape[1]):
@@ -36,10 +36,10 @@ class ToConstantExprMixin(ExpressionBaseMixin):
                     continue
 
                 if tuple() in polynomial:
-                    terms[row, col][tuple()] = polynomial[tuple()]
+                    poly_matrix_data[row, col][tuple()] = polynomial[tuple()]
         
         poly_matrix = init_poly_matrix(
-            data=dict(terms),
+            data=dict(poly_matrix_data),
             shape=underlying.shape,
         )
 
diff --git a/polymatrix/expression/op.py b/polymatrix/expression/op.py
index f19d310..612a2db 100644
--- a/polymatrix/expression/op.py
+++ b/polymatrix/expression/op.py
@@ -40,6 +40,27 @@ def filter_(
     )
 
 
+def integrate(
+    expression: ExpressionBaseMixin, 
+    variables: ExpressionBaseMixin,
+    from_: tuple[float, ...],
+    to: tuple[float, ...],
+) -> ExpressionBaseMixin:
+            
+    if not isinstance(variables, ExpressionBaseMixin):
+        variables=polymatrix.expression.init.init_from_expr(variables)
+
+    assert len(from_) == len(to)
+
+    return polymatrix.expression.impl.IntegrateExprImpl(
+		underlying=expression,
+		variables=variables,
+		from_=from_,
+        to=to,
+		stack=get_stack_lines(),
+	)
+
+
 def legendre(
     expression: ExpressionBaseMixin,
     degrees: tuple[int, ...] = None,
diff --git a/polymatrix/expression/utils/getderivativemonomials.py b/polymatrix/expression/utils/getderivativemonomials.py
index 83083aa..8fb4729 100644
--- a/polymatrix/expression/utils/getderivativemonomials.py
+++ b/polymatrix/expression/utils/getderivativemonomials.py
@@ -74,7 +74,7 @@ def differentiate_polynomial(
         confirmed_variables = tuple()
 
 
-    derivation_terms = collections.defaultdict(float)
+    diff_polynomial = collections.defaultdict(float)
 
     for monomial, value in polynomial.items():
 
@@ -104,7 +104,7 @@ def differentiate_polynomial(
 
         if diff_wrt_variable in monomial_cnt:
             diff_monomial, value = differentiate_monomial(diff_wrt_variable)
-            derivation_terms[diff_monomial] += value
+            diff_polynomial[diff_monomial] += value
 
         # only used if introduce_derivatives == True
         if introduce_derivatives:
@@ -120,7 +120,7 @@ def differentiate_polynomial(
                         dependent_variable=candidate_variable, 
                         derivation_variable=derivation_variable,
                     )
-                    derivation_terms[diff_monomial] += value
+                    diff_polynomial[diff_monomial] += value
 
-    return state, dict(derivation_terms)
+    return state, dict(diff_polynomial)
     # 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 dcb7932..238f80a 100644
--- a/polymatrix/expression/utils/getmonomialindices.py
+++ b/polymatrix/expression/utils/getmonomialindices.py
@@ -4,16 +4,16 @@ from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 
 def get_monomial_indices(
         state: ExpressionState, 
-        monomials: ExpressionBaseMixin,
-    ) -> tuple[ExpressionState, tuple[int, ...]]:
+        expression: ExpressionBaseMixin,
+) -> tuple[ExpressionState, tuple[int, ...]]:
 
-    state, monomials_obj = monomials.apply(state)
+    state, poly_matrix = expression.apply(state)
 
-    assert monomials_obj.shape[1] == 1, f'{monomials_obj.shape=}'
+    assert poly_matrix.shape[1] == 1, f'{poly_matrix.shape=}'
 
     def gen_indices():
-        for row in range(monomials_obj.shape[0]):
-            row_terms = monomials_obj.get_poly(row, 0)
+        for row in range(poly_matrix.shape[0]):
+            row_terms = poly_matrix.get_poly(row, 0)
 
             assert len(row_terms) == 1, f'{row_terms} contains more than one term'
         
diff --git a/polymatrix/polymatrix/typing.py b/polymatrix/polymatrix/typing.py
index 0540054..2518d9e 100644
--- a/polymatrix/polymatrix/typing.py
+++ b/polymatrix/polymatrix/typing.py
@@ -4,5 +4,5 @@
 # is represented as ((0, 2), (1, 1))
 MonomialData = tuple[tuple[int, int], ...]
 
-PolynomialData = dict[MonomialData, float]
-PolynomialMatrixData = dict[tuple[int, int], dict[MonomialData, float]]
\ No newline at end of file
+PolynomialData = dict[MonomialData, int | float]
+PolynomialMatrixData = dict[tuple[int, int], PolynomialData]
\ No newline at end of file