diff --git a/dynamics_learning/dynamics_learning/preprocessing/dataset_analysis.py b/dynamics_learning/dynamics_learning/preprocessing/dataset_analysis.py
index ff50461f4aead6dd61349b26eb61280e848a49c7..c45524fc2c805a27d27c9015a01086f56067f4e8 100644
--- a/dynamics_learning/dynamics_learning/preprocessing/dataset_analysis.py
+++ b/dynamics_learning/dynamics_learning/preprocessing/dataset_analysis.py
@@ -1,154 +1,177 @@
-
-# Import and variables
+# Import necessary libraries and modules
 import sys
-sys.path.append("src")
-sys.path.append("..")
+sys.path.append("src")  # Add 'src' directory to system path
+sys.path.append("..")   # Add parent directory to system path
 import numpy as np
 import matplotlib.pyplot as plt
 import pandas as pd
 import rwth_style
-import panda_limits 
+import panda_limits
 import os
 from scipy.optimize import curve_fit
 from sklearn.metrics import r2_score
 from rich import print
 from rich.progress import track
 
-def gauss(x, *p):
+# Define a Gaussian function for curve fitting
+def gauss(x: np.ndarray, *p: float) -> np.ndarray:
+    """
+    Gaussian function for curve fitting.
+    
+    Args:
+        x (np.ndarray): The input data.
+        p (float): Parameters A, mu, sigma.
+
+    Returns:
+        np.ndarray: The Gaussian function applied to the input data.
+    """
     A, mu, sigma = p
-    return A*np.exp(-(x-mu)**2/(2.*sigma**2))
+    return A * np.exp(-(x - mu)**2 / (2. * sigma**2))
 
+# Initial guess for Gaussian parameters
 p0 = [1, 0, 1]
 
-def analize():
-    for directory in ["dataset_v1", "dataset_v2", "dataset_v3"]:
-        for dataset in ["train", "test"]:
-            attained_freqs = np.empty([0,1])
+def analyze() -> None:
+    """
+    Analyzes datasets in specified directories, generates statistics, and creates plots.
+
+    This function processes datasets, performs statistical analysis, and saves the results
+    and plots for each dataset.
+
+    Returns:
+        None
+    """
+    directories = ["dataset_v1", "dataset_v2", "dataset_v3"]
+    datasets = ["train", "test"]
+    
+    for directory in directories:
+        for dataset in datasets:
+            # Initialize empty arrays for data collection
+            attained_freqs = np.empty([0, 1])
             duration_meas = np.empty(0, np.float32)
             num_samples_meas = np.empty(0, np.float32)
-            q_command = np.empty((0,7), np.float32)
-            qd_command = np.empty((0,7), np.float32)
-            qdd_command = np.empty((0,7), np.float32)
-            q_meas = np.empty((0,7), np.float32)
-            qd_meas = np.empty((0,7), np.float32)
-            tau_meas = np.empty((0,7), np.float32)
-            directory_in = "data/{}/{}".format(directory, dataset)
+            q_command = np.empty((0, 7), np.float32)
+            qd_command = np.empty((0, 7), np.float32)
+            qdd_command = np.empty((0, 7), np.float32)
+            q_meas = np.empty((0, 7), np.float32)
+            qd_meas = np.empty((0, 7), np.float32)
+            tau_meas = np.empty((0, 7), np.float32)
+            
+            directory_in = f"data/{directory}/{dataset}"
             file_list = list(filter(lambda k: 'meas.csv' in k, sorted(os.listdir(directory_in))))
             
             for filename in track(file_list):
