Add comments for code review

28 files changed, 184 insertions, 33 deletions

diff --git a/polymatrix/denserepr/from_.py b/polymatrix/denserepr/from_.py
index 4af345a..aa2e015 100644
--- a/polymatrix/denserepr/from_.py
+++ b/polymatrix/denserepr/from_.py
@@ -12,7 +12,10 @@ from polymatrix.expression.utils.getvariableindices import (
 from polymatrix.denserepr.utils.monomialtoindex import variable_indices_to_column_index
 from polymatrix.denserepr.impl import DenseReprBufferImpl, DenseReprImpl
 
-
+# NP: create a dense representation from a polymatrix expression
+# FIXME:
+#  - typing problems
+#  - create custom exception for error (instead of AssertionError & Exception)
 def from_polymatrix(
     expressions: Expression | tuple[Expression],
     variables: Expression = None,
diff --git a/polymatrix/denserepr/impl.py b/polymatrix/denserepr/impl.py
index 9007c2f..294a396 100644
--- a/polymatrix/denserepr/impl.py
+++ b/polymatrix/denserepr/impl.py
@@ -12,7 +12,7 @@ from polymatrix.expression.utils.getvariableindices import (
 
 @dataclasses.dataclass
 class DenseReprBufferImpl:
-    data: dict[int, np.ndarray]
+    data: dict[int, np.ndarray] # FIXME: use np.typing.NDArray
     n_row: int
     n_param: int
 
@@ -43,6 +43,12 @@ class DenseReprBufferImpl:
         return self.data[key]
 
 
+# NP: as discussed in meeting #6, this is not actually dense
+# NP: also representation here does not mean much, this is actually a series of
+# NP: matrices from a cone, consider changing name to reflet this fact.
+
+# NP: also this really needs documentation on how it is indexed, it is _not_
+# NP: intuitive from the method names
 @dataclasses.dataclass
 class DenseReprImpl:
     data: tuple[DenseReprBufferImpl, ...]
@@ -50,7 +56,10 @@ class DenseReprImpl:
     variable_mapping: tuple[int, ...]
     state: ExpressionState
 
+    # NP: as discussed in meeting #5 the "equations" are not equations
+    # NP: rename to something more meaningful
     def merge_matrix_equations(self):
+        # FIXME: this function is never used and shadowed by the other gen_matrices below
         def gen_matrices(index: int):
             for equations in self.data:
                 if index < len(equations):
@@ -74,7 +83,8 @@ class DenseReprImpl:
 
         return dict(gen_matrices())
 
-    def get_value(self, variable, value):
+    # NP: get variable value out of solution vector
+    def get_value(self, variable, value): # NP: variable index (state), value result from optimiz
         variable_indices = get_variable_indices_from_variable(self.state, variable)[1]
 
         def gen_value_index():
@@ -103,10 +113,13 @@ class DenseReprImpl:
     # def get_matrix(self, eq_idx: int):
     #     equations = self.data[eq_idx].data
 
-    def get_func(self, eq_idx: int):
+    # NP: Get a function that evaluates the polynomials
+    # NP: rename to something more meaningful 
+    def get_func(self, eq_idx: int): # -> Callable[[npt.NDArray], [npt.NDArray]
         equations = self.data[eq_idx].data
         max_idx = max(equations.keys())
 
+        # NP: applies slight optimization for small sizes
         if 2 <= max_idx:
 
             def func(x: np.ndarray) -> np.ndarray:
diff --git a/polymatrix/denserepr/utils/monomialtoindex.py b/polymatrix/denserepr/utils/monomialtoindex.py
index a40f2ee..a41b9fc 100644
--- a/polymatrix/denserepr/utils/monomialtoindex.py
+++ b/polymatrix/denserepr/utils/monomialtoindex.py
@@ -1,6 +1,7 @@
 import itertools
 
 
+# NP: document this function, especially magic return line
 def variable_indices_to_column_index(
     n_var: int,
     variable_indices: tuple[int, ...],
diff --git a/polymatrix/expression/expression.py b/polymatrix/expression/expression.py
index 6a84807..fe3957c 100644
--- a/polymatrix/expression/expression.py
+++ b/polymatrix/expression/expression.py
@@ -7,6 +7,7 @@ import numpy as np
 import polymatrix.expression.init
 
 from polymatrix.utils.getstacklines import get_stack_lines
+
 from polymatrix.polymatrix.abc import PolyMatrix
 from polymatrix.expressionstate.abc import ExpressionState
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
@@ -28,7 +29,10 @@ class Expression(
 ):
     @property
     @abc.abstractmethod
-    def underlying(self) -> ExpressionBaseMixin: ...
+    def underlying(self) -> ExpressionBaseMixin:
+        # NP: I know what it is but it was very confusing in the beginning.
+        # FIXME: provide documentation on how underlying works / what it means
+        ...
 
     # overwrites the abstract method of `ExpressionBaseMixin`
     def apply(self, state: ExpressionState) -> tuple[ExpressionState, PolyMatrix]:
@@ -73,6 +77,8 @@ class Expression(
     def __pow__(self, num):
         curr = 1
 
+        # FIXME: this only works for positive integral powers, consider raising
+        # an error if the power is not a positive integer
         for _ in range(num):
             curr = curr * self
 
@@ -112,6 +118,9 @@ class Expression(
             right = polymatrix.expression.init.init_from_expr_or_none(right)
 
             # delegate to upper level
+
+            # NP: what is upper level? Class inherits from ABC and base class,
+            # NP: neither has a notion of left / right
             if right is None:
                 return NotImplemented
 
@@ -486,6 +495,8 @@ class Expression(
         )
 
 
+# NP: why is this impl class here?
+# FIXME: move to impl.py
 @dataclassabc.dataclassabc(frozen=True, repr=False)
 class ExpressionImpl(Expression):
     underlying: ExpressionBaseMixin
@@ -500,6 +511,9 @@ class ExpressionImpl(Expression):
         )
 
 
+
+# NP: why is this here and not in init.py?
+# FIXME: move to init.py
 def init_expression(
     underlying: ExpressionBaseMixin,
 ):
diff --git a/polymatrix/expression/from_.py b/polymatrix/expression/from_.py
index 00f484e..ac43069 100644
--- a/polymatrix/expression/from_.py
+++ b/polymatrix/expression/from_.py
@@ -7,7 +7,8 @@ from polymatrix.expression.expression import init_expression, Expression
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 from polymatrix.statemonad.abc import StateMonad
 
-
+# NP: how is this related to expression.init.DATA_TYPE?
+# FIXME: move to typing module
 FromDataTypes = (
     str
     | np.ndarray
@@ -19,6 +20,9 @@ FromDataTypes = (
 )
 
 
+# NP: this function name makes no sense to me, 
+# NP: also why does this allow None return type? init_expression always return
+# NP: something so None object is in the underlying, not the returned object, which is confusing.
 def from_expr_or_none(
     data: FromDataTypes,
 ) -> Expression | None:
@@ -29,6 +33,8 @@ def from_expr_or_none(
     )
 
 
+# NP: from is not an ideal name because it is a keyword
+# NP: consider differnt name like make_from?
 def from_(
     data: FromDataTypes,
 ) -> Expression:
diff --git a/polymatrix/expression/impl.py b/polymatrix/expression/impl.py
index b902c22..32f2552 100644
--- a/polymatrix/expression/impl.py
+++ b/polymatrix/expression/impl.py
@@ -296,7 +296,7 @@ class ProductExprImpl(ProductExprMixin):
 class QuadraticInExprImpl(QuadraticInExprMixin):
     underlying: ExpressionBaseMixin
     monomials: ExpressionBaseMixin
-    variables: tuple
+    variables: tuple # FIXME: typing
     stack: tuple[FrameSummary]
 
     # implement custom __repr__ method that returns a representation without the stack
diff --git a/polymatrix/expression/init.py b/polymatrix/expression/init.py
index 681c0d7..17e7ae7 100644
--- a/polymatrix/expression/init.py
+++ b/polymatrix/expression/init.py
@@ -104,6 +104,7 @@ def init_elem_mult_expr(
     )
 
 
+# FIXME: typing
 def init_eval_expr(
     underlying: ExpressionBaseMixin,
     variables: typing.Union[typing.Any, tuple, dict],
@@ -161,11 +162,12 @@ def init_from_expr_or_none(
 ) -> ExpressionBaseMixin | None:
     if isinstance(data, str):
         return init_parametrize_expr(
-            underlying=init_from_expr_or_none(1),
+            underlying=init_from_expr_or_none(1), # FIXME: typing
             name=data,
         )
 
     elif isinstance(data, StateMonad):
+        # FIXME: typing
         return data.flat_map(lambda inner_data: init_from_expr_or_none(inner_data))
 
     elif isinstance(data, np.ndarray):
@@ -218,12 +220,13 @@ def init_from_expr(data: DATA_TYPE):
     return expr
 
 
+# FIXME: typing
 def init_from_terms_expr(
     data: PolyMatrixMixin | PolynomialMatrixData,
     shape: tuple[int, int] = None,
 ):
     if isinstance(data, PolyMatrixMixin):
-        shape = data.shape
+        shape = data.shape # FIXME: typing, also cf. PolyMatrixMixin.shape() comment
         poly_matrix_data = data.gen_data()
 
     else:
@@ -233,7 +236,7 @@ def init_from_terms_expr(
             poly_matrix_data = data
 
         elif isinstance(data, dict):
-            poly_matrix_data = data.items()
+            poly_matrix_data = data.items() # FIXME: typing
 
         else:
             raise Exception(f"{data=}")
@@ -283,6 +286,8 @@ def init_linear_matrix_in_expr(
     underlying: ExpressionBaseMixin,
     variable: int,
 ):
+    # FIXME: LinearMatrixInExprImpl has an abstract method variables and cannot
+    # be instantiated
     return polymatrix.expression.impl.LinearMatrixInExprImpl(
         underlying=underlying,
         variable=variable,
@@ -311,7 +316,7 @@ def init_max_expr(
 
 def init_parametrize_expr(
     underlying: ExpressionBaseMixin,
-    name: str = None,
+    name: str = None, # FIXME: typing
 ):
     if name is None:
         name = "undefined"
@@ -343,7 +348,7 @@ def init_quadratic_in_expr(
     return polymatrix.expression.impl.QuadraticInExprImpl(
         underlying=underlying,
         monomials=monomials,
-        variables=variables,
+        variables=variables, # FIXME: typing
         stack=stack,
     )
 
@@ -356,7 +361,7 @@ def init_quadratic_monomials_expr(
 
     return polymatrix.expression.impl.QuadraticMonomialsExprImpl(
         underlying=underlying,
-        variables=variables,
+        variables=variables, # FIXME: typing
     )
 
 
@@ -510,6 +515,7 @@ def init_truncate_expr(
     if inverse is None:
         inverse = False
 
+    # FIME: typing
     return polymatrix.expression.impl.TruncateExprImpl(
         underlying=underlying,
         variables=variables,
diff --git a/polymatrix/expression/op.py b/polymatrix/expression/op.py
index 0ca55ee..ae58935 100644
--- a/polymatrix/expression/op.py
+++ b/polymatrix/expression/op.py
@@ -3,6 +3,8 @@ import polymatrix.expression.init
 from polymatrix.utils.getstacklines import get_stack_lines
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 
+# NP: IIRC this file is redundant and will be removed, soon~ish?
+# TODO: remove this file in favour of operations directly via Expression objects
 
 def diff(
     expression: ExpressionBaseMixin,
diff --git a/polymatrix/expression/utils/broadcastpolymatrix.py b/polymatrix/expression/utils/broadcastpolymatrix.py
index 56f7568..c24e31e 100644
--- a/polymatrix/expression/utils/broadcastpolymatrix.py
+++ b/polymatrix/expression/utils/broadcastpolymatrix.py
@@ -4,6 +4,11 @@ from polymatrix.utils.tooperatorexception import to_operator_exception
 from polymatrix.polymatrix.abc import PolyMatrix
 
 
+# NP: "broadcast" does not mean scalar to me cf. BroadcastPolyMatrix
+# NP: broadcast means only to change into something else cf. numpy.broadcast_to
+# NP: consider renaming to broadcast_to_scalar or just get rid of broadcast
+# NP: altogether and call it to_scalar
+# FIXME: typing problems
 def broadcast_poly_matrix(
     left: PolyMatrix,
     right: PolyMatrix,
diff --git a/polymatrix/expression/utils/getderivativemonomials.py b/polymatrix/expression/utils/getderivativemonomials.py
index 9f2c48e..d671976 100644
--- a/polymatrix/expression/utils/getderivativemonomials.py
+++ b/polymatrix/expression/utils/getderivativemonomials.py
@@ -6,6 +6,25 @@ from polymatrix.expressionstate.abc import ExpressionState
 from polymatrix.polymatrix.typing import PolynomialData
 
 
+# NP: why does this return a dict?
+# NP: why doe we need considered_variables and introduce_derivatives? Wouldn't
+# NP: be better to do partial NP: derivatives and the define total derivative in
+# NP: terms of partial derivative?
+# NP:
+# NP: this is a recursive function, consider making it more clear by using the
+# NP: structure (=> Where is the code that makes sure that the recursion stops?)
+# NP:
+# NP:   def recursive_fn(args):
+# NP:       if base_case:
+# NP:           return base_case_result
+# NP:
+# NP:       # compute new args for recursion
+# NP:       result = recursive_fn(new_args)
+# NP:       # compute complete result after recursion or do tail call recursion
+# NP:       return complete_result
+# NP:
+# NP: if recursion is more complex split into smaller functions
+# FIXME: typing
 def differentiate_polynomial(
     polynomial: PolynomialData,
     diff_wrt_variable: int,
diff --git a/polymatrix/expression/utils/getvariableindices.py b/polymatrix/expression/utils/getvariableindices.py
index 9799da3..743f6b9 100644
--- a/polymatrix/expression/utils/getvariableindices.py
+++ b/polymatrix/expression/utils/getvariableindices.py
@@ -5,6 +5,7 @@ from polymatrix.expressionstate.abc import ExpressionState
 from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
 
 
+# FIXME: typing
 def get_variable_indices_from_variable(
     state: ExpressionState,
     variable: ExpressionBaseMixin | int | typing.Any,
diff --git a/polymatrix/expressionstate/abc.py b/polymatrix/expressionstate/abc.py
index 0c88dd4..bdc4484 100644
--- a/polymatrix/expressionstate/abc.py
+++ b/polymatrix/expressionstate/abc.py
@@ -1,5 +1,6 @@
 from polymatrix.expressionstate.mixins import ExpressionStateMixin
 
 
+# NP: "state" of an expression that maps indices to variable / parameter objects
 class ExpressionState(ExpressionStateMixin):
     pass
diff --git a/polymatrix/expressionstate/init.py b/polymatrix/expressionstate/init.py
index 31b40f0..8fdbabf 100644
--- a/polymatrix/expressionstate/init.py
+++ b/polymatrix/expressionstate/init.py
@@ -5,6 +5,8 @@ def init_expression_state(
     n_param: int = None,
     offset_dict: dict = None,
 ):
+    # FIXME: just set the defaults above instead of None which btw is not
+    # allowed by the type checker ("implicit none is not allowed")
     if n_param is None:
         n_param = 0
 
diff --git a/polymatrix/expressionstate/mixins.py b/polymatrix/expressionstate/mixins.py
index ae8d0fb..7122e80 100644
--- a/polymatrix/expressionstate/mixins.py
+++ b/polymatrix/expressionstate/mixins.py
@@ -5,6 +5,7 @@ import typing
 from polymatrix.statemonad.mixins import StateCacheMixin
 
 
+# NP: "state" of an expression that maps indices to variable / parameter objects
 class ExpressionStateMixin(
     StateCacheMixin,
 ):
@@ -14,7 +15,6 @@ class ExpressionStateMixin(
         """
         number of parameters used in polynomial matrix expressions
         """
-
         ...
 
     @property
@@ -24,27 +24,38 @@ class ExpressionStateMixin(
         a variable consists of one or more parameters indexed by a start
         and an end index
         """
-
+        # NP: I call a thing (start, end) a _range_ or _interval_ of indices to index multiple varilables
+        # NP: offset_dict is confusing IMHO, consider renaming
         ...
 
     @property
     @abc.abstractmethod
-    def auxillary_equations(self) -> dict[int, dict[tuple[int], float]]: ...
+    def auxillary_equations(self) -> dict[int, dict[tuple[int], float]]:
+        # NP: TODO explanation of how auxiliary equaitons work
+        ...
 
+    # NP: get a variable name from the offset, hovever you can only ever
+    # NP: get names using offsets, so maybe rename to something more intutive like
+    # FIXME: rename to get_variable_name() or get_parameter_name()
     def get_name_from_offset(self, offset: int):
         for variable, (start, end) in self.offset_dict.items():
             if start <= offset < end:
                 return f"{str(variable)}_{offset-start}"
 
+    # NP: key does not mean anything for someone who does not know how this class works inside
+    # FIXME: rename to just get_variable() or get_parameter()
     def get_key_from_offset(self, offset: int):
         for variable, (start, end) in self.offset_dict.items():
             if start <= offset < end:
                 return variable
 
+    # NP: register a variable / parameter into the state object
+    # NP: why are you allowed to not give a key (use case)? also, rename key to variable / parameter
+    # NP: register() is good, other good names are index(), index_variable(), index_parameter()
     def register(
         self,
         n_param: int,
-        key: typing.Any = None,
+        key: typing.Any = None, # NP: Any is close to useless, specify type
     ) -> "ExpressionStateMixin":
         if key is None:
             updated_state = dataclasses.replace(
diff --git a/polymatrix/polymatrix/abc.py b/polymatrix/polymatrix/abc.py
index ba909c4..8b64caf 100644
--- a/polymatrix/polymatrix/abc.py
+++ b/polymatrix/polymatrix/abc.py
@@ -4,4 +4,5 @@ from polymatrix.polymatrix.mixins import PolyMatrixAsDictMixin
 
 
 class PolyMatrix(PolyMatrixAsDictMixin, abc.ABC):
-    pass
+    # NP: Interface for matrix whose entries are multivariate polynomials
+	pass
diff --git a/polymatrix/polymatrix/impl.py b/polymatrix/polymatrix/impl.py
index 1777ea8..27b6260 100644
--- a/polymatrix/polymatrix/impl.py
+++ b/polymatrix/polymatrix/impl.py
@@ -7,11 +7,11 @@ from polymatrix.polymatrix.typing import PolynomialData
 
 @dataclassabc.dataclassabc(frozen=True)
 class PolyMatrixImpl(PolyMatrix):
-    data: dict[tuple[int, int], PolynomialData]
+    data: dict[tuple[int, int], PolynomialData] # FIXME: use polymatrix.typing
     shape: tuple[int, int]
 
 
 @dataclassabc.dataclassabc(frozen=True)
 class BroadcastPolyMatrixImpl(BroadcastPolyMatrixMixin):
-    data: tuple[tuple[int], float]
+    data: tuple[tuple[int], float] # FIXME: use polymatrix.typing
     shape: tuple[int, int]
diff --git a/polymatrix/polymatrix/init.py b/polymatrix/polymatrix/init.py
index 55f941b..c4912c7 100644
--- a/polymatrix/polymatrix/init.py
+++ b/polymatrix/polymatrix/init.py
@@ -2,6 +2,7 @@ from polymatrix.polymatrix.impl import BroadcastPolyMatrixImpl, PolyMatrixImpl
 from polymatrix.polymatrix.typing import PolynomialData
 
 
+# FIXME: use polymatrix.typing
 def init_poly_matrix(
     data: dict[tuple[int, int], PolynomialData],
     shape: tuple[int, int],
@@ -12,6 +13,8 @@ def init_poly_matrix(
     )
 
 
+# NP: cosider renaming to scalar cf. comment on broadcast_poly_matrix and BroadcastPolymatrix
+# FIXME: use polymatrix.typing
 def init_broadcast_poly_matrix(
     data: PolynomialData,
     shape: tuple[int, int],
diff --git a/polymatrix/polymatrix/mixins.py b/polymatrix/polymatrix/mixins.py
index 997ff9f..d4381e7 100644
--- a/polymatrix/polymatrix/mixins.py
+++ b/polymatrix/polymatrix/mixins.py
@@ -5,13 +5,18 @@ from polymatrix.polymatrix.typing import PolynomialData
 
 
 class PolyMatrixMixin(abc.ABC):
+    # NP: is this still in use? cf. Filter operation, that breaks the
+    # NP: assumption that shape can be computed in advance
+    # FIXME: classmethod as properties are no longer allowed, change to abstractmethod
     @property
     @abc.abstractclassmethod
     def shape(self) -> tuple[int, int]: ...
 
-    def gen_data(
-        self,
-    ) -> typing.Generator[tuple[tuple[int, int], PolynomialData], None, None]:
+    # NP: the "data" is an iteration over the entries of the matrix
+    # FIXME: rename to something meaningful like entries()
+    # NP: the fact that it returns the inner dictionary directly leaks the abstraction,
+    # NP: as discussed in in meeting #6 IMHO it is better to return indices and query like in mdpoly
+    def gen_data(self) -> typing.Generator[tuple[tuple[int, int], PolynomialData], None, None]:
         for row in range(self.shape[0]):
             for col in range(self.shape[1]):
                 polynomial = self.get_poly(row, col)
@@ -20,10 +25,17 @@ class PolyMatrixMixin(abc.ABC):
 
                 yield (row, col), polynomial
 
+    
+    # NP: get a polynomial at a given entry
+    # NP: why is this an abstract class method if it uses self?
+    # FIXME: change to abstractmethod
     @abc.abstractclassmethod
     def get_poly(self, row: int, col: int) -> PolynomialData | None: ...
 
 
+# NP: As far as I can tell the dict is the only way you ever store a polymatrix
+# NP: (cf. PolynomialData), why separate interface for PolyMatrixAsDictMixin?
+# NP: Also applies to BroadcastPolyMatrixMixin class below.
 class PolyMatrixAsDictMixin(
     PolyMatrixMixin,
     abc.ABC,
@@ -33,11 +45,16 @@ class PolyMatrixAsDictMixin(
     def data(self) -> dict[tuple[int, int], PolynomialData]: ...
 
     # overwrites the abstract method of `PolyMatrixMixin`
+    # NP: Returning none is not great design, from math perspective it would
+    # NP: make more sense to return zero
     def get_poly(self, row: int, col: int) -> PolynomialData | None:
         if (row, col) in self.data:
             return self.data[row, col]
 
 
+# NP: What does broadcast mean here? Is this a scalar multivariate polynomial
+# NP: that does not require the (row, col) indexing?
+# NP: If yes, consider new name ScalarPolynomial, ScalarPolyMatrix?
 class BroadcastPolyMatrixMixin(
     PolyMatrixMixin,
     abc.ABC,
diff --git a/polymatrix/polymatrix/typing.py b/polymatrix/polymatrix/typing.py
index f8fb9f6..e304e5c 100644
--- a/polymatrix/polymatrix/typing.py
+++ b/polymatrix/polymatrix/typing.py
@@ -1,7 +1,9 @@
 # monomial x1**2 x2
 # with indices {x1: 0, x2: 1}
 # is represented as ((0, 2), (1, 1))
-MonomialData = tuple[tuple[int, int], ...]
 
+# NP: consider renaming Data to Index (data does not mean anything to me)
+# NP: introduce index classes?
+MonomialData = tuple[tuple[int, int], ...]
 PolynomialData = dict[MonomialData, int | float]
 PolynomialMatrixData = dict[tuple[int, int], PolynomialData]
diff --git a/polymatrix/polymatrix/utils/mergemonomialindices.py b/polymatrix/polymatrix/utils/mergemonomialindices.py
index 5a05265..ec8f309 100644
--- a/polymatrix/polymatrix/utils/mergemonomialindices.py
+++ b/polymatrix/polymatrix/utils/mergemonomialindices.py
@@ -2,6 +2,8 @@ from polymatrix.polymatrix.typing import MonomialData
 from polymatrix.polymatrix.utils.sortmonomialindices import sort_monomial_indices
 
 
+# NP: this is a werid case of the product
+# NP: merge is not intutive IMHO, find better name
 def merge_monomial_indices(
     monomials: tuple[MonomialData, ...],
 ) -> MonomialData:
diff --git a/polymatrix/polymatrix/utils/multiplypolynomial.py b/polymatrix/polymatrix/utils/multiplypolynomial.py
index fe8b13b..4b54813 100644
--- a/polymatrix/polymatrix/utils/multiplypolynomial.py
+++ b/polymatrix/polymatrix/utils/multiplypolynomial.py
@@ -5,6 +5,7 @@ from polymatrix.polymatrix.utils.mergemonomialindices import merge_monomial_indi
 from polymatrix.polymatrix.typing import PolynomialData
 
 
+# NP: compute index of product of polynomials
 def multiply_polynomial(
     left: PolynomialData,
     right: PolynomialData,
diff --git a/polymatrix/polymatrix/utils/sortmonomialindices.py b/polymatrix/polymatrix/utils/sortmonomialindices.py
index cf6fa61..ff4ce59 100644
--- a/polymatrix/polymatrix/utils/sortmonomialindices.py
+++ b/polymatrix/polymatrix/utils/sortmonomialindices.py
@@ -1,6 +1,7 @@
 from polymatrix.polymatrix.typing import MonomialData
 
 
+# NP: sort mononial indices with respect to variable index in state object
 def sort_monomial_indices(
     monomial: MonomialData,
 ) -> MonomialData:
diff --git a/polymatrix/polymatrix/utils/sortmonomials.py b/polymatrix/polymatrix/utils/sortmonomials.py
index 810ef9b..276a521 100644
--- a/polymatrix/polymatrix/utils/sortmonomials.py
+++ b/polymatrix/polymatrix/utils/sortmonomials.py
@@ -1,6 +1,6 @@
 from polymatrix.polymatrix.typing import MonomialData
 
-
+# NP: Sort list of monomials according to ... what?
 def sort_monomials(
     monomials: tuple[MonomialData],
 ) -> tuple[MonomialData]:
diff --git a/polymatrix/polymatrix/utils/splitmonomialindices.py b/polymatrix/polymatrix/utils/splitmonomialindices.py
index 9f15670..29890b2 100644
--- a/polymatrix/polymatrix/utils/splitmonomialindices.py
+++ b/polymatrix/polymatrix/utils/splitmonomialindices.py
@@ -1,6 +1,7 @@
 from polymatrix.polymatrix.typing import MonomialData
 
 
+# NP: what does this function do? split according to what?
 def split_monomial_indices(
     monomial: MonomialData,
 ) -> tuple[MonomialData, MonomialData]:
diff --git a/polymatrix/polymatrix/utils/subtractmonomialindices.py b/polymatrix/polymatrix/utils/subtractmonomialindices.py
index 3a671c5..f79eda9 100644
--- a/polymatrix/polymatrix/utils/subtractmonomialindices.py
+++ b/polymatrix/polymatrix/utils/subtractmonomialindices.py
@@ -3,10 +3,13 @@ from polymatrix.polymatrix.typing import MonomialData
 from polymatrix.polymatrix.utils.sortmonomialindices import sort_monomial_indices
 
 
+# NP: consider making a module containing all exceptions
 class SubtractError(Exception):
     pass
 
 
+# NP: set difference / division of polynomials?
+# NP: name is very confusing, find better name once is clear what it does
 def subtract_monomial_indices(
     m1: MonomialData,
     m2: MonomialData,
diff --git a/polymatrix/statemonad/__init__.py b/polymatrix/statemonad/__init__.py
index bff288a..48e369e 100644
--- a/polymatrix/statemonad/__init__.py
+++ b/polymatrix/statemonad/__init__.py
@@ -1,7 +1,8 @@
 from polymatrix.statemonad.init import init_state_monad
 from polymatrix.statemonad.abc import StateMonad
 
-
+# NP: this is the unit operation for the monad, why not move it inside the
+# NP: monad class?
 def from_(val):
     def func(state):
         return state, val
@@ -9,6 +10,8 @@ def from_(val):
     return init_state_monad(func)
 
 
+# NP: duplicate-ish with StateMonadMixin.zip, this one is more generic
+# NP: consider moving this into the mixin and deleting this function
 def zip(monads: tuple[StateMonad]):
     def zip_func(state):
         values = tuple()
diff --git a/polymatrix/statemonad/mixins.py b/polymatrix/statemonad/mixins.py
index 4b65353..e2bfda9 100644
--- a/polymatrix/statemonad/mixins.py
+++ b/polymatrix/statemonad/mixins.py
@@ -14,31 +14,57 @@ U = TypeVar("U")
 V = TypeVar("V")
 
 
+# NP: Monadic type, use shorthand M for StateMonadMixin.
+# NP: (will use haskell-like notation in this comment)
+#
+# NP: typical operations for monad (do you agree wit this? If not please explain your conventions)
+#
+# NP: - unit operation (aka return, bad name) take a value `u :: U` and make a `m :: M[U]`
+# NP:   you often call operation "from"
+#
+# NP: - map operation (aka lift) take a function (U -> V) and make a new function (M[U] -> M[V])
+#
+# NP: - bind operation (aka flat map) take a function (U -> M[V])
+# NP:   and make a new function (M[U] -> M[V])
+# NP    you call this operation flat_map
+#
+# NP: - apply operation take a M[U -> V] and make (M[U] -> M[V])
+#
+# NP: - zip operation take a function (U -> V -> W)
+# NP:   and make a new function (M[U] -> M[V] -> M[W])
+#
+# NP: TODO: text comparing the above to implementation below
 class StateMonadMixin(
     Generic[State, U],
     abc.ABC,
 ):
     @property
     @abc.abstractmethod
-    def apply_func(self) -> Callable[[State], tuple[State, U]]: ...
+    def apply_func(self) -> Callable[[State], tuple[State, U]]:
+        # NP: TODO comment
+        ...
 
-    def map(self, fn: Callable[[U], V]) -> "StateMonadMixin[State, V]":
+    # NP: typing, use from __future__ import annotations
+    def map(self, fn: Callable[[U], V]) -> 'StateMonadMixin[State, V]':
+        # NP: add functools.wrap(fn) decorator to copy docstrings etc.
         def internal_map(state: State) -> Tuple[State, U]:
             n_state, val = self.apply(state)
             return n_state, fn(val)
 
         return dataclasses.replace(self, apply_func=internal_map)
 
-    def flat_map(
-        self, fn: Callable[[U], "StateMonadMixin"]
-    ) -> "StateMonadMixin[State, V]":
+    # NP: shouldn't typing be
+    # NP: flat_map(self, fn: Callable[[U], StateMonadMixin[State, V]]) -> StateMonadMixin[State, V]
+    def flat_map(self, fn: Callable[[U], 'StateMonadMixin']) -> 'StateMonadMixin[State, V]':
+        # NP: add functools.wrap(fn) decorator
         def internal_map(state: State) -> Tuple[State, V]:
             n_state, val = self.apply(state)
             return fn(val).apply(n_state)
 
         return dataclasses.replace(self, apply_func=internal_map)
 
-    def zip(self, other: "StateMonadMixin") -> "StateMonadMixin":
+    # FIXME: typing
+    def zip(self, other: 'StateMonadMixin') -> 'StateMonadMixin':
         def internal_map(state: State) -> Tuple[State, V]:
             state, val1 = self.apply(state)
             state, val2 = other.apply(state)
@@ -46,7 +72,8 @@ class StateMonadMixin(
 
         return dataclasses.replace(self, apply_func=internal_map)
 
-    def cache(self) -> "StateMonadMixin":
+    # FIXME: typing
+    def cache(self) -> 'StateMonadMixin':
         def internal_map(state: State) -> Tuple[State, V]:
             if self in state.cache:
                 return state, state.cache[self]
@@ -62,8 +89,12 @@ class StateMonadMixin(
 
         return dataclasses.replace(self, apply_func=internal_map)
 
+    # NP: Need to find consistent naming and explain somewhere naming convention
+    # NP: of monad operations (is this from scala conventions? I have never used scala)
     def apply(self, state: State) -> Tuple[State, U]:
         return self.apply_func(state)
 
+    # NP: find better name or add explaination somewhere
+    # NP: (I know what it does but the name is very vague)
     def read(self, state: State) -> U:
         return self.apply_func(state)[1]
diff --git a/polymatrix/utils/getstacklines.py b/polymatrix/utils/getstacklines.py
index 6fdd69f..aed2915 100644
--- a/polymatrix/utils/getstacklines.py
+++ b/polymatrix/utils/getstacklines.py
@@ -12,6 +12,8 @@ class FrameSummary:
 
 def get_stack_lines(index: int = 2) -> tuple[FrameSummary]:
     def gen_stack_lines():
+        # NP: according to the stdlib traceback documentation these objects
+        # NP: should not be left around indefintely as they can cause memory leaks
         for obj in traceback.extract_stack()[:-index]:
             yield FrameSummary(
                 filename=obj.filename,