diff --git a/n3fit/src/n3fit/backends/keras_backend/MetaLayer.py b/n3fit/src/n3fit/backends/keras_backend/MetaLayer.py
index 156d22a87b..53a28c7aec 100644
--- a/n3fit/src/n3fit/backends/keras_backend/MetaLayer.py
+++ b/n3fit/src/n3fit/backends/keras_backend/MetaLayer.py
@@ -114,80 +114,3 @@ def select_initializer(ini_name, seed=None, **kwargs):
             if key in ini_args.keys():
                 ini_args[key] = value
         return ini_class(**ini_args)
-
-    # Make numpy array into a tensor
-    def np_to_tensor(self, np_array, **kwargs):
-        """
-        Given a numpy array, returns a constant tensor
-        """
-        return K.constant(np_array, **kwargs)
-
-    # Common tensor operations
-    def tensor_ones(self, shape, **kwargs):
-        """
-        Generates a tensor of ones of the given shape
-        """
-        return K.ones(shape, **kwargs)
-
-    def tensor_ones_like(self, tensor, **kwargs):
-        """
-        Generates a tensor of ones of the same shape as the input tensor
-        """
-        return K.ones_like(tensor, **kwargs)
-
-    def many_replication(self, grid, replications, axis=0, **kwargs):
-        """
-            Generates a tensor with one extra dimension:
-                a repetition of "grid" n times along the given axis
-            from keras documentation:
-            If x has shape (s1, s2, s3) and axis is 1, the output will have shape (s1, s2 * rep, s3)
-        """
-        return K.repeat_elements(grid, rep=replications, axis=axis, **kwargs)
-
-    def sum(self, tensor, axis=None, **kwargs):
-        """
-        Computes the sum of the elements of the tensor
-        """
-        return K.sum(tensor, axis=axis, **kwargs)
-
-    def tensor_product(self, tensor_x, tensor_y, axes, **kwargs):
-        """
-        Computes the tensordot product between tensor_x and tensor_y
-        """
-        return tf.tensordot(tensor_x, tensor_y, axes=axes, **kwargs)
-
-    def transpose(self, tensor, **kwargs):
-        """
-        Transpose a tensor
-        """
-        return K.transpose(tensor, **kwargs)
-
-    def boolean_mask(self, tensor, mask, axis=None, **kwargs):
-        """
-        Applies boolean mask to a tensor
-        """
-        return tf.boolean_mask(tensor, mask, axis=axis, **kwargs)
-
-    def concatenate(self, tensor_list, axis=-1, target_shape=None):
-        """
-        Concatenates a list of numbers or tenosr into a bigger tensor
-        If the target shape is given, the output is reshaped to said shape
-        """
-        concatenated_tensor = K.concatenate(tensor_list, axis=axis)
-        if target_shape:
-            return K.reshape(concatenated_tensor, target_shape)
-        else:
-            return concatenated_tensor
-
-    def flatten(self, x):
-        """ Flatten tensor x """
-        return tf.reshape(x, (-1,))
-
-    def permute_dimensions(self, tensor, permutation, **kwargs):
-        """
-        Receives a tensor and a tuple and permutes the axes of the tensor according to it.
-        i.e.
-        if permutation = (1,0,2)
-        does the permutation: axis_0 -> axis_1, axis_1 -> axis_0, axis_2 -> axis_2
-        """
-        return K.permute_dimensions(tensor, permutation, **kwargs)
diff --git a/n3fit/src/n3fit/backends/keras_backend/MetaModel.py b/n3fit/src/n3fit/backends/keras_backend/MetaModel.py
index 37f7003c39..254a46c548 100644
--- a/n3fit/src/n3fit/backends/keras_backend/MetaModel.py
+++ b/n3fit/src/n3fit/backends/keras_backend/MetaModel.py
@@ -8,6 +8,18 @@
 from tensorflow.keras.models import Model
 from tensorflow.keras import optimizers as Kopt
 