-                    data2 = np.genfromtxt(os.path.join(directory_in, filename), dtype=float, delimiter=',') 
-                    t_meas = data2[:,0]
-                    time_diffs = np.diff(t_meas)
-                    attained_freq = np.mean(1/time_diffs)
-                    attained_freqs = np.append(attained_freqs, attained_freq)
-                    num_samples_meas = np.concatenate((num_samples_meas, [data2.shape[0]]))
-                    duration_meas = np.concatenate((duration_meas, [data2[-1, 0]]), axis=0)
-                    q_meas = np.concatenate((q_meas, data2[:, 1:8]), axis=0)
-                    qd_meas = np.concatenate((qd_meas, data2[:, 8:15]), axis=0)
-                    tau_meas = np.concatenate((tau_meas, data2[:, 15:22]), axis=0)
-                    filename2 = filename.replace("meas", "com")
-                    data1 = np.genfromtxt(os.path.join(directory_in, filename2), dtype=float, delimiter=',') 
-                    q_command = np.concatenate((q_command, data1[:,1:8]), axis=0)
-                    qd_command = np.concatenate((qd_command, data1[:, 8:15]), axis=0)
-                    qdd_command = np.concatenate((qdd_command, data1[:, 15:22]), axis=0)
+                # Read measurement data
+                data2 = np.genfromtxt(os.path.join(directory_in, filename), dtype=float, delimiter=',')
+                t_meas = data2[:, 0]
+                time_diffs = np.diff(t_meas)
+                attained_freq = np.mean(1 / time_diffs)
+                attained_freqs = np.append(attained_freqs, attained_freq)
+                num_samples_meas = np.concatenate((num_samples_meas, [data2.shape[0]]))
+                duration_meas = np.concatenate((duration_meas, [data2[-1, 0]]), axis=0)
+                q_meas = np.concatenate((q_meas, data2[:, 1:8]), axis=0)
+                qd_meas = np.concatenate((qd_meas, data2[:, 8:15]), axis=0)
+                tau_meas = np.concatenate((tau_meas, data2[:, 15:22]), axis=0)
+                
+                # Read command data
+                filename2 = filename.replace("meas", "com")
+                data1 = np.genfromtxt(os.path.join(directory_in, filename2), dtype=float, delimiter=',')
+                q_command = np.concatenate((q_command, data1[:, 1:8]), axis=0)
+                qd_command = np.concatenate((qd_command, data1[:, 8:15]), axis=0)
+                qdd_command = np.concatenate((qdd_command, data1[:, 15:22]), axis=0)
             
-            df = pd.DataFrame({"# trajecotries":[duration_meas.size], 
-                            "Duration Sum [s]": np.sum(duration_meas),
-                            "Duration Min [s]": np.min(duration_meas),
-                            "Duration Max [s]":np.max(duration_meas),
-                            "Duration Mean [s]":np.mean(duration_meas),
-                            "# Samples Sum":np.sum(num_samples_meas),
-                            "# Samples Min":np.min(num_samples_meas),
-                            "# Samples Max":np.max(num_samples_meas),
-                            "# Samples Mean":np.mean(num_samples_meas)
-                            })
+            # Create DataFrame for trajectory analysis
+            df = pd.DataFrame({
+                "# trajectories": [duration_meas.size],
+                "Duration Sum [s]": np.sum(duration_meas),
+                "Duration Min [s]": np.min(duration_meas),
+                "Duration Max [s]": np.max(duration_meas),
+                "Duration Mean [s]": np.mean(duration_meas),
+                "# Samples Sum": np.sum(num_samples_meas),
+                "# Samples Min": np.min(num_samples_meas),
+                "# Samples Max": np.max(num_samples_meas),
+                "# Samples Mean": np.mean(num_samples_meas)
+            })
             df.index = ['Value']
-            df1 = df.T
-            df1.to_csv('data/{}/analysis/trajectory_analysis_{}.csv'.format(directory, dataset), float_format='%.3f')  
-
+            df.T.to_csv(f'data/{directory}/analysis/trajectory_analysis_{dataset}.csv', float_format='%.3f')
+            
             # Attained frequency analysis
             freq_mu = np.mean(attained_freqs)
             freq_sigma = np.std(attained_freqs)
             freq_worst_case = np.max(freq_mu - attained_freqs)
+            
             with plt.style.context("default"):
-                num_col = 1
-                num_row = 1
-                fig, axs = plt.subplots(num_row, num_col, figsize=(10,4))
-                count, bins, ignored = axs.hist(attained_freqs, density=True, bins=50, color=rwth_style.blue)
-                plt.plot(bins, 1/(freq_sigma * np.sqrt(2 * np.pi)) *
-                        np.exp( - (bins - freq_mu)**2 / (2 * freq_sigma**2) ),linewidth=2, color=rwth_style.red)
-                axs.grid()
-                axs.set_title(r'$\mathrm{Histogram\ of\ Frequency:}\ \mu=%.3f Hz,\ \sigma=%.3f$ Hz' %(freq_mu, freq_sigma))
-                axs.set_ylabel("Probability")
-                axs.set_xlabel("Frequency")
+                fig, ax = plt.subplots(figsize=(10, 4))
+                ax.hist(attained_freqs, density=True, bins=50, color='blue')
+                bins = np.linspace(min(attained_freqs), max(attained_freqs), 50)
+                ax.plot(bins, 1 / (freq_sigma * np.sqrt(2 * np.pi)) * np.exp(-(bins - freq_mu)**2 / (2 * freq_sigma**2)), linewidth=2, color='red')
+                ax.grid()
+                ax.set_title(f'Histogram of Frequency: μ={freq_mu:.3f} Hz, σ={freq_sigma:.3f} Hz')
+                ax.set_ylabel("Probability")
+                ax.set_xlabel("Frequency")
                 fig.tight_layout()
