diff --git a/src/odatse/_info.py b/src/odatse/_info.py
index 37b0caab..03e3663f 100644
--- a/src/odatse/_info.py
+++ b/src/odatse/_info.py
@@ -26,18 +26,37 @@
 
 
 class Info:
+    """
+    A class to represent the information structure for the data-analysis software.
+    """
+
     base: dict
     algorithm: dict
     solver: dict
     runner: dict
 
     def __init__(self, d: Optional[MutableMapping] = None):
+        """
+        Initialize the Info object.
+
+        Parameters:
+        d (Optional[MutableMapping]): A dictionary to initialize the Info object.
+        """
         if d is not None:
             self.from_dict(d)
         else:
             self._cleanup()
 
     def from_dict(self, d: MutableMapping) -> None:
+        """
+        Initialize the Info object from a dictionary.
+
+        Parameters:
+        d (MutableMapping): A dictionary containing the information to initialize the Info object.
+
+        Raises:
+        exception.InputError: If any required section is missing in the input dictionary.
+        """
         for section in ["base", "algorithm", "solver"]:
             if section not in d:
                 raise exception.InputError(
@@ -58,6 +77,9 @@ def from_dict(self, d: MutableMapping) -> None:
         )
 
     def _cleanup(self) -> None:
+        """
+        Reset the Info object to its default state.
+        """
         self.base = {}
         self.base["root_dir"] = Path(".").absolute()
         self.base["output_dir"] = self.base["root_dir"]
@@ -67,6 +89,20 @@ def _cleanup(self) -> None:
 
     @classmethod
     def from_file(cls, file_name, fmt="", **kwargs):
+        """
+        Create an Info object from a file.
+
+        Parameters:
+        file_name (str): The name of the file to load the information from.
+        fmt (str): The format of the file (default is "").
+        **kwargs: Additional keyword arguments.
+
+        Returns:
+        Info: An Info object initialized with the data from the file.
+
+        Raises:
+        ValueError: If the file format is unsupported.
+        """
         if fmt == "toml" or fnmatch(file_name.lower(), "*.toml"):
             inp = {}
             if mpi.rank() == 0:
@@ -75,4 +111,4 @@ def from_file(cls, file_name, fmt="", **kwargs):
                 inp = mpi.comm().bcast(inp, root=0)
             return cls(inp)
         else:
-            raise ValueError("unsupported file format: {}".format(file_name))
+            raise ValueError("unsupported file format: {}".format(file_name))
\ No newline at end of file
diff --git a/src/odatse/_main.py b/src/odatse/_main.py
index 9517798d..0bb62619 100644
--- a/src/odatse/_main.py
+++ b/src/odatse/_main.py
@@ -22,6 +22,11 @@
 
 
 def main():
