diff --git a/README.rst b/README.rst
index 97d3920d..d909b678 100644
--- a/README.rst
+++ b/README.rst
@@ -34,20 +34,18 @@ you are using Python 3.10 and below, you can use ``cpctools`` to access
 
   $ pip install cpctools
 
-If you want to use the ``dynsight.tica`` module you will need to install the
-deeptime package. This can be done with with pip::
-
-  $ pip install deeptime
-
-or with conda::
-
-  $ conda install -c conda-forge deeptime
-
 If you want to use the ``dynsight.vision`` and ``dynsight.track`` modules
 you will need to install a series of packages. This can be done with with pip::
 
   $ pip install ultralytics PyYAML
 
+How to get started
+------------------
+
+We suggest you give a read to the ``dynsight.trajectory`` module documentation,
+which offers a compact and easy way of using most of the ``dynsight`` tools. 
+Also, the documentation offers some copiable Recipes and Examples for the most
+common analyses.
 
 How to contribute
 -----------------
diff --git a/docs/source/_static/info_plot.png b/docs/source/_static/info_plot.png
new file mode 100644
index 00000000..96c9d0af
Binary files /dev/null and b/docs/source/_static/info_plot.png differ
diff --git a/docs/source/_static/recipes/__init__.py b/docs/source/_static/recipes/__init__.py
new file mode 100644
index 00000000..4121ce73
--- /dev/null
+++ b/docs/source/_static/recipes/__init__.py
@@ -0,0 +1 @@
+"""Copyable code for the recipes."""
diff --git a/docs/source/_static/recipes/descr_from_trj.py b/docs/source/_static/recipes/descr_from_trj.py
new file mode 100644
index 00000000..01cba368
--- /dev/null
+++ b/docs/source/_static/recipes/descr_from_trj.py
@@ -0,0 +1,50 @@
+"""Copiable code from Recipe #1."""
+
+from pathlib import Path
+
+from dynsight.trajectory import Trj
+
+
+def main() -> None:
+    """Code from the Recipe #1."""
+    # Loading an example trajectory
+    files_path = Path("../tests/systems/")
+    trj = Trj.init_from_xtc(
+        traj_file=files_path / "balls_7_nvt.xtc",
+        topo_file=files_path / "balls_7_nvt.gro",
+    )
+
+    _ = trj.get_soap(
+        r_cut=2.0,  # cutoff radius for neighbors list
+        n_max=4,  # n_max SOAP parameter
+        l_max=4,  # l_max SOAP parameter
+        selection="all",  # compute on a selection of particles
+        centers="all",  # compute for a selection of centers
+        respect_pbc=False,  # consider PBC
+        n_core=1,  # use multiprocessing to speed up calculations
+    )
+
+    # Loading an example trajectory
+    files_path = Path("../tests/systems/")
+    trj = Trj.init_from_xtc(
+        traj_file=files_path / "balls_7_nvt.xtc",
+        topo_file=files_path / "balls_7_nvt.gro",
+    )
+
+    # Computing number of neighbors from scratch
+    neigcounts, n_neig = trj.get_coord_number(
+        r_cut=2.0,  # cutoff radius for neighbors list
+        selection="all",  # compute on a selection of particles
+        neigcounts=None,  # it will be computed and returned
+    )
+
+    # Now for LENS we already have neigcounts
+    _, lens = trj.get_lens(
+        r_cut=2.0,  # cutoff radius for neighbors list
+        selection="all",  # compute on a selection of particles
+        neigcounts=neigcounts,  # no need to compute it again
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/docs/source/_static/recipes/info_gain.py b/docs/source/_static/recipes/info_gain.py
new file mode 100644
index 00000000..873839a0
--- /dev/null
+++ b/docs/source/_static/recipes/info_gain.py
@@ -0,0 +1,179 @@
+"""Copiable code from Recipe #3."""
+
+from __future__ import annotations
+
+from pathlib import Path
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from numpy.typing import NDArray
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import dynsight
+
+
+def info_gain_with_onion(
+    delta_t_list: NDArray[np.int64] | list[int],
+    data: NDArray[np.float64],
+    n_bins: int = 40,
+) -> tuple[
+    NDArray[np.int64],
+    NDArray[np.float64],
+    NDArray[np.float64],
+    NDArray[np.float64],
+    float,
+]:
+    """Performs full information gain analysis with Onion clustering."""
+    data_range = (np.min(data), np.max(data))
+
+    n_clusters = np.zeros(len(delta_t_list), dtype=int)
+    clusters_frac = []
+    info_gain = np.zeros(len(delta_t_list))
+    clusters_entr = []
+
+    for i, delta_t in enumerate(delta_t_list):
+        state_list, labels = dynsight.onion.onion_uni_smooth(
+            data,
+            delta_t=delta_t,
+        )
+
+        n_clusters[i] = len(state_list)
+        tmp_frac = [0.0] + [state.perc for state in state_list]
+        tmp_frac[0] = 1.0 - np.sum(tmp_frac)
+        clusters_frac.append(tmp_frac)
+
+        flat_data = data.flatten()
+        flat_labels = labels.flatten()
+        info_gain[i], _, h_0, _ = dynsight.analysis.compute_entropy_gain(
+            flat_data,
+            flat_labels,
+            n_bins=n_bins,
+        )
+
+        tmp_entr = []
+        label_list = np.unique(labels)
+        if label_list[0] != -1:
+            tmp_entr.append(-1.0)
+
+        for _, lab in enumerate(label_list):
+            mask = labels == lab
+            selected_points = data[mask]
+            tmp_entr.append(
+                dynsight.analysis.compute_shannon(
+                    selected_points,
+                    data_range,
+                    n_bins=n_bins,
+                )
+            )
+        clusters_entr.append(tmp_entr)
+
+    max_n_envs = np.max([len(elem) for elem in clusters_entr])
+    for i, elem in enumerate(clusters_entr):
+        while len(elem) < max_n_envs:
+            elem.append(-1.0)
+            clusters_frac[i].append(0.0)
+
+    cl_frac = np.array(clusters_frac, dtype=float)
+    cl_entr = np.array(clusters_entr)
+
+    return n_clusters, cl_frac, info_gain, cl_entr, h_0
+
+
+def plot_info_results(
+    delta_t_list: NDArray[np.int64] | list[int],
+    cl_frac: NDArray[np.float64],
+    cl_entr: NDArray[np.float64],
+    h_0: float,
+    file_path: Path,
+) -> None:
+    """Plot information gain as a function of ∆t."""
+    frac = cl_frac.T
+    entr = cl_entr.T
+    s_list = []
+    for i, st_fr in enumerate(frac):
+        s_list.append(st_fr * entr[i])
+    s_cumul = [s_list[0]]
+    for _, tmp_s in enumerate(s_list[1:]):
+        s_cumul.append(s_cumul[-1] + tmp_s)
+
+    fig, ax = plt.subplots()
+
+    i_0 = (1 - h_0) * np.ones(len(delta_t_list))
+    ax.plot(delta_t_list, i_0, ls="--", c="black", marker="")  # I_0
+    ax.fill_between(
+        delta_t_list,
+        1,
+        1 - s_cumul[0],
+        alpha=0.5,
+    )
+    for i, tmp_s in enumerate(s_cumul[1:]):
+        ax.fill_between(
+            delta_t_list,
+            1 - s_cumul[i],
+            1 - tmp_s,
+            alpha=0.5,
+        )
+    ax.fill_between(
+        delta_t_list,
+        1 - s_cumul[-1],
+        1 - h_0,
+        color="gainsboro",
+    )
+    ax.plot(
+        delta_t_list,
+        1 - s_cumul[-1],
+        c="black",
+        marker="",
+    )  # I_clust
+
+    ax.set_ylim(0.0, 1.0)
+    ax.set_xlabel(r"Time resolution $\Delta t$")
+    ax.set_ylabel(r"Information $I$")
+    ax.set_xscale("log")
+
+    fig.savefig(file_path, dpi=600)
+    plt.close()
+
+
+def main() -> None:
+    """Copiable code from Recipe #3."""
+    rng = np.random.default_rng(1234)
+
+    n_atoms = 10
+    num_blocks = 10
+    block_size = 100
+    sigma = 0.1
+
+    # Generate the array
+    tmp_data = [
+        np.concatenate(
+            [
+                rng.normal(loc=(i % 2), scale=sigma, size=block_size)
+                for i in range(num_blocks)
+            ]
+        )
+        for _ in range(n_atoms)
+    ]
+    data = np.array(tmp_data)
+
+    _, n_frames = data.shape
+    delta_t_list = np.unique(np.geomspace(2, n_frames, 10, dtype=int))
+
+    n_cl, cl_frac, info_gain, cl_entr, h_0 = info_gain_with_onion(
+        delta_t_list,
+        data,
+    )
+
+    plot_info_results(
+        delta_t_list,
+        cl_frac,
+        cl_entr,
+        h_0,
+        Path("./source/_static/info_plot.png"),
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/docs/source/_static/recipes/soap_dim_red.py b/docs/source/_static/recipes/soap_dim_red.py
new file mode 100644
index 00000000..b93d6314
--- /dev/null
+++ b/docs/source/_static/recipes/soap_dim_red.py
@@ -0,0 +1,179 @@
+"""Copiable code from Recipe #2."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+from sklearn.decomposition import PCA
+
+import dynsight
+from dynsight.trajectory import Insight, Trj
+from dynsight.utilities import load_or_compute_soap
+
+
+def compute_soap_pca(
+    trj: Trj,
+    r_cut: float,
+    n_max: int,
+    l_max: int,
+    n_components: int,
+    soap_path: Path | None = None,
+    pca_path: Path | None = None,
+    selection: str = "all",
+    centers: str = "all",
+    respect_pbc: bool = True,
+    n_core: int = 1,
+) -> Insight:
+    """Computes the PCA of SOAP on a Trajectory."""
+    if pca_path is not None and pca_path.exists():
+        return Insight.load_from_json(pca_path)
+
+    soap = load_or_compute_soap(
+        trj=trj,
+        r_cut=r_cut,
+        n_max=n_max,
+        l_max=l_max,
+        selection=selection,
+        centers=centers,
+        respect_pbc=respect_pbc,
+        n_core=n_core,
+        soap_path=soap_path,
+    )
+
+    n_atom, n_frames, n_dims = soap.dataset.shape
+    reshaped_soap = soap.dataset.reshape(n_atom * n_frames, n_dims)
+    pca = PCA(n_components=n_components)
+    transformed_soap = pca.fit_transform(reshaped_soap)
+    pca_ds = transformed_soap.reshape(n_atom, n_frames, -1)
+
+    soap_pca = Insight(pca_ds, meta=soap.meta.copy())
+
+    if pca_path is not None:
+        soap_pca.dump_to_json(pca_path)
+
+    return soap_pca
+
+
+def compute_soap_tica(
+    trj: Trj,
+    r_cut: float,
+    n_max: int,
+    l_max: int,
+    lag_time: int,
+    tica_dim: int,
+    soap_path: Path | None = None,
+    tica_path: Path | None = None,
+    selection: str = "all",
+    centers: str = "all",
+    respect_pbc: bool = True,
+    n_core: int = 1,
+) -> Insight:
+    """Computes the tICA of SOAP on a Trajectory."""
+    if tica_path is not None and tica_path.exists():
+        return Insight.load_from_json(tica_path)
+
+    soap = load_or_compute_soap(
+        trj=trj,
+        r_cut=r_cut,
+        n_max=n_max,
+        l_max=l_max,
+        selection=selection,
+        centers=centers,
+        respect_pbc=respect_pbc,
+        n_core=n_core,
+        soap_path=soap_path,
+    )
+
+    rel_times, _, tica_ds = dynsight.tica.many_body_tica(
+        soap.dataset,
+        lag_time=lag_time,
+        tica_dim=tica_dim,
+    )
+
+    meta = soap.meta.copy()
+    meta.update(
+        {
+            "lag_time": lag_time,
+            "rel_times": rel_times,
+        }
+    )
+    soap_tica = Insight(tica_ds, meta=meta)
+
+    if tica_path is not None:
+        soap_tica.dump_to_json(tica_path)
+
+    return soap_tica
+
+
+def compute_timesoap(
+    trj: Trj,
+    r_cut: float,
+    n_max: int,
+    l_max: int,
+    delay: int = 1,
+    soap_path: Path | None = None,
+    tsoap_path: Path | None = None,
+    selection: str = "all",
+    centers: str = "all",
+    respect_pbc: bool = True,
+    n_core: int = 1,
+) -> Insight:
+    """Computes timeSOAP on a Trajectory."""
+    if tsoap_path is not None and tsoap_path.exists():
+        return Insight.load_from_json(tsoap_path)
+
+    soap = load_or_compute_soap(
+        trj=trj,
+        r_cut=r_cut,
+        n_max=n_max,
+        l_max=l_max,
+        selection=selection,
+        centers=centers,
+        respect_pbc=respect_pbc,
+        n_core=n_core,
+        soap_path=soap_path,
+    )
+
+    tsoap = soap.get_angular_velocity(delay=delay)
+
+    if tsoap_path is not None:
+        tsoap.dump_to_json(tsoap_path)
+
+    return tsoap
+
+
+def main() -> None:
+    """Copiable code from Recipe #2."""
+    files_path = Path("../tests/systems/")
+    trj = Trj.init_from_xtc(
+        traj_file=files_path / "balls_7_nvt.xtc",
+        topo_file=files_path / "balls_7_nvt.gro",
+    )
+
+    _ = compute_soap_pca(
+        trj=trj,
+        r_cut=2.0,
+        n_max=4,
+        l_max=4,
+        n_components=1,
+    )
+
+    _ = compute_soap_tica(
+        trj=trj,
+        r_cut=10.0,
+        n_max=4,
+        l_max=4,
+        lag_time=10,
+        tica_dim=1,
+    )
+
+    _ = compute_timesoap(
+        trj=trj,
+        r_cut=10.0,
+        n_max=4,
+        l_max=4,
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/docs/source/descr_from_trj.rst b/docs/source/descr_from_trj.rst
new file mode 100644
index 00000000..b7d21716
--- /dev/null
+++ b/docs/source/descr_from_trj.rst
@@ -0,0 +1,104 @@
+Descriptors from a :class:`.trajectory.Trj` 
+===========================================
+
+This recipe explains how to compute descriptors directly from a 
+:class:`.trajectory.Trj` object. 
+
+.. warning::
+
+    This code works when run from the ``/docs`` directory of the ``dynsight`` 
+    repo. To use it elsewhere, you have to change the ``Path`` variables
+    accordingly.
+
+First of all, we import all the packages and objects we'll need:
+
+.. testcode:: recipe1-test
+
+    from pathlib import Path
+    from dynsight.trajectory import Trj
+
+SOAP
+----
+
+Computing SOAP for every particle and frame of the trajectory is easy, since
+it's directly calculated by the :class:`.trajectory.Trj.get_soap()` method.
+
+.. warning::
+
+    Please consider that the SOAP dataset can be very large, due to the high
+    dimensionality, thus calculations can be expensive, and saving to/loading 
+    from file quite slow.
+
+.. testcode:: recipe1-test
+
+    # Loading an example trajectory
+    files_path = Path("../tests/systems/")
+    trj = Trj.init_from_xtc(
+        traj_file=files_path / "balls_7_nvt.xtc",
+        topo_file=files_path / "balls_7_nvt.gro",
+    )
+
+    soap = trj.get_soap(
+        r_cut=2.0,          # cutoff radius for neighbors list
+        n_max=4,            # n_max SOAP parameter
+        l_max=4,            # l_max SOAP parameter
+        selection="all",    # compute on a selection of particles
+        centers="all",      # compute for a selection of centers
+        respect_pbc=False,  # consider PBC
+        n_core=1,           # use multiprocessing to speed up calculations
+    )
+
+Number of neighbors and LENS
+----------------------------
+
+Similarly to SOAP, computing number of neighbors or LENS can be done with the
+respective methods of the :class:`.trajectory.Trj` class.
+
+Since both calculations need to compute the list of neighbors for each
+particle at each frame, this list is also returned, and can be passed as an
+optional parameter, so that when computing both quantities the second
+calculation can be sped up significantly.
+
+.. testcode:: recipe1-test
+
+    # Loading an example trajectory
+    files_path = Path("../tests/systems/")
+    trj = Trj.init_from_xtc(
+        traj_file=files_path / "balls_7_nvt.xtc",
+        topo_file=files_path / "balls_7_nvt.gro",
+    )
+
+    # Computing number of neighbors from scratch
+    neigcounts, n_neig = trj.get_coord_number(
+        r_cut=2.0,          # cutoff radius for neighbors list
+        selection="all",    # compute on a selection of particles
+        neigcounts=None,    # it will be computed and returned
+    )
+
+    # Now for LENS we already have neigcounts
+    _, lens = trj.get_lens(
+        r_cut=2.0,               # cutoff radius for neighbors list
+        selection="all",         # compute on a selection of particles
+        neigcounts=neigcounts,   # no need to compute it again
+    )
+
+Notice that, differently from SOAP - which is computed for every frame, LENS
+is computed for every pair of frames. Thus, the LENS dataset has shape 
+``(n_particles, n_frames - 1)``. Consequently, if you need to match the LENS
+values with the particles along the trajectory, you will need to use a sliced
+trajectory (removing the last frame). The easiest way to do this is:
+
+.. testcode:: recipe1-test
+
+    trajslice = slice(0, -1, 1)
+    shorter_trj = trj.with_slice(trajslice=trajslice)
+
+.. raw:: html
+
+    <a class="btn-download" href="../_static/recipes/descr_from_trj.py" download>⬇️ Download Python Script</a>
+
+.. testcode:: recipe1-test
+    :hide:
+
+    assert soap.dataset.shape == (7, 201, 50)
+    assert lens.dataset.shape == (7, 200)
diff --git a/docs/source/analysis_workflow.rst b/docs/source/example_analysis_workflow.rst
similarity index 100%
rename from docs/source/analysis_workflow.rst
rename to docs/source/example_analysis_workflow.rst
diff --git a/docs/source/example_info_gain.rst b/docs/source/example_info_gain.rst
new file mode 100644
index 00000000..884f9c6d
--- /dev/null
+++ b/docs/source/example_info_gain.rst
@@ -0,0 +1,257 @@
+Information gain computation
+============================
+
+For the theoretical aspects of this work, see see https://doi.org/10.48550/arXiv.2504.12990.
+
+Here we show how to compute the information gain through clustering. To this end, we create two datasets by simulating the Langevin Dynamics of particles moving in two different bidimensional potential energy landscapes, one with two and one with four minima. We compare the information gain after clustering the particles' trajectories, using either one or both variables. 
+
+To start, let's import the packages we will need and create a folder in the cwd to save the results in.
+
+.. code-block:: python
+
+    from pathlib import Path
+    from typing import Callable
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    cwd = Path.cwd()
+    folder_name = "info_gain"
+    folder_path = cwd / folder_name
+    if not folder_path.exists():
+        folder_path.mkdir()
+
+
+Now, we want to simulate the Langevin Dynamics. We start defining the potential energy landscapes, with two and four minima. 
+
+.. code-block:: python
+
+    def energy_landscape_1(x: float, y: float) -> float:
+        """A potential energy landscape with 2 minima."""
+        sigma = 0.12  # Width of the Gaussian wells
+        gauss1 = np.exp(-(x**2 + y**2) / (2 * sigma**2))
+        gauss2 = np.exp(-(x**2 + (y - 1) ** 2) / (2 * sigma**2))
+        return -np.log(gauss1 + gauss2 + 1e-6)
+
+
+    def energy_landscape_2(x: float, y: float) -> float:
+        """A potential energy landscape with 4 minima."""
+        sigma = 0.12  # Width of the Gaussian wells
+        gauss1 = np.exp(-(x**2 + y**2) / (2 * sigma**2))
+        gauss2 = np.exp(-((x - 1) ** 2 + y**2) / (2 * sigma**2))
+        gauss3 = np.exp(-(x**2 + (y - 1) ** 2) / (2 * sigma**2))
+        gauss4 = np.exp(-((x - 1) ** 2 + (y - 1) ** 2) / (2 * sigma**2))
+        return -np.log(gauss1 + gauss2 + gauss3 + gauss4 + 1e-6)
+
+
+To compute the force acting on the particle, we need to compute the potential energy gradient. 
+
+.. code-block:: python
+
+    def numerical_gradient(
+        f: Callable[[float, float], float], x: float, y: float, h: float = 1e-5
+    ) -> tuple[float, float]:
+        """Compute numerical gradient using finite differences."""
+        grad_x = (f(x + h, y) - f(x - h, y)) / (2 * h)
+        grad_y = (f(x, y + h) - f(x, y - h)) / (2 * h)
+        return -grad_x, -grad_y
+
+
+This function simulates, for both energy landscapes, the dynamics of 100 particles for 10000 timesteps. Particles are initialized close to the minima. 
+
+.. code-block:: python
+
+    def create_trajectory(
+        energy_landscape: Callable[[float, float], float], file_name: Path
+    ) -> np.ndarray:
+        """Simulate Langevin Dynamics on a given energy landscape."""
+        rng = np.random.default_rng(0)
+        n_atoms = 100
+        time_steps = 10000
+        dt = 0.01  # Time step
+        diffusion_coeff = 0.6  # Diffusion coefficient (random noise strength)
+
+        if energy_landscape == energy_landscape_1:
+            particles = rng.standard_normal((n_atoms, 2)) * 0.2
+            particles[n_atoms // 2 :, 1] += 1.0
+        else:
+            particles = rng.standard_normal((n_atoms, 2)) * 0.2
+            n_group = n_atoms // 4
+            particles[n_group : 2 * n_group, 1] += 1  # (0, 1)
+            particles[2 * n_group : 3 * n_group, 0] += 1  # (1, 0)
+            particles[3 * n_group :, 0] += 1  # (1, 1)
+            particles[3 * n_group :, 1] += 1
+
+        trajectory = np.zeros((time_steps, n_atoms, 2))
+        for t in range(time_steps):
+            for i in range(n_atoms):
+                x, y = particles[i]
+                fx, fy = numerical_gradient(energy_landscape, x, y)
+                noise_x = np.sqrt(2 * diffusion_coeff * dt) * rng.standard_normal()
+                noise_y = np.sqrt(2 * diffusion_coeff * dt) * rng.standard_normal()
+
+                # Update position with deterministic force and stochastic term
+                particles[i, 0] += fx * dt + noise_x
+                particles[i, 1] += fy * dt + noise_y
+
+                trajectory[t, i] = particles[i]
+
+        plt.figure()
+        plt.plot(trajectory[:, :, 0], trajectory[:, :, 1])
+        plt.show()
+
+        dataset = np.transpose(trajectory, (1, 0, 2))
+        np.save(filename, dataset)
+        return dataset
+
+
+Let's simulate the trajectories and store them in two variables. We also save them as .npy files so that we don't have to simulate them every time. 
+
+.. code-block:: python
+
+    file_1 = folder_path / "trj_2.npy"  #  With 2 minima
+    file_2 = folder_path / "trj_4.npy"  #  With 4 minima
+
+    if not file_1.exists():
+        dataset_1 = create_trajectory(energy_landscape_1, file_1)
+
+    dataset_1 = np.load(file_1)
+
+    if not file_2.exists():
+        dataset_2 = create_trajectory(energy_landscape_2, file_2)
+
+    dataset_2 = np.load(file_2)
+
+
+Now that we have the trajectories, we can, for each one, perform Onion clustering and computing the corresponding information gain, using either 
+
+- only the y coordinate, or
+- both x and y coordinates
+
+For each case, we do the analysis for a range of values of the Onion clustering time resolution ∆t. The information gain for each analysis is saved in the variable "results". 
+To check if the clustering is working in a meaningful way, we also plot the results of Onion clustering for one specific value of ∆t. 
+
+
+.. code-block:: python
+
+    delta_t_list = np.unique(np.geomspace(2, 1000, 45, dtype=int))
+    results = np.zeros((4, delta_t_list.size))
+    example_delta_t = 4  #  Choosing a ∆t which works well to plot results
+
+    for i, dataset in enumerate([dataset_1, dataset_2]):
+        n_atoms, n_frames, n_dims = dataset.shape
+
+        # We can do clustering using only the y variable:
+        y_positions = dataset[:, :, 1]
+        info_gain_y = np.zeros(delta_t_list.size)
+
+        for j, delta_t in enumerate(delta_t_list):
+            reshaped_data = dynsight.onion.helpers.reshape_from_nt(
+                y_positions, delta_t
+            )
+            state_list, labels = dynsight.onion.onion_uni(reshaped_data)
+
+            if j == example_delta_t:
+                dynsight.onion.plot.plot_output_uni(
+                    f"Example_{i}_1D.png",
+                    reshaped_data,
+                    n_atoms,
+                    state_list,
+                )
+
+            # and compute the information gain:
+            info_gain_y[j], *_ = dynsight.analysis.compute_entropy_gain(
+                reshaped_data, labels, n_bins=40
+            )
+        results.append(info_gain_y)
+
+        # Or we can do clustering using both (x, y) variables:
+        info_gain_xy = np.zeros(delta_t_list.size)
+        tmp1_dataset = np.transpose(dataset, (2, 0, 1))
+        for j, delta_t in enumerate(delta_t_list):
+            reshaped_data = dynsight.onion.helpers.reshape_from_dnt(
+                tmp1_dataset, delta_t
+            )
+            state_list, labels = dynsight.onion.onion_multi(reshaped_data)
+
+            if j == example_delta_t:
+                dynsight.onion.plot.plot_output_multi(
+                    f"Example_{i}_2D.png",
+                    tmp1_dataset,
+                    state_list,
+                    labels,
+                    delta_t,
+                )
+
+            # and compute the information gain:
+            # We need an array (n_samples, n_dims), and labels (n_samples,)
+            n_sequences = int(n_frames / delta_t)
+            long_labels = np.repeat(labels, delta_t)
+            tmp = dataset[:, : n_sequences * delta_t, :]
+            ds_reshaped = tmp.reshape((-1, n_dims))
+
+            info_gain_xy[j], *_ = dynsight.analysis.compute_entropy_gain_multi(
+                ds_reshaped, long_labels, n_bins=[40, 40]
+            )
+        # Need to multiply by two because it's 2 dimensional, and the output
+        # of the info_gain functions is normalized by the log volume of the
+        # phase space, which is 2D is doubled
+        info_gain_xy *= 2
+        results.append(info_gain_xy)
+
+
+Here are the plots of the two datasets, with the different clusters identified when clustering the full, bi-dimensional data, using ∆t = 4 frames:
+
+.. list-table::
+   :widths: auto
+   :align: center
+
+   * - .. image:: _static/info_gain_clusters_1d.png
+     - .. image:: _static/info_gain_clusters_2d.png
+
+
+As can be seen, all the clusters are correctly identified at this time resolution ∆t. When we are using only the y-coordinate instead, as expected in both cases just two clusters can be identified (the two plots look the same but they are actually from the two different systems):
+
+.. list-table::
+   :widths: auto
+   :align: center
+
+   * - .. image:: _static/info_gain_clusters_1d_y.png
+     - .. image:: _static/info_gain_clusters_2d_y.png
+
+
+We can now plot, for every case and for every choice of ∆t, the corresponding information gain. 
+
+.. code-block:: python
+
+    colorlist = ["C0", "C2", "C1", "C3"]
+    markerlist = ["s", "o", "d", "o"]
+    labellist = [
+        "2 peaks - 1D clustering",
+        "2 peaks - 2D clustering",
+        "4 peaks - 1D clustering",
+        "4 peaks - 2D clustering",
+    ]
+
+    fig, ax = plt.subplots()
+    for i, system in enumerate(results):
+        ax.plot(
+            delta_t_list,
+            system,
+            label=labellist[i],
+            c=colorlist[i],
+            marker=markerlist[i],
+        )
+
+    ax.set_xlabel(r"Time resolution $\Delta t$ [frame]")
+    ax.set_ylabel(r"Information gain $\Delta H$ [bit]")
+    ax.set_xscale("log")
+    ax.legend()
+    plt.show()
+
+As can be seen in the plot below, clustering both datasets using only the y coordinate gives the same information gain, because only two clusters can be distinguished. 
+
+Clustering the trajectories in the energy potential with two minima using both variables gives once again the same information gain for small values of ∆t; then, the clustering performance degrades because the fraction of classifiable data points starts to decreases. 
+
+Finally, clustering the trajectories in the energy potential with four minima using both variables gives an information gain which is double the previous ones (at least for small ∆t), which makes sense, because 4 clusters are discovered instead of 2. For larger ∆t, we see the same degrading in performance that always affects clustering on multivariate distributions. 
+
+.. image:: _static/Information_gains.png
diff --git a/docs/source/sample_entropy.rst b/docs/source/example_sample_entropy.rst
similarity index 100%
rename from docs/source/sample_entropy.rst
rename to docs/source/example_sample_entropy.rst
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 31c71c99..e68865ab 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -3,6 +3,7 @@
    :caption: dynsight
    :maxdepth: 2
 
+   trajectory <trajectory>
    SOAP <soap>
    timeSOAP <time_soap>
    LENS <lens>
@@ -18,9 +19,18 @@
   :maxdepth: 2
   :caption: Examples:
 
-  Typical analysis workflow <analysis_workflow>
-  Information gain <info_gain>
-  Sample Entropy <sample_entropy>
+  Typical analysis workflow <example_analysis_workflow>
+  Information gain <example_info_gain>
+  Sample Entropy <example_sample_entropy>
+
+.. toctree::
+  :hidden:
+  :maxdepth: 2
+  :caption: Recipes:
+
+  Descriptors from a Trj <descr_from_trj>
+  Dimensionality reduction methods <soap_dim_red>
+  Information gain analysis <info_gain>
 
 .. toctree::
   :hidden:
@@ -63,20 +73,18 @@ you are using Python 3.10 and below, you can use :mod:`cpctools` to access
 
   $ pip install cpctools
 
-If you want to use the :mod:`dynsight.tica` module you will need to install the
-deeptime package. This can be done with with pip::
-
-  $ pip install deeptime
-
-or with conda::
-
-  $ conda install -c conda-forge deeptime
-
 If you want to use the :mod:`dynsight.vision` and :mod:`dynsight.track` modules
 you will need to install a series of packages. This can be done with with pip::
 
   $ pip install ultralytics PyYAML
 
+How to get started
+------------------
+
+We suggest you give a read to the ``dynsight.trajectory`` module documentation,
+which offers a compact and easy way of using most of the ``dynsight`` tools. 
+Also, the documentation offers some copiable Recipes and Examples for the most
+common analyses.
 
 How to contribute
 -----------------
diff --git a/docs/source/info_gain.rst b/docs/source/info_gain.rst
index 9fac3187..44b1d2f2 100644
--- a/docs/source/info_gain.rst
+++ b/docs/source/info_gain.rst
@@ -1,257 +1,230 @@
-Information gain computation
-============================
+Information gain analysis
+=========================
 
-For the theoretical aspects of this work, see INSERT REF.
+For the theoretical aspects of this work, see https://doi.org/10.48550/arXiv.2504.12990.
 
-Here we show how to compute the information gain through clustering. To this end, we create two datasets by simulating the Langevin Dynamics of particles moving in two different bidimensional potential energy landscapes, one with two and one with four minima. We compare the information gain after clustering the particles' trajectories, using either one or both variables. 
+This recipe explains how to compute the information gain through clustering 
+analysis. We use a synthetic dataset containing a signal that oscillates
+between 0 and 1, with Gaussian noise. Onion clustering is run on a broad
+range of time resolutions ∆t. The information gain and the Shannon entropy of
+the environments is computed for each value of ∆t. The analysis is implemented 
+using onion 2.0.0 ("onion smooth").
 
-To start, let's import the packages we will need and create a folder in the cwd to save the results in.
+.. warning::
 
-.. code-block:: python
+    This code works when run from the ``/docs`` directory of the ``dynsight`` 
+    repo. To use it elsewhere, you have to change the ``Path`` variables
+    accordingly.
+
+First of all, we import all the packages and objects we'll need:
+
+.. testcode:: recipe3-test
 
     from pathlib import Path
-    from typing import Callable
     import numpy as np
+    import dynsight
+    from dynsight.trajectory import Insight
     import matplotlib.pyplot as plt
 
-    cwd = Path.cwd()
-    folder_name = "info_gain"
-    folder_path = cwd / folder_name
-    if not folder_path.exists():
-        folder_path.mkdir()
-
-
-Now, we want to simulate the Langevin Dynamics. We start defining the potential energy landscapes, with two and four minima. 
-
-.. code-block:: python
-
-    def energy_landscape_1(x: float, y: float) -> float:
-        """A potential energy landscape with 2 minima."""
-        sigma = 0.12  # Width of the Gaussian wells
-        gauss1 = np.exp(-(x**2 + y**2) / (2 * sigma**2))
-        gauss2 = np.exp(-(x**2 + (y - 1) ** 2) / (2 * sigma**2))
-        return -np.log(gauss1 + gauss2 + 1e-6)
-
-
-    def energy_landscape_2(x: float, y: float) -> float:
-        """A potential energy landscape with 4 minima."""
-        sigma = 0.12  # Width of the Gaussian wells
-        gauss1 = np.exp(-(x**2 + y**2) / (2 * sigma**2))
-        gauss2 = np.exp(-((x - 1) ** 2 + y**2) / (2 * sigma**2))
-        gauss3 = np.exp(-(x**2 + (y - 1) ** 2) / (2 * sigma**2))
-        gauss4 = np.exp(-((x - 1) ** 2 + (y - 1) ** 2) / (2 * sigma**2))
-        return -np.log(gauss1 + gauss2 + gauss3 + gauss4 + 1e-6)
-
-
-To compute the force acting on the particle, we need to compute the potential energy gradient. 
-
-.. code-block:: python
-
-    def numerical_gradient(
-        f: Callable[[float, float], float], x: float, y: float, h: float = 1e-5
-    ) -> tuple[float, float]:
-        """Compute numerical gradient using finite differences."""
-        grad_x = (f(x + h, y) - f(x - h, y)) / (2 * h)
-        grad_y = (f(x, y + h) - f(x, y - h)) / (2 * h)
-        return -grad_x, -grad_y
-
-
-This function simulates, for both energy landscapes, the dynamics of 100 particles for 10000 timesteps. Particles are initialized close to the minima. 
-
-.. code-block:: python
-
-    def create_trajectory(
-        energy_landscape: Callable[[float, float], float], file_name: Path
-    ) -> np.ndarray:
-        """Simulate Langevin Dynamics on a given energy landscape."""
-        rng = np.random.default_rng(0)
-        n_atoms = 100
-        time_steps = 10000
-        dt = 0.01  # Time step
-        diffusion_coeff = 0.6  # Diffusion coefficient (random noise strength)
-
-        if energy_landscape == energy_landscape_1:
-            particles = rng.standard_normal((n_atoms, 2)) * 0.2
-            particles[n_atoms // 2 :, 1] += 1.0
-        else:
-            particles = rng.standard_normal((n_atoms, 2)) * 0.2
-            n_group = n_atoms // 4
-            particles[n_group : 2 * n_group, 1] += 1  # (0, 1)
-            particles[2 * n_group : 3 * n_group, 0] += 1  # (1, 0)
-            particles[3 * n_group :, 0] += 1  # (1, 1)
-            particles[3 * n_group :, 1] += 1
-
-        trajectory = np.zeros((time_steps, n_atoms, 2))
-        for t in range(time_steps):
-            for i in range(n_atoms):
-                x, y = particles[i]
-                fx, fy = numerical_gradient(energy_landscape, x, y)
-                noise_x = np.sqrt(2 * diffusion_coeff * dt) * rng.standard_normal()
-                noise_y = np.sqrt(2 * diffusion_coeff * dt) * rng.standard_normal()
-
-                # Update position with deterministic force and stochastic term
-                particles[i, 0] += fx * dt + noise_x
-                particles[i, 1] += fy * dt + noise_y
-
-                trajectory[t, i] = particles[i]
-
-        plt.figure()
-        plt.plot(trajectory[:, :, 0], trajectory[:, :, 1])
-        plt.show()
-
-        dataset = np.transpose(trajectory, (1, 0, 2))
-        np.save(filename, dataset)
-        return dataset
 
+Let's start by creating a the synthetic dataset:
 
-Let's simulate the trajectories and store them in two variables. We also save them as .npy files so that we don't have to simulate them every time. 
+.. testcode:: recipe3-test
 
-.. code-block:: python
+    rng = np.random.default_rng(1234)
 
-    file_1 = folder_path / "trj_2.npy"  #  With 2 minima
-    file_2 = folder_path / "trj_4.npy"  #  With 4 minima
+    # Parameters
+    n_atoms = 10
+    num_blocks = 10
+    block_size = 100
+    sigma = 0.1
 
-    if not file_1.exists():
-        dataset_1 = create_trajectory(energy_landscape_1, file_1)
+    # Generate the array
+    tmp_data = [
+        np.concatenate(
+            [
+                rng.normal(loc=(i % 2), scale=sigma, size=block_size)
+                for i in range(num_blocks)
+            ]
+        )
+        for _ in range(n_atoms)
+    ]
+    data = np.array(tmp_data)
 
-    dataset_1 = np.load(file_1)
 
-    if not file_2.exists():
-        dataset_2 = create_trajectory(energy_landscape_2, file_2)
+The following function takes as input the dataset, and a list of values
+of time resolutions ∆t, and for each of these it performs Onion clustering, and
+computes the information gain achieved through clustering with that ∆t. 
 
-    dataset_2 = np.load(file_2)
+.. warning::
 
+    For now, this only works with univariate datasets.
 
-Now that we have the trajectories, we can, for each one, perform Onion clustering and computing the corresponding information gain, using either 
+The function's output is a tuple of ``np.ndarray``, which for each value of ∆t
+contain:
 
-- only the y coordinate, or
-- both x and y coordinates
+* the number of identified clusters - shape (delta_t_list.size,);
+* the population fraction of each cluster - shape (delta_t_list.size, n_clust);
+* the information gain - shape (delta_t_list.size,);
+* the Shannon entropy of each cluster - shape (delta_t_list.size, n_clust).
 
-For each case, we do the analysis for a range of values of the Onion clustering time resolution ∆t. The information gain for each analysis is saved in the variable "results". 
-To check if the clustering is working in a meaningful way, we also plot the results of Onion clustering for one specific value of ∆t. 
+Additionally, the (float) dataset Shannon entropy is returned.
 
+.. testcode:: recipe3-test
 
-.. code-block:: python
+    def info_gain_with_onion(
+        delta_t_list: np.ndarray | list[int],
+        data: np.array,
+        n_bins: int = 40,
+    ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, float]:
+        """Performs full information gain analysis with Onion clustering."""
+        data_range = (np.min(data), np.max(data))
 
-    delta_t_list = np.unique(np.geomspace(2, 1000, 45, dtype=int))
-    results = np.zeros((4, delta_t_list.size))
-    example_delta_t = 4  #  Choosing a ∆t which works well to plot results
+        n_clusters = np.zeros(delta_t_list.size)
+        clusters_frac = []
+        info_gain = np.zeros(delta_t_list.size)
+        clusters_entr = []
 
-    for i, dataset in enumerate([dataset_1, dataset_2]):
-        n_atoms, n_frames, n_dims = dataset.shape
+        for i, delta_t in enumerate(delta_t_list):
+            state_list, labels = dynsight.onion.onion_uni_smooth(
+                data,
+                delta_t=delta_t,
+            )
 
-        # We can do clustering using only the y variable:
-        y_positions = dataset[:, :, 1]
-        info_gain_y = np.zeros(delta_t_list.size)
+            n_clusters[i] = len(state_list)
+            tmp_frac = [0.0] + [state.perc for state in state_list]
+            tmp_frac[0] = 1.0 - np.sum(tmp_frac)
+            clusters_frac.append(tmp_frac)
 
-        for j, delta_t in enumerate(delta_t_list):
-            reshaped_data = dynsight.onion.helpers.reshape_from_nt(
-                y_positions, delta_t
+            flat_data = data.flatten()
+            flat_labels = labels.flatten()
+            info_gain[i], _, h_0, _ = dynsight.analysis.compute_entropy_gain(
+                flat_data, flat_labels, n_bins=n_bins,
             )
-            state_list, labels = dynsight.onion.onion_uni(reshaped_data)
-
-            if j == example_delta_t:
-                dynsight.onion.plot.plot_output_uni(
-                    f"Example_{i}_1D.png",
-                    reshaped_data,
-                    n_atoms,
-                    state_list,
-                )
 
-            # and compute the information gain:
-            info_gain_y[j], *_ = dynsight.analysis.compute_entropy_gain(
-                reshaped_data, labels, n_bins=40
-            )
-        results.append(info_gain_y)
-
-        # Or we can do clustering using both (x, y) variables:
-        info_gain_xy = np.zeros(delta_t_list.size)
-        tmp1_dataset = np.transpose(dataset, (2, 0, 1))
-        for j, delta_t in enumerate(delta_t_list):
-            reshaped_data = dynsight.onion.helpers.reshape_from_dnt(
-                tmp1_dataset, delta_t
-            )
-            state_list, labels = dynsight.onion.onion_multi(reshaped_data)
-
-            if j == example_delta_t:
-                dynsight.onion.plot.plot_output_multi(
-                    f"Example_{i}_2D.png",
-                    tmp1_dataset,
-                    state_list,
-                    labels,
-                    delta_t,
+            tmp_entr = []
+            label_list = np.unique(labels)
+            if label_list[0] != -1:
+                tmp_entr.append(-1.0)
+
+            for _, lab in enumerate(label_list):
+                mask = labels == lab
+                selected_points = data[mask]
+                tmp_entr.append(
+                    dynsight.analysis.compute_shannon(
+                        selected_points,
+                        data_range,
+                        n_bins=n_bins,
+                    )
                 )
-
-            # and compute the information gain:
-            # We need an array (n_samples, n_dims), and labels (n_samples,)
-            n_sequences = int(n_frames / delta_t)
-            long_labels = np.repeat(labels, delta_t)
-            tmp = dataset[:, : n_sequences * delta_t, :]
-            ds_reshaped = tmp.reshape((-1, n_dims))
-
-            info_gain_xy[j], *_ = dynsight.analysis.compute_entropy_gain_multi(
-                ds_reshaped, long_labels, n_bins=[40, 40]
+            clusters_entr.append(tmp_entr)
+
+        max_n_envs = np.max([len(elem) for elem in clusters_entr])
+        for i, elem in enumerate(clusters_entr):
+            while len(elem) < max_n_envs:
+                elem.append(-1.0)
+                clusters_frac[i].append(0.0)
+
+        cl_frac = np.array(clusters_frac)
+        cl_entr = np.array(clusters_entr)
+
+        return n_clusters, cl_frac, info_gain, cl_entr, h_0
+
+    # Example usage
+    _, n_frames = data.shape
+    delta_t_list = np.unique(np.geomspace(2, n_frames, 10, dtype=int))
+
+    n_cl, cl_frac, info_gain, cl_entr, h_0 = info_gain_with_onion(
+        delta_t_list,
+        data,
+    )
+
+
+A default visualization of the results of this analysis can be obtained with
+the following function. Be aware that this could require some tweaking to ensure
+that clusters identified at different ∆t are matched in the way the user want
+them to, expecially when different number of clusters are found for different ∆t
+values.
+
+.. testcode:: recipe3-test
+
+    def plot_info_results(
+        delta_t_list: np.ndarray | list[int],
+        cl_frac: np.ndarray,
+        cl_entr: np.ndarray,
+        h_0: float,
+        file_path: Path,
+    ) -> None:
+        frac = cl_frac.T
+        entr = cl_entr.T
+        s_list = []
+        for i, st_fr in enumerate(frac):
+            s_list.append(st_fr * entr[i])
+        s_cumul = [s_list[0]]
+        for _, tmp_s in enumerate(s_list[1:]):
+            s_cumul.append(s_cumul[-1] + tmp_s)
+
+        fig, ax = plt.subplots()
+
+        i_0 = (1 - h_0) * np.ones(len(delta_t_list))
+        ax.plot(delta_t_list, i_0, ls="--", c="black", marker="")  # I_0
+        ax.fill_between(
+            delta_t_list,
+            1,
+            1 - s_cumul[0],
+            alpha=0.5,
+        )
+        for i, tmp_s in enumerate(s_cumul[1:]):
+            ax.fill_between(
+                delta_t_list,
+                1 - s_cumul[i],
+                1 - tmp_s,
+                alpha=0.5,
             )
-        # Need to multiply by two because it's 2 dimensional, and the output
-        # of the info_gain functions is normalized by the log volume of the
-        # phase space, which is 2D is doubled
-        info_gain_xy *= 2
-        results.append(info_gain_xy)
-
-
-Here are the plots of the two datasets, with the different clusters identified when clustering the full, bi-dimensional data, using ∆t = 4 frames:
-
-.. list-table::
-   :widths: auto
-   :align: center
-
-   * - .. image:: _static/info_gain_clusters_1d.png
-     - .. image:: _static/info_gain_clusters_2d.png
-
-
-As can be seen, all the clusters are correctly identified at this time resolution ∆t. When we are using only the y-coordinate instead, as expected in both cases just two clusters can be identified (the two plots look the same but they are actually from the two different systems):
+        ax.fill_between(
+            delta_t_list, 1 - s_cumul[-1], 1 - h_0, color="gainsboro",
+        )
+        ax.plot(
+            delta_t_list, 1 - s_cumul[-1], c="black", marker="",
+        )  # I_clust
 
-.. list-table::
-   :widths: auto
-   :align: center
+        ax.set_ylim(0.0, 1.0)
+        ax.set_xlabel(r"Time resolution $\Delta t$")
+        ax.set_ylabel(r"Information $I$")
+        ax.set_xscale("log")
 
-   * - .. image:: _static/info_gain_clusters_1d_y.png
-     - .. image:: _static/info_gain_clusters_2d_y.png
+        fig.savefig(file_path, dpi=600)
+        plt.close()
 
+    # Example usage
+    plot_info_results(
+        delta_t_list,
+        cl_frac,
+        cl_entr,
+        h_0,
+        Path("./source/_static/info_plot.png"),
+    )
 
-We can now plot, for every case and for every choice of ∆t, the corresponding information gain. 
+The figure obtained (see below) shows, for each value of ∆t:
 
-.. code-block:: python
+* The initial information (1 - H) of the entire dataset: dashed line;
+* The information after clustering: solid line;
+* The information gained through clustering ∆I: gray area;
+* The Shannon entropy of each of the discovered clusters: colored bands.
 
-    colorlist = ["C0", "C2", "C1", "C3"]
-    markerlist = ["s", "o", "d", "o"]
-    labellist = [
-        "2 peaks - 1D clustering",
-        "2 peaks - 2D clustering",
-        "4 peaks - 1D clustering",
-        "4 peaks - 2D clustering",
-    ]
+In this case, 2 states are correctly identified for ∆t <= 100 (green and orange),
+with an information gain of around 0.2.
+For ∆t > 100 all the data points remain unclassified (blue), and the information
+gain goes to 0.
 
-    fig, ax = plt.subplots()
-    for i, system in enumerate(results):
-        ax.plot(
-            delta_t_list,
-            system,
-            label=labellist[i],
-            c=colorlist[i],
-            marker=markerlist[i],
-        )
+.. image:: _static/info_plot.png
 
-    ax.set_xlabel(r"Time resolution $\Delta t$ [frame]")
-    ax.set_ylabel(r"Information gain $\Delta H$ [bit]")
-    ax.set_xscale("log")
-    ax.legend()
-    plt.show()
 
-As can be seen in the plot below, clustering both datasets using only the y coordinate gives the same information gain, because only two clusters can be distinguished. 
+.. raw:: html
 
-Clustering the trajectories in the energy potential with two minima using both variables gives once again the same information gain for small values of ∆t; then, the clustering performance degrades because the fraction of classifiable data points starts to decreases. 
+    <a class="btn-download" href="../_static/recipes/info_gain.py" download>⬇️ Download Python Script</a>
 
-Finally, clustering the trajectories in the energy potential with four minima using both variables gives an information gain which is double the previous ones (at least for small ∆t), which makes sense, because 4 clusters are discovered instead of 2. For larger ∆t, we see the same degrading in performance that always affects clustering on multivariate distributions. 
+.. testcode:: recipe3-test
+    :hide:
 
-.. image:: _static/Information_gains.png
+    assert np.isclose(info_gain[0], 0.19043795503255656)
diff --git a/docs/source/soap_dim_red.rst b/docs/source/soap_dim_red.rst
new file mode 100644
index 00000000..c230b678
--- /dev/null
+++ b/docs/source/soap_dim_red.rst
@@ -0,0 +1,279 @@
+Dimensionality reduction methods 
+================================
+
+This recipe explains how to compute descriptors via dimensionality reduction
+of a multivariated descriptor. These example use SOAP, but the same approaches
+can be applied to a variety of other quantities. The newly computed
+descriptors are always stored in an :class:`.trajectory.Insight` variable.
+
+All the functions take as optional parameters one or more file paths; if these
+paths are passed, the function, before computing some quantity from scratch,
+tries to load it from file, in case it has already been previously computed
+and saved. For SOAP, which is required in all the examples, we use a function
+from dynsight.utilities.
+
+.. warning::
+
+    Please consider that the SOAP dataset can be very large, due to the high
+    dimensionality, thus calculations can be expensive, and saving to/loading 
+    from file quite slow.
+
+.. warning::
+
+    This code works when run from the ``/docs`` directory of the ``dynsight`` 
+    repo. To use it elsewhere, you have to change the ``Path`` variables
+    accordingly.
+
+First of all, we import all the packages and objects we'll need:
+
+.. testcode:: recipe2-test
+
+    from pathlib import Path
+    import dynsight
+    from dynsight.trajectory import Trj, Insight
+    from dynsight.utilities import load_or_compute_soap
+    from sklearn.decomposition import PCA
+
+Let's start by creating a :class:`.trajectory.Trj` object to use as a
+starting point for the examples:
+
+.. testcode:: recipe2-test
+
+    # Loading an example trajectory
+    files_path = Path("../tests/systems/")
+    trj = Trj.init_from_xtc(
+        traj_file=files_path / "balls_7_nvt.xtc",
+        topo_file=files_path / "balls_7_nvt.gro",
+    )
+
+Principal Component Analysis (PCA)
+----------------------------------
+
+Principal Component Analysis is a dimensionality reduction method that finds
+the combinations of components that maximize the data variance. More details
+on the algorithm here https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.
+
+This function takes as input a :class:`.trajectory.Trj` and all the relevant
+parameters, and performs the PCA of the corresponding SOAP dataset.
+
+``n_components`` is the number of PCs that the function stores in the output.
+
+.. testcode:: recipe2-test
+
+    def compute_soap_pca(
+        trj: Trj,
+        r_cut: float,
+        n_max: int,
+        l_max: int,
+        n_components: int,
+        soap_path: Path | None = None,
+        pca_path: Path | None = None,
+        selection: str = "all",
+        centers: str = "all",
+        respect_pbc: bool = True,
+        n_core: int = 1,
+    ) -> Insight:
+        if pca_path is not None and pca_path.exists():
+            return Insight.load_from_json(pca_path)
+
+        soap = load_or_compute_soap(
+            trj=trj,
+            r_cut=r_cut,
+            n_max=n_max,
+            l_max=l_max,
+            selection=selection,
+            centers=centers,
+            respect_pbc=respect_pbc,
+            n_core=n_core,
+            soap_path=soap_path,
+        )
+
+        n_atom, n_frames, n_dims = soap.dataset.shape
+        reshaped_soap = soap.dataset.reshape(n_atom * n_frames, n_dims)
+        pca = PCA(n_components=n_components)
+        transformed_soap = pca.fit_transform(reshaped_soap)
+        pca_ds = transformed_soap.reshape(n_atom, n_frames, -1)
+
+        soap_pca = Insight(pca_ds, meta=soap.meta.copy())
+
+        if pca_path is not None:
+            soap_pca.dump_to_json(pca_path)
+
+        return soap_pca
+
+    # Example of how to use
+    soap_pc1 = compute_soap_pca(
+        trj=trj,
+        r_cut=2.0,
+        n_max=4,
+        l_max=4,
+        n_components=1,
+    )
+
+The output :class:`.trajectory.Insight` stores the SOAP information in its
+"meta" attribute.
+
+Time-lagged Independent Component Analysis (TICA)
+-------------------------------------------------
+
+More details on the algorithm here:
+
+.. toctree::
+  :maxdepth: 1
+
+  many_body_tica <_autosummary/dynsight.tica.many_body_tica>
+
+This function takes as input a :class:`.trajectory.Trj` and all the relevant
+parameters, and performs the TICA of the corresponding SOAP dataset.
+
+``lag_time`` is the time lag used to perform TICA.
+``tica_dim`` is the number of TICs that the function stores in the output.
+
+.. testcode:: recipe2-test
+
+    def compute_soap_tica(
+        trj: Trj,
+        r_cut: float,
+        n_max: int,
+        l_max: int,
+        lag_time: int,
+        tica_dim: int,
+        soap_path: Path | None = None,
+        tica_path: Path | None = None,
+        selection: str = "all",
+        centers: str = "all",
+        respect_pbc: bool = True,
+        n_core: int = 1,
+    ) -> Insight:
+        if tica_path is not None and tica_path.exists():
+            return Insight.load_from_json(tica_path)
+
+        soap = load_or_compute_soap(
+            trj=trj,
+            r_cut=r_cut,
+            n_max=n_max,
+            l_max=l_max,
+            selection=selection,
+            centers=centers,
+            respect_pbc=respect_pbc,
+            n_core=n_core,
+            soap_path=soap_path,
+        )
+
+        rel_times, _, tica_ds = dynsight.tica.many_body_tica(
+            soap.dataset,
+            lag_time=lag_time,
+            tica_dim=tica_dim,
+        )
+
+        meta = soap.meta.copy()
+        meta.update({
+            "lag_time": lag_time,
+            "rel_times": rel_times,
+        })
+        soap_tica = Insight(tica_ds, meta=meta)
+
+        if tica_path is not None:
+            soap_tica.dump_to_json(tica_path)
+
+        return soap_tica
+
+    # Example of how to use
+    soap_tic1 = compute_soap_tica(
+        trj=trj,
+        r_cut=10.0,
+        n_max=4,
+        l_max=4,
+        lag_time=10,
+        tica_dim=1,
+    )
+
+The output :class:`.trajectory.Insight` stores the SOAP information in its
+"meta" attribute, together with the ``lag_time`` parameter and ``rel_times``, 
+the relaxation times of the computed TICs.
+
+
+timeSOAP (tSOAP)
+----------------
+
+More details on the algorithm here:
+
+.. toctree::
+  :maxdepth: 1
+
+  timesoap <_autosummary/dynsight.soap.timesoap>
+
+This function takes as input a :class:`.trajectory.Trj` and all the relevant
+parameters, and computes the corresponding timeSOAP dataset.
+
+``delay`` is the time lag used to perform timeSOAP.
+
+.. testcode:: recipe2-test
+
+    def compute_timesoap(
+        trj: Trj,
+        r_cut: float,
+        n_max: int,
+        l_max: int,
+        delay: int = 1,
+        soap_path: Path | None = None,
+        tsoap_path: Path | None = None,
+        selection: str = "all",
+        centers: str = "all",
+        respect_pbc: bool = True,
+        n_core: int = 1,
+    ) -> Insight:
+        if tsoap_path is not None and tsoap_path.exists():
+            return Insight.load_from_json(tsoap_path)
+
+        soap = load_or_compute_soap(
+            trj=trj,
+            r_cut=r_cut,
+            n_max=n_max,
+            l_max=l_max,
+            selection=selection,
+            centers=centers,
+            respect_pbc=respect_pbc,
+            n_core=n_core,
+            soap_path=soap_path,
+        )
+
+        tsoap = soap.get_angular_velocity(delay=delay)
+
+        if tsoap_path is not None:
+            tsoap.dump_to_json(tsoap_path)
+
+        return tsoap
+
+    # Example of how to use
+    tsoap = compute_timesoap(
+        trj=trj,
+        r_cut=10.0,
+        n_max=4,
+        l_max=4,
+    )
+
+The output :class:`.trajectory.Insight` stores the SOAP information in its
+"meta" attribute, together with the ``delay`` parameter.
+
+Notice that, differently from SOAP - which is computed for every frame, tSOAP
+is computed for every pair of frames. Thus, the tSOAP dataset has shape 
+``(n_particles, n_frames - 1)``. Consequently, if you need to match the tSOAP
+values with the particles along the trajectory, you will need to use a sliced
+trajectory (removing the last frame). The easiest way to do this is:
+
+.. testcode:: recipe2-test
+
+    trajslice = slice(0, -1, 1)
+    shorter_trj = trj.with_slice(trajslice=trajslice)
+
+.. raw:: html
+
+    <a class="btn-download" href="../_static/recipes/soap_dim_red.py" download>⬇️ Download Python Script</a>
+
+.. testcode:: recipe2-test
+    :hide:
+
+    assert soap_pc1.dataset.shape == (7, 201, 1)
+    assert soap_tic1.dataset.shape == (7, 201, 1)
+    assert tsoap.dataset.shape == (7, 200)
diff --git a/docs/source/trajectory.rst b/docs/source/trajectory.rst
new file mode 100644
index 00000000..07a3f4df
--- /dev/null
+++ b/docs/source/trajectory.rst
@@ -0,0 +1,29 @@
+trajectory
+==========
+
+This module offers an object-oriented implementation of the ``dynsight``
+utilities, to facilitate the workflow of any analysis.
+Moreover, the ``trajectory`` classes make it easier to save and share results,
+and have a built-in logging module.
+
+The :class:`.trajectory.Trj` class is a container for a MDAnalysis.Universe,
+which allows for the computation of all the descriptors implemented in ``dynsight``.
+
+These descriptors, as well as the output of subsequent analyses, are stored in
+:class:`.trajectory.Insight` or :class:`.trajectory.ClusterInsight` objects.
+
+We recommend the users, when possible, to write code using this module.
+
+Complete examples can be found in the recipes section of this documentation. 
+
+Classes
+-------
+
+.. toctree::
+  :maxdepth: 1
+
+  Trj <_autosummary/dynsight.trajectory.Trj>
+  Insight <_autosummary/dynsight.trajectory.Insight>
+  ClusterInsight <_autosummary/dynsight.trajectory.ClusterInsight>
+  OnionInsight <_autosummary/dynsight.trajectory.OnionInsight>
+  OnionSmoothInsight <_autosummary/dynsight.trajectory.OnionSmoothInsight>
diff --git a/examples/info_gain.py b/examples/info_gain_example.py
similarity index 100%
rename from examples/info_gain.py
rename to examples/info_gain_example.py
diff --git a/pyproject.toml b/pyproject.toml
index b1949aee..801d3f0d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,7 +10,7 @@ maintainers = [
     { name = "Simone Martino", email = "s.martino0898@gmail.com" },
 ]
 
-dependencies = ["numpy", "dscribe", "tropea-clustering", "MDAnalysis"]
+dependencies = ["numpy", "dscribe", "tropea-clustering", "MDAnalysis", "deeptime"]
 # Set by cpctools.
 requires-python = ">=3.8"
 dynamic = ["version"]
@@ -92,6 +92,7 @@ check_untyped_defs = true
 disallow_untyped_decorators = true
 warn_unreachable = true
 disallow_any_generics = true
+exclude = 'docs/build/html/_static'
 
 [[tool.mypy.overrides]]
 module = [
@@ -111,5 +112,7 @@ module = [
     'ultralytics.*',
     'opencv-python.*',
     'cv2.*',
+    'deeptime.*',
+    'sklearn.decomposition.*'
 ]
 ignore_missing_imports = true
diff --git a/src/dynsight/__init__.py b/src/dynsight/__init__.py
index 36bc1232..0d970042 100644
--- a/src/dynsight/__init__.py
+++ b/src/dynsight/__init__.py
@@ -5,6 +5,7 @@
     data_processing,
     hdf5er,
     lens,
+    logs,
     onion,
     soap,
     tica,
@@ -18,6 +19,7 @@
     "data_processing",
     "hdf5er",
     "lens",
+    "logs",
     "onion",
     "soap",
     "tica",
diff --git a/src/dynsight/_internal/logs.py b/src/dynsight/_internal/logs.py
new file mode 100644
index 00000000..6bf9afe6
--- /dev/null
+++ b/src/dynsight/_internal/logs.py
@@ -0,0 +1,29 @@
+"""logging package."""
+
+from datetime import datetime, timezone
+from pathlib import Path
+
+
+class Logger:
+    """Creates and save human-readible log."""
+
+    def __init__(self) -> None:
+        self._log: list[str] = []
+
+    def log(self, msg: str) -> None:
+        timestamp = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M:%S")
+        entry = f"[{timestamp}] {msg}"
+        self._log.append(entry)
+
+    def save(self, filename: Path) -> None:
+        with filename.open("w", encoding="utf-8") as f:
+            f.write("\n".join(self._log))
+
+    def clear(self) -> None:
+        self._log = []
+
+    def get(self) -> str:
+        return "\n".join(self._log)
+
+
+logger = Logger()
diff --git a/src/dynsight/_internal/trajectory/cluster_insight.py b/src/dynsight/_internal/trajectory/cluster_insight.py
index db4cce55..70d641cb 100644
--- a/src/dynsight/_internal/trajectory/cluster_insight.py
+++ b/src/dynsight/_internal/trajectory/cluster_insight.py
@@ -16,6 +16,7 @@
 
 
 import dynsight
+from dynsight.logs import logger
 
 UNIVAR_DIM = 2
 
@@ -36,6 +37,7 @@ def dump_to_json(self, file_path: Path) -> None:
         data["labels"] = data["labels"].tolist()
         with file_path.open("w") as file:
             json.dump(data, file, indent=4)
+        logger.log(f"ClusterInsight saved to {file_path}.")
 
     @classmethod
     def load_from_json(cls, file_path: Path) -> ClusterInsight:
@@ -48,7 +50,10 @@ def load_from_json(cls, file_path: Path) -> ClusterInsight:
             data = json.load(file)
         if "labels" not in data:
             msg = "'labels' key not found in JSON file."
+            logger.log(msg)
             raise ValueError(msg)
+
+        logger.log(f"ClusterInsight loaded from {file_path}.")
         return cls(labels=np.array(data.get("labels"), dtype=np.int64))
 
 
@@ -89,6 +94,7 @@ def dump_to_json(self, file_path: Path) -> None:
         data["state_list"] = new_state_list
         with file_path.open("w") as file:
             json.dump(data, file, indent=4)
+        logger.log(f"OnionInsight saved to {file_path}.")
 
     @classmethod
     def load_from_json(cls, file_path: Path) -> OnionInsight:
@@ -101,7 +107,10 @@ def load_from_json(cls, file_path: Path) -> OnionInsight:
             data = json.load(file)
         if "state_list" not in data:
             msg = "'state_list' key not found in JSON file."
+            logger.log(msg)
             raise ValueError(msg)
+
+        logger.log(f"OnionInsight loaded from {file_path}.")
         return cls(
             labels=np.array(data.get("labels")),
             state_list=data.get("state_list"),
@@ -126,6 +135,8 @@ def plot_output(self, file_path: Path, data_insight: Insight) -> None:
                 self.labels,
                 self.meta["delta_t"],
             )
+        attr_dict = {"file_path": file_path}
+        logger.log(f"plot_output() with args {attr_dict}.")
 
     def plot_one_trj(
         self,
@@ -151,6 +162,9 @@ def plot_one_trj(
                 self.labels,
             )
 
+        attr_dict = {"file_path": file_path, "particle_id": particle_id}
+        logger.log(f"plot_one_trj() with args {attr_dict}.")
+
     def plot_medoids(self, file_path: Path, data_insight: Insight) -> None:
         """Plot the average sequence of each onion cluster."""
         if data_insight.dataset.ndim == UNIVAR_DIM:
@@ -166,6 +180,8 @@ def plot_medoids(self, file_path: Path, data_insight: Insight) -> None:
                 data_insight.dataset,
                 self.labels,
             )
+        attr_dict = {"file_path": file_path}
+        logger.log(f"plot_medoids() with args {attr_dict}.")
 
     def plot_state_populations(
         self,
@@ -179,6 +195,8 @@ def plot_state_populations(
             self.meta["delta_t"],
             self.labels,
         )
+        attr_dict = {"file_path": file_path}
+        logger.log(f"plot_state_populations() with args {attr_dict}.")
 
     def plot_sankey(
         self,
@@ -194,6 +212,9 @@ def plot_sankey(
             frame_list,
         )
 
+        attr_dict = {"file_path": file_path, "frame_list": frame_list}
+        logger.log(f"plot_state_populations() with args {attr_dict}.")
+
 
 @dataclass(frozen=True)
 class OnionSmoothInsight(ClusterInsight):
@@ -229,6 +250,7 @@ def dump_to_json(self, file_path: Path) -> None:
         data["state_list"] = new_state_list
         with file_path.open("w") as file:
             json.dump(data, file, indent=4)
+        logger.log(f"OnionSmoothInsight saved to {file_path}.")
 
     @classmethod
     def load_from_json(cls, file_path: Path) -> OnionSmoothInsight:
@@ -241,7 +263,10 @@ def load_from_json(cls, file_path: Path) -> OnionSmoothInsight:
             data = json.load(file)
         if "state_list" not in data:
             msg = "'state_list' key not found in JSON file."
+            logger.log(msg)
             raise ValueError(msg)
+
+        logger.log(f"OnionSmoothInsight loaded from {file_path}.")
         return cls(
             labels=np.array(data.get("labels")),
             state_list=data.get("state_list"),
@@ -263,6 +288,8 @@ def plot_output(self, file_path: Path, data_insight: Insight) -> None:
                 self.state_list,
                 self.labels,
             )
+        attr_dict = {"file_path": file_path}
+        logger.log(f"plot_output() with args {attr_dict}.")
 
     def plot_one_trj(
         self,
@@ -285,6 +312,8 @@ def plot_one_trj(
                 data_insight.dataset,
                 self.labels,
             )
+        attr_dict = {"file_path": file_path, "particle_id": particle_id}
+        logger.log(f"plot_one_trj() with args {attr_dict}.")
 
     def plot_state_populations(
         self,
@@ -295,6 +324,8 @@ def plot_state_populations(
             file_path,
             self.labels,
         )
+        attr_dict = {"file_path": file_path}
+        logger.log(f"plot_state_populations() with args {attr_dict}.")
 
     def plot_sankey(
         self,
@@ -307,6 +338,8 @@ def plot_sankey(
             self.labels,
             frame_list,
         )
+        attr_dict = {"file_path": file_path, "frame_list": frame_list}
+        logger.log(f"plot_state_populations() with args {attr_dict}.")
 
     def dump_colored_trj(self, trj: Trj, file_path: Path) -> None:
         """Save an .xyz file with the clustering labels for each atom."""
@@ -322,6 +355,7 @@ def dump_colored_trj(self, trj: Trj, file_path: Path) -> None:
                 f"atoms, {self.labels.shape[0]} frames, but has {n_atoms} "
                 f"atoms, {n_frames} frames."
             )
+            logger.log(msg)
             raise ValueError(msg)
 
         with file_path.open("w") as f:
@@ -335,3 +369,4 @@ def dump_colored_trj(self, trj: Trj, file_path: Path) -> None:
                         f"{trj.universe.atoms[atom_idx].name} {x:.5f}"
                         f" {y:.5f} {z:.5f} {label}\n"
                     )
+        logger.log(f"Colored trj saved to {file_path}.")
diff --git a/src/dynsight/_internal/trajectory/insight.py b/src/dynsight/_internal/trajectory/insight.py
index 27bbbb3f..181b30d9 100644
--- a/src/dynsight/_internal/trajectory/insight.py
+++ b/src/dynsight/_internal/trajectory/insight.py
@@ -14,6 +14,7 @@
     from dynsight.trajectory import Trj
 
 import dynsight
+from dynsight.logs import logger
 from dynsight.trajectory import OnionInsight, OnionSmoothInsight
 
 UNIVAR_DIM = 2
@@ -37,6 +38,7 @@ def dump_to_json(self, file_path: Path) -> None:
         data["dataset"] = data["dataset"].tolist()
         with file_path.open("w") as file:
             json.dump(data, file, indent=4)
+        logger.log(f"Insight saved to {file_path}.")
 
     @classmethod
     def load_from_json(cls, file_path: Path) -> Insight:
@@ -50,8 +52,10 @@ def load_from_json(cls, file_path: Path) -> Insight:
 
         if "dataset" not in data:
             msg = "'dataset' key not found in JSON file."
+            logger.log(msg)
             raise ValueError(msg)
 
+        logger.log(f"Insight loaded from {file_path}.")
         return cls(
             dataset=np.array(data.get("dataset"), dtype=np.float64),
             meta=data.get("meta"),
@@ -77,9 +81,12 @@ def spatial_average(
             trajslice=trj.trajslice,
             num_processes=num_processes,
         )
+        attr_dict = {"sp_av_r_cut": r_cut, "selection": selection}
+
+        logger.log(f"Computed spatial average with args {attr_dict}.")
         return Insight(
             dataset=averaged_dataset,
-            meta={"sp_av_r_cut": r_cut, "selection": selection},
+            meta=attr_dict,
         )
 
     def get_time_correlation(
@@ -87,6 +94,10 @@ def get_time_correlation(
         max_delay: int | None = None,
     ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
         """Self time correlation function of the time-series signal."""
+        attr_dict = {"max_delay": max_delay}
+        logger.log(
+            f"Computed time corrleation function with args {attr_dict}."
+        )
         return dynsight.analysis.self_time_correlation(
             self.dataset,
             max_delay,
@@ -100,10 +111,13 @@ def get_angular_velocity(self, delay: int = 1) -> Insight:
         """
         if self.dataset.ndim != UNIVAR_DIM + 1:
             msg = "dataset.ndim != 3."
+            logger.log(msg)
             raise ValueError(msg)
         theta = dynsight.soap.timesoap(self.dataset, delay=delay)
         attr_dict = self.meta.copy()
         attr_dict.update({"delay": delay})
+
+        logger.log(f"Computed angular velocity with args {attr_dict}.")
         return Insight(dataset=theta, meta=attr_dict)
 
     def get_onion(
@@ -136,16 +150,18 @@ def get_onion(
             )
 
         onion_clust.fit(reshaped_data)
+        attr_dict = {
+            "delta_t": delta_t,
+            "bins": bins,
+            "number_of_sigmas": number_of_sigmas,
+        }
 
+        logger.log(f"Performed onion clustering with args {attr_dict}.")
         return OnionInsight(
             labels=onion_clust.labels_,
             state_list=onion_clust.state_list_,
             reshaped_data=reshaped_data,
-            meta={
-                "delta_t": delta_t,
-                "bins": bins,
-                "number_of_sigmas": number_of_sigmas,
-            },
+            meta=attr_dict,
         )
 
     def get_onion_smooth(
@@ -175,16 +191,18 @@ def get_onion_smooth(
             )
 
         onion_clust.fit(self.dataset)
-
+        attr_dict = {
+            "delta_t": delta_t,
+            "bins": bins,
+            "number_of_sigmas": number_of_sigmas,
+            "max_area_overlap": max_area_overlap,
+        }
+
+        logger.log(f"Performed onion clustering smooth with args {attr_dict}.")
         return OnionSmoothInsight(
             labels=onion_clust.labels_,
             state_list=onion_clust.state_list_,
-            meta={
-                "delta_t": delta_t,
-                "bins": bins,
-                "number_of_sigmas": number_of_sigmas,
-                "max_area_overlap": max_area_overlap,
-            },
+            meta=attr_dict,
         )
 
     def get_onion_analysis(
@@ -264,4 +282,18 @@ def get_onion_analysis(
                 fig2_path, list_of_pop, tra
             )
 
+        attr_dict = {
+            "delta_t_min": delta_t_min,
+            "delta_t_max": delta_t_max,
+            "delta_t_num": delta_t_num,
+            "fig1_path": fig1_path,
+            "fig2_path": fig2_path,
+            "bins": bins,
+            "number_of_sigmas": number_of_sigmas,
+            "max_area_overlap": max_area_overlap,
+        }
+
+        logger.log(
+            f"Performed full onion clustering analysis with args {attr_dict}."
+        )
         return delta_t_list, n_clust, unclass_frac
diff --git a/src/dynsight/_internal/trajectory/trajectory.py b/src/dynsight/_internal/trajectory/trajectory.py
index 58156fcf..bc80e888 100644
--- a/src/dynsight/_internal/trajectory/trajectory.py
+++ b/src/dynsight/_internal/trajectory/trajectory.py
@@ -15,6 +15,7 @@
 from MDAnalysis.coordinates.memory import MemoryReader
 
 import dynsight
+from dynsight.logs import logger
 from dynsight.trajectory import Insight
 
 UNIVAR_DIM = 2
@@ -35,6 +36,17 @@ class Trj:
 
     universe: MDAnalysis.Universe = field()
     trajslice: slice | None = None
+    n_atoms: int = field(init=False)
+    n_frames: int = field(init=False)
+
+    def __post_init__(self) -> None:
+        n_atoms = len(self.universe.atoms)
+        if self.trajslice is None:
+            n_frames = len(self.universe.trajectory)
+        else:
+            n_frames = sum(1 for _ in self.universe.trajectory[self.trajslice])
+        object.__setattr__(self, "n_atoms", n_atoms)
+        object.__setattr__(self, "n_frames", n_frames)
 
     @classmethod
     def init_from_universe(cls, universe: MDAnalysis.Universe) -> Trj:
@@ -42,6 +54,7 @@ def init_from_universe(cls, universe: MDAnalysis.Universe) -> Trj:
 
         See https://docs.mdanalysis.org/2.9.0/documentation_pages/core/universe.html#MDAnalysis.core.universe.Universe.
         """
+        logger.log("Created Trj from MDAnalysis.Universe.")
         return Trj(universe)
 
     @classmethod
@@ -53,6 +66,7 @@ def init_from_xyz(cls, traj_file: Path, dt: float) -> Trj:
         Parameters:
         dt: the trajectory's time-step.
         """
+        logger.log(f"Created Trj from {traj_file} with dt = {dt}.")
         universe = MDAnalysis.Universe(traj_file, dt=dt)
         return Trj(universe)
 
@@ -62,6 +76,7 @@ def init_from_xtc(cls, traj_file: Path, topo_file: Path) -> Trj:
 
         See https://docs.mdanalysis.org/2.9.0/documentation_pages/core/universe.html#MDAnalysis.core.universe.Universe.
         """
+        logger.log(f"Created Trj from {traj_file}, {topo_file}.")
         universe = MDAnalysis.Universe(topo_file, traj_file)
         return Trj(universe)
 
@@ -73,6 +88,8 @@ def get_coordinates(self, selection: str) -> NDArray[np.float64]:
         atoms = self.universe.select_atoms(selection)
         trajslice = slice(None) if self.trajslice is None else self.trajslice
 
+        attr_dict = {"selection": selection}
+        logger.log(f"Extracted coordinates array with args {attr_dict}.")
         return np.array(
             [
                 atoms.positions.copy()
@@ -82,6 +99,8 @@ def get_coordinates(self, selection: str) -> NDArray[np.float64]:
 
     def with_slice(self, trajslice: slice | None) -> Trj:
         """Returns a Trj with a different frames' slice."""
+        attr_dict = {"trajslice": trajslice}
+        logger.log(f"Created a sliced Trj with args {attr_dict}.")
         return Trj(self.universe, trajslice=trajslice)
 
     def get_slice(self, start: int, stop: int, step: int) -> Trj:
@@ -104,6 +123,8 @@ def get_slice(self, start: int, stop: int, step: int) -> Trj:
         u_new = MDAnalysis.Universe(topology=self.universe._topology)  # noqa: SLF001
         u_new.trajectory = mem_reader
 
+        attr_dict = {"start": start, "stop": stop, "step": step}
+        logger.log(f"Created a sliced Trj with args {attr_dict}.")
         return Trj(u_new)
 
     def get_coord_number(
@@ -128,9 +149,12 @@ def get_coord_number(
                 trajslice=self.trajslice,
             )
         _, nn, *_ = dynsight.lens.neighbour_change_in_time(neigcounts)
+
+        attr_dict = {"r_cut": r_cut, "selection": selection}
+        logger.log(f"Computed coord_number using args {attr_dict}.")
         return neigcounts, Insight(
             dataset=nn.astype(np.float64),
-            meta={"r_cut": r_cut, "selection": selection},
+            meta=attr_dict,
         )
 
     def get_lens(
@@ -155,9 +179,13 @@ def get_lens(
                 trajslice=self.trajslice,
             )
         lens, *_ = dynsight.lens.neighbour_change_in_time(neigcounts)
+
+        attr_dict = {"r_cut": r_cut, "selection": selection}
+        logger.log(f"Computed LENS using args {attr_dict}.")
+
         return neigcounts, Insight(
             dataset=lens[:, 1:],
-            meta={"r_cut": r_cut, "selection": selection},
+            meta=attr_dict,
         )
 
     def get_soap(
@@ -194,6 +222,7 @@ def get_soap(
             "selection": selection,
             "centers": centers,
         }
+        logger.log(f"Computed SOAP with args {attr_dict}.")
         return Insight(dataset=soap, meta=attr_dict)
 
     def get_rdf(
@@ -234,4 +263,5 @@ def get_rdf(
             "nbins": nbins,
             "norm": norm,
         }
+        logger.log(f"Computed g(r) with args {attr_dict}.")
         return Insight(dataset=dataset, meta=attr_dict)
diff --git a/src/dynsight/_internal/utilities/utilities.py b/src/dynsight/_internal/utilities/utilities.py
index 442f5125..8d521bfc 100644
--- a/src/dynsight/_internal/utilities/utilities.py
+++ b/src/dynsight/_internal/utilities/utilities.py
@@ -1,9 +1,16 @@
-from typing import Literal
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Literal
+
+if TYPE_CHECKING:
+    from pathlib import Path
 
 import numpy as np
 import numpy.typing as npt
 from scipy.signal import find_peaks
 
+from dynsight.trajectory import Insight, Trj
+
 
 def normalize_array(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
     """Normalizes the further axis of the given array.
@@ -56,3 +63,42 @@ def find_extrema_points(
     minima_xcoord = x_axis[peaks]
     minima_ycoord = y_axis[peaks]
     return np.column_stack((minima_xcoord, minima_ycoord))
+
+
+def load_or_compute_soap(
+    trj: Trj,
+    r_cut: float,
+    n_max: int,
+    l_max: int,
+    selection: str = "all",
+    centers: str = "all",
+    respect_pbc: bool = True,
+    n_core: int = 1,
+    soap_path: Path | None = None,
+) -> Insight:
+    """Load or compute SOAP.
+
+    If a valid path to a .json file with a SOAP Insight is provided, that
+    Insight is loaded and returned. Otherwise, the Insight is computed from
+    the Trj.
+
+    Returns:
+        Insight object containing SOAP descriptors.
+    """
+    if soap_path is not None and soap_path.exists():
+        return Insight.load_from_json(soap_path)
+
+    soap = trj.get_soap(
+        r_cut=r_cut,
+        n_max=n_max,
+        l_max=l_max,
+        selection=selection,
+        centers=centers,
+        respect_pbc=respect_pbc,
+        n_core=n_core,
+    )
+
+    if soap_path is not None:
+        soap.dump_to_json(soap_path)
+
+    return soap
diff --git a/src/dynsight/logs.py b/src/dynsight/logs.py
new file mode 100644
index 00000000..00128318
--- /dev/null
+++ b/src/dynsight/logs.py
@@ -0,0 +1,7 @@
+"""logging package."""
+
+from dynsight._internal.logs import logger
+
+__all__ = [
+    "logger",
+]
diff --git a/src/dynsight/utilities.py b/src/dynsight/utilities.py
index 081bdea1..b1315855 100644
--- a/src/dynsight/utilities.py
+++ b/src/dynsight/utilities.py
@@ -2,10 +2,12 @@
 
 from dynsight._internal.utilities.utilities import (
     find_extrema_points,
+    load_or_compute_soap,
     normalize_array,
 )
 
 __all__ = [
     "find_extrema_points",
+    "load_or_compute_soap",
     "normalize_array",
 ]
diff --git a/tests/trajectory/test_trj.py b/tests/trajectory/test_trj.py
index 1b504037..a845fd7e 100644
--- a/tests/trajectory/test_trj.py
+++ b/tests/trajectory/test_trj.py
@@ -8,6 +8,7 @@
 import numpy as np
 import pytest
 
+from dynsight.logs import logger
 from dynsight.trajectory import (
     ClusterInsight,
     Insight,
@@ -16,6 +17,8 @@
     Trj,
 )
 
+TRJ_SHAPE = (2, 21)
+
 # ---------------- Fixtures ----------------
 
 
@@ -61,6 +64,10 @@ def test_trj_inits(
     trj_4 = Trj.init_from_xtc(file_paths["xtc"], file_paths["gro"])
     assert len(trj_4.universe.trajectory) == n_frames_xtc
 
+    assert trj_1.n_atoms, trj_1.n_frames == TRJ_SHAPE
+
+    logger.get()
+
 
 def test_insight(
     tmp_path: Path, file_paths: dict[str, Path], universe: MDAnalysis.Universe
@@ -103,6 +110,8 @@ def test_insight(
     loaded_onion_smooth = OnionSmoothInsight.load_from_json(onion_smooth_json)
     assert loaded_onion_smooth is not None
 
+    logger.get()
+
 
 def test_onion_analysis(universe: MDAnalysis.Universe) -> None:
     """Test the onion clustering complete analysis tool."""
@@ -132,3 +141,5 @@ def test_insight_load_errors(file_paths: dict[str, Path]) -> None:
         _ = OnionInsight.load_from_json(
             file_paths["files_dir"] / "cl_ins_test.json"
         )
+
+    logger.get()