-                fig.savefig("data/{}/analysis/hist_freq_{}.png".format(directory, dataset))
+                fig.savefig(f"data/{directory}/analysis/hist_freq_{dataset}.png")
                 plt.close(fig)
-
+            
+            # Feature analysis
             values = [q_command, qd_command, qdd_command, tau_meas]
             names = ["q", "qd", "qdd", "tau"]
-            labels_y = ["Position [rad]", "Velocity [rad-s]", "Acceleration [rad-s^2]", "Torque [Nm]"]
-            units = ["[rad]", "[rad-s]", "[rad-s^2]", "[Nm]"]
+            labels_y = ["Position [rad]", "Velocity [rad/s]", "Acceleration [rad/s^2]", "Torque [Nm]"]
+            units = ["[rad]", "[rad/s]", "[rad/s^2]", "[Nm]"]
             titles = ["Command Position", "Command Velocity", "Command Acceleration", "Measured Torque"]
             labels = ['Axis 1', 'Axis 2', 'Axis 3', 'Axis 4', 'Axis 5', 'Axis 6', 'Axis 7']
-            # phyiscal limits
+            
             min_limits_phys = [panda_limits.q_lim_min_phys, panda_limits.qd_lim_min_phys, panda_limits.qdd_lim_min_phys, panda_limits.tau_lim_min_phys]
             max_limits_phys = [panda_limits.q_lim_max_phys, panda_limits.qd_lim_max_phys, panda_limits.qdd_lim_max_phys, panda_limits.tau_lim_max_phys]
-            # Moveit limits
             min_limits_moveit = [panda_limits.q_lim_min_moveit, panda_limits.qd_lim_min_moveit, panda_limits.qdd_lim_min_moveit, panda_limits.tau_lim_min_moveit]
             max_limits_moveit = [panda_limits.q_lim_max_moveit, panda_limits.qd_lim_max_moveit, panda_limits.qdd_lim_max_moveit, panda_limits.tau_lim_max_moveit]
-            length= len(values)
-
-            for i in range(length):        
+            
+            for i in range(len(values)):
                 value = values[i]
                 std = np.std(value, axis=0)
                 mean = np.mean(value, axis=0)
-                min = np.min(value, axis=0)
-                max = np.max(value, axis=0)
-                phy_mean = (max_limits_phys[i]+min_limits_phys[i])/2
-                cov_max_phys  = (max-min) / (max_limits_phys[i]-min_limits_phys[i]) * 100 # Percentage of used working space
-                cov_std_phys  = (2*std) / (max_limits_phys[i]-min_limits_phys[i]) * 100 # Percentage of used working space based on standard deviaton
-                cov_max_moveit  = (max-min) / (max_limits_moveit[i]-min_limits_moveit[i]) * 100 # Percentage of used working space
-                cov_std_moveit  = (2*std) / (max_limits_moveit[i]-min_limits_moveit[i]) * 100 # Percentage of used working space based on standard deviaton
+                min_val = np.min(value, axis=0)
+                max_val = np.max(value, axis=0)
+                phy_mean = (max_limits_phys[i] + min_limits_phys[i]) / 2
+                cov_max_phys = (max_val - min_val) / (max_limits_phys[i] - min_limits_phys[i]) * 100
+                cov_std_phys = (2 * std) / (max_limits_phys[i] - min_limits_phys[i]) * 100
+                cov_max_moveit = (max_val - min_val) / (max_limits_moveit[i] - min_limits_moveit[i]) * 100
+                cov_std_moveit = (2 * std) / (max_limits_moveit[i] - min_limits_moveit[i]) * 100
                 