+    """
+        Main function to run the data-analysis software for quantum beam diffraction experiments
+        on 2D material structures. It parses command-line arguments, loads the input file,
+        selects the appropriate algorithm and solver, and executes the analysis.
+    """
     import argparse
 
     parser = argparse.ArgumentParser(
diff --git a/src/odatse/_runner.py b/src/odatse/_runner.py
index 968cd489..a84c9732 100644
--- a/src/odatse/_runner.py
+++ b/src/odatse/_runner.py
@@ -36,14 +36,16 @@
 class Run(metaclass=ABCMeta):
     def __init__(self, nprocs=None, nthreads=None, comm=None):
         """
+        Initialize the Run class.
+
         Parameters
         ----------
         nprocs : int
-            Number of process which one solver uses
+            Number of processes which one solver uses.
         nthreads : int
-            Number of threads which one solver process uses
+            Number of threads which one solver process uses.
         comm : MPI.Comm
-            MPI Communicator
+            MPI Communicator.
         """
         self.nprocs = nprocs
         self.nthreads = nthreads
@@ -51,6 +53,14 @@ def __init__(self, nprocs=None, nthreads=None, comm=None):
 
     @abstractmethod
     def submit(self, solver):
+        """
+        Abstract method to submit a solver.
+
+        Parameters
+        ----------
+        solver : object
+            Solver object to be submitted.
+        """
         pass
 
 
@@ -64,10 +74,18 @@ def __init__(self,
                  mapping = None,
                  limitation = None) -> None:
         """
+        Initialize the Runner class.
 
         Parameters
         ----------
-        Solver: odatse.solver.SolverBase object
+        solver : odatse.solver.SolverBase
+            Solver object.
+        info : Optional[odatse.Info]
+            Information object.
+        mapping : object, optional
+            Mapping object.
+        limitation : object, optional
+            Limitation object.
         """
         self.solver = solver
         self.solver_name = solver.name
@@ -82,7 +100,7 @@ def __init__(self,
         else:
             # trivial mapping
             self.mapping = odatse.util.mapping.TrivialMapping()
-        
+
         if limitation is not None:
             self.limitation = limitation
         elif "limitation" in info.runner:
@@ -92,11 +110,38 @@ def __init__(self,
             self.limitation = odatse.util.limitation.Unlimited()
 
     def prepare(self, proc_dir: Path):
+        """
+        Prepare the logger with the given process directory.
+
+        Parameters
+        ----------
+        proc_dir : Path
+            Path to the process directory.
+        """
         self.logger.prepare(proc_dir)
 
     def submit(
             self, x: np.ndarray, args = (), nprocs: int = 1, nthreads: int = 1
     ) -> float:
+        """
+        Submit the solver with the given parameters.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Input array.
+        args : tuple, optional
+            Additional arguments.
+        nprocs : int, optional
+            Number of processes.
+        nthreads : int, optional
+            Number of threads.
+
+        Returns
+        -------
+        float
+            Result of the solver evaluation.
+        """
         if self.limitation.judge(x):
             xp = self.mapping(x)
             result = self.solver.evaluate(xp, args)
@@ -106,4 +151,7 @@ def submit(
         return result
 
     def post(self) -> None:
+        """
+        Write the logger data.
+        """
         self.logger.write()
diff --git a/src/odatse/algorithm/_algorithm.py b/src/odatse/algorithm/_algorithm.py
index 674a5328..5f8ce014 100644
--- a/src/odatse/algorithm/_algorithm.py
+++ b/src/odatse/algorithm/_algorithm.py
@@ -41,11 +41,14 @@
     from mpi4py import MPI
 
 class AlgorithmStatus(IntEnum):
+    """Enumeration for the status of the algorithm."""
     INIT = 1
     PREPARE = 2
     RUN = 3
 
 class AlgorithmBase(metaclass=ABCMeta):
+    """Base class for algorithms, providing common functionality and structure."""
+
     mpicomm: Optional["MPI.Comm"]
     mpisize: int
     mpirank: int
@@ -70,6 +73,13 @@ def __init__(
             runner: Optional[odatse.Runner] = None,
             run_mode: str = "initial"
     ) -> None:
+        """
+        Initialize the algorithm with the given information and runner.
+
+        :param info: Information object containing algorithm and base parameters.
+        :param runner: Optional runner object to execute the algorithm.
+        :param run_mode: Mode in which the algorithm should run.
+        """
         self.mpicomm = mpi.comm()
         self.mpisize = mpi.size()
         self.mpirank = mpi.rank()
@@ -119,6 +129,11 @@ def __init__(
             self.set_runner(runner)
 
     def __init_rng(self, info: odatse.Info) -> None:
+        """
+        Initialize the random number generator.
+
+        :param info: Information object containing algorithm parameters.
+        """
         seed = info.algorithm.get("seed", None)
         seed_delta = info.algorithm.get("seed_delta", 314159)
 
@@ -128,9 +143,17 @@ def __init_rng(self, info: odatse.Info) -> None:
             self.rng = np.random.RandomState(seed + self.mpirank * seed_delta)
 
     def set_runner(self, runner: odatse.Runner) -> None:
+        """
+        Set the runner for the algorithm.
+
+        :param runner: Runner object to execute the algorithm.
+        """
         self.runner = runner
 
     def prepare(self) -> None:
+        """
+        Prepare the algorithm for execution.
+        """
         if self.runner is None:
             msg = "Runner is not assigned"
             raise RuntimeError(msg)
@@ -139,9 +162,13 @@ def prepare(self) -> None:
 
     @abstractmethod
     def _prepare(self) -> None:
+        """Abstract method to be implemented by subclasses for preparation steps."""
         pass
 
     def run(self) -> None:
+        """
+        Run the algorithm.
+        """
         if self.status < AlgorithmStatus.PREPARE:
             msg = "algorithm has not prepared yet"
             raise RuntimeError(msg)
@@ -155,9 +182,15 @@ def run(self) -> None:
 
     @abstractmethod
     def _run(self) -> None:
+        """Abstract method to be implemented by subclasses for running steps."""
         pass
 
     def post(self) -> Dict:
+        """
+        Perform post-processing after the algorithm has run.
+
+        :return: Dictionary containing post-processing results.
+        """
         if self.status < AlgorithmStatus.RUN:
             msg = "algorithm has not run yet"
             raise RuntimeError(msg)
@@ -169,9 +202,13 @@ def post(self) -> Dict:
 
     @abstractmethod
     def _post(self) -> Dict:
+        """Abstract method to be implemented by subclasses for post-processing steps."""
         pass
 
     def main(self):
+        """
+        Main method to execute the algorithm.
+        """
         time_sta = time.perf_counter()
         self.prepare()
         time_end = time.perf_counter()
@@ -197,6 +234,11 @@ def main(self):
         return result
 
     def write_timer(self, filename: Path):
+        """
+        Write the timing information to a file.
+
+        :param filename: Path to the file where timing information will be written.
+        """
         with open(filename, "w") as fw:
             fw.write("#in units of seconds\n")
 
@@ -214,6 +256,13 @@ def output_file(type):
             output_file("post")
 
     def _save_data(self, data, filename="state.pickle", ngen=3) -> None:
+        """
+        Save data to a file with versioning.
+
+        :param data: Data to be saved.
+        :param filename: Name of the file to save the data.
+        :param ngen: Number of generations for versioning.
+        """
         try:
             fn = Path(filename + ".tmp")
             with open(fn, "wb") as f:
@@ -235,6 +284,12 @@ def _save_data(self, data, filename="state.pickle", ngen=3) -> None:
         print("save_state: write to {}".format(filename))
 
     def _load_data(self, filename="state.pickle") -> Dict:
+        """
+        Load data from a file.
+
+        :param filename: Name of the file to load the data from.
+        :return: Dictionary containing the loaded data.
+        """
         if Path(filename).exists():
             try:
                 fn = Path(filename)
@@ -250,12 +305,20 @@ def _load_data(self, filename="state.pickle") -> Dict:
         return data
 
     def _show_parameters(self):
+        """
+        Show the parameters of the algorithm.
+        """
         if self.mpirank == 0:
             info = flatten_dict(self.info)
             for k, v in info.items():
                 print("{:16s}: {}".format(k, v))
 
     def _check_parameters(self, param=None):
+        """
+        Check the parameters of the algorithm against previous parameters.
+
+        :param param: Previous parameters to check against.
+        """
         info = flatten_dict(self.info)
         info_prev = flatten_dict(param)
 
@@ -269,6 +332,14 @@ def _check_parameters(self, param=None):
 
 # utility
 def flatten_dict(d, parent_key="", separator="."):
+    """
+    Flatten a nested dictionary.
+
+    :param d: Dictionary to flatten.
+    :param parent_key: Key for the parent dictionary.
+    :param separator: Separator to use between keys.
+    :return: Flattened dictionary.
+    """
     items = []
     if d:
         for key_, val in d.items():
@@ -277,4 +348,4 @@ def flatten_dict(d, parent_key="", separator="."):
                 items.extend(flatten_dict(val, key, separator=separator).items())
             else:
                 items.append((key, val))
-    return dict(items)
+    return dict(items)
\ No newline at end of file
diff --git a/src/odatse/algorithm/bayes.py b/src/odatse/algorithm/bayes.py
index 32f25bdd..3a819f06 100644
--- a/src/odatse/algorithm/bayes.py
+++ b/src/odatse/algorithm/bayes.py
@@ -27,30 +27,52 @@
 import odatse.domain
 
 class Algorithm(odatse.algorithm.AlgorithmBase):
-
-    # inputs
-    mesh_list: np.ndarray
-    label_list: List[str]
-
-    # hyperparameters of Bayesian optimization
-    random_max_num_probes: int
-    bayes_max_num_probes: int
-    score: str
-    interval: int
-    num_rand_basis: int
-
-    # results
-    xopt: np.ndarray
-    best_fx: List[float]
-    best_action: List[int]
-    fx_list: List[float]
-    param_list: List[np.ndarray]
-
-    def __init__(self, info: odatse.Info,
-                 runner: odatse.Runner = None,
-                 domain = None,
-                 run_mode: str = "initial"
-    ) -> None:
+    """
+    A class to represent the Bayesian optimization algorithm.
+
+    Attributes
+    ----------
+    mesh_list : np.ndarray
+        The mesh grid list.
+    label_list : List[str]
+        The list of labels.
+    random_max_num_probes : int
+        The maximum number of random probes.
+    bayes_max_num_probes : int
+        The maximum number of Bayesian probes.
+    score : str
+        The scoring method.
+    interval : int
+        The interval for Bayesian optimization.
+    num_rand_basis : int
+        The number of random basis.
+    xopt : np.ndarray
+        The optimal solution.
+    best_fx : List[float]
+        The list of best function values.
+    best_action : List[int]
+        The list of best actions.
+    fx_list : List[float]
+        The list of function values.
+    param_list : List[np.ndarray]
+        The list of parameters.
+    """
+
+    def __init__(self, info: odatse.Info, runner: odatse.Runner = None, domain = None, run_mode: str = "initial") -> None:
+        """
+        Constructs all the necessary attributes for the Algorithm object.
+
+        Parameters
+        ----------
+        info : odatse.Info
+            The information object.
+        runner : odatse.Runner, optional
+            The runner object (default is None).
+        domain : optional
+            The domain object (default is None).
+        run_mode : str, optional
+            The run mode (default is "initial").
+        """
         super().__init__(info=info, runner=runner, run_mode=run_mode)
 
         info_param = info.algorithm.get("param", {})
@@ -75,7 +97,6 @@ def __init__(self, info: odatse.Info,
             print(f"interval = {self.interval}")
             print(f"num_rand_basis = {self.num_rand_basis}")
 
-        #self.mesh_list, actions = self._meshgrid(info, split=False)
         if domain and isinstance(domain, odatse.domain.MeshGrid):
             self.domain = domain
         else:
@@ -90,12 +111,14 @@ def __init__(self, info: odatse.Info,
             seed = info.algorithm["seed"]
             self.policy.set_seed(seed)
 
-        # store state
         self.file_history = "history.npz"
         self.file_training = "training.npz"
         self.file_predictor = "predictor.dump"
 
     def _initialize(self):
+        """
+        Initializes the algorithm parameters and timers.
+        """
         self.istep = 0
         self.param_list = []
         self.fx_list = []
@@ -105,13 +128,27 @@ def _initialize(self):
         self._show_parameters()
 
     def _run(self) -> None:
+        """
+        Runs the Bayesian optimization process.
+        """
         runner = self.runner
         mesh_list = self.mesh_list
-        # fx_list = []
-        # param_list = []
 
         class simulator:
             def __call__(self, action: np.ndarray) -> float:
+                """
+                Simulates the function evaluation for a given action.
+
+                Parameters
+                ----------
+                action : np.ndarray
+                    The action to be evaluated.
+
+                Returns
+                -------
+                float
+                    The negative function value.
+                """
                 a = int(action[0])
                 args = (a, 0)
                 x = mesh_list[a, 1:]
@@ -145,9 +182,7 @@ def __call__(self, action: np.ndarray) -> float:
             time_end = time.perf_counter()
             self.timer["run"]["random_search"] = time_end - time_sta
 
-            # store initial state
             if self.checkpoint:
-                # print(">>> store initial state")
                 self._save_state(self.checkpoint_file)
         else:
             if self.istep >= self.bayes_max_num_probes:
@@ -157,7 +192,6 @@ def __call__(self, action: np.ndarray) -> float:
         next_checkpoint_time = time.time() + self.checkpoint_interval
 
         while self.istep < self.bayes_max_num_probes:
-            # print(">>> step {}".format(self.istep+1))
             intv = 0 if self.istep % self.interval == 0 else -1
 
             time_sta = time.perf_counter()
@@ -176,8 +210,6 @@ def __call__(self, action: np.ndarray) -> float:
             if self.checkpoint:
                 time_now = time.time()
                 if self.istep >= next_checkpoint_step or time_now >= next_checkpoint_time:
-                    # print(">>> checkpointing")
-
                     self.fx_list = fx_list
                     self.param_list = param_list
 
@@ -190,18 +222,19 @@ def __call__(self, action: np.ndarray) -> float:
         self.fx_list = fx_list
         self.param_list = param_list
 
-        # physbo.search.utility.show_search_results(self.policy.history, 20)
-
-        # store final state for continuation
         if self.checkpoint:
-            # print(">>> store final state")
             self._save_state(self.checkpoint_file)
 
     def _prepare(self) -> None:
-        # do nothing
+        """
+        Prepares the algorithm for execution.
+        """
         pass
 
     def _post(self) -> None:
+        """
+        Finalizes the algorithm execution and writes the results to a file.
+        """
         label_list = self.label_list
         if self.mpirank == 0:
             with open("BayesData.txt", "w") as file_BD:
@@ -229,14 +262,20 @@ def _post(self) -> None:
         return {"x": self.xopt, "fx": self.best_fx}
 
     def _save_state(self, filename):
+        """
+        Saves the current state of the algorithm to a file.
+
+        Parameters
+        ----------
+        filename : str
+            The name of the file to save the state.
+        """
         data = {
-            #-- _algorithm
             "mpisize": self.mpisize,
             "mpirank": self.mpirank,
             "rng": self.rng.get_state(),
             "timer": self.timer,
             "info": self.info,
-            #-- bayes
             "istep": self.istep,
             "param_list": self.param_list,
             "fx_list": self.fx_list,
@@ -247,19 +286,28 @@ def _save_state(self, filename):
         }
         self._save_data(data, filename)
 
-        #-- bayes
         self.policy.save(file_history=Path(self.output_dir, self.file_history),
                          file_training=Path(self.output_dir, self.file_training),
                          file_predictor=Path(self.output_dir, self.file_predictor))
 
-
     def _load_state(self, filename, mode="resume", restore_rng=True):
+        """
+        Loads the state of the algorithm from a file.
+
+        Parameters
+        ----------
+        filename : str
+            The name of the file to load the state from.
+        mode : str, optional
+            The mode to load the state (default is "resume").
+        restore_rng : bool, optional
+            Whether to restore the random number generator state (default is True).
+        """
         data = self._load_data(filename)
         if not data:
             print("ERROR: Load status file failed")
             sys.exit(1)
 
-        #-- _algorithm
         assert self.mpisize == data["mpisize"]
         assert self.mpirank == data["mpirank"]
 
@@ -272,11 +320,10 @@ def _load_state(self, filename, mode="resume", restore_rng=True):
         info = data["info"]
         self._check_parameters(info)
 
-        #-- bayes
         self.istep = data["istep"]
         self.param_list = data["param_list"]
         self.fx_list = data["fx_list"]
 
         self.policy.load(file_history=Path(self.output_dir, self.file_history),
                          file_training=Path(self.output_dir, self.file_training),
-                         file_predictor=Path(self.output_dir, self.file_predictor))
+                         file_predictor=Path(self.output_dir, self.file_predictor))
\ No newline at end of file
diff --git a/src/odatse/algorithm/exchange.py b/src/odatse/algorithm/exchange.py
index 6888bcac..c6831f66 100644
--- a/src/odatse/algorithm/exchange.py
+++ b/src/odatse/algorithm/exchange.py
@@ -81,6 +81,18 @@ def __init__(self,
                  runner: odatse.Runner = None,
                  run_mode: str = "initial"
     ) -> None:
+        """
+        Initialize the Algorithm class.
+
+        Parameters
+        ----------
+        info : odatse.Info
+            Information object containing algorithm parameters.
+        runner : odatse.Runner, optional
+            Runner object for executing the algorithm.
+        run_mode : str, optional
+            Mode to run the algorithm in, by default "initial".
+        """
         time_sta = time.perf_counter()
 
         info_exchange = info.algorithm["exchange"]
@@ -97,12 +109,18 @@ def __init__(self,
         self.timer["init"]["total"] = time_end - time_sta
 
     def _print_info(self) -> None:
+        """
+        Print information about the algorithm.
+        """
         if self.mpirank == 0:
             pass
         if self.mpisize > 1:
             self.mpicomm.barrier()
 
     def _initialize(self) -> None:
+        """
+        Initialize the algorithm parameters and state.
+        """
         super()._initialize()
 
         self.Tindex = np.arange(
@@ -117,6 +135,9 @@ def _initialize(self) -> None:
         self._show_parameters()
 
     def _run(self) -> None:
+        """
+        Run the algorithm.
+        """
         # print(">>> _run")
 
         if self.mode is None:
@@ -196,13 +217,28 @@ def _run(self) -> None:
             self._save_state(self.checkpoint_file)
 
     def _exchange(self, direction: bool) -> None:
-        """try to exchange temperatures"""
+        """
+        Try to exchange temperatures.
+
+        Parameters
+        ----------
+        direction : bool
+            Direction of the exchange.
+        """
         if self.nwalkers == 1:
             self.__exchange_single_walker(direction)
         else:
             self.__exchange_multi_walker(direction)
 
     def __exchange_single_walker(self, direction: bool) -> None:
+        """
+        Exchange temperatures for a single walker.
+
+        Parameters
+        ----------
+        direction : bool
+            Direction of the exchange.
+        """
         if self.mpisize > 1:
             self.mpicomm.Barrier()
         if direction:
@@ -256,6 +292,14 @@ def __exchange_single_walker(self, direction: bool) -> None:
         self.mpicomm.Bcast(self.T2rep, root=0)
 
     def __exchange_multi_walker(self, direction: bool) -> None:
+        """
+        Exchange temperatures for multiple walkers.
+
+        Parameters
+        ----------
+        direction : bool
+            Direction of the exchange.
+        """
         comm = self.mpicomm
         if self.mpisize > 1:
             fx_all = comm.allgather(self.fx)
@@ -300,10 +344,16 @@ def __exchange_multi_walker(self, direction: bool) -> None:
         ]
 
     def _prepare(self) -> None:
+        """
+        Prepare the algorithm for execution.
+        """
         self.timer["run"]["submit"] = 0.0
         self.timer["run"]["exchange"] = 0.0
 
     def _post(self) -> None:
+        """
+        Post-process the results of the algorithm.
+        """
         Ts = self.betas if self.input_as_beta else 1.0 / self.betas
         if self.mpirank == 0:
             print(f"start separateT {self.mpirank}")
@@ -360,6 +410,14 @@ def _post(self) -> None:
         }
 
     def _save_state(self, filename) -> None:
+        """
+        Save the current state of the algorithm.
+
+        Parameters
+        ----------
+        filename : str
+            The name of the file to save the state to.
+        """
         data = {
             #-- _algorithm
             "mpisize": self.mpisize,
@@ -388,6 +446,18 @@ def _save_state(self, filename) -> None:
         self._save_data(data, filename)
 
     def _load_state(self, filename, mode="resume", restore_rng=True):
+        """
+        Load the state of the algorithm from a file.
+
+        Parameters
+        ----------
+        filename : str
+            The name of the file to load the state from.
+        mode : str, optional
+            The mode to load the state in, by default "resume".
+        restore_rng : bool, optional
+            Whether to restore the random number generator state, by default True.
+        """
         data = self._load_data(filename)
         if not data:
             print("ERROR: Load status file failed")
@@ -426,5 +496,4 @@ def _load_state(self, filename, mode="resume", restore_rng=True):
         self.Tindex = data["Tindex"]
         self.rep2T = data["rep2T"]
         self.T2rep = data["T2rep"]
-        self.exchange_direction = data["exchange_direction"]
-
+        self.exchange_direction = data["exchange_direction"]
\ No newline at end of file
diff --git a/src/odatse/algorithm/mapper_mpi.py b/src/odatse/algorithm/mapper_mpi.py
index be6f3e68..bc0a0022 100644
--- a/src/odatse/algorithm/mapper_mpi.py
+++ b/src/odatse/algorithm/mapper_mpi.py
@@ -26,6 +26,10 @@
 import odatse.domain
 
 class Algorithm(odatse.algorithm.AlgorithmBase):
+    """
+    Algorithm class for data analysis of quantum beam diffraction experiments.
+    Inherits from odatse.algorithm.AlgorithmBase.
+    """
     mesh_list: List[Union[int, float]]
 
     def __init__(self, info: odatse.Info,
@@ -33,6 +37,14 @@ def __init__(self, info: odatse.Info,
                  domain = None,
                  run_mode: str = "initial"
     ) -> None:
+        """
+        Initialize the Algorithm instance.
+
+        :param info: Information object containing algorithm parameters.
+        :param runner: Optional runner object for submitting tasks.
+        :param domain: Optional domain object, defaults to MeshGrid.
+        :param run_mode: Mode to run the algorithm, defaults to "initial".
+        """
         super().__init__(info=info, runner=runner, run_mode=run_mode)
 
         if domain and isinstance(domain, odatse.domain.MeshGrid):
@@ -47,21 +59,25 @@ def __init__(self, info: odatse.Info,
         self.local_colormap_file = Path(self.colormap_file).name + ".tmp"
 
     def _initialize(self) -> None:
+        """
+        Initialize the algorithm parameters and timer.
+        """
         self.fx_list = []
         self.timer["run"]["submit"] = 0.0
         self._show_parameters()
 
     def _run(self) -> None:
+        """
+        Execute the main algorithm process.
+        """
         # Make ColorMap
 
         if self.mode is None:
             raise RuntimeError("mode unset")
 
         if self.mode.startswith("init"):
-            # print(">>> initialize")
             self._initialize()
         elif self.mode.startswith("resume"):
-            # print(">>> resume")
             self._load_state(self.checkpoint_file)
         else:
             raise RuntimeError("unknown mode {}".format(self.mode))
@@ -74,8 +90,6 @@ def _run(self) -> None:
         iterations = len(self.mesh_list)
         istart = len(self.fx_list)
 
-        # print(">>> iterations={}, istart={}".format(iterations, istart))
-
         next_checkpoint_step = istart + self.checkpoint_steps
         next_checkpoint_time = time.time() + self.checkpoint_interval
 
@@ -111,8 +125,6 @@ def _run(self) -> None:
             opt_id, opt_fx = self.fx_list[opt_index]
             opt_mesh = self.mesh_list[opt_index]
 
-            # assert opt_id == opt_mesh[0]
-
             self.opt_fx = opt_fx
             self.opt_mesh = opt_mesh
 
@@ -121,12 +133,14 @@ def _run(self) -> None:
         self._output_results()
 
         if Path(self.local_colormap_file).exists():
-            # print(">>> remove local colormap file {}".format(self.local_colormap_file))
             os.remove(Path(self.local_colormap_file))
 
         print("complete main process : rank {:08d}/{:08d}".format(self.mpirank, self.mpisize))
 
     def _output_results(self):
+        """
+        Output the results to the colormap file.
+        """
         print("Make ColorMap")
         time_sta = time.perf_counter()
 
@@ -151,10 +165,17 @@ def _output_results(self):
         self.timer["run"]["file_CM"] = time_end - time_sta
 
     def _prepare(self) -> None:
-        # do nothing
+        """
+        Prepare the algorithm (no operation).
+        """
         pass
 
     def _post(self) -> Dict:
+        """
+        Post-process the results and gather data from all MPI ranks.
+
+        :return: Dictionary of results.
+        """
         if self.mpisize > 1:
             fx_lists = self.mpicomm.allgather(self.fx_list)
             results = [v for vs in fx_lists for v in vs]
@@ -163,7 +184,6 @@ def _post(self) -> Dict:
 
         if self.mpirank == 0:
             with open(self.colormap_file, "w") as fp:
-                # fp.write("#" + " ".join(self.label_list) + " fval\n")
                 for x, (idx, fx) in zip(self.domain.grid, results):
                     assert x[0] == idx
                     fp.write(" ".join(
@@ -173,38 +193,41 @@ def _post(self) -> Dict:
         return {}
 
     def _save_state(self, filename) -> None:
+        """
+        Save the current state of the algorithm to a file.
+
+        :param filename: The name of the file to save the state to.
+        """
         data = {
-            #-- _algorithm
             "mpisize": self.mpisize,
             "mpirank": self.mpirank,
-            #"rng": self.rng.get_state(),
             "timer": self.timer,
             "info": self.info,
-            #-- mapper
             "fx_list": self.fx_list,
             "mesh_size": len(self.mesh_list),
         }
         self._save_data(data, filename)
 
     def _load_state(self, filename, restore_rng=True):
+        """
+        Load the state of the algorithm from a file.
+
+        :param filename: The name of the file to load the state from.
+        :param restore_rng: Whether to restore the random number generator state.
+        """
         data = self._load_data(filename)
         if not data:
             print("ERROR: Load status file failed")
             sys.exit(1)
 
-        #-- _algorithm
         assert self.mpisize == data["mpisize"]
         assert self.mpirank == data["mpirank"]
 
-        # if restore_rng:
-        #     self.rng = np.random.RandomState()
-        #     self.rng.set_state(data["rng"])
         self.timer = data["timer"]
 
         info = data["info"]
         self._check_parameters(info)
 
-        #-- mapper
         self.fx_list = data["fx_list"]
 
-        assert len(self.mesh_list) == data["mesh_size"]
+        assert len(self.mesh_list) == data["mesh_size"]
\ No newline at end of file
diff --git a/src/odatse/algorithm/min_search.py b/src/odatse/algorithm/min_search.py
index 91fd18f0..62bbb399 100644
--- a/src/odatse/algorithm/min_search.py
+++ b/src/odatse/algorithm/min_search.py
@@ -26,6 +26,9 @@
 
 
 class Algorithm(odatse.algorithm.AlgorithmBase):
+    """
+    Algorithm class for performing minimization using the Nelder-Mead method.
+    """
 
     # inputs
     label_list: np.ndarray
@@ -54,6 +57,14 @@ def __init__(self, info: odatse.Info,
                  domain = None,
                  run_mode: str = "initial"
     ) -> None:
+        """
+        Initialize the Algorithm class.
+
+        :param info: Information object containing algorithm settings.
+        :param runner: Runner object for submitting jobs.
+        :param domain: Domain object defining the search space.
+        :param run_mode: Mode of running the algorithm.
+        """
         super().__init__(info=info, runner=runner, run_mode=run_mode)
 
         if domain and isinstance(domain, odatse.domain.Region):
@@ -80,6 +91,9 @@ def __init__(self, info: odatse.Info,
         self._show_parameters()
 
     def _run(self) -> None:
+        """
+        Run the minimization algorithm.
+        """
         run = self.runner
 
         min_list = self.min_list
@@ -98,17 +112,30 @@ def _run(self) -> None:
 
         if scipy_version[0] >= 1 and scipy_version[1] >= 11:
             def _cb(intermediate_result):
+                """
+                Callback function for scipy.optimize.minimize.
+                """
                 x = intermediate_result.x
                 fun = intermediate_result.fun
                 print("eval: x={}, fun={}".format(x, fun))
                 iter_history.append([*x, fun])
         else:
             def _cb(x):
+                """
+                Callback function for scipy.optimize.minimize.
+                """
                 fun = _f_calc(x, 1)
                 print("eval: x={}, fun={}".format(x, fun))
                 iter_history.append([*x, fun])
 
         def _f_calc(x_list: np.ndarray, iset) -> float:
+            """
+            Calculate the objective function value.
+
+            :param x_list: List of variables.
+            :param iset: Set index.
+            :return: Objective function value.
+            """
             # check if within region -> boundary option in minimize
             # note: 'bounds' option supported in scipy >= 1.7.0
             in_range = np.all((min_list < x_list) & (x_list < max_list))
@@ -167,6 +194,9 @@ def _f_calc(x_list: np.ndarray, iset) -> float:
             self.mpicomm.barrier()
 
     def _prepare(self):
+        """
+        Prepare the initial simplex for the Nelder-Mead algorithm.
+        """
         # make initial simplex
         #   [ v0, v0+a_1*e_1, v0+a_2*e_2, ... v0+a_d*e_d ]
         # where a = ( a_1 a_2 a_3 ... a_d ) and e_k is a unit vector along k-axis
@@ -175,6 +205,9 @@ def _prepare(self):
         self.initial_simplex_list = np.vstack((v, v + np.diag(a)))
 
     def _output_results(self):
+        """
+        Output the results of the minimization to files.
+        """
         label_list = self.label_list
 
         with open("SimplexData.txt", "w") as fp:
@@ -195,6 +228,9 @@ def _output_results(self):
             fp.write(f"function_evaluations = {self.funcalls}\n")
 
     def _post(self):
+        """
+        Post-process the results after minimization.
+        """
         result = {
             "x": self.xopt,
             "fx": self.fopt,
@@ -223,4 +259,4 @@ def _post(self):
                     for x, y in zip(label_list, x0s[idx]):
                         fp.write(f"initial {x} = {y}\n")
 
-        return {"x": xs[idx], "fx": fxs[idx], "x0": x0s[idx]}
+        return {"x": xs[idx], "fx": fxs[idx], "x0": x0s[idx]}
\ No newline at end of file
diff --git a/src/odatse/algorithm/montecarlo.py b/src/odatse/algorithm/montecarlo.py
index 451d3dba..b3db2846 100644
--- a/src/odatse/algorithm/montecarlo.py
+++ b/src/odatse/algorithm/montecarlo.py
@@ -93,11 +93,26 @@ class AlgorithmBase(odatse.algorithm.AlgorithmBase):
     naccepted: int
 
     def __init__(self, info: odatse.Info,
-                 runner: odatse.Runner = None,
-                 domain = None,
-                 nwalkers: int = 1,
-                 run_mode: str = "initial"
-    ) -> None:
+             runner: odatse.Runner = None,
+             domain = None,
+             nwalkers: int = 1,
+             run_mode: str = "initial") -> None:
+        """
+        Initialize the AlgorithmBase class.
+
+        Parameters
+        ----------
+        info : odatse.Info
+            Information object containing algorithm parameters.
+        runner : odatse.Runner, optional
+            Runner object for executing the algorithm (default is None).
+        domain : optional
+            Domain object defining the problem space (default is None).
+        nwalkers : int, optional
+            Number of walkers (default is 1).
+        run_mode : str, optional
+            Mode of the run, e.g., "initial" (default is "initial").
+        """
         time_sta = time.perf_counter()
         super().__init__(info=info, runner=runner, run_mode=run_mode)
         self.nwalkers = nwalkers
@@ -146,6 +161,13 @@ def __init__(self, info: odatse.Info,
         self.input_as_beta = False
 
     def _initialize(self):
+        """
+            Initialize the algorithm state.
+
+            This method sets up the initial state of the algorithm, including the
+            positions and energies of the walkers, and resets the counters for
+            accepted and trial steps.
+        """
         if self.iscontinuous:
             self.domain.initialize(rng=self.rng, limitation=self.runner.limitation, num_walkers=self.nwalkers)
             self.x = self.domain.initial_list
@@ -162,6 +184,21 @@ def _initialize(self):
         self.ntrial = 0
 
     def _setup_neighbour(self, info_param):
+        """
+        Set up the neighbor list for the discrete problem.
+
+        Parameters
+        ----------
+        info_param : dict
+            Dictionary containing algorithm parameters, including the path to the neighbor list file.
+
+        Raises
+        ------
+        ValueError
+            If the neighbor list path is not specified in the parameters.
+        RuntimeError
+            If the transition graph made from the neighbor list is not connected or not bidirectional.
+        """
         if "neighborlist_path" in info_param:
             nn_path = self.root_dir / Path(info_param["neighborlist_path"]).expanduser()
             self.neighbor_list = load_neighbor_list(nn_path, nnodes=self.nnodes)
@@ -184,17 +221,23 @@ def _setup_neighbour(self, info_param):
             )
             # otherwise find neighbourlist
 
+
     def _evaluate(self, in_range: np.ndarray = None) -> np.ndarray:
-        """evaluate current "Energy"s
+        """
+            Evaluate the current "Energy"s.
 
-        ``self.fx`` will be overwritten with the result
+            This method overwrites `self.fx` with the result.
 
-        Parameters
-        ==========
-        run_info: dict
-            Parameter set.
-            Some parameters will be overwritten.
-        """
+            Parameters
+            ----------
+            in_range : np.ndarray, optional
+                Array indicating whether each walker is within the valid range (default is None).
+
+            Returns
+            -------
+            np.ndarray
+                Array of evaluated energies for the current configurations.
+            """
         # print(">>> _evaluate")
         for iwalker in range(self.nwalkers):
             x = self.x[iwalker, :]
@@ -210,18 +253,19 @@ def _evaluate(self, in_range: np.ndarray = None) -> np.ndarray:
         return self.fx
 
     def propose(self, current: np.ndarray) -> np.ndarray:
-        """propose next candidate
-
-        Parameters
-        ==========
-        current: np.ndarray
-            current position
-
-        Returns
-        =======
-        proposed: np.ndarray
-            proposal
         """
+            Propose the next candidate positions for the walkers.
+
+            Parameters
+            ----------
+            current : np.ndarray
+                Current positions of the walkers.
+
+            Returns
+            -------
+            proposed : np.ndarray
+                Proposed new positions for the walkers.
+            """
         if self.iscontinuous:
             dx = self.rng.normal(size=(self.nwalkers, self.dimension)) * self.xstep
             proposed = current + dx
@@ -313,6 +357,16 @@ def local_update(
         self._write_result(file_result, extra_info_to_write=extra_info_to_write)
 
     def _write_result_header(self, fp, extra_names=None) -> None:
+        """
+        Write the header for the result file.
+
+        Parameters
+        ----------
+        fp : TextIO
+            File pointer to the result file.
+        extra_names : list of str, optional
+            Additional column names to include in the header.
+        """
         if self.input_as_beta:
             fp.write("# step walker beta fx")
         else:
@@ -325,35 +379,60 @@ def _write_result_header(self, fp, extra_names=None) -> None:
         fp.write("\n")
 
     def _write_result(self, fp, extra_info_to_write: Union[List, Tuple] = None) -> None:
-        for iwalker in range(self.nwalkers):
-            if isinstance(self.Tindex, int):
-                beta = self.betas[self.Tindex]
-            else:
-                beta = self.betas[self.Tindex[iwalker]]
-            fp.write(f"{self.istep}")
-            fp.write(f" {iwalker}")
-            if self.input_as_beta:
-                fp.write(f" {beta}")
-            else:
-                fp.write(f" {1.0/beta}")
-            fp.write(f" {self.fx[iwalker]}")
-            for x in self.x[iwalker, :]:
-                fp.write(f" {x}")
-            if extra_info_to_write is not None:
-                for ex in extra_info_to_write:
-                    fp.write(f" {ex[iwalker]}")
-            fp.write("\n")
-        fp.flush()
+        """
+        Write the result of the current step to the file.
 
+        Parameters
+        ----------
+        fp : TextIO
+            File pointer to the result file.
+        extra_info_to_write : Union[List, Tuple], optional
+            Additional information to write for each walker (default is None).
+        """
+    for iwalker in range(self.nwalkers):
+        if isinstance(self.Tindex, int):
+            beta = self.betas[self.Tindex]
+        else:
+            beta = self.betas[self.Tindex[iwalker]]
+        fp.write(f"{self.istep}")
+        fp.write(f" {iwalker}")
+        if self.input_as_beta:
+            fp.write(f" {beta}")
+        else:
+            fp.write(f" {1.0/beta}")
+        fp.write(f" {self.fx[iwalker]}")
+        for x in self.x[iwalker, :]:
+            fp.write(f" {x}")
+        if extra_info_to_write is not None:
+            for ex in extra_info_to_write:
+                fp.write(f" {ex[iwalker]}")
+        fp.write("\n")
+    fp.flush()
 def read_Ts(info: dict, numT: int = None) -> Tuple[bool, np.ndarray]:
     """
+    Read temperature or inverse-temperature values from the provided info dictionary.
+
+    Parameters
+    ----------
+    info : dict
+        Dictionary containing temperature or inverse-temperature parameters.
+    numT : int, optional
+        Number of temperature or inverse-temperature values to generate (default is None).
 
     Returns
     -------
-    as_beta: bool
-        true when using inverse-temperature
-    betas: np.ndarray
-        sequence of inverse-temperature
+    as_beta : bool
+        True if using inverse-temperature, False if using temperature.
+    betas : np.ndarray
+        Sequence of inverse-temperature values.
+
+    Raises
+    ------
+    ValueError
+        If numT is not specified, or if both Tmin/Tmax and bmin/bmax are defined, or if neither are defined,
+        or if bmin/bmax or Tmin/Tmax values are invalid.
+    RuntimeError
+        If the mode is unknown (neither set_T nor set_b).
     """
     if numT is None:
         raise ValueError("read_Ts: numT is not specified")
diff --git a/src/odatse/algorithm/pamc.py b/src/odatse/algorithm/pamc.py
index e6aaacbb..c71477db 100644
--- a/src/odatse/algorithm/pamc.py
+++ b/src/odatse/algorithm/pamc.py
@@ -100,6 +100,18 @@ def __init__(self,
                  runner: odatse.Runner = None,
                  run_mode: str = "initial"
     ) -> None:
+        """
+        Initialize the Algorithm class.
+
+        Parameters
+        ----------
+        info : odatse.Info
+            Information object containing algorithm parameters.
+        runner : odatse.Runner, optional
+            Runner object for executing the algorithm, by default None.
+        run_mode : str, optional
+            Mode in which to run the algorithm, by default "initial".
+        """
         time_sta = time.perf_counter()
 
         info_pamc = info.algorithm["pamc"]
@@ -147,6 +159,19 @@ def _initialize(self) -> None:
         self._show_parameters()
 
     def _find_scheduling(self, info_pamc) -> int:
+        """
+        Determine the scheduling for the algorithm based on the provided parameters.
+
+        Parameters
+        ----------
+        info_pamc : dict
+            Dictionary containing the parameters for the PAMC algorithm.
+
+        Returns
+        -------
+        int
+            The number of temperature steps (numT) determined from the input parameters.
+        """
         numsteps = info_pamc.get("numsteps", 0)
         numsteps_annealing = info_pamc.get("numsteps_annealing", 0)
         numT = info_pamc.get("Tnum", 0)
@@ -381,6 +406,19 @@ def _gather_information(self, numT: int = None) -> Dict[str, np.ndarray]:
         return res
 
     def _save_stats(self, info: Dict[str, np.ndarray]) -> None:
+        """
+        Save statistical information from the algorithm run.
+
+        Parameters
+        ----------
+        info : Dict[str, np.ndarray]
+            Dictionary containing the following keys:
+            - fxs: Objective function of each walker over all processes.
+            - logweights: Logarithm of weights.
+            - ns: Number of walkers in each process.
+            - ancestors: Ancestor (origin) of each walker.
+            - acceptance ratio: Acceptance ratio for each temperature.
+        """
         fxs = info["fxs"]
         numT, nreplicas = fxs.shape
         endTindex = self.Tindex + 1
@@ -391,7 +429,7 @@ def _save_stats(self, info: Dict[str, np.ndarray]) -> None:
             logweights - logweights.max(axis=1).reshape(-1, 1)
         )  # to avoid overflow
 
-        # bias-corrected jackknife resampling method
+        # Bias-corrected jackknife resampling method
         fs = np.zeros((numT, nreplicas))
         fw_sum = (fxs * weights).sum(axis=1)
         w_sum = weights.sum(axis=1)
@@ -412,7 +450,7 @@ def _save_stats(self, info: Dict[str, np.ndarray]) -> None:
         logz = np.log(np.mean(weights, axis=1))
         self.logZs[startTindex:endTindex] = self.logZ + logz
         if endTindex < len(self.betas):
-            # calculate the next weight before reset and evalute dF
+            # Calculate the next weight before reset and evaluate dF
             bdiff = self.betas[endTindex] - self.betas[endTindex - 1]
             w = np.exp(logweights[-1, :] - bdiff * fxs[-1, :])
             self.logZ = self.logZs[startTindex] + np.log(w.mean())
@@ -429,6 +467,13 @@ def _save_stats(self, info: Dict[str, np.ndarray]) -> None:
                 ])))
 
     def _resample(self) -> None:
+        """
+            Perform the resampling of walkers.
+
+            This method gathers information, saves statistical data, and performs resampling
+            using either fixed or varied weights. The method ensures that the algorithm
+            maintains a balanced set of walkers across different temperature steps.
+        """
         res = self._gather_information()
         self._save_stats(res)
 
@@ -448,6 +493,19 @@ def _resample(self) -> None:
             self.logweights = np.zeros(self.nwalkers)
 
     def _resample_varied(self, weights: np.ndarray, offset: int) -> None:
+        """
+        Perform resampling with varied weights.
+
+        This method resamples the walkers based on the provided weights and updates
+        the state of the algorithm accordingly.
+
+        Parameters
+        ----------
+        weights : np.ndarray
+            Array of weights for resampling.
+        offset : int
+            Offset for the weights array.
+        """
         weights_sum = np.sum(weights)
         expected_numbers = (self.nreplicas[0] / weights_sum) * weights[
             offset : offset + self.nwalkers
@@ -482,6 +540,17 @@ def _resample_varied(self, weights: np.ndarray, offset: int) -> None:
         self.nwalkers = np.sum(next_numbers)
 
     def _resample_fixed(self, weights: np.ndarray) -> None:
+        """
+        Perform resampling with fixed weights.
+
+        This method resamples the walkers based on the provided weights and updates
+        the state of the algorithm accordingly.
+
+        Parameters
+        ----------
+        weights : np.ndarray
+            Array of weights for resampling.
+        """
         resampler = odatse.util.resampling.WalkerTable(weights)
         new_index = resampler.sample(self.rng, self.nwalkers)
 
@@ -512,10 +581,22 @@ def _resample_fixed(self, weights: np.ndarray) -> None:
             self.x = self.node_coordinates[self.inode, :]
 
     def _prepare(self) -> None:
-        self.timer["run"]["submit"] = 0.0
-        self.timer["run"]["resampling"] = 0.0
+        """
+        Prepare the algorithm for execution.
 
+        This method initializes the timers for the 'submit' and 'resampling' phases
+        of the algorithm run.
+        """
+    self.timer["run"]["submit"] = 0.0
+    self.timer["run"]["resampling"] = 0.0
     def _post(self) -> None:
+        """
+        Post-processing after the algorithm execution.
+
+        This method consolidates the results from different temperature steps
+        into single files for 'result' and 'trial'. It also gathers the best
+        results from all processes and writes them to 'best_result.txt'.
+        """
         for name in ("result", "trial"):
             with open(self.proc_dir / f"{name}.txt", "w") as fout:
                 self._write_result_header(fout, ["weight", "ancestor"])
@@ -581,6 +662,14 @@ def _post(self) -> None:
         }
 
     def _save_state(self, filename) -> None:
+        """
+        Save the current state of the algorithm to a file.
+
+        Parameters
+        ----------
+        filename : str
+            The name of the file where the state will be saved.
+        """
         data = {
             #-- _algorithm
             "mpisize": self.mpisize,
@@ -622,6 +711,18 @@ def _save_state(self, filename) -> None:
         self._save_data(data, filename)
 
     def _load_state(self, filename, mode="resume", restore_rng=True):
+        """
+        Load the saved state of the algorithm from a file.
+
+        Parameters
+        ----------
+        filename : str
+            The name of the file from which the state will be loaded.
+        mode : str, optional
+            The mode in which to load the state. Can be "resume" or "continue", by default "resume".
+        restore_rng : bool, optional
+            Whether to restore the random number generator state, by default True.
+        """
         data = self._load_data(filename)
         if not data:
             print("ERROR: Load status file failed")
diff --git a/src/odatse/domain/_domain.py b/src/odatse/domain/_domain.py
index 7980097e..347b26c5 100644
--- a/src/odatse/domain/_domain.py
+++ b/src/odatse/domain/_domain.py
@@ -22,7 +22,22 @@
 import odatse
 
 class DomainBase:
+    """
+    Base class for domain management in the 2DMAT software.
+
+    Attributes:
+        root_dir (Path): The root directory for the domain.
+        output_dir (Path): The output directory for the domain.
+        mpisize (int): The size of the MPI communicator.
+        mpirank (int): The rank of the MPI process.
+    """
     def __init__(self, info: odatse.Info = None):
+        """
+        Initializes the DomainBase instance.
+
+        Args:
+            info (odatse.Info, optional): An instance of odatse.Info containing base directory information.
+        """
         if info:
             self.root_dir = info.base["root_dir"]
             self.output_dir = info.base["output_dir"]
@@ -32,4 +47,3 @@ def __init__(self, info: odatse.Info = None):
 
         self.mpisize = odatse.mpi.size()
         self.mpirank = odatse.mpi.rank()
-            
diff --git a/src/odatse/domain/meshgrid.py b/src/odatse/domain/meshgrid.py
index 123f4e9c..a438093f 100644
--- a/src/odatse/domain/meshgrid.py
+++ b/src/odatse/domain/meshgrid.py
@@ -23,13 +23,21 @@
 from ._domain import DomainBase
 
 class MeshGrid(DomainBase):
+    """
+    MeshGrid class for handling grid data for quantum beam diffraction experiments.
+    """
+
     grid: List[Union[int, float]] = []
     grid_local: List[Union[int, float]] = []
     candicates: int
-    
-    def __init__(self, info: odatse.Info = None,
-                 *,
-                 param: Dict[str, Any] = None):
+
+    def __init__(self, info: odatse.Info = None, *, param: Dict[str, Any] = None):
+        """
+        Initialize the MeshGrid object.
+
+        :param info: Information object containing algorithm parameters.
+        :param param: Dictionary containing parameters for setting up the grid.
+        """
         super().__init__(info)
 
         if info:
@@ -42,26 +50,36 @@ def __init__(self, info: odatse.Info = None,
         else:
             pass
 
-
     def do_split(self):
+        """
+        Split the grid data among MPI processes.
+        """
         if self.mpisize > 1:
             index = [idx for idx, *v in self.grid]
             index_local = np.array_split(index, self.mpisize)[self.mpirank]
             self.grid_local = [[idx, *v] for idx, *v in self.grid if idx in index_local]
         else:
             self.grid_local = self.grid
-            
 
     def _setup(self, info_param):
+        """
+        Setup the grid based on provided parameters.
+
+        :param info_param: Dictionary containing parameters for setting up the grid.
+        """
         if "mesh_path" in info_param:
             self._setup_from_file(info_param)
         else:
             self._setup_grid(info_param)
 
         self.ncandicates = len(self.grid)
-            
 
     def _setup_from_file(self, info_param):
+        """
+        Setup the grid from a file.
+
+        :param info_param: Dictionary containing parameters for setting up the grid.
+        """
         if "mesh_path" not in info_param:
             raise ValueError("ERROR: mesh_path not defined")
         mesh_path = self.root_dir / Path(info_param["mesh_path"]).expanduser()
@@ -79,7 +97,7 @@ def _setup_from_file(self, info_param):
                 data = data.reshape(1, -1)
 
             # old format: index x1 x2 ... -> omit index
-            data = data[:,1:]
+            data = data[:, 1:]
         else:
             data = None
 
@@ -88,8 +106,12 @@ def _setup_from_file(self, info_param):
 
         self.grid = [[idx, *v] for idx, v in enumerate(data)]
 
-
     def _setup_grid(self, info_param):
+        """
+        Setup the grid based on min, max, and num lists.
+
+        :param info_param: Dictionary containing parameters for setting up the grid.
+        """
         if "min_list" not in info_param:
             raise ValueError("ERROR: algorithm.param.min_list is not defined in the input")
         min_list = np.array(info_param["min_list"], dtype=float)
@@ -104,7 +126,7 @@ def _setup_grid(self, info_param):
 
         if len(min_list) != len(max_list) or len(min_list) != len(num_list):
             raise ValueError("ERROR: lengths of min_list, max_list, num_list do not match")
-        
+
         xs = [
             np.linspace(mn, mx, num=nm)
             for mn, mx, nm in zip(min_list, max_list, num_list)
@@ -118,21 +140,35 @@ def _setup_grid(self, info_param):
             )
         ]
 
-
     def store_file(self, store_path, *, header=""):
+        """
+        Store the grid data to a file.
+
+        :param store_path: Path to the file where the grid data will be stored.
+        :param header: Header to be included in the file.
+        """
         if self.mpirank == 0:
             np.savetxt(store_path, [[*v] for idx, *v in self.grid], header=header)
 
-
     @classmethod
     def from_file(cls, mesh_path):
-        return cls(param={"mesh_path": mesh_path})
+        """
+        Create a MeshGrid object from a file.
 
+        :param mesh_path: Path to the file containing the grid data.
+        :return: MeshGrid object.
+        """
+        return cls(param={"mesh_path": mesh_path})
 
     @classmethod
     def from_dict(cls, param):
-        return cls(param=param)
+        """
+        Create a MeshGrid object from a dictionary of parameters.
 
+        :param param: Dictionary containing parameters for setting up the grid.
+        :return: MeshGrid object.
+        """
+        return cls(param=param)
 
 
 if __name__ == "__main__":
diff --git a/src/odatse/domain/region.py b/src/odatse/domain/region.py
index 547f01c6..f4bb362a 100644
--- a/src/odatse/domain/region.py
+++ b/src/odatse/domain/region.py
@@ -23,14 +23,37 @@
 from ._domain import DomainBase
 
 class Region(DomainBase):
+    """
+    A class to represent a region in the domain.
+
+    Attributes
+    ----------
+    min_list : np.array
+        Minimum values for each dimension.
+    max_list : np.array
+        Maximum values for each dimension.
+    unit_list : np.array
+        Unit values for each dimension.
+    initial_list : np.array
+        Initial values for each dimension.
+    """
+
     min_list: np.array
     max_list: np.array
     unit_list: np.array
     initial_list: np.array
 
-    def __init__(self, info: odatse.Info = None,
-                 *,
-                 param: Dict[str, Any] = None):
+    def __init__(self, info: odatse.Info = None, *, param: Dict[str, Any] = None):
+        """
+        Initialize the Region object.
+
+        Parameters
+        ----------
+        info : odatse.Info, optional
+            Information object containing algorithm parameters.
+        param : dict, optional
+            Dictionary containing algorithm parameters.
+        """
         super().__init__(info)
 
         if info:
@@ -42,9 +65,16 @@ def __init__(self, info: odatse.Info = None,
             self._setup(param)
         else:
             pass
-    
 
     def _setup(self, info_param):
+        """
+        Setup the region with the given parameters.
+
+        Parameters
+        ----------
+        info_param : dict
+            Dictionary containing the parameters for the region.
+        """
         if "min_list" not in info_param:
             raise ValueError("ERROR: algorithm.param.min_list is not defined in the input")
         min_list = np.array(info_param["min_list"])
@@ -57,7 +87,7 @@ def _setup(self, info_param):
             raise ValueError("ERROR: lengths of min_list and max_list do not match")
 
         self.dimension = len(min_list)
-        
+
         unit_list = np.array(info_param.get("unit_list", [1.0] * self.dimension))
 
         self.min_list = min_list
@@ -77,10 +107,19 @@ def _setup(self, info_param):
 
         self.initial_list = initial_list
 
-    def initialize(self,
-                   rng=np.random,
-                   limitation=odatse.util.limitation.Unlimited(),
-                   num_walkers: int = 1):
+    def initialize(self, rng=np.random, limitation=odatse.util.limitation.Unlimited(), num_walkers: int = 1):
+        """
+        Initialize the region with random values or predefined initial values.
+
+        Parameters
+        ----------
+        rng : numpy.random, optional
+            Random number generator.
+        limitation : odatse.util.limitation, optional
+            Limitation object to judge the validity of the values.
+        num_walkers : int, optional
+            Number of walkers to initialize.
+        """
         if num_walkers > self.num_walkers:
             self.num_walkers = num_walkers
 
@@ -89,10 +128,19 @@ def initialize(self,
         else:
             self._init_random(rng=rng, limitation=limitation)
 
-    def _init_random(self,
-                     rng=np.random,
-                     limitation=odatse.util.limitation.Unlimited(),
-                     max_count=100):
+    def _init_random(self, rng=np.random, limitation=odatse.util.limitation.Unlimited(), max_count=100):
+        """
+        Initialize the region with random values within the specified limits.
+
+        Parameters
+        ----------
+        rng : numpy.random, optional
+            Random number generator.
+        limitation : odatse.util.limitation, optional
+            Limitation object to judge the validity of the values.
+        max_count : int, optional
+            Maximum number of trials to generate valid values.
+        """
         initial_list = np.zeros((self.num_walkers, self.dimension), dtype=float)
         is_ok = np.full(self.num_walkers, False)
 
@@ -110,7 +158,6 @@ def _init_random(self,
                 raise RuntimeError("ERROR: init_random: trial count exceeds {}".format(max_count))
         self.initial_list = initial_list
 
-
 if __name__ == "__main__":
     reg = Region(param={
         "min_list": [0.0, 0.0, 0.0],
diff --git a/src/odatse/solver/_solver.py b/src/odatse/solver/_solver.py
index d4bd7259..efaf628a 100644
--- a/src/odatse/solver/_solver.py
+++ b/src/odatse/solver/_solver.py
@@ -30,6 +30,10 @@
 
 
 class SolverBase(object, metaclass=ABCMeta):
+    """
+    Abstract base class for solvers in the 2DMAT software.
+    """
+
     root_dir: Path
     output_dir: Path
     proc_dir: Path
@@ -40,6 +44,12 @@ class SolverBase(object, metaclass=ABCMeta):
 
     @abstractmethod
     def __init__(self, info: odatse.Info) -> None:
+        """
+        Initialize the solver with the given information.
+
+        Args:
+            info (odatse.Info): Information object containing configuration details.
+        """
         self.root_dir = info.base["root_dir"]
         self.output_dir = info.base["output_dir"]
         self.proc_dir = self.output_dir / str(odatse.mpi.rank())
@@ -53,8 +63,26 @@ def __init__(self, info: odatse.Info) -> None:
 
     @property
     def name(self) -> str:
+        """
+        Get the name of the solver.
+
+        Returns:
+            str: The name of the solver.
+        """
         return self._name
 
     @abstractmethod
     def evaluate(self, x: np.ndarray, arg: Tuple = (), nprocs: int = 1, nthreads: int = 1) -> None:
-        raise NotImplementedError()
+        """
+        Evaluate the solver with the given parameters.
+
+        Args:
+            x (np.ndarray): Input data array.
+            arg (Tuple, optional): Additional arguments for evaluation. Defaults to ().
+            nprocs (int, optional): Number of processes to use. Defaults to 1.
+            nthreads (int, optional): Number of threads to use. Defaults to 1.
+
+        Raises:
+            NotImplementedError: This method should be implemented by subclasses.
+        """
+        raise NotImplementedError()
\ No newline at end of file
diff --git a/src/odatse/solver/analytical.py b/src/odatse/solver/analytical.py
index 74568c37..8b52af3a 100644
--- a/src/odatse/solver/analytical.py
+++ b/src/odatse/solver/analytical.py
@@ -19,28 +19,64 @@
 import odatse
 import odatse.solver.function
 
-
 def quadratics(xs: np.ndarray) -> float:
-    """quadratic (sphear) function
-
-    It has one global miminum f(xs)=0 at xs = [0,0,...,0].
     """
-    return np.sum(xs * xs)
+    Quadratic (sphere) function.
 
+    Parameters
+    ----------
+    xs : np.ndarray
+        Input array.
 
-def quartics(xs: np.ndarray) -> float:
-    """quartic function with two minimum
+    Returns
+    -------
+    float
+        The calculated value of the quadratic function.
 
-    It has two global minimum f(xs)=0 at xs = [1,1,...,1] and [0,0,...,0].
-    It has one suddle point f(0,0,...,0) = 1.0.
+    Notes
+    -----
+    It has one global minimum f(xs)=0 at xs = [0,0,...,0].
     """
+    return np.sum(xs * xs)
 
+def quartics(xs: np.ndarray) -> float:
+    """
+    Quartic function with two global minima.
+
+    Parameters
+    ----------
+    xs : np.ndarray
+        Input array.
+
+    Returns
+    -------
+    float
+        The calculated value of the quartic function.
+
+    Notes
+    -----
+    It has two global minima f(xs)=0 at xs = [1,1,...,1] and [0,0,...,0].
+    It has one saddle point f(0,0,...,0) = 1.0.
+    """
     return np.mean((xs - 1.0) ** 2) * np.mean((xs + 1.0) ** 2)
 
 
 def ackley(xs: np.ndarray) -> float:
-    """Ackley's function in arbitrary dimension
+    """
+    Ackley's function in arbitrary dimension
 
+    Parameters
+    ----------
+    xs : np.ndarray
+        Input array.
+
+    Returns
+    -------
+    float
+        The calculated value of Ackley's function.
+
+    Notes
+    -----
     It has one global minimum f(xs)=0 at xs=[0,0,...,0].
     It has many local minima.
     """
@@ -50,18 +86,42 @@ def ackley(xs: np.ndarray) -> float:
     b = np.exp(0.5 * np.sum(b))
     return 20.0 + np.exp(1.0) - a - b
 
-
 def rosenbrock(xs: np.ndarray) -> float:
-    """Rosenbrock's function
+    """
+    Rosenbrock's function.
+
+    Parameters
+    ----------
+    xs : np.ndarray
+        Input array.
+
+    Returns
+    -------
+    float
+        The calculated value of Rosenbrock's function.
 
+    Notes
+    -----
     It has one global minimum f(xs) = 0 at xs=[1,1,...,1].
     """
     return np.sum(100.0 * (xs[1:] - xs[:-1] ** 2) ** 2 + (1.0 - xs[:-1]) ** 2)
 
-
 def himmelblau(xs: np.ndarray) -> float:
-    """Himmelblau's function
+    """
+    Himmelblau's function.
+
+    Parameters
+    ----------
+    xs : np.ndarray
+        Input array of shape (2,).
 
+    Returns
+    -------
+    float
+        The calculated value of Himmelblau's function.
+
+    Notes
+    -----
     It has four global minima f(xs) = 0 at
     xs=[3,2], [-2.805118..., 3.131312...], [-3.779310..., -3.2831860], and [3.584428..., -1.848126...].
     """
@@ -71,9 +131,9 @@ def himmelblau(xs: np.ndarray) -> float:
         )
     return (xs[0] ** 2 + xs[1] - 11.0) ** 2 + (xs[0] + xs[1] ** 2 - 7.0) ** 2
 
-
 def linear_regression_test(xs: np.ndarray) -> float:
-    """ Negative log likelihood of linear regression with Gaussian noise N(0,sigma)
+    """
+    Negative log likelihood of linear regression with Gaussian noise N(0,sigma)
 
     y = ax + b
 
@@ -83,7 +143,17 @@ def linear_regression_test(xs: np.ndarray) -> float:
     a = xs[0], b = xs[1], log(sigma**2) = xs[2]
 
     It has a global minimum f(xs) = 1.005071.. at
-    xs = [0.628571..., 0.8, -0.664976...].    
+    xs = [0.628571..., 0.8, -0.664976...].
+
+    Parameters
+    ----------
+    xs : np.ndarray
+        Input array of model parameters.
+
+    Returns
+    -------
+    float
+        The negative log likelihood of the linear regression model.
     """
     if xs.shape[0] != 3:
         raise RuntimeError(
@@ -98,7 +168,6 @@ def linear_regression_test(xs: np.ndarray) -> float:
         n * xs[2] + np.sum((xs[0] * xdata + xs[1] - ydata) ** 2) / np.exp(xs[2])
     )
 
-
 class Solver(odatse.solver.function.Solver):
     """Function Solver with pre-defined benchmark functions"""
 
@@ -112,6 +181,7 @@ def __init__(self, info: odatse.Info) -> None:
         Parameters
         ----------
         info: Info
+            Information object containing solver configuration.
         """
         super().__init__(info)
         self._name = "analytical"
diff --git a/src/odatse/solver/function.py b/src/odatse/solver/function.py
index 0d9088c1..fcf58220 100644
--- a/src/odatse/solver/function.py
+++ b/src/odatse/solver/function.py
@@ -25,6 +25,9 @@
 
 
 class Solver(odatse.solver.SolverBase):
+    """
+    Solver class for evaluating functions with given parameters.
+    """
     x: np.ndarray
     fx: float
     _func: Optional[Callable[[np.ndarray], float]]
@@ -36,6 +39,7 @@ def __init__(self, info: odatse.Info) -> None:
         Parameters
         ----------
         info: Info
+            Information object containing solver configuration.
         """
         super().__init__(info)
         self._name = "function"
@@ -45,6 +49,25 @@ def __init__(self, info: odatse.Info) -> None:
         self.delay = info.solver.get("delay", 0.0)
 
     def evaluate(self, x: np.ndarray, args: Tuple = (), nprocs: int = 1, nthreads: int = 1) -> float:
+        """
+        Evaluate the function with given parameters.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Input array for the function.
+        args : Tuple, optional
+            Additional arguments for the function.
+        nprocs : int, optional
+            Number of processes to use.
+        nthreads : int, optional
+            Number of threads to use.
+
+        Returns
+        -------
+        float
+            Result of the function evaluation.
+        """
         self.prepare(x, args)
         cwd = os.getcwd()
         os.chdir(self.work_dir)
@@ -54,9 +77,34 @@ def evaluate(self, x: np.ndarray, args: Tuple = (), nprocs: int = 1, nthreads: i
         return result
 
     def prepare(self, x: np.ndarray, args = ()) -> None:
+        """
+        Prepare the solver with the given parameters.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Input array for the function.
+        args : tuple, optional
+            Additional arguments for the function.
+        """
         self.x = x
 
     def run(self, nprocs: int = 1, nthreads: int = 1) -> None:
+        """
+        Run the function evaluation.
+
+        Parameters
+        ----------
+        nprocs : int, optional
+            Number of processes to use.
+        nthreads : int, optional
+            Number of threads to use.
+
+        Raises
+        ------
+        RuntimeError
+            If the function is not set.
+        """
         if self._func is None:
             raise RuntimeError(
                 "ERROR: function is not set. Make sure that `set_function` is called."
@@ -67,7 +115,23 @@ def run(self, nprocs: int = 1, nthreads: int = 1) -> None:
             time.sleep(self.delay)
 
     def get_results(self) -> float:
+        """
+        Get the results of the function evaluation.
+
+        Returns
+        -------
+        float
+            Result of the function evaluation.
+        """
         return self.fx
 
     def set_function(self, f: Callable[[np.ndarray], float]) -> None:
+        """
+        Set the function to be evaluated.
+
+        Parameters
+        ----------
+        f : Callable[[np.ndarray], float]
+            Function to be evaluated.
+        """
         self._func = f
diff --git a/src/odatse/util/graph.py b/src/odatse/util/graph.py
index e6afb69c..ac7a0f48 100644
--- a/src/odatse/util/graph.py
+++ b/src/odatse/util/graph.py
@@ -21,6 +21,15 @@
 
 
 def is_connected(nnlist: List[List[int]]) -> bool:
+    """
+    Check if the graph represented by the neighbor list is connected.
+
+    Parameters:
+    nnlist (List[List[int]]): A list of lists where each sublist represents the neighbors of a node.
+
+    Returns:
+    bool: True if the graph is connected, False otherwise.
+    """
     nnodes = len(nnlist)
     visited = np.full(nnodes, False)
     nvisited = 1
@@ -37,13 +46,21 @@ def is_connected(nnlist: List[List[int]]) -> bool:
 
 
 def is_bidirectional(nnlist: List[List[int]]) -> bool:
+    """
+    Check if the graph represented by the neighbor list is bidirectional.
+
+    Parameters:
+    nnlist (List[List[int]]): A list of lists where each sublist represents the neighbors of a node.
+
+    Returns:
+    bool: True if the graph is bidirectional, False otherwise.
+    """
     for i in range(len(nnlist)):
         for j in nnlist[i]:
             if i not in nnlist[j]:
                 return False
     return True
 
-
 if __name__ == "__main__":
     filename = "./neighborlist.txt"
     nnlist = []
diff --git a/src/odatse/util/limitation.py b/src/odatse/util/limitation.py
index 8a26c067..6663dd2d 100644
--- a/src/odatse/util/limitation.py
+++ b/src/odatse/util/limitation.py
@@ -6,22 +6,62 @@
 
 
 class LimitationBase(metaclass=ABCMeta):
+    """
+    Abstract base class for limitations.
+    """
+
     @abstractmethod
     def __init__(self, is_limitary: bool):
+        """
+        Initialize the limitation.
+
+        :param is_limitary: Boolean indicating if the limitation is active.
+        """
         self.is_limitary = is_limitary
 
     @abstractmethod
     def judge(self, x: np.ndarray) -> bool:
+        """
+        Abstract method to judge if the limitation is satisfied.
+
+        :param x: Input array to be judged.
+        :return: Boolean indicating if the limitation is satisfied.
+        """
         raise NotImplementedError
 
 class Unlimited(LimitationBase):
+    """
+    Class representing an unlimited (no limitation) condition.
+    """
+
     def __init__(self):
+        """
+        Initialize the unlimited condition.
+        """
         super().__init__(False)
+
     def judge(self, x: np.ndarray) -> bool:
+        """
+        Always returns True as there is no limitation.
+
+        :param x: Input array to be judged.
+        :return: Always True.
+        """
         return True
 
 class Inequality(LimitationBase):
+    """
+    Class representing an inequality limitation.
+    """
+
     def __init__(self, a: np.ndarray, b: np.ndarray, is_limitary: bool):
+        """
+        Initialize the inequality limitation.
+
+        :param a: Coefficient matrix.
+        :param b: Constant vector.
+        :param is_limitary: Boolean indicating if the limitation is active.
+        """
         super().__init__(is_limitary)
         if self.is_limitary:
             self.a = np.array(a)
@@ -31,6 +71,12 @@ def __init__(self, a: np.ndarray, b: np.ndarray, is_limitary: bool):
             self.ndim = a.shape[1]
 
     def judge(self, x: np.ndarray) -> bool:
+        """
+        Judge if the inequality limitation is satisfied.
+
+        :param x: Input array to be judged.
+        :return: Boolean indicating if the limitation is satisfied.
+        """
         if self.is_limitary:
             Ax_b = np.dot(self.a, x) + self.b
             judge_result = np.all(Ax_b > 0)
@@ -40,12 +86,15 @@ def judge(self, x: np.ndarray) -> bool:
 
     @classmethod
     def from_dict(cls, d):
+        """
+        Create an Inequality instance from a dictionary.
+
+        :param d: Dictionary containing 'co_a' and 'co_b' keys.
+        :return: Inequality instance.
+        """
         co_a: np.ndarray = read_matrix(d.get("co_a", []))
         co_b: np.ndarray = read_matrix(d.get("co_b", []))
 
-        # is_set_co_a = (co_a.size > 0 and co_a.ndim == 2 and co_a.shape[1] == dimension)
-        # is_set_co_b = (co_b.size > 0 and co_b.ndim == 2 and co_b.shape == (co_a.shape[0], 1))
-
         if co_a.size == 0:
             is_set_co_a = False
         else:
diff --git a/src/odatse/util/logger.py b/src/odatse/util/logger.py
index 6aeb8006..ca173bc4 100644
--- a/src/odatse/util/logger.py
+++ b/src/odatse/util/logger.py
@@ -33,6 +33,10 @@
 #   write_result
 
 class Logger:
+    """
+    Logger class to handle logging of calls, elapsed time, and optionally input and result data.
+    """
+
     logfile: Path
     buffer_size: int
     buffer: List[str]
@@ -50,7 +54,18 @@ def __init__(self, info: Optional[odatse.Info] = None,
                  write_result: bool = False,
                  params: Optional[Dict[str,Any]] = None,
                  **rest) -> None:
-
+        """
+        Initialize the Logger.
+
+        Parameters:
+        info (Optional[odatse.Info]): Information object containing logging parameters.
+        buffer_size (int): Size of the buffer before writing to the log file.
+        filename (str): Name of the log file.
+        write_input (bool): Flag to indicate if input should be logged.
+        write_result (bool): Flag to indicate if result should be logged.
+        params (Optional[Dict[str,Any]]): Additional parameters for logging.
+        **rest: Additional keyword arguments.
+        """
         if info is not None:
             info_log = info.runner.get("log", {})
         else:
@@ -67,9 +82,21 @@ def __init__(self, info: Optional[odatse.Info] = None,
         self.buffer = []
 
     def is_active(self) -> bool:
+        """
+        Check if logging is active.
+
+        Returns:
+        bool: True if logging is active, False otherwise.
+        """
         return self.buffer_size > 0
 
     def prepare(self, proc_dir: Path) -> None:
+        """
+        Prepare the log file for writing.
+
+        Parameters:
+        proc_dir (Path): Directory where the log file will be created.
+        """
         if not self.is_active():
             return
 
@@ -88,6 +115,14 @@ def prepare(self, proc_dir: Path) -> None:
             f.write("\n")
 
     def count(self, x: np.ndarray, args, result: float) -> None:
+        """
+        Log a call with input and result data.
+
+        Parameters:
+        x (np.ndarray): Input data.
+        args: Additional arguments.
+        result (float): Result data.
+        """
         if not self.is_active():
             return
 
@@ -112,9 +147,12 @@ def count(self, x: np.ndarray, args, result: float) -> None:
             self.write()
 
     def write(self) -> None:
+        """
+        Write the buffered log entries to the log file.
+        """
         if not self.is_active():
             return
         with open(self.logfile, "a") as f:
             for w in self.buffer:
                 f.write(w)
-        self.buffer.clear()
+        self.buffer.clear()
\ No newline at end of file
diff --git a/src/odatse/util/mapping.py b/src/odatse/util/mapping.py
index 059f404b..d8cb0719 100644
--- a/src/odatse/util/mapping.py
+++ b/src/odatse/util/mapping.py
@@ -19,32 +19,69 @@
 import copy
 import numpy as np
 
-from .read_matrix import read_matrix, read_vector
+from .read_matrix import read_matrix
 
 # type hints
 from typing import Optional
 
 class MappingBase:
+    """
+    Base class for mapping operations.
+    """
+
     def __init__(self):
         pass
 
     def __call__(self, x: np.ndarray) -> np.ndarray:
-        raise NotImplemented
+        """
+        Apply the mapping to the input array.
+
+        Parameters:
+        x (np.ndarray): Input array.
+
+        Returns:
+        np.ndarray: Mapped array.
+        """
+        raise NotImplementedError
 
 
 class TrivialMapping(MappingBase):
+    """
+    A trivial mapping that returns the input array unchanged.
+    """
+
     def __init__(self):
         super().__init__()
 
     def __call__(self, x: np.ndarray) -> np.ndarray:
+        """
+        Return the input array unchanged.
+
+        Parameters:
+        x (np.ndarray): Input array.
+
+        Returns:
+        np.ndarray: The same input array.
+        """
         return x
 
 
 class Affine(MappingBase):
+    """
+    An affine mapping defined by a matrix A and a vector b.
+    """
+
     A: Optional[np.ndarray]
     b: Optional[np.ndarray]
 
     def __init__(self, A: Optional[np.ndarray] = None, b: Optional[np.ndarray] = None):
+        """
+        Initialize the affine mapping.
+
+        Parameters:
+        A (Optional[np.ndarray]): Transformation matrix.
+        b (Optional[np.ndarray]): Translation vector.
+        """
         # copy arguments
         self.A = np.array(A) if A is not None else None
         self.b = np.array(b) if b is not None else None
@@ -60,8 +97,16 @@ def __init__(self, A: Optional[np.ndarray] = None, b: Optional[np.ndarray] = Non
             if not self.A.shape[0] == self.b.shape[0]:
                 raise ValueError("shape of A and b mismatch")
 
-
     def __call__(self, x: np.ndarray) -> np.ndarray:
+        """
+        Apply the affine mapping to the input array.
+
+        Parameters:
+        x (np.ndarray): Input array.
+
+        Returns:
+        np.ndarray: Mapped array.
+        """
         if self.A is None:
             ret = copy.copy(x)
         else:
@@ -73,6 +118,15 @@ def __call__(self, x: np.ndarray) -> np.ndarray:
 
     @classmethod
     def from_dict(cls, d):
+        """
+        Create an Affine instance from a dictionary.
+
+        Parameters:
+        d (dict): Dictionary containing 'A' and 'b' keys.
+
+        Returns:
+        Affine: An instance of the Affine class.
+        """
         A: Optional[np.ndarray] = read_matrix(d.get("A", []))
         b: Optional[np.ndarray] = read_matrix(d.get("b", []))
 
@@ -92,7 +146,7 @@ def from_dict(cls, d):
             if not (b.ndim == 2 and b.shape[1] == 1):
                 raise ValueError("b should be a column vector")
             if not (A is not None and b.shape[0] == A.shape[0]):
-                raise ValueError("shape of A and b does not match")
+                raise ValueError("shape of A and b mismatch")
             b = b.reshape(-1)
 
-        return cls(A, b)
+        return cls(A, b)
\ No newline at end of file
diff --git a/src/odatse/util/neighborlist.py b/src/odatse/util/neighborlist.py
index 762286a3..0a82f9ea 100644
--- a/src/odatse/util/neighborlist.py
+++ b/src/odatse/util/neighborlist.py
@@ -34,6 +34,10 @@
 
 
 class Cells:
+    """
+    A class to represent a grid of cells for spatial partitioning.
+    """
+
     cells: List[Set[int]]
     dimension: int
     mins: np.ndarray
@@ -43,6 +47,14 @@ class Cells:
     cellsize: float
 
     def __init__(self, mins: np.ndarray, maxs: np.ndarray, cellsize: float):
+        """
+        Initialize the Cells object.
+
+        Parameters:
+        mins (np.ndarray): The minimum coordinates of the grid.
+        maxs (np.ndarray): The maximum coordinates of the grid.
+        cellsize (float): The size of each cell.
+        """
         self.dimension = len(mins)
         self.mins = mins
         Ls = (maxs - mins) * 1.001
@@ -53,12 +65,39 @@ def __init__(self, mins: np.ndarray, maxs: np.ndarray, cellsize: float):
         self.cells = [set() for _ in range(self.ncell)]
 
     def coord2cellindex(self, x: np.ndarray) -> int:
+        """
+        Convert coordinates to a cell index.
+
+        Parameters:
+        x (np.ndarray): The coordinates to convert.
+
+        Returns:
+        int: The index of the cell.
+        """
         return self.cellcoord2cellindex(self.coord2cellcoord(x))
 
     def coord2cellcoord(self, x: np.ndarray) -> np.ndarray:
+        """
+        Convert coordinates to cell coordinates.
+
+        Parameters:
+        x (np.ndarray): The coordinates to convert.
+
+        Returns:
+        np.ndarray: The cell coordinates.
+        """
         return np.floor((x - self.mins) / self.cellsize).astype(np.int64)
 
     def cellcoord2cellindex(self, ns: np.ndarray) -> int:
+        """
+        Convert cell coordinates to a cell index.
+
+        Parameters:
+        ns (np.ndarray): The cell coordinates to convert.
+
+        Returns:
+        int: The index of the cell.
+        """
         index = 0
         oldN = 1
         for n, N in zip(ns, self.Ns):
@@ -68,6 +107,15 @@ def cellcoord2cellindex(self, ns: np.ndarray) -> int:
         return index
 
     def cellindex2cellcoord(self, index: int) -> np.ndarray:
+        """
+        Convert a cell index to cell coordinates.
+
+        Parameters:
+        index (int): The index of the cell.
+
+        Returns:
+        np.ndarray: The cell coordinates.
+        """
         ns = np.zeros(self.dimension, dtype=np.int64)
         for d in range(self.dimension):
             d = self.dimension - d - 1
@@ -77,6 +125,15 @@ def cellindex2cellcoord(self, index: int) -> np.ndarray:
         return ns
 
     def out_of_bound(self, ns: np.ndarray) -> bool:
+        """
+        Check if cell coordinates are out of bounds.
+
+        Parameters:
+        ns (np.ndarray): The cell coordinates to check.
+
+        Returns:
+        bool: True if out of bounds, False otherwise.
+        """
         if np.any(ns < 0):
             return True
         if np.any(ns >= self.Ns):
@@ -84,6 +141,15 @@ def out_of_bound(self, ns: np.ndarray) -> bool:
         return False
 
     def neighborcells(self, index: int) -> List[int]:
+        """
+        Get the indices of neighboring cells.
+
+        Parameters:
+        index (int): The index of the cell.
+
+        Returns:
+        List[int]: The indices of the neighboring cells.
+        """
         neighbors: List[int] = []
         center_coord = self.cellindex2cellcoord(index)
         for diff in itertools.product([-1, 0, 1], repeat=self.dimension):
@@ -94,7 +160,6 @@ def neighborcells(self, index: int) -> List[int]:
             neighbors.append(other_coord_index)
         return neighbors
 
-
 def make_neighbor_list_cell(
     X: np.ndarray,
     radius: float,
@@ -102,6 +167,19 @@ def make_neighbor_list_cell(
     show_progress: bool,
     comm: mpi.Comm = None,
 ) -> List[List[int]]:
+    """
+    Create a neighbor list using cell-based spatial partitioning.
+
+    Parameters:
+    X (np.ndarray): The coordinates of the points.
+    radius (float): The radius within which neighbors are considered.
+    allow_selfloop (bool): Whether to allow self-loops in the neighbor list.
+    show_progress (bool): Whether to show a progress bar.
+    comm (mpi.Comm, optional): The MPI communicator.
+
+    Returns:
+    List[List[int]]: The neighbor list.
+    """
     if comm is None:
         mpisize = 1
         mpirank = 0
@@ -153,6 +231,19 @@ def make_neighbor_list_naive(
     show_progress: bool,
     comm: mpi.Comm = None,
 ) -> List[List[int]]:
+    """
+    Create a neighbor list using a naive all-pairs approach.
+
+    Parameters:
+    X (np.ndarray): The coordinates of the points.
+    radius (float): The radius within which neighbors are considered.
+    allow_selfloop (bool): Whether to allow self-loops in the neighbor list.
+    show_progress (bool): Whether to show a progress bar.
+    comm (mpi.Comm, optional): The MPI communicator.
+
+    Returns:
+    List[List[int]]: The neighbor list.
+    """
     if comm is None:
         mpisize = 1
         mpirank = 0
@@ -195,6 +286,20 @@ def make_neighbor_list(
     show_progress: bool = False,
     comm: mpi.Comm = None,
 ) -> List[List[int]]:
+    """
+    Create a neighbor list for given points.
+
+    Parameters:
+    X (np.ndarray): The coordinates of the points.
+    radius (float): The radius within which neighbors are considered.
+    allow_selfloop (bool): Whether to allow self-loops in the neighbor list.
+    check_allpairs (bool): Whether to use the naive all-pairs approach.
+    show_progress (bool): Whether to show a progress bar.
+    comm (mpi.Comm, optional): The MPI communicator.
+
+    Returns:
+    List[List[int]]: The neighbor list.
+    """
     if check_allpairs:
         return make_neighbor_list_naive(
             X,
@@ -214,6 +319,16 @@ def make_neighbor_list(
 
 
 def load_neighbor_list(filename: PathLike, nnodes: int = None) -> List[List[int]]:
+    """
+    Load a neighbor list from a file.
+
+    Parameters:
+    filename (PathLike): The path to the file containing the neighbor list.
+    nnodes (int, optional): The number of nodes. If None, it will be determined from the file.
+
+    Returns:
+    List[List[int]]: The neighbor list.
+    """
     if nnodes is None:
         nnodes = 0
         with open(filename) as f:
@@ -222,6 +337,7 @@ def load_neighbor_list(filename: PathLike, nnodes: int = None) -> List[List[int]
                 if len(line) == 0:
                     continue
                 nnodes += 1
+
     neighbor_list: List[List[int]] = [[] for _ in range(nnodes)]
     with open(filename) as f:
         for line in f:
@@ -241,6 +357,15 @@ def write_neighbor_list(
     radius: float = None,
     unit: np.ndarray = None,
 ):
+    """
+    Write the neighbor list to a file.
+
+    Parameters:
+    filename (str): The path to the output file.
+    nnlist (List[List[int]]): The neighbor list to write.
+    radius (float, optional): The neighborhood radius. Defaults to None.
+    unit (np.ndarray, optional): The unit for each coordinate. Defaults to None.
+    """
     with open(filename, "w") as f:
         if radius is not None:
             f.write(f"# radius = {radius}\n")
@@ -255,7 +380,6 @@ def write_neighbor_list(
                 f.write(f" {o}")
             f.write("\n")
 
-
 def main():
     import argparse
 
diff --git a/src/odatse/util/read_matrix.py b/src/odatse/util/read_matrix.py
index 819cfe4f..115fd947 100644
--- a/src/odatse/util/read_matrix.py
+++ b/src/odatse/util/read_matrix.py
@@ -20,6 +20,18 @@
 
 
 def read_vector(inp: Union[str, List[float]]) -> np.ndarray:
+    """
+    Converts an input string or list of floats into a numpy array vector.
+
+    Parameters:
+    inp (Union[str, List[float]]): Input data, either as a space-separated string of numbers or a list of floats.
+
+    Returns:
+    np.ndarray: A numpy array representing the vector.
+
+    Raises:
+    RuntimeError: If the input is not a vector.
+    """
     if isinstance(inp, str):
         vlist = [float(w) for w in inp.split()]
     else:
@@ -30,8 +42,19 @@ def read_vector(inp: Union[str, List[float]]) -> np.ndarray:
         raise RuntimeError(msg)
     return v
 
-
 def read_matrix(inp: Union[str, List[List[float]]]) -> np.ndarray:
+    """
+    Converts an input string or list of lists of floats into a numpy array matrix.
+
+    Parameters:
+    inp (Union[str, List[List[float]]]): Input data, either as a string with rows of space-separated numbers or a list of lists of floats.
+
+    Returns:
+    np.ndarray: A numpy array representing the matrix.
+
+    Raises:
+    RuntimeError: If the input is not a matrix.
+    """
     if isinstance(inp, str):
         Alist: List[List[float]] = []
         for line in inp.split("\n"):
diff --git a/src/odatse/util/resampling.py b/src/odatse/util/resampling.py
index 356a4a83..0a8d1a5d 100644
--- a/src/odatse/util/resampling.py
+++ b/src/odatse/util/resampling.py
@@ -35,43 +35,95 @@ def sample(self, rs: np.random.RandomState, size=None) -> Union[int, np.ndarray]
 
 
 class BinarySearch(Resampler):
+    """
+    A resampler that uses binary search to sample based on given weights.
+    """
     weights_accumulated: List[float]
     wmax: float
 
     def __init__(self, weights: Iterable):
+        """
+        Initialize the BinarySearch resampler with the given weights.
+
+        :param weights: An iterable of weights.
+        """
         self.reset(weights)
 
     def reset(self, weights: Iterable):
+        """
+        Reset the resampler with new weights.
+
+        :param weights: An iterable of weights.
+        """
         self.weights_accumulated = list(itertools.accumulate(weights))
         self.wmax = self.weights_accumulated[-1]
 
     @typing.overload
     def sample(self, rs: np.random.RandomState) -> int:
+        """
+        Sample a single index based on the weights.
+
+        :param rs: A random state for generating random numbers.
+        :return: A single sampled index.
+        """
         ...
 
     @typing.overload
     def sample(self, rs: np.random.RandomState, size) -> np.ndarray:
+        """
+        Sample multiple indices based on the weights.
+
+        :param rs: A random state for generating random numbers.
+        :param size: The number of samples to generate.
+        :return: An array of sampled indices.
+        """
         ...
 
     def sample(self, rs: np.random.RandomState, size=None) -> Union[int, np.ndarray]:
+        """
+        Sample indices based on the weights.
+
+        :param rs: A random state for generating random numbers.
+        :param size: The number of samples to generate. If None, a single sample is generated.
+        :return: A single sampled index or an array of sampled indices.
+        """
         if size is None:
             return self._sample(self.wmax * rs.rand())
         else:
             return np.array([self._sample(r) for r in self.wmax * rs.rand(size)])
 
     def _sample(self, r: float) -> int:
+        """
+        Perform a binary search to find the index corresponding to the given random number.
+
+        :param r: A random number scaled by the maximum weight.
+        :return: The index corresponding to the random number.
+        """
         return typing.cast(int, np.searchsorted(self.weights_accumulated, r))
 
 
 class WalkerTable(Resampler):
+    """
+    A resampler that uses Walker's alias method to sample based on given weights.
+    """
     N: int
     itable: np.ndarray
     ptable: np.ndarray
 
     def __init__(self, weights: Iterable):
+        """
+        Initialize the WalkerTable resampler with the given weights.
+
+        :param weights: An iterable of weights.
+        """
         self.reset(weights)
 
     def reset(self, weights: Iterable):
+        """
+        Reset the resampler with new weights.
+
+        :param weights: An iterable of weights.
+        """
         self.ptable = np.array(weights).astype(np.float64).flatten()
         self.N = len(self.ptable)
         self.itable = np.full(self.N, -1)
@@ -93,13 +145,33 @@ def reset(self, weights: Iterable):
 
     @typing.overload
     def sample(self, rs: np.random.RandomState) -> int:
+        """
+        Sample a single index based on the weights.
+
+        :param rs: A random state for generating random numbers.
+        :return: A single sampled index.
+        """
         ...
 
     @typing.overload
     def sample(self, rs: np.random.RandomState, size) -> np.ndarray:
+        """
+        Sample multiple indices based on the weights.
+
+        :param rs: A random state for generating random numbers.
+        :param size: The number of samples to generate.
+        :return: An array of sampled indices.
+        """
         ...
 
     def sample(self, rs: np.random.RandomState, size=None) -> Union[int, np.ndarray]:
+        """
+        Sample indices based on the weights.
+
+        :param rs: A random state for generating random numbers.
+        :param size: The number of samples to generate. If None, a single sample is generated.
+        :return: A single sampled index or an array of sampled indices.
+        """
         if size is None:
             r = rs.rand() * self.N
             return self._sample(r)
@@ -107,10 +179,16 @@ def sample(self, rs: np.random.RandomState, size=None) -> Union[int, np.ndarray]
             r = rs.rand(size) * self.N
             i = np.floor(r).astype(np.int64)
             p = r - i
-            ret =  np.where(p < self.ptable[i], i, self.itable[i])
+            ret = np.where(p < self.ptable[i], i, self.itable[i])
             return ret
 
     def _sample(self, r: float) -> int:
+        """
+        Perform a sampling operation based on the given random number.
+
+        :param r: A random number scaled by the number of weights.
+        :return: The index corresponding to the random number.
+        """
         i = int(np.floor(r))
         p = r - i
         if p < self.ptable[i]:
@@ -118,7 +196,6 @@ def _sample(self, r: float) -> int:
         else:
             return self.itable[i]
 
-
 if __name__ == "__main__":
     import argparse
 
diff --git a/src/odatse/util/separateT.py b/src/odatse/util/separateT.py
index 4951c96a..38e8cf58 100644
--- a/src/odatse/util/separateT.py
+++ b/src/odatse/util/separateT.py
@@ -34,6 +34,20 @@ def separateT(
     use_beta: bool,
     buffer_size: int = 10000,
 ) -> None:
+    """
+    Separates and processes temperature data for quantum beam diffraction experiments.
+
+    Parameters:
+    Ts (np.ndarray): Array of temperature values.
+    nwalkers (int): Number of walkers.
+    output_dir (PathLike): Directory to store the output files.
+    comm (Optional[mpi.Comm]): MPI communicator for parallel processing.
+    use_beta (bool): Flag to determine if beta values are used instead of temperature.
+    buffer_size (int, optional): Size of the buffer for reading input data. Default is 10000.
+
+    Returns:
+    None
+    """
     if comm is None:
         mpisize = 1
         mpirank = 0
diff --git a/tests/bayes/do.sh b/tests/bayes/do.sh
index d1244560..c25a0b92 100644
--- a/tests/bayes/do.sh
+++ b/tests/bayes/do.sh
@@ -1,14 +1,20 @@
 #!/bin/sh
 
+# Remove the output directory if it exists
 rm -rf output
 
+# Run the Python script with the input file and measure the time taken
 time python3 ../../src/odatse_main.py input.toml
 
+# Define the result file path
 resfile=output/BayesData.txt
 
+# Compare the result file with the reference file
 echo diff $resfile ref.txt
 res=0
 diff $resfile ref.txt || res=$?
+
+# Check the result of the diff command
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
@@ -16,4 +22,3 @@ else
   echo TEST FAILED: $resfile and ref.txt differ
   false
 fi
-
diff --git a/tests/bayes_continue/do.sh b/tests/bayes_continue/do.sh
index f488bc82..1685c3c3 100644
--- a/tests/bayes_continue/do.sh
+++ b/tests/bayes_continue/do.sh
@@ -1,23 +1,31 @@
 #!/bin/sh
 
+# Command to run the main Python script
 CMD="python3 -u ../../src/odatse_main.py"
 
+# Remove the output1 directory if it exists
 rm -rf output1
 
+# Run the main script with input1a.toml and input1b.toml
 time $CMD input1a.toml
 time $CMD --cont input1b.toml
 
+# Remove the output2 directory if it exists
 rm -rf output2
 
+# Run the main script with input2.toml
 time $CMD input2.toml
 
-
+# Define the result and reference files
 resfile=output1/BayesData.txt
 reffile=output2/BayesData.txt
 
+# Compare the result and reference files
 echo diff $resfile $reffile
 res=0
 diff $resfile $reffile || res=$?
+
+# Check the result of the comparison
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
diff --git a/tests/exchange/do.sh b/tests/exchange/do.sh
index 4f5e6447..0f5efac2 100644
--- a/tests/exchange/do.sh
+++ b/tests/exchange/do.sh
@@ -1,19 +1,26 @@
 #!/bin/sh
 
+# Remove the output directory if it exists
 rm -rf output
 
+# Run the Python script using mpiexec with 4 processes and measure the time taken
 time mpiexec --oversubscribe -np 4 python3 ../../src/odatse_main.py input.toml
 
+# Define the result file path
 resfile=output/best_result.txt
 
+# Compare the result file with the reference file
 echo diff $resfile ref.txt
 res=0
 diff $resfile ref.txt || res=$?
+
+# Check the result of the diff command
 if [ $res -eq 0 ]; then
+  # If the files are the same, print TEST PASS
   echo TEST PASS
   true
 else
+  # If the files differ, print TEST FAILED with the result file path
   echo TEST FAILED: $resfile and ref.txt differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/exchange_continue/do.sh b/tests/exchange_continue/do.sh
index 451e9fa0..d84dfc2b 100644
--- a/tests/exchange_continue/do.sh
+++ b/tests/exchange_continue/do.sh
@@ -1,30 +1,38 @@
 #!/bin/sh
 
+# Command to run the Python script with MPI
 CMD="mpiexec --oversubscribe -np 4 python3 -u ../../src/odatse_main.py"
 
+# Remove the output1 directory if it exists
 rm -rf output1
 
+# Run the command with the first input file and measure the time taken
 time $CMD input1a.toml
+# Run the command with the continuation input file and measure the time taken
 time $CMD --cont input1b.toml
 
+# Remove the output2 directory if it exists
 rm -rf output2
 
+# Run the command with the second input file and measure the time taken
 time $CMD input2.toml
 
-
-#resfile=output1/best_result.txt
-#reffile=output2/best_result.txt
+# Define the result files to compare
+# resfile=output1/best_result.txt
+# reffile=output2/best_result.txt
 resfile=output1/result_T0.txt
 reffile=output2/result_T0.txt
 
+# Compare the result files and store the result of the comparison
 echo diff $resfile $reffile
 res=0
 diff $resfile $reffile || res=$?
+
+# Check the result of the comparison and print the appropriate message
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
 else
   echo TEST FAILED: $resfile and $reffile differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/exchange_mesh/do.sh b/tests/exchange_mesh/do.sh
index 09bf2322..40761310 100644
--- a/tests/exchange_mesh/do.sh
+++ b/tests/exchange_mesh/do.sh
@@ -1,28 +1,35 @@
 #!/bin/sh
 
+# Remove the output directory if it exists
 rm -rf output
 
+# Generate MeshData.txt using makemesh.py
 echo generate MeshData.txt
 time python3 ./makemesh.py > MeshData.txt
 
 echo
+# Generate neighborlist.txt using odatse_neighborlist.py with a radius of 0.11
 echo generate neighborlist.txt
 time python3 ../../src/odatse_neighborlist.py -r 0.11 MeshData.txt
 
 echo
+# Perform exchange Monte Carlo simulation using odatse_main.py with input.toml
 echo perform exchange mc
 time python3 ../../src/odatse_main.py input.toml
 
+# Define the result file path
 resfile=output/best_result.txt
 
+# Compare the result file with the reference file
 echo diff $resfile ref.txt
 res=0
 diff $resfile ref.txt || res=$?
 if [ $res -eq 0 ]; then
+  # Output TEST PASS if files are identical
   echo TEST PASS
   true
 else
+  # Output TEST FAILED if files differ
   echo TEST FAILED: $resfile and ref.txt differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/mapper/do.sh b/tests/mapper/do.sh
index 2f192d47..14aa65a7 100644
--- a/tests/mapper/do.sh
+++ b/tests/mapper/do.sh
@@ -1,19 +1,24 @@
 #!/bin/sh
 
+# Remove the output directory if it exists
 rm -rf output
 
+# Run the Python script with MPI, using 2 processes
 time mpiexec --oversubscribe -np 2 python3 ../../src/odatse_main.py input.toml
 
+# Define the result file path
 resfile=output/ColorMap.txt
 
+# Compare the result file with the reference file
 echo diff $resfile ref.txt
 res=0
 diff $resfile ref.txt || res=$?
+
+# Check the result of the diff command
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
 else
   echo TEST FAILED: $resfile and ref.txt differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/mapper_resume/do.sh b/tests/mapper_resume/do.sh
index 92cf7727..6a59618a 100644
--- a/tests/mapper_resume/do.sh
+++ b/tests/mapper_resume/do.sh
@@ -1,30 +1,43 @@
 #!/bin/sh
 
+# Command to run the main Python script
 CMD="python3 -u ../../src/odatse_main.py"
-#CMD="mpiexec -np 2 python3 -u ../../src/odatse_main.py"
+# Uncomment the following line to run with MPI
+# CMD="mpiexec -np 2 python3 -u ../../src/odatse_main.py"
 
+# Remove the output1 directory if it exists
 rm -rf output1
 
+# Run the command with a timeout of 12 seconds using input1.toml
 time timeout 12s $CMD input1.toml
 
+# Run the command with the --resume option using input1.toml
 time $CMD --resume input1.toml
 
+# Remove the output2 directory if it exists
 rm -rf output2
 
+# Run the command using input2.toml
 time $CMD input2.toml
 
-
+# Define the result and reference files
 resfile=output1/ColorMap.txt
 reffile=output2/ColorMap.txt
 
+# Print the diff command to be executed
 echo diff $resfile $reffile
+
+# Initialize the result variable
 res=0
+
+# Compare the result and reference files, update res if they differ
 diff $resfile $reffile || res=$?
+
+# Check if the files are identical
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
 else
   echo TEST FAILED: $resfile and $reffile differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/minsearch/do.sh b/tests/minsearch/do.sh
index bd4bf163..40982da5 100644
--- a/tests/minsearch/do.sh
+++ b/tests/minsearch/do.sh
@@ -1,19 +1,30 @@
 #!/bin/sh
 
+# Remove the output directory if it exists
 rm -rf output
 
+# Run the Python script with the input file and measure the time taken
 time python3 ../../src/odatse_main.py input.toml
 
+# Define the result file path
 resfile=output/res.txt
 
+# Display the diff command being executed
 echo diff $resfile ref.txt
+
+# Initialize the result variable
 res=0
+
+# Compare the result file with the reference file, update the result variable if they differ
 diff $resfile ref.txt || res=$?
+
+# Check the result of the diff command
 if [ $res -eq 0 ]; then
+  # If the files are the same, print TEST PASS
   echo TEST PASS
   true
 else
+  # If the files differ, print TEST FAILED with the file names
   echo TEST FAILED: $resfile and ref.txt differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/pamc/do.sh b/tests/pamc/do.sh
index 01838abb..29336be1 100644
--- a/tests/pamc/do.sh
+++ b/tests/pamc/do.sh
@@ -1,19 +1,24 @@
 #!/bin/sh
 
+# Remove the output directory if it exists
 rm -rf output
 
+# Run the Python script using MPI with 2 processes
 time mpiexec --oversubscribe -np 2 python3 -m mpi4py ../../src/odatse_main.py input.toml
 
+# Define the result file path
 resfile=output/best_result.txt
 
+# Compare the result file with the reference file
 echo diff $resfile ref.txt
 res=0
 diff $resfile ref.txt || res=$?
+
+# Check the result of the diff command
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
 else
   echo TEST FAILED: $resfile and ref.txt differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/pamc_continue/do.sh b/tests/pamc_continue/do.sh
index afd70103..c9968b37 100644
--- a/tests/pamc_continue/do.sh
+++ b/tests/pamc_continue/do.sh
@@ -1,30 +1,43 @@
 #!/bin/sh
 
+# Command to run the Python script with MPI
 CMD="mpiexec --oversubscribe -np 2 python3 -u ../../src/odatse_main.py"
 
+# Remove the output1 directory if it exists
 rm -rf output1
 
+# Run the command with input1a.toml and log the output
 time $CMD input1a.toml 2>&1 | tee run.log.1a
+
+# Run the command with input1b.toml in continuation mode and log the output
 time $CMD --cont input1b.toml 2>&1 | tee run.log.1b
 
+# Remove the output2 directory if it exists
 rm -rf output2
 
+# Run the command with input2.toml and log the output
 time $CMD input2.toml 2>&1 | tee run.log.2
 
-
-#resfile=output1/best_result.txt
-#reffile=output2/best_result.txt
+# Define the result and reference files for comparison
+# resfile=output1/best_result.txt
+# reffile=output2/best_result.txt
 resfile=output1/fx.txt
 reffile=output2/fx.txt
 
+# Print the diff command to be executed
 echo diff $resfile $reffile
+
+# Initialize the result variable
 res=0
+
+# Compare the result and reference files, update the result variable if they differ
 diff $resfile $reffile || res=$?
+
+# Check the result of the diff command and print the appropriate message
 if [ $res -eq 0 ]; then
   echo TEST PASS
   true
 else
   echo TEST FAILED: $resfile and $reffile differ
   false
-fi
-
+fi
\ No newline at end of file
diff --git a/tests/transform/do.sh b/tests/transform/do.sh
index db8cc6fb..a6ceae14 100644
--- a/tests/transform/do.sh
+++ b/tests/transform/do.sh
@@ -1,19 +1,28 @@
 #!/bin/sh
 
+# Remove the existing ColorMap.txt file from the output_transform directory
 rm -f output_transform/ColorMap.txt
+
+# Run the odatse_main.py script with the input_transform.toml configuration file
 python3 ../../src/odatse_main.py input_transform.toml
 
+# Remove the existing ColorMap.txt file from the output_meshlist directory
 rm -f output_meshlist/ColorMap.txt
+
+# Run the odatse_main.py script with the input_meshlist.toml configuration file
 python3 ../../src/odatse_main.py input_meshlist.toml
 
+# Calculate the difference between the ColorMap.txt files from both outputs
 res=$(
 paste output_transform/ColorMap.txt output_meshlist/ColorMap.txt \
   | awk 'BEGIN {diff = 0.0} {diff += ($2 - $(NF))**2} END {print diff/NR}'
 )
+
+# Check if the difference is zero
 if [ $res = 0 ]; then
   echo TEST PASS
   true
 else
   echo "TEST FAILED (diff = $res)"
   false
-fi
+fi
\ No newline at end of file