diff --git a/app/function_approximation.py b/app/function_approximation.py
new file mode 100644
index 0000000000000000000000000000000000000000..d2cd790ac9e5759b79b754fe17c837109579ab04
--- /dev/null
+++ b/app/function_approximation.py
@@ -0,0 +1,21 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import torch
+from src.misc import time_stamp, timeit
+from src import target_function
+
+
+def main() -> None:
+    print(time_stamp(), "Initialize main file")
+    with timeit("create x_train data ({:4.2f} s)"):
+        x_train = np.random.uniform(0, 1, (100000, 2))
+    with timeit("create y_train data ({:4.2f} s)"):
+        y_train = target_function.sin2d(x_train)
+
+    plt.figure()
+    plt.hexbin(x_train[:, 0], x_train[:, 1], y_train, gridsize=50)
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/misc.py b/src/misc.py
new file mode 100644
index 0000000000000000000000000000000000000000..5e1f8a2b3e56aa8c66a967fd79f3e779557542e9
--- /dev/null
+++ b/src/misc.py
@@ -0,0 +1,45 @@
+from typing import List, Union
+import time
+from contextlib import contextmanager
+import numpy as np
+
+
+def time_stamp():
+    """Time stamp for logging."""
+    return time.strftime("[%Y-%m-%d %H:%M:%S]", time.localtime())
+
+
+@contextmanager
+def timeit(title: str, verbose: bool = True) -> None:
+    """Time functions and print time."""
+    start = time.perf_counter()
+    try:
+        yield
+    finally:
+        if verbose:
+            print(time_stamp(), title.format(time.perf_counter()-start))
+
+
+def cart_prod(array_list: Union[List[np.ndarray], np.ndarray]) -> np.ndarray:
+    """Compute the outer product of two or more arrays.
+
+    Assemble an array containing all possible combinations of the elements
+    of the input vectors :math:`v_1,\\dots,v_n`.
+
+    Parameters
+    ----------
+    array_list : list of array_like
+        List of vectors :math:`v_1,\\dots,v_n`.
+
+    Returns
+    -------
+    :
+        Cartesian product array.
+    """
+    dim = len(array_list)
+    if dim == 1:
+        return np.array(array_list).T
+    x = np.hstack((np.meshgrid(*array_list))).swapaxes(0, 1).reshape(dim, -1).T
+    return x
+
+
diff --git a/src/target_function.py b/src/target_function.py
new file mode 100644
index 0000000000000000000000000000000000000000..37c9851b0a287c81bdcda510c1ac75a5fe35bd3c
--- /dev/null
+++ b/src/target_function.py
@@ -0,0 +1,23 @@
+import numpy as np
+
+
+def sin2d(xs: np.ndarray) -> np.ndarray:
+    """Benchmark function for approximation.
+
+    The benchmark function is given by
+
+    .. math::
+        f(x_1, x_2) = \\sin(4\\pi x_1) \\sin(3\\pi x_2)
+
+    Parameters
+    ----------
+    xs : np.ndarray
+        Evaluation points. Shape is (n, 2).
+
+    Returns
+    -------
+    :
+        Function evaluations.
+    """
+    assert xs.ndim == 2 and xs.shape[1] == 2
+    return np.sin(4*np.pi*xs[:, 0]) * np.sin(3*np.pi*xs[:, 1])