-                df = pd.DataFrame({"Att. Mean {}".format(units[i]):mean, 
-                                "Phys. Mean {}".format(units[i]):phy_mean,
-                                "MoveIt Mean {}".format(units[i]):phy_mean,
-                                
-                                "Att. Min {}".format(units[i]):min,
-                                "Phys. Min {}".format(units[i]):min_limits_phys[i],
-                                "MoveIt Min {}".format(units[i]):min_limits_moveit[i],
-                                                            
-                                "Att. Max {}".format(units[i]):max,
-                                "Phys. Max {}".format(units[i]):max_limits_phys[i],
-                                "MoveIt Max {}".format(units[i]):max_limits_moveit[i],
-                                
-                                "Phy. Coverage Min/Max [\\%]":cov_max_phys ,
-                                "MoveIt Coverage Min/Max [\\%]":cov_max_moveit ,
-                                
-                                "Att. Std {}".format(units[i]):std,
-                                "Phy. Coverage Std [\\%]":cov_std_phys ,                           
-                                "MoveIt Coverage Std [\\%]":cov_std_moveit ,
-                                })
-                df.index = ['Axis_1', 'Axis_2', 'Axis_3', 'Axis_4', 'Axis_5', 'Axis_6', 'Axis_7']
-                df1 = df.T
-                df1.to_csv('data/{}/analysis/feature_analysis_{}_{}.csv'.format(directory, dataset, names[i]), float_format='%.3f')  
+                df = pd.DataFrame({
+                    f"Att. Mean {units[i]}": mean,
+                    f"Phys. Mean {units[i]}": phy_mean,
+                    f"MoveIt Mean {units[i]}": phy_mean,
+                    f"Att. Min {units[i]}": min_val,
+                    f"Phys. Min {units[i]}": min_limits_phys[i],
+                    f"MoveIt Min {units[i]}": min_limits_moveit[i],
+                    f"Att. Max {units[i]}": max_val,
+                    f"Phys. Max {units[i]}": max_limits_phys[i],
+                    f"MoveIt Max {units[i]}": max_limits_moveit[i],
+                    "Phy. Coverage Min/Max [%]": cov_max_phys,
+                    "MoveIt Coverage Min/Max [%]": cov_max_moveit,
+                    f"Att. Std {units[i]}": std,
+                    "Phy. Coverage Std [%]": cov_std_phys,
+                    "MoveIt Coverage Std [%]": cov_std_moveit,
+                })
+                df.index = ['Axis 1', 'Axis 2', 'Axis 3', 'Axis 4', 'Axis 5', 'Axis 6', 'Axis 7']
+                df.T.to_csv(f'data/{directory}/analysis/feature_analysis_{dataset}_{names[i]}.csv', float_format='%.3f')
                 
+                # Create bar plots for limits and measured values
                 x_pos = np.arange(len(labels))
                 with plt.style.context("default"):
-                    f4, ax = plt.subplots(figsize=(10,5))
+                    fig, ax = plt.subplots(figsize=(10, 5))
                     bar1_alpha = 1
                     bar1_hatch = ''
-                    plt.rcParams.update({'hatch.color': rwth_style.blue_light})
-                    if (i==2 or i==0):
-                        ax.bar(x_pos, max_limits_moveit[i] - min_limits_moveit[i], 0.4, bottom = min_limits_moveit[i], color=rwth_style.blue_light, label='MoveIt Limits', zorder=3)
+                    if i == 2 or i == 0:
+                        ax.bar(x_pos, max_limits_moveit[i] - min_limits_moveit[i], 0.4, bottom=min_limits_moveit[i], color='lightblue', label='MoveIt Limits', zorder=3)
                         bar1_alpha = 0.5
                         bar1_hatch = '//'
-                    bar1 = ax.bar(x_pos, max_limits_phys[i] - min_limits_phys[i], 0.4, bottom = min_limits_phys[i], color=rwth_style.blue_light, alpha=bar1_alpha, hatch=bar1_hatch, label='Physical Limits', zorder=2)
-                    ax.scatter(x_pos, (max_limits_phys[i] + min_limits_phys[i])/2, c=[rwth_style.blue], s=400, marker='_', label='Physical Mean', zorder=4)
-                    errorbar_1 = ax.errorbar(x_pos, mean, std, fmt='o', capsize=5, label='Mean/Std Attained {}'.format(titles[i]),  c=rwth_style.blue, zorder=5)
-                    scatter_min = ax.scatter(x_pos, min, c=[rwth_style.blue], s=20, marker='x', label='Min/Max Attained {}'.format(titles[i]), zorder=6)
-                    ax.scatter(x_pos, max, c=[rwth_style.blue], s=20, marker='x', zorder=6)
+                    ax.bar(x_pos, max_limits_phys[i] - min_limits_phys[i], 0.4, bottom=min_limits_phys[i], color='lightblue', alpha=bar1_alpha, hatch=bar1_hatch, label='Physical Limits', zorder=2)
+                    ax.scatter(x_pos, (max_limits_phys[i] + min_limits_phys[i]) / 2, c='blue', s=400, marker='_', label='Physical Mean', zorder=4)
+                    ax.errorbar(x_pos, mean, std, fmt='o', capsize=5, label=f'Mean/Std Attained {titles[i]}', c='blue', zorder=5)
+                    ax.scatter(x_pos, min_val, c='blue', s=20, marker='x', label=f'Min/Max Attained {titles[i]}', zorder=6)
+                    ax.scatter(x_pos, max_val, c='blue', s=20, marker='x', zorder=6)
                     ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.08), ncol=2)
                     box = ax.get_position()
                     ax.set_position([box.x0, box.y0 + box.height * 0.15, box.width, box.height])
@@ -158,37 +181,34 @@ def analize():
                     ax.set_xticklabels(labels)
                     ax.yaxis.grid(True)
                     ax.set_axisbelow(True)
-                    f4.savefig("data/{}/analysis/errorbar_{}_{}.png".format(directory, dataset, names[i]))
-                    plt.close(f4)
-
-                with plt.style.context("default"):
-                    num_col = 3
-                    num_row = 3
-                    fig, axs = plt.subplots(num_row, num_col, figsize=(20,10))
+                    fig.savefig(f"data/{directory}/analysis/errorbar_{dataset}_{names[i]}.png")
+                    plt.close(fig)
                 
+                # Create histograms for the values
+                with plt.style.context("default"):
+                    fig, axs = plt.subplots(3, 3, figsize=(20, 10))
                     for k, label in enumerate(labels):
-                        count, bins, ignored = axs[k//num_col, k%num_col].hist(value[:,k], density=True, bins=40, color=rwth_style.blue)
+                        count, bins, _ = axs[k // 3, k % 3].hist(value[:, k], density=True, bins=40, color='blue')
                         bins = bins[:-1]
                         try:
-                            coeff, var_matrix = curve_fit(gauss, bins, count, p0=p0, maxfev = 2000)
-                            # Get the fitted curve
+                            coeff, _ = curve_fit(gauss, bins, count, maxfev=2000)
                             hist_fit = gauss(bins, *coeff)
                             r_squared = r2_score(count, hist_fit)
-                            axs[k//num_col, k%num_col].plot(bins,hist_fit, color=rwth_style.red)
-                            t = axs[k//num_col, k%num_col].text(0.15,0.85, f"mean: {round(coeff[1],3)}\nstd: {round(coeff[2],3)}\nr squared: {round(round(r_squared,4)*100,2)}%", horizontalalignment='center', verticalalignment='center', transform=axs[k//num_col, k%num_col].transAxes)
+                            axs[k // 3, k % 3].plot(bins, hist_fit, color='red')
+                            t = axs[k // 3, k % 3].text(0.15, 0.85, f"mean: {round(coeff[1], 3)}\nstd: {round(coeff[2], 3)}\nr squared: {round(r_squared * 100, 2)}%", horizontalalignment='center', verticalalignment='center', transform=axs[k // 3, k % 3].transAxes)
                             t.set_bbox(dict(facecolor='white', alpha=0.5, edgecolor='black'))
                         except Exception as e:
                             print(e)
-                        axs[k//num_col, k%num_col].grid()
-                        axs[k//num_col, k%num_col].set_title(labels[k])
-                        axs[k//num_col, k%num_col].set_ylabel("Probability Density")
-                        axs[k//num_col, k%num_col].set_xlabel(labels_y[i]) 
-                        axs[k//num_col, k%num_col].set_xlim([min_limits_phys[i][k],max_limits_phys[i][k]])
-                        
-                    fig.delaxes(axs[-1,-1])
-                    fig.delaxes(axs[-1,-2])
+                        axs[k // 3, k % 3].grid()
+                        axs[k // 3, k % 3].set_title(label)
+                        axs[k // 3, k % 3].set_ylabel("Probability Density")
+                        axs[k // 3, k % 3].set_xlabel(labels_y[i])
+                        axs[k // 3, k % 3].set_xlim([min_limits_phys[i][k], max_limits_phys[i][k]])
+                    
+                    fig.delaxes(axs[-1, -1])
+                    fig.delaxes(axs[-1, -2])
                     fig.tight_layout()
-                    fig.savefig("data/{}/analysis/hist_{}_{}.png".format(directory, dataset, names[i]))
+                    fig.savefig(f"data/{directory}/analysis/hist_{dataset}_{names[i]}.png")
                     plt.close(fig)