+# Define in this dictionary new optimizers as well as the arguments they accept
+# (with default values if needed be)
+optimizers = {
+    "RMSprop": (Kopt.RMSprop, {"lr": 0.01}),
+    "Adam": (Kopt.Adam, {"lr": 0.01}),
+    "Adagrad": (Kopt.Adagrad, {}),
+    "Adadelta": (Kopt.Adadelta, {"lr": 1.0}),
+    "Adamax": (Kopt.Adamax, {}),
+    "Nadam": (Kopt.Nadam, {}),
+    "Amsgrad": (Kopt.Adam, {"lr": 0.01, "amsgrad": True}),
+}
+
 
 def _fill_placeholders(original_input, new_input=None):
     """
@@ -57,18 +69,6 @@ class MetaModel(Model):
             keyword arguments to pass directly to Model
     """
 
-    # Define in this dictionary new optimizers as well as the arguments they accept
-    # (with default values if needed be)
-    optimizers = {
-        "RMSprop": (Kopt.RMSprop, {"lr": 0.01}),
-        "Adam": (Kopt.Adam, {"lr": 0.01}),
-        "Adagrad": (Kopt.Adagrad, {}),
-        "Adadelta": (Kopt.Adadelta, {"lr": 1.0}),
-        "Adamax": (Kopt.Adamax, {}),
-        "Nadam": (Kopt.Nadam, {}),
-        "Amsgrad": (Kopt.Adam, {"lr": 0.01, "amsgrad": True}),
-    }
-
     def __init__(self, input_tensors, output_tensors, **kwargs):
         self.has_dataset = False
 
@@ -227,7 +227,7 @@ def compile(
                 if given further calls to fit/evaluate must be done with y = None.
         """
         try:
-            opt_tuple = self.optimizers[optimizer_name]
+            opt_tuple = optimizers[optimizer_name]
         except KeyError as e:
             raise NotImplementedError(
                 f"[MetaModel.select_initializer] optimizer not implemented: {optimizer_name}"
diff --git a/n3fit/src/n3fit/backends/keras_backend/operations.py b/n3fit/src/n3fit/backends/keras_backend/operations.py
index 1c2e4996d8..1673e1ed32 100644
--- a/n3fit/src/n3fit/backends/keras_backend/operations.py
+++ b/n3fit/src/n3fit/backends/keras_backend/operations.py
@@ -1,12 +1,22 @@
 """
-    This module containg a list of useful operations translated in the keras language
-
-    All operations accept as input an iterable of keras layers or tensors
-    and (when necessary) keyword arguments.
-    The return operation is always a keras layer (or tensor)
+    This module contains the list of operations that can be used within the
+    ``call`` method of the ``n3fit`` layers as well as operations that can
+    act on layers.
 
     This includes an implementation of the NNPDF operations on fktable in the keras
-    language (hence the mapping `c_to_py_fun`)
+    language (with the mapping ``c_to_py_fun``) into Keras ``Lambda`` layers.
+
+    The rest of the operations in this module are divided into three categories:
+    numpy to tensor:
+        Operations that take a numpy array and return a tensorflow tensor
+    tensor to tensor:
+        Operations that take a tensor and return a tensor
+    layer generation:
+        Instanciate a layer to be applied by the calling function
+
+    Some of these are just aliases to the backend (tensorflow or Keras) operations
+    Note that tensor operations can also be applied to layers as the output of a layer is a tensor
+    equally operations are automatically converted to layers when used as such.
 """
 
 import numpy as np
@@ -15,12 +25,39 @@
 from tensorflow.keras.layers import multiply as keras_multiply
 from tensorflow.keras.layers import Concatenate as keras_concatenate
 
-from tensorflow.keras.layers import Input, Layer
+from tensorflow.keras.layers import Input
 from tensorflow.keras import backend as K
 
 from validphys.convolution import OP
 
 
+def evaluate(tensor):
+    """ Evaluate input tensor using the backend """
+    return K.eval(tensor)
+
+
+# NNPDF operations
+def c_to_py_fun(op_name, name="dataset"):
+    """
+    Map the NNPDF operations to Keras layers
+    NNPDF operations are defined in :py:func:`validphys.convolution.OP`
+
+    Parameters
+    ----------
+        op_name: str
+            A string defining the operation name
+    """
+    try:
+        operation = OP[op_name]
+    except KeyError as e:
+        raise ValueError(f"Operation {op_name} not recognised") from e
+
+    # Convert the operation into a lambda layer
+    operation_layer = keras_Lambda(lambda x: operation(*x), name=f"op_{name}_{op_name}")
+    return operation_layer
+
+
+# f(x: numpy) -> y: tensor
 def numpy_to_tensor(ival):
     """
         Make the input into a tensor
@@ -28,11 +65,13 @@ def numpy_to_tensor(ival):
     return K.constant(ival)
 
 
+# f(x: tensor) -> y: tensor
 def batchit(x, batch_dimension=0):
     """ Add a batch dimension to tensor x """
     return tf.expand_dims(x, batch_dimension)
 
 
+# layer generation
 def concatenate_split(splitting_sizes, axis=1):
     """ Generate a pair of concatention and splitting layer
     so that they invert each other
@@ -49,6 +88,7 @@ def concatenate_split(splitting_sizes, axis=1):
     return concatenation_layer, splitting_layer
 
 
+# layer generation
 def numpy_to_input(numpy_array, no_reshape=False, name=None):
     """
     Takes a numpy array and generates a Input layer.
@@ -77,33 +117,76 @@ def numpy_to_input(numpy_array, no_reshape=False, name=None):
     return input_layer
 
 
-def evaluate(tensor):
-    """ Evaluate input tensor using the backend """
-    return K.eval(tensor)
+#
+# Tensor operations
+# f(x: tensor[s]) -> y: tensor
+#
+
+# Generation operations
+# generate tensors of given shape/content
+def tensor_ones_like(*args, **kwargs):
+    """
+    Generates a tensor of ones of the same shape as the input tensor
+    See full `docs <https://www.tensorflow.org/api_docs/python/tf/keras/backend/ones_like>`_
+    """
+    return K.ones_like(*args, **kwargs)
 
 
-def c_to_py_fun(op_name, name = "dataset"):
+def many_replication(grid, replications, axis=0, **kwargs):
     """
-    Map the NNPDF operations to Keras layers
-    NNPDF operations are defined in :py:func:`validphys.convolution.OP
+    Generates a tensor with one extra dimension:
+        a repetition of "grid" n times along the given axis
+    from keras documentation:
+    If x has shape (s1, s2, s3) and axis is 1, the output will have shape (s1, s2 * rep, s3)
+    see full `docs <https://www.tensorflow.org/api_docs/python/tf/keras/backend/repeat_elements>`_
+    """
+    return K.repeat_elements(grid, rep=replications, axis=axis, **kwargs)
 
-    Parameters
-    ----------
-        op_name: str
-            A string defining the operation name
+
+# Property operations
+# modify properties of the tensor like the shape or elements it has
+def flatten(x):
+    """ Flatten tensor x """
+    return tf.reshape(x, (-1,))
+
+
+def boolean_mask(*args, **kwargs):
     """
-    try:
-        operation = OP[op_name]
-    except KeyError as e:
-        raise ValueError(f"Operation {op_name} not recognised") from e
+    Applies a boolean mask to a tensor
 
-    # Convert the operation into a lambda layer
-    operation_layer = keras_Lambda(lambda x: operation(*x), name=f"op_{name}_{op_name}")
-    return operation_layer
+    Relevant parameters: (tensor, mask, axis=None)
+    see full `docs <https://www.tensorflow.org/api_docs/python/tf/boolean_mask>`_.
+    """
+    return tf.boolean_mask(*args, **kwargs)
 
-def op_subtract(o_list): #TODO to be removed, not used once other PRs are merged
-    from tensorflow.keras.layers import subtract
-    return subtract(o_list)
+
+def transpose(tensor, **kwargs):
+    """
+    Transpose a layer,
+    see full `docs <https://www.tensorflow.org/api_docs/python/tf/keras/backend/transpose>`_
+    """
+    return K.transpose(tensor, **kwargs)
+
+
+def concatenate(tensor_list, axis=-1, target_shape=None):
+    """
+    Concatenates a list of numbers or tenosr into a bigger tensor
+    If the target shape is given, the output is reshaped to said shape
+    """
+    concatenated_tensor = K.concatenate(tensor_list, axis=axis)
+    if target_shape:
+        return K.reshape(concatenated_tensor, target_shape)
+    else:
+        return concatenated_tensor
+
+
+# Mathematical operations
+def tensor_product(*args, **kwargs):
+    """
+    Computes the tensordot product between tensor_x and tensor_y
+    See full `docs <https://www.tensorflow.org/api_docs/python/tf/tensordot>`_
+    """
+    return tf.tensordot(*args, **kwargs)
 
 
 def op_multiply(o_list, **kwargs):
@@ -134,3 +217,11 @@ def op_log(o_tensor, **kwargs):
     Computes the logarithm of the input
     """
     return K.log(o_tensor)
+
+
+def sum(*args, **kwargs):
+    """
+    Computes the sum of the elements of the tensor
+    see full `docs <https://www.tensorflow.org/api_docs/python/tf/keras/backend/sum>`_
+    """
+    return K.sum(*args, **kwargs)
diff --git a/n3fit/src/n3fit/layers/DIS.py b/n3fit/src/n3fit/layers/DIS.py
index c800c6b09e..50710d8694 100644
--- a/n3fit/src/n3fit/layers/DIS.py
+++ b/n3fit/src/n3fit/layers/DIS.py
@@ -1,5 +1,6 @@
 import numpy as np
 from n3fit.layers.Observable import Observable
+from n3fit.backends import operations as op
 
 
 class DIS(Observable):
@@ -36,8 +37,8 @@ def call(self, pdf_in):
                 - `result`: rank 1 tensor (ndata)
         """
         pdf_in = self.digest_pdf(pdf_in)
-        pdf_masked = self.boolean_mask(pdf_in, self.basis, axis=1)
+        pdf_masked = op.boolean_mask(pdf_in, self.basis, axis=1)
 
-        pdfT = self.transpose(pdf_masked)
-        result = self.tensor_product(self.fktable, pdfT, axes=2)
+        pdfT = op.transpose(pdf_masked)
+        result = op.tensor_product(self.fktable, pdfT, axes=2)
         return result
diff --git a/n3fit/src/n3fit/layers/DY.py b/n3fit/src/n3fit/layers/DY.py
index f1a2b0cc01..c1166d07bd 100644
--- a/n3fit/src/n3fit/layers/DY.py
+++ b/n3fit/src/n3fit/layers/DY.py
@@ -1,5 +1,6 @@
 import numpy as np
 from n3fit.layers.Observable import Observable
+from n3fit.backends import operations as op
 
 
 class DY(Observable):
@@ -39,31 +40,31 @@ def call(self, pdf_in):
         # This is a convoluted way of applying a mask, but it is faster
         # mask-version below
         lumi_fun = []
-        pdfT = self.transpose(pdf_in)
+        pdfT = op.transpose(pdf_in)
 
         for i, j in self.basis:
-            lumi_fun.append(self.tensor_product(pdfT[i], pdfT[j], axes=0))
+            lumi_fun.append(op.tensor_product(pdfT[i], pdfT[j], axes=0))
 
-        pdf_X_pdf = self.concatenate(lumi_fun, axis=0, target_shape=(self.basis_size, self.xgrid_size, self.xgrid_size))
+        pdf_X_pdf = op.concatenate(lumi_fun, axis=0, target_shape=(self.basis_size, self.xgrid_size, self.xgrid_size))
 
-        result = self.tensor_product(self.fktable, pdf_X_pdf, axes=3)
+        result = op.tensor_product(self.fktable, pdf_X_pdf, axes=3)
         return result
 
 
 # Another example on how to performt the DY convolution
 # this code is equivalent to the previos one, with a slightly greater cost
-class DY_mask(Observable):
-    def gen_basis(self, basis):
-        if basis is not None:
-            self.basis = np.zeros((self.nfl, self.nfl), dtype=bool)
-            for i, j in basis.reshape(-1, 2):
-                self.basis[i, j] = True
-        else:
-            self.basis = np.ones((self.nfl, self.nfl), dtype=bool)
-
-    def call(self, pdf_in):
-        lfun = self.tensor_product(pdf_in, pdf_in, axes=0)
-        lfunT = self.permute_dimensions(lfun, (3, 1, 2, 0))
-        x = self.boolean_mask(lfunT, self.basis, axis=0)
-        result = self.tensor_product(self.fktable, x, axes=3)
-        return result
+# class DY_mask(Observable):
+#     def gen_basis(self, basis):
+#         if basis is not None:
+#             self.basis = np.zeros((self.nfl, self.nfl), dtype=bool)
+#             for i, j in basis.reshape(-1, 2):
+#                 self.basis[i, j] = True
+#         else:
+#             self.basis = np.ones((self.nfl, self.nfl), dtype=bool)
+# 
+#     def call(self, pdf_in):
+#         lfun = op.tensor_product(pdf_in, pdf_in, axes=0)
+#         lfunT = tensorflow.keras.backend.permute_dimensions(lfun, (3, 1, 2, 0))
+#         x = op.boolean_mask(lfunT, self.basis, axis=0)
+#         result = op.tensor_product(self.fktable, x, axes=3)
+#         return result
diff --git a/n3fit/src/n3fit/layers/MSR_Normalization.py b/n3fit/src/n3fit/layers/MSR_Normalization.py
index 49f3bfec45..1b68948ea4 100644
--- a/n3fit/src/n3fit/layers/MSR_Normalization.py
+++ b/n3fit/src/n3fit/layers/MSR_Normalization.py
@@ -1,4 +1,5 @@
 from n3fit.backends import MetaLayer
+from n3fit.backends import operations as op
 
 
 class MSR_Normalization(MetaLayer):
@@ -8,8 +9,8 @@ class MSR_Normalization(MetaLayer):
 
     def __init__(self, output_dim=14, **kwargs):
         self.output_dim = output_dim
-        self.one = self.tensor_ones((1, 1))
-        self.three = 3 * self.tensor_ones((1, 1))
+        self.one = op.numpy_to_tensor([[1.0]])
+        self.three = op.numpy_to_tensor([[3.0]])
         super(MSR_Normalization, self).__init__(**kwargs, name="normalizer")
 
     def call(self, xgrid):
@@ -17,8 +18,8 @@ def call(self, xgrid):
             Receives as input a tensor with the value of the MSR for each PDF
             and returns a rank-1 tensor with the normalization factor A_i of each flavour
         """
-        x = self.flatten(xgrid)
-        pdf_sr = self.concatenate(
+        x = op.flatten(xgrid)
+        pdf_sr = op.concatenate(
             [
                 self.one,  # photon
                 self.one,  # sigma
diff --git a/n3fit/src/n3fit/layers/Mask.py b/n3fit/src/n3fit/layers/Mask.py
index 65a9f4145d..f33b13c06c 100644
--- a/n3fit/src/n3fit/layers/Mask.py
+++ b/n3fit/src/n3fit/layers/Mask.py
@@ -1,6 +1,6 @@
 import numpy as np
 from n3fit.backends import MetaLayer
-from n3fit.backends import operations
+from n3fit.backends import operations as op
 
 
 class Mask(MetaLayer):
@@ -33,8 +33,8 @@ def build(self, input_shape):
         super(Mask, self).build(input_shape)
 
     def call(self, prediction_in):
-        ret = self.boolean_mask(self.kernel * prediction_in, self.mask)
+        ret = op.boolean_mask(self.kernel * prediction_in, self.mask)
         if self.batch_it:
-            return operations.batchit(ret)
+            return op.batchit(ret)
         else:
             return ret
diff --git a/n3fit/src/n3fit/layers/Observable.py b/n3fit/src/n3fit/layers/Observable.py
index 6958bada2b..375663585d 100644
--- a/n3fit/src/n3fit/layers/Observable.py
+++ b/n3fit/src/n3fit/layers/Observable.py
@@ -1,6 +1,7 @@
 from n3fit.backends import MetaLayer
 import tensorflow as tf
 from abc import abstractmethod, ABC
+from n3fit.backends import operations as op
 
 
 class Observable(MetaLayer, ABC):
@@ -26,7 +27,7 @@ def __init__(self, output_dim, fktable, basis=None, nfl=14, **kwargs):
         self.nfl = nfl
 
         self.output_dim = output_dim
-        self.fktable = self.np_to_tensor(fktable)
+        self.fktable = op.numpy_to_tensor(fktable)
         self.xgrid_size = self.fktable.shape[-1]
 
         self.gen_basis(basis)
diff --git a/n3fit/src/n3fit/layers/Obsrotation.py b/n3fit/src/n3fit/layers/Obsrotation.py
index 4d2b996f12..8c3772773a 100644
--- a/n3fit/src/n3fit/layers/Obsrotation.py
+++ b/n3fit/src/n3fit/layers/Obsrotation.py
@@ -1,4 +1,5 @@
 from n3fit.backends import MetaLayer
+from n3fit.backends import operations as op
 
 class ObsRotation(MetaLayer):
     """
@@ -9,9 +10,9 @@ class ObsRotation(MetaLayer):
     """
 
     def __init__(self, transform_matrix, **kwargs):
-        self.rotation = self.np_to_tensor(transform_matrix)
+        self.rotation = op.numpy_to_tensor(transform_matrix)
         super(MetaLayer, self).__init__(**kwargs)
 
     def call(self, prediction_in):
-        pinT = self.transpose(prediction_in)
-        return self.tensor_product(self.rotation, pinT, axes=1)
+        pinT = op.transpose(prediction_in)
+        return op.tensor_product(self.rotation, pinT, axes=1)
diff --git a/n3fit/src/n3fit/layers/Preprocessing.py b/n3fit/src/n3fit/layers/Preprocessing.py
index 6a3bb29dc7..9833fa1697 100644
--- a/n3fit/src/n3fit/layers/Preprocessing.py
+++ b/n3fit/src/n3fit/layers/Preprocessing.py
@@ -1,5 +1,6 @@
 from n3fit.backends import MetaLayer
 from n3fit.backends import constraints
+from n3fit.backends import operations as op
 
 BASIS_SIZE = 8
 
@@ -99,7 +100,7 @@ def build(self, input_shape):
 
     def call(self, inputs, **kwargs):
         x = inputs
-        pdf_raw = self.concatenate(
+        pdf_raw = op.concatenate(
             [
                 x ** (1 - self.kernel[0][0]) * (1 - x) ** self.kernel[1][0],  # sigma
                 x ** (1 - self.kernel[2][0]) * (1 - x) ** self.kernel[3][0],  # g
diff --git a/n3fit/src/n3fit/layers/Rotation.py b/n3fit/src/n3fit/layers/Rotation.py
index c5f071606b..65b7b5f8af 100644
--- a/n3fit/src/n3fit/layers/Rotation.py
+++ b/n3fit/src/n3fit/layers/Rotation.py
@@ -3,7 +3,7 @@
 """
 
 from n3fit.backends import MetaLayer
-from n3fit.backends import operations
+from n3fit.backends import operations as op
 
 
 class Rotation(MetaLayer):
@@ -22,7 +22,7 @@ def __init__(self, output_dim=14, **kwargs):
 
     def call(self, pdf_raw):
         # Transpose the PDF so that the flavour index is the first one
-        x = self.transpose(pdf_raw)
+        x = op.transpose(pdf_raw)
         pdf_raw_list = [
             0 * x[0],  # photon
             x[0],  # sigma
@@ -39,6 +39,6 @@ def call(self, pdf_raw):
             x[0],  # t24
             x[0],  # t35
         ]
-        ret = self.concatenate(pdf_raw_list)
+        ret = op.concatenate(pdf_raw_list)
         # Concatenating destroys the batch index so we have to regenerate it
-        return operations.batchit(ret)
+        return op.batchit(ret)
diff --git a/n3fit/src/n3fit/layers/x_operations.py b/n3fit/src/n3fit/layers/x_operations.py
index 3bc83cd099..a2197bb5ab 100644
--- a/n3fit/src/n3fit/layers/x_operations.py
+++ b/n3fit/src/n3fit/layers/x_operations.py
@@ -12,6 +12,7 @@
 """
 
 from n3fit.backends import MetaLayer
+from n3fit.backends import operations as op
 
 
 class xDivide(MetaLayer):
@@ -36,14 +37,14 @@ def __init__(self, output_dim=8, div_list=None, **kwargs):
 
     def call(self, x):
         out_array = []
-        one = self.tensor_ones_like(x)
+        one = op.tensor_ones_like(x)
         for i in range(self.output_dim):
             if i in self.div_list:
                 res = one / x
             else:
                 res = one
             out_array.append(res)
-        out_tensor = self.concatenate(out_array, axis=-1)
+        out_tensor = op.concatenate(out_array, axis=-1)
         return out_tensor
 
 
@@ -62,11 +63,11 @@ class xIntegrator(MetaLayer):
     """
 
     def __init__(self, grid_weights, **kwargs):
-        grid_weights_tensor = self.np_to_tensor(grid_weights)
+        grid_weights_tensor = op.numpy_to_tensor(grid_weights)
         # Open up the grid weights
-        self.grid_weights = self.many_replication(grid_weights_tensor, replications=8, axis=1)
+        self.grid_weights = op.many_replication(grid_weights_tensor, 8, axis=1)
         super(MetaLayer, self).__init__(**kwargs)
 
     def call(self, x):
         xx = x * self.grid_weights
-        return self.sum(xx, axis=-2)
+        return op.sum(xx, axis=-2)
diff --git a/n3fit/src/n3fit/tests/test_backend.py b/n3fit/src/n3fit/tests/test_backend.py
index db0b1713a8..2693fd435a 100644
--- a/n3fit/src/n3fit/tests/test_backend.py
+++ b/n3fit/src/n3fit/tests/test_backend.py
@@ -3,6 +3,7 @@
     and ensures they do the same thing as their numpy counterparts
 """
 import operator
+import functools
 import numpy as np
 from n3fit.backends import operations as op
 from n3fit.backends import losses
@@ -40,6 +41,8 @@ def numpy_check(backend_op, python_op, mode="same"):
      - `same` (default): two arrays of the same dimensionality
      - `diff`: first array has one extra dimension that second array
      - `single`: only one array enters the operation
+     - (tensor, array): if passed a tuple (backend tensor, numpy array), uses these
+        values as tensor and array inputs for the operations
     """
     if mode == "same":
         tensors = [T1, T2]
@@ -53,6 +56,10 @@ def numpy_check(backend_op, python_op, mode="same"):
     elif mode == "single":
         tensors = [T1]
         arrays = [ARR1]
+    elif isinstance(mode, tuple):
+        tensors = mode[0]
+        arrays = mode[1]
+
     result = backend_op(tensors)
     reference = python_op(*arrays)
     are_equal(result, reference)
@@ -89,7 +96,27 @@ def test_op_multiply_dim():
 
 
 def test_op_log():
-    numpy_check(op.op_log, np.log, mode="single")
+    numpy_check(op.op_log, np.log, mode='single')
+
+
+def test_flatten():
+    numpy_check(op.flatten, np.ndarray.flatten, mode=(T3, [ARR3]))
+
+
+def test_boolean_mask():
+    bools = np.random.randint(0, 2, DIM, dtype=bool)
+    np_result = ARR1[bools]
+    tf_bools = op.numpy_to_tensor(bools)
+    tf_result = op.boolean_mask(T1, tf_bools, axis=0)
+    are_equal(np_result, tf_result)
+
+def test_tensor_product():
+    np_result = np.tensordot(ARR3, ARR1, axes=1)
+    tf_result = op.tensor_product(T3, T1, axes=1)
+    are_equal(np_result, tf_result)
+
+def test_sum():
+    numpy_check(op.sum, np.sum, mode='single')
 
 
 # Tests loss functions