added more detailed comments to bond_analysis

analysis/src/oxDNA_analysis_tools/bond_analysis.py

@@ -19,26 +19,36 @@
                                               "designed_pairs",
                                               "input_file"])
 
+# This is the function which gets parallelized
 def compute(ctx:ComputeContext, chunk_size:int, chunk_id:int) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    # Create the oxpy context
     with oxpy.Context():
-        inp = oxpy.InputFile()
-        inp.init_from_filename(ctx.input_file)
-        inp["list_type"] = "cells"
+        inp = oxpy.InputFile()                      # Create an empty input file
+        inp.init_from_filename(ctx.input_file)      # Fill in input file with the provided input
+
+        # Modify the input file to analyze the target trajectory and select a chunk to analyze
+        inp["list_type"] = "cells"                   
         inp["trajectory_file"] = ctx.traj_info.path
         inp["analysis_bytes_to_skip"] = str(ctx.traj_info.idxs[chunk_id*chunk_size].offset)
         inp["confs_to_analyse"] = str(chunk_size)
+
+        # Use oxDNA's observables interface to define what parameter we want to measure
         inp["analysis_data_output_1"] = '{ \n name = stdout \n print_every = 1e10 \n col_1 = { \n id = my_obs \n type = hb_list \n } \n }'
 
         if (not inp["use_average_seq"] or inp.get_bool("use_average_seq")) and "RNA" in inp["interaction_type"]:
             log("Sequence dependence not set for RNA model, wobble base pairs will be ignored", level="warning")
 
+        # Start up an analysis backend with the input file we just made
         backend = oxpy.analysis.AnalysisBackend(inp)
 
+        # The parameters we're going to read from the trajectory
         i = 0
         count_correct_bonds = []
         count_incorrect_bonds = []
         tot_bonds = []
         out_array = np.zeros(ctx.top_info.nbases, dtype=int)
+
+        # Iterate over configurations in the target chunk and process the hydrogen bonds at each step
         while backend.read_next_configuration():
             conf_corr_bonds = 0
             conf_incorr_bonds = 0
@@ -62,7 +72,8 @@ def compute(ctx:ComputeContext, chunk_size:int, chunk_id:int) -> Tuple[np.ndarra
             count_incorrect_bonds.append(conf_incorr_bonds)
             tot_bonds.append(conf_tot_bonds)
 
-        return(np.array(tot_bonds), np.array(count_correct_bonds), np.array(count_incorrect_bonds), out_array)
+    # Send the processed data back to the main process
+    return(np.array(tot_bonds), np.array(count_correct_bonds), np.array(count_incorrect_bonds), out_array)
 
 
 def bond_analysis(traj_info:TrajInfo, top_info:TopInfo, pairs:Dict[int, int], inputfile:str, ncpus:int=1) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
@@ -83,13 +94,16 @@ def bond_analysis(traj_info:TrajInfo, top_info:TopInfo, pairs:Dict[int, int], in
             | Number of correct bonds among specified nucleotides at each step in the simulation
             | Per-nucleotide correct bond occupancy
     '''
+    # oat_multiprocessor requires all arguments to compute be passed as a single variable.  We use a namedtuple for this.
     ctx = ComputeContext(traj_info, top_info, pairs, inputfile)
 
+    # Pre-allocate the arrays where the results will be stored
     total_bonds = np.empty(traj_info.nconfs)
     correct_bonds = np.empty(traj_info.nconfs)
     incorrect_bonds = np.empty(traj_info.nconfs)
     nt_array = np.zeros(ctx.top_info.nbases, dtype=int)
 
+    # The callback function defines what to do with the returned data when a chunk is finished processing
     chunk_size = get_chunk_size()
     def callback(i, r):
         nonlocal total_bonds, correct_bonds, incorrect_bonds, nt_array
@@ -98,6 +112,7 @@ def callback(i, r):
         incorrect_bonds[i*chunk_size:i*chunk_size+len(r[2])] = r[2]
         nt_array += r[3]
 
+    # Run the <compute> function on <ncpus> over <nconfs> with arguments <ctx> and catch the return with <callback>
     oat_multiprocesser(traj_info.nconfs, ncpus, compute, callback, ctx)
 
     nt_array  = nt_array / traj_info.nconfs