diff --git a/include/disk_utils.h b/include/disk_utils.h
index dabeeeb0c..6aaa02426 100644
--- a/include/disk_utils.h
+++ b/include/disk_utils.h
@@ -75,7 +75,7 @@ namespace diskann {
                                      const _u64 nshards, unsigned max_degree,
                                      const std::string &output_vamana,
                                      const std::string &medoids_file, bool use_filters = false,
-      const std::string &labels_to_medoids_file = std::string(""));
+                                     const std::string &labels_to_medoids_file = std::string(""));
 
 DISKANN_DLLEXPORT void extract_shard_labels(
       const std::string &in_label_file, const std::string &shard_ids_bin,
diff --git a/include/index.h b/include/index.h
index 7ea309e34..a866591bf 100644
--- a/include/index.h
+++ b/include/index.h
@@ -192,6 +192,8 @@ namespace diskann {
     DISKANN_DLLEXPORT consolidation_report
     consolidate_deletes(const Parameters &parameters);
 
+    DISKANN_DLLEXPORT void prune_all_nbrs(const Parameters &parameters);
+
     DISKANN_DLLEXPORT bool is_index_saved();
 
     // repositions frozen points to the end of _data - if they have been moved
@@ -238,7 +240,7 @@ namespace diskann {
     // determines navigating node of the graph by calculating medoid of datafopt
     unsigned calculate_entry_point();
 
-    void parse_label_file(const std::string &map_file);
+    size_t parse_label_file(const std::string &map_file);
 
     std::pair<uint32_t, uint32_t> iterate_to_fixed_point(
         const T *node_coords, const unsigned Lindex,
diff --git a/src/disk_utils.cpp b/src/disk_utils.cpp
index 9b7c3ca24..e6aa02fdb 100644
--- a/src/disk_utils.cpp
+++ b/src/disk_utils.cpp
@@ -323,6 +323,7 @@ namespace diskann {
       }
 
       std::ofstream mapping_writer(labels_to_medoids_file);
+      assert(mapping_writer.is_open());
       for (auto iter : global_label_to_medoids) {
         mapping_writer << iter.first << ", ";
         auto &vec = iter.second;
@@ -534,6 +535,7 @@ namespace diskann {
       diskann::cout << labels_per_point.size() << " is the new number of points"
                     << std::endl;
       std::ofstream label_writer(out_labels_file);
+      assert(label_writer.is_open());
       for (_u32 i = 0; i < labels_per_point.size(); i++) {
         for (_u32 j = 0; j < (labels_per_point[i].size() - 1); j++) {
           label_writer << labels_per_point[i][j] << ",";
@@ -551,6 +553,7 @@ namespace diskann {
       data = (T *) std::realloc((void *) data,
                                 labels_per_point.size() * ndims * sizeof(T));
       std::ofstream dummy_writer(out_metadata_file);
+      assert(dummy_writer.is_open());
       for (auto i = dummy_pt_ids.begin(); i != dummy_pt_ids.end(); i++) {
         dummy_writer << i->first << "," << i->second << std::endl;
         std::memcpy(data + i->first * ndims, data + i->second * ndims,
@@ -566,7 +569,7 @@ namespace diskann {
       const std::string &in_label_file, const std::string &shard_ids_bin,
       const std::string &shard_label_file) {  // assumes ith row is for ith
                                               // point in labels file
-    std::cout << "Extracting labels for shard" << std::endl;
+    diskann::cout << "Extracting labels for shard" << std::endl;
 
     _u32 *ids = nullptr;
     _u64  num_ids, tmp_dim;
@@ -578,6 +581,8 @@ namespace diskann {
     std::ifstream label_reader(in_label_file);
     std::ofstream label_writer(shard_label_file);
 
+    assert(label_reader.is_open());
+    assert(label_writer.is_open());
     if (label_reader && label_writer) {
       while (std::getline(label_reader, cur_line)) {
         if (shard_counter >= num_ids) {
diff --git a/src/index.cpp b/src/index.cpp
index 35eab43c3..344ca86d6 100644
--- a/src/index.cpp
+++ b/src/index.cpp
@@ -279,9 +279,9 @@ namespace diskann {
         if (_filter_to_medoid_id.size() > 0) {
           std::ofstream medoid_writer(std::string(filename) +
                                       "_labels_to_medoids.txt");
+          assert(medoid_writer.is_open());
           for (auto iter : _filter_to_medoid_id) {
             medoid_writer << iter.first << ", " << iter.second << std::endl;
-            // std::cout << iter.first << ", " << iter.second << std::endl;
           }
           medoid_writer.close();
         }
@@ -289,12 +289,14 @@ namespace diskann {
         if (_use_universal_label) {
           std::ofstream universal_label_writer(std::string(filename) +
                                                "_universal_label.txt");
+          assert(universal_label_writer.is_open());
           universal_label_writer << _universal_label << std::endl;
           universal_label_writer.close();
         }
 
         if (_pts_to_labels.size() > 0) {
           std::ofstream label_writer(std::string(filename) + "_labels.txt");
+          assert(label_writer.is_open());
           for (_u32 i = 0; i < _pts_to_labels.size(); i++) {
             for (_u32 j = 0; j < (_pts_to_labels[i].size() - 1); j++) {
               label_writer << _pts_to_labels[i][j] << ",";
@@ -477,52 +479,12 @@ namespace diskann {
 
     _has_built = true;
 
-    size_t tags_file_num_pts = 0, graph_num_pts = 0, data_file_num_pts = 0;
+    size_t tags_file_num_pts = 0, graph_num_pts = 0, data_file_num_pts = 0, label_num_pts = 0;
 
     std::string mem_index_file(filename);
     std::string labels_file = mem_index_file + "_labels.txt";
     std::string labels_to_medoids = mem_index_file + "_labels_to_medoids.txt";
-    if (file_exists(labels_file)) {
-      parse_label_file(labels_file);
-      if (file_exists(labels_to_medoids)) {
-        std::ifstream medoid_stream(labels_to_medoids);
-
-        std::string line, token;
-        unsigned    line_cnt = 0;
-
-        _filter_to_medoid_id.clear();
-
-        while (std::getline(medoid_stream, line)) {
-          std::istringstream iss(line);
-          _u32               cnt = 0;
-          _u32               medoid = 0;
-          label              label;
-          while (std::getline(iss, token, ',')) {
-            token.erase(std::remove(token.begin(), token.end(), '\n'),
-                        token.end());
-            token.erase(std::remove(token.begin(), token.end(), '\r'),
-                        token.end());
-            unsigned token_as_num = std::stoul(token);
-            if (cnt == 0)
-              label = token_as_num;
-            else
-              medoid = token_as_num;
-            cnt++;
-          }
-          _filter_to_medoid_id[label] = medoid;
-          line_cnt++;
-        }
-      }
-
-      std::string universal_label_file(filename);
-      universal_label_file += "_universal_label.txt";
-      if (file_exists(universal_label_file)) {
-        std::ifstream universal_label_reader(universal_label_file);
-        universal_label_reader >> _universal_label;
-        _use_universal_label = true;
-        universal_label_reader.close();
-      }
-    }
+    
 
     if (!_save_as_one_file) {
       // For DLVS Store, we will not support saving the index in multiple files.
@@ -561,6 +523,53 @@ namespace diskann {
       throw diskann::ANNException(stream.str(), -1, __FUNCSIG__, __FILE__,
                                   __LINE__);
     }
+    if (file_exists(labels_file)) {
+      label_num_pts = parse_label_file(labels_file);
+      assert(label_num_pts == data_file_num_pts);
+      if (file_exists(labels_to_medoids)) {
+        std::ifstream medoid_stream(labels_to_medoids);
+        assert(medoid_stream.is_open());        
+        std::string line, token;
+        unsigned    line_cnt = 0;
+
+        _filter_to_medoid_id.clear();
+        try {
+          while (std::getline(medoid_stream, line)) {
+            std::istringstream iss(line);
+            _u32               cnt = 0;
+            _u32               medoid = 0;
+            label              label;
+            while (std::getline(iss, token, ',')) {
+              token.erase(std::remove(token.begin(), token.end(), '\n'),
+                          token.end());
+              token.erase(std::remove(token.begin(), token.end(), '\r'),
+                          token.end());
+              unsigned token_as_num = std::stoul(token);
+              if (cnt == 0)
+                label = token_as_num;
+              else
+                medoid = token_as_num;
+              cnt++;
+            }
+            _filter_to_medoid_id[label] = medoid;
+            line_cnt++;
+          }
+        } catch (std::system_error &e) {
+          throw FileException(labels_to_medoids, e, __FUNCSIG__, __FILE__,
+                              __LINE__);
+        }
+      }
+
+      std::string universal_label_file(filename);
+      universal_label_file += "_universal_label.txt";
+      if (file_exists(universal_label_file)) {
+        std::ifstream universal_label_reader(universal_label_file);
+        assert(universal_label_reader.is_open());
+        universal_label_reader >> _universal_label;
+        _use_universal_label = true;
+        universal_label_reader.close();
+      }
+    }
 
     _nd = data_file_num_pts - _num_frozen_pts;
     _empty_slots.clear();
@@ -688,13 +697,13 @@ namespace diskann {
       }
     }
 #else
-
     size_t   bytes_read = vamana_metadata_size;
     size_t   cc = 0;
     unsigned nodes_read = 0;
     while (bytes_read != expected_file_size) {
       unsigned k;
       in.read((char *) &k, sizeof(unsigned));
+
       if (k == 0) {
         diskann::cerr << "ERROR: Point found with no out-neighbors, point#"
                       << nodes_read << std::endl;
@@ -1378,6 +1387,65 @@ namespace diskann {
     }
   }
 
+  template<typename T, typename TagT>
+  void Index<T, TagT>::prune_all_nbrs(const Parameters &parameters) {
+    const unsigned range = parameters.Get<unsigned>("R");
+    const unsigned maxc = parameters.Get<unsigned>("C");
+    const float alpha = parameters.Get<unsigned>("alpha");
+    _filtered_index = true;
+
+    diskann::Timer timer;
+#pragma omp parallel for
+    for (_s64 node = 0; node < (_s64)(_max_points + _num_frozen_pts); node++) {
+      if ((size_t) node < _nd || (size_t) node == _max_points) {
+        if (_final_graph[node].size() > range) {
+          tsl::robin_set<unsigned> dummy_visited(0);
+          std::vector<Neighbor>    dummy_pool(0);
+          std::vector<unsigned>    new_out_neighbors;
+
+					ScratchStoreManager<InMemQueryScratch<T>> manager(_query_scratch);
+      		auto scratch = manager.scratch_space();
+
+          for (auto cur_nbr : _final_graph[node]) {
+            if (dummy_visited.find(cur_nbr) == dummy_visited.end() &&
+                cur_nbr != node) {
+              float dist =
+                  _distance->compare(_data + _aligned_dim * (size_t) node,
+                                     _data + _aligned_dim * (size_t) cur_nbr,
+                                     (unsigned) _aligned_dim);
+              dummy_pool.emplace_back(Neighbor(cur_nbr, dist));
+              dummy_visited.insert(cur_nbr);
+            }
+          }
+          
+          prune_neighbors((_u32) node, dummy_pool, range, maxc, alpha, new_out_neighbors, scratch);
+          _final_graph[node].clear();
+          for (auto id : new_out_neighbors)
+            _final_graph[node].emplace_back(id);
+        }
+      }
+    }
+
+    diskann::cout << "Prune time : " << timer.elapsed() / 1000 << "ms"
+                  << std::endl;
+    size_t max = 0, min = 1 << 30, total = 0, cnt = 0;
+    for (size_t i = 0; i < (_nd + _num_frozen_pts); i++) {
+      std::vector<unsigned> pool = _final_graph[i];
+      max = (std::max)(max, pool.size());
+      min = (std::min)(min, pool.size());
+      total += pool.size();
+      if (pool.size() < 2)
+        cnt++;
+    }
+    if (min > max)
+      min = max;
+    if (_nd > 0) {
+      diskann::cout << "Index built with degree: max:" << max << "  avg:"
+                    << (float) total / (float) (_nd + _num_frozen_pts)
+                    << "  min:" << min << "  count(deg<2):" << cnt << std::endl;
+    }
+  }
+
   template<typename T, typename TagT>
   void Index<T, TagT>::set_start_point(T *data) {
     std::unique_lock<std::shared_timed_mutex> ul(_update_lock);
@@ -1648,13 +1716,13 @@ namespace diskann {
   }
 
   template<typename T, typename TagT>
-  void Index<T, TagT>::parse_label_file(const std::string &map_file) {
+  size_t Index<T, TagT>::parse_label_file(const std::string &label_file) {
     // Format of Label txt file: filters with comma separators
-
-    std::ifstream infile(map_file);
+    std::ifstream infile(label_file);
+    assert(infile.is_open());
     std::string   line, token;
     unsigned      line_cnt = 0;
-
+    
     while (std::getline(infile, line)) {
       line_cnt++;
     }
@@ -1663,28 +1731,35 @@ namespace diskann {
     infile.clear();
     infile.seekg(0, std::ios::beg);
     line_cnt = 0;
-    while (std::getline(infile, line)) {
-      std::istringstream iss(line);
-      std::vector<label> lbls(0);
-      getline(iss, token, '\t');
-      std::istringstream new_iss(token);
-      while (getline(new_iss, token, ',')) {
-        token.erase(std::remove(token.begin(), token.end(), '\n'), token.end());
-        token.erase(std::remove(token.begin(), token.end(), '\r'), token.end());
-        unsigned token_as_num = std::stoul(token);
-        lbls.push_back(token_as_num);
-        _labels.insert(token_as_num);
-      }
-      if (lbls.size() <= 0) {
-        std::cout << "No label found";
-        exit(-1);
+    try {
+      while (std::getline(infile, line)) {
+        std::istringstream iss(line);
+        std::vector<label> lbls(0);
+        getline(iss, token, '\t');
+        std::istringstream new_iss(token);
+        while (getline(new_iss, token, ',')) {
+          token.erase(std::remove(token.begin(), token.end(), '\n'),
+                      token.end());
+          token.erase(std::remove(token.begin(), token.end(), '\r'),
+                      token.end());
+          unsigned token_as_num = std::stoul(token);
+          lbls.push_back(token_as_num);
+          _labels.insert(token_as_num);
+        }
+        if (lbls.size() <= 0) {
+          diskann::cout << "No label found";
+          exit(-1);
+        }
+        std::sort(lbls.begin(), lbls.end());
+        _pts_to_labels[line_cnt] = lbls;
+        line_cnt++;
       }
-      std::sort(lbls.begin(), lbls.end());
-      _pts_to_labels[line_cnt] = lbls;
-      line_cnt++;
+      diskann::cout << "Identified " << _labels.size() << " distinct label(s)"
+                    << std::endl;
+    } catch (std::system_error &e) {
+      throw FileException(label_file, e, __FUNCSIG__, __FILE__, __LINE__);
     }
-    std::cout << "Identified " << _labels.size() << " distinct label(s)"
-              << std::endl;
+    return line_cnt;
   }
 
   template<typename T, typename TagT>
@@ -1702,9 +1777,9 @@ namespace diskann {
     _labels_file = label_file;
     _filtered_index = true;
     _filter_to_medoid_id.clear();
-    parse_label_file(label_file);  // determines medoid for each label and
+    size_t label_num_pts = parse_label_file(label_file);  // determines medoid for each label and
                                    // identifies the points to label mapping
-
+    assert(label_num_pts != num_points_to_load);
     _u32 counter = 0;
 #pragma omp parallel for schedule(dynamic, 1)
     for (int lbl = 0; lbl < _labels.size(); lbl++) {
@@ -1744,15 +1819,13 @@ namespace diskann {
           _filter_to_medoid_id[x] = best_medoid;
           _medoid_counts[best_medoid]++;
           std::stringstream a;
-          // a << "Medoid of " << x << " is " << best_medoid <<
-          // std::endl; std::cout << a.str();
         }
       }
 #pragma omp critical
       counter++;
       std::stringstream a;
       a << ((100.0 * counter) / _labels.size()) << "\% processed \r";
-      std::cout << a.str() << std::flush;
+      diskann::cout << a.str() << std::flush;
     }
 
     this->build(filename, num_points_to_load, parameters, tags);
diff --git a/src/pq_flash_index.cpp b/src/pq_flash_index.cpp
index fe326a10e..a86c14869 100644
--- a/src/pq_flash_index.cpp
+++ b/src/pq_flash_index.cpp
@@ -464,7 +464,7 @@ namespace diskann {
       num_pts++;
     }
 
-    std::cout << "Labels file metadata: num_points: " << num_pts
+    diskann::cout << "Labels file metadata: num_points: " << num_pts
               << ", #total_labels: " << num_total_labels << std::endl;
     infile.close();
   }
@@ -484,54 +484,64 @@ namespace diskann {
     return ret_val;
   }
   template<typename T>
-  void PQFlashIndex<T>::parse_label_file(const std::string &map_file) {
+  void PQFlashIndex<T>::parse_label_file(const std::string &label_file) {
     //_filtered_ann = 1;
 
-    std::ifstream infile(map_file);
+    std::ifstream infile(label_file);
+    
     std::string   line, token;
     _u32          line_cnt = 0;
 
     _u32 num_pts_in_label_file;
     _u32 num_total_labels;
-    get_label_file_metadata(map_file, num_pts_in_label_file, num_total_labels);
-
+    get_label_file_metadata(label_file, num_pts_in_label_file, num_total_labels);
+    assert(num_pts_in_label_file == this->num_points);
+    
     _pts_to_label_offsets = new _u32[num_pts_in_label_file];
     _pts_to_labels = new _u32[num_pts_in_label_file + num_total_labels];
     _u32 counter = 0;
-
-    while (std::getline(infile, line)) {
-      std::istringstream iss(line);
-      std::vector<_u32>  lbls(0);
-
-      _pts_to_label_offsets[line_cnt] = counter;
-      _u32 &num_lbls_in_cur_pt = _pts_to_labels[counter];
-      num_lbls_in_cur_pt = 0;
-      counter++;
-      getline(iss, token, '\t');
-      std::istringstream new_iss(token);
-      while (getline(new_iss, token, ',')) {
-        token.erase(std::remove(token.begin(), token.end(), '\n'), token.end());
-        token.erase(std::remove(token.begin(), token.end(), '\r'), token.end());
-        if (_labels.find(token) == _labels.end()) {
-          _filter_list.emplace_back(token);
+    infile.clear();
+    infile.seekg(0, std::ios::beg);
+    assert(infile.is_open());
+    try {
+      while (std::getline(infile, line)) {
+        std::istringstream iss(line);
+        std::vector<_u32>  lbls(0);
+
+        _pts_to_label_offsets[line_cnt] = counter;
+        _u32 &num_lbls_in_cur_pt = _pts_to_labels[counter];
+        num_lbls_in_cur_pt = 0;
+        counter++;
+        getline(iss, token, '\t');
+        std::istringstream new_iss(token);
+        while (getline(new_iss, token, ',')) {
+          token.erase(std::remove(token.begin(), token.end(), '\n'),
+                      token.end());
+          token.erase(std::remove(token.begin(), token.end(), '\r'),
+                      token.end());
+          if (_labels.find(token) == _labels.end()) {
+            _filter_list.emplace_back(token);
+          }
+          int32_t filter_num = get_filter_number(token);
+          if (filter_num == -1) {
+            diskann::cout << "Error!! " << std::endl;
+            exit(-1);
+          }
+          _pts_to_labels[counter++] = filter_num;
+          num_lbls_in_cur_pt++;
+          _labels.insert(token);
         }
-        int32_t filter_num = get_filter_number(token);
-        if (filter_num == -1) {
-          std::cout << "Error!! " << std::endl;
+
+        if (num_lbls_in_cur_pt == 0) {
+          diskann::cout << "No label found for point " << line_cnt << std::endl;
           exit(-1);
         }
-        _pts_to_labels[counter++] = filter_num;
-        num_lbls_in_cur_pt++;
-        _labels.insert(token);
+        line_cnt++;
       }
-
-      if (num_lbls_in_cur_pt == 0) {
-        std::cout << "No label found for point " << line_cnt << std::endl;
-        exit(-1);
-      }
-      line_cnt++;
+    }catch (std::system_error &e) {
+      throw FileException(label_file, e, __FUNCSIG__, __FILE__,
+                          __LINE__);
     }
-    // MallocExtension::instance()->ReleaseFreeMemory();
     infile.close();
   }
 
@@ -539,7 +549,7 @@ namespace diskann {
   void PQFlashIndex<T>::set_universal_label(const std::string &label) {
     int32_t temp_filter_num = get_filter_number(label);
     if (temp_filter_num == -1) {
-      std::cout << "Error, could not find universal label. Exitting."
+      diskann::cout << "Error, could not find universal label. Exitting."
                 << std::endl;
       exit(-1);
     } else {
@@ -595,29 +605,70 @@ namespace diskann {
         std ::string(disk_index_file) + "_labels_to_medoids.txt";
     std::string dummy_map_file =
         std ::string(disk_index_file) + "_dummy_map.txt";
+
+    size_t pq_file_dim, pq_file_num_centroids;
+#ifdef EXEC_ENV_OLS
+    get_bin_metadata(files, pq_table_bin, pq_file_num_centroids, pq_file_dim,
+                     METADATA_SIZE);
+#else
+    get_bin_metadata(pq_table_bin, pq_file_num_centroids, pq_file_dim,
+                     METADATA_SIZE);
+#endif
+
+    this->disk_index_file = disk_index_file;
+
+    if (pq_file_num_centroids != 256) {
+      diskann::cout << "Error. Number of PQ centroids is not 256. Exiting."
+                    << std::endl;
+      return -1;
+    }
+
+    this->data_dim = pq_file_dim;
+    // will reset later if we use PQ on disk
+    this->disk_data_dim = this->data_dim;
+    // will change later if we use PQ on disk or if we are using
+    // inner product without PQ
+    this->disk_bytes_per_point = this->data_dim * sizeof(T);
+    this->aligned_dim = ROUND_UP(pq_file_dim, 8);
+
+    size_t npts_u64, nchunks_u64;
+#ifdef EXEC_ENV_OLS
+    diskann::load_bin<_u8>(files, pq_compressed_vectors, this->data, npts_u64,
+                           nchunks_u64);
+#else
+    diskann::load_bin<_u8>(pq_compressed_vectors, this->data, npts_u64,
+                           nchunks_u64);
+#endif
+
+    this->num_points = npts_u64;
+    this->n_chunks = nchunks_u64;
+
     if (file_exists(labels_file)) {
       parse_label_file(labels_file);
       if (file_exists(labels_to_medoids)) {
         std::ifstream medoid_stream(labels_to_medoids);
-
+        assert(medoid_stream.is_open());
         std::string line, token;
-        //        unsigned    line_cnt = 0;
-
         _filter_to_medoid_id.clear();
 
-        while (std::getline(medoid_stream, line)) {
-          std::istringstream iss(line);
-          _u32               cnt = 0;
-          _u32               medoid;
-          std::string        label;
-          while (std::getline(iss, token, ',')) {
-            if (cnt == 0)
-              label = token;
-            else
-              medoid = (_u32) stoul(token);
-            cnt++;
+        try {
+          while (std::getline(medoid_stream, line)) {
+            std::istringstream iss(line);
+            _u32               cnt = 0;
+            _u32               medoid;
+            std::string        label;
+            while (std::getline(iss, token, ',')) {
+              if (cnt == 0)
+                label = token;
+              else
+                medoid = (_u32) stoul(token);
+              cnt++;
+            }
+            _filter_to_medoid_id[label] = medoid;
           }
-          _filter_to_medoid_id[label] = medoid;
+        } catch (std::system_error &e) {
+          throw FileException(labels_to_medoids, e, __FUNCSIG__, __FILE__,
+                              __LINE__);
         }
       }
       std::string univ_label_file =
@@ -633,8 +684,6 @@ namespace diskann {
         std::ifstream dummy_map_stream(dummy_map_file);
 
         std::string line, token;
-        //        unsigned    line_cnt = 0;
-
         while (std::getline(dummy_map_stream, line)) {
           std::istringstream iss(line);
           _u32               cnt = 0;
@@ -657,47 +706,10 @@ namespace diskann {
           _real_to_dummy_map[real_id].emplace_back(dummy_id);
         }
         dummy_map_stream.close();
-        std::cout << "Loaded dummy map" << std::endl;
+        diskann::cout << "Loaded dummy map" << std::endl;
       }
     }
 
-    size_t pq_file_dim, pq_file_num_centroids;
-#ifdef EXEC_ENV_OLS
-    get_bin_metadata(files, pq_table_bin, pq_file_num_centroids, pq_file_dim,
-                     METADATA_SIZE);
-#else
-    get_bin_metadata(pq_table_bin, pq_file_num_centroids, pq_file_dim,
-                     METADATA_SIZE);
-#endif
-
-    this->disk_index_file = disk_index_file;
-
-    if (pq_file_num_centroids != 256) {
-      diskann::cout << "Error. Number of PQ centroids is not 256. Exiting."
-                    << std::endl;
-      return -1;
-    }
-
-    this->data_dim = pq_file_dim;
-    // will reset later if we use PQ on disk
-    this->disk_data_dim = this->data_dim;
-    // will change later if we use PQ on disk or if we are using
-    // inner product without PQ
-    this->disk_bytes_per_point = this->data_dim * sizeof(T);
-    this->aligned_dim = ROUND_UP(pq_file_dim, 8);
-
-    size_t npts_u64, nchunks_u64;
-#ifdef EXEC_ENV_OLS
-    diskann::load_bin<_u8>(files, pq_compressed_vectors, this->data, npts_u64,
-                           nchunks_u64);
-#else
-    diskann::load_bin<_u8>(pq_compressed_vectors, this->data, npts_u64,
-                           nchunks_u64);
-#endif
-
-    this->num_points = npts_u64;
-    this->n_chunks = nchunks_u64;
-
 #ifdef EXEC_ENV_OLS
     pq_table.load_pq_centroid_bin(files, pq_table_bin.c_str(), nchunks_u64);
 #else
@@ -735,7 +747,7 @@ namespace diskann {
       disk_bytes_per_point =
           disk_pq_n_chunks *
           sizeof(_u8);  // revising disk_bytes_per_point since DISK PQ is used.
-      std::cout << "Disk index uses PQ data compressed down to "
+      diskann::cout << "Disk index uses PQ data compressed down to "
                 << disk_pq_n_chunks << " bytes per point." << std::endl;
     }
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index c34be5690..6aa5532ef 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -6,6 +6,9 @@ set(CMAKE_CXX_STANDARD 14)
 add_executable(build_memory_index build_memory_index.cpp)
 target_link_libraries(build_memory_index ${PROJECT_NAME} ${DISKANN_TOOLS_TCMALLOC_LINK_OPTIONS} Boost::program_options)
 
+add_executable(build_stitched_index build_stitched_index.cpp)
+target_link_libraries(build_stitched_index ${PROJECT_NAME} ${DISKANN_TOOLS_TCMALLOC_LINK_OPTIONS} Boost::program_options)
+
 add_executable(search_memory_index search_memory_index.cpp)
 target_link_libraries(search_memory_index ${PROJECT_NAME} ${DISKANN_ASYNC_LIB} ${DISKANN_TOOLS_TCMALLOC_LINK_OPTIONS} Boost::program_options)
 
diff --git a/tests/build_stitched_index.cpp b/tests/build_stitched_index.cpp
new file mode 100644
index 000000000..a714a51f5
--- /dev/null
+++ b/tests/build_stitched_index.cpp
@@ -0,0 +1,857 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT license.
+
+#include <boost/program_options.hpp>
+#include <chrono>
+#include <cstdio>
+#include <cstring>
+#include <random>
+#include <string>
+#include <tuple>
+
+#include <omp.h>
+#ifndef _WINDOWS
+#include <sys/uio.h>
+#endif
+
+#include "index.h"
+#include "memory_mapper.h"
+#include "parameters.h"
+#include "utils.h"
+
+namespace po = boost::program_options;
+
+// macros
+#define PBSTR "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||"
+#define PBWIDTH 60
+
+// custom types (for readability)
+typedef unsigned              label;
+typedef tsl::robin_set<label> label_set;
+typedef std::string           path;
+
+// structs for returning multiple items from a function
+typedef std::tuple<std::vector<label_set>, tsl::robin_map<label, _u32>,
+                   label_set>
+    parse_label_file_return_values;
+typedef std::tuple<std::vector<std::vector<_u32>>, _u64>
+    load_label_index_return_values;
+typedef std::tuple<std::vector<std::vector<_u32>>, _u64>
+    stitch_indices_return_values;
+
+/*
+ * Inline function to display progress bar.
+ */
+inline void print_progress(double percentage) {
+  int val = (int) (percentage * 100);
+  int lpad = (int) (percentage * PBWIDTH);
+  int rpad = PBWIDTH - lpad;
+  printf("\r%3d%% [%.*s%*s]", val, lpad, PBSTR, rpad, "");
+  fflush(stdout);
+}
+
+/*
+ * Inline function to generate a random integer in a range.
+ */
+inline size_t random(size_t range_from, size_t range_to) {
+  std::random_device                    rand_dev;
+  std::mt19937                          generator(rand_dev());
+  std::uniform_int_distribution<size_t> distr(range_from, range_to);
+  return distr(generator);
+}
+
+/*
+ * function to handle command line parsing.
+ *
+ * Arguments are merely the inputs from the command line.
+ */
+size_t handle_args(int argc, char **argv, std::string &data_type,
+                   path &input_data_path, path &final_index_path_prefix,
+                   path &label_data_path, label &universal_label,
+                   unsigned &num_threads, unsigned &R, unsigned &L,
+                   unsigned &stitched_R, float &alpha) {
+  po::options_description desc{"Arguments"};
+  try {
+    desc.add_options()("help,h", "Print information on arguments");
+    desc.add_options()("data_type",
+                       po::value<std::string>(&data_type)->required(),
+                       "data type <int8/uint8/float>");
+    desc.add_options()("data_path",
+                       po::value<path>(&input_data_path)->required(),
+                       "Input data file in bin format");
+    desc.add_options()("index_path_prefix",
+                       po::value<path>(&final_index_path_prefix)->required(),
+                       "Path prefix for saving index file components");
+    desc.add_options()("max_degree,R",
+                       po::value<uint32_t>(&R)->default_value(64),
+                       "Maximum graph degree");
+    desc.add_options()(
+        "Lbuild,L", po::value<uint32_t>(&L)->default_value(100),
+        "Build complexity, higher value results in better graphs");
+    desc.add_options()("stitched_R",
+                       po::value<uint32_t>(&stitched_R)->default_value(100),
+                       "Degree to prune final graph down to");
+    desc.add_options()(
+        "alpha", po::value<float>(&alpha)->default_value(1.2f),
+        "alpha controls density and diameter of graph, set 1 for sparse graph, "
+        "1.2 or 1.4 for denser graphs with lower diameter");
+    desc.add_options()(
+        "num_threads,T",
+        po::value<uint32_t>(&num_threads)->default_value(omp_get_num_procs()),
+        "Number of threads used for building index (defaults to "
+        "omp_get_num_procs())");
+    desc.add_options()("label_file",
+                       po::value<path>(&label_data_path)->default_value(""),
+                       "Input label file in txt format if present");
+    desc.add_options()(
+        "universal_label",
+        po::value<label>(&universal_label)->default_value(UINT32_MAX),
+        "If a point comes with the specified universal label (and only the "
+        "univ. "
+        "label), then the point is considered to have every possible label");
+
+    po::variables_map vm;
+    po::store(po::parse_command_line(argc, argv, desc), vm);
+    if (vm.count("help")) {
+      std::cout << desc;
+      exit(0);
+    }
+    po::notify(vm);
+  } catch (const std::exception &ex) {
+    std::cerr << ex.what() << '\n';
+    throw;
+  }
+  return 1;
+}
+
+/*
+ * Parses the label datafile, which has comma-separated labels on
+ * each line. Line i corresponds to point id i.
+ *
+ * Returns three objects via std::tuple:
+ * 1. map: key is point id, value is vector of labels said point has
+ * 2. map: key is label, value is number of points with the label
+ * 3. the label universe as a set
+ */
+parse_label_file_return_values parse_label_file(path  label_data_path,
+                                                label universal_label) {
+  std::ifstream label_data_stream(label_data_path);
+  std::string   line, token;
+  unsigned      line_cnt = 0;
+
+  // allows us to reserve space for the points_to_labels vector
+  while (std::getline(label_data_stream, line))
+    line_cnt++;
+  label_data_stream.clear();
+  label_data_stream.seekg(0, std::ios::beg);
+
+  // values to return
+  std::vector<label_set>      point_ids_to_labels(line_cnt);
+  tsl::robin_map<label, _u32> labels_to_number_of_points;
+  label_set                   all_labels;
+
+  std::vector<_u32> points_with_universal_label;
+  line_cnt = 0;
+  while (std::getline(label_data_stream, line)) {
+    std::istringstream current_labels_comma_separated(line);
+    label_set          current_labels;
+
+    // get point id
+    _u32 point_id = line_cnt;
+
+    // parse comma separated labels
+    bool current_universal_label_check = false;
+    while (getline(current_labels_comma_separated, token, ',')) {
+      token.erase(std::remove(token.begin(), token.end(), '\n'), token.end());
+      token.erase(std::remove(token.begin(), token.end(), '\r'), token.end());
+
+      // if token is empty, there's no labels for the point
+      unsigned token_as_num = std::stoul(token);
+      if (token_as_num == universal_label) {
+        points_with_universal_label.push_back(point_id);
+        current_universal_label_check = true;
+      } else {
+        all_labels.insert(token_as_num);
+        current_labels.insert(token_as_num);
+        labels_to_number_of_points[token_as_num]++;
+      }
+    }
+
+    if (current_labels.size() <= 0 && !current_universal_label_check) {
+      std::cerr << "Error: " << point_id << " has no labels." << std::endl;
+      exit(-1);
+    }
+    point_ids_to_labels[point_id] = current_labels;
+    line_cnt++;
+  }
+
+  // for every point with universal label, set its label set to all labels
+  // also, increment the count for number of points a label has
+  for (const auto &point_id : points_with_universal_label) {
+    point_ids_to_labels[point_id] = all_labels;
+    for (const auto &label : all_labels)
+      labels_to_number_of_points[label]++;
+  }
+
+  std::cout << "Identified " << all_labels.size() << " distinct label(s) for "
+            << point_ids_to_labels.size() << " points\n"
+            << std::endl;
+
+  return std::make_tuple(point_ids_to_labels, labels_to_number_of_points,
+                         all_labels);
+}
+
+/*
+ * For each label, generates a file containing all vectors that have said label.
+ * Also copies data from original bin file to new dimension-aligned file.
+ *
+ * Utilizes POSIX functions mmap and writev in order to minimize memory
+ * overhead, so we include an STL version as well.
+ *
+ * Each data file is saved under the following format:
+ *    input_data_path + "_" + label
+ */
+template<typename T>
+tsl::robin_map<label, std::vector<_u32>> generate_label_specific_vector_files(
+    path                        input_data_path,
+    tsl::robin_map<label, _u32> labels_to_number_of_points,
+    std::vector<label_set> point_ids_to_labels, label_set all_labels) {
+  auto file_writing_timer = std::chrono::high_resolution_clock::now();
+  diskann::MemoryMapper input_data(input_data_path);
+  char *                input_start = input_data.getBuf();
+
+  _u32 number_of_points, dimension;
+  std::memcpy(&number_of_points, input_start, sizeof(_u32));
+  std::memcpy(&dimension, input_start + sizeof(_u32), sizeof(_u32));
+  const _u32   VECTOR_SIZE = dimension * sizeof(T);
+  const size_t METADATA = 2 * sizeof(_u32);
+  if (number_of_points != point_ids_to_labels.size()) {
+    std::cerr << "Error: number of points in labels file and data file differ."
+              << std::endl;
+    throw;
+  }
+
+  tsl::robin_map<label, iovec *>           label_to_iovec_map;
+  tsl::robin_map<label, _u32>              label_to_curr_iovec;
+  tsl::robin_map<label, std::vector<_u32>> label_id_to_orig_id;
+
+  // setup iovec list for each label
+  for (const auto &label : all_labels) {
+    iovec *label_iovecs =
+        (iovec *) malloc(labels_to_number_of_points[label] * sizeof(iovec));
+    if (label_iovecs == nullptr) {
+      throw;
+    }
+    label_to_iovec_map[label] = label_iovecs;
+    label_to_curr_iovec[label] = 0;
+    label_id_to_orig_id[label].reserve(labels_to_number_of_points[label]);
+  }
+
+  // each point added to corresponding per-label iovec list
+  for (_u32 point_id = 0; point_id < number_of_points; point_id++) {
+    char *curr_point = input_start + METADATA + (VECTOR_SIZE * point_id);
+    iovec curr_iovec;
+
+    curr_iovec.iov_base = curr_point;
+    curr_iovec.iov_len = VECTOR_SIZE;
+    for (const auto &label : point_ids_to_labels[point_id]) {
+      *(label_to_iovec_map[label] + label_to_curr_iovec[label]) = curr_iovec;
+      label_to_curr_iovec[label]++;
+      label_id_to_orig_id[label].push_back(point_id);
+    }
+  }
+
+  // write each label iovec to resp. file
+  for (const auto &label : all_labels) {
+    int  label_input_data_fd;
+    path curr_label_input_data_path(input_data_path + "_" +
+                                    std::to_string(label));
+    _u32 curr_num_pts = labels_to_number_of_points[label];
+
+    label_input_data_fd =
+        open(curr_label_input_data_path.c_str(),
+             O_CREAT | O_WRONLY | O_TRUNC | O_APPEND, (mode_t) 0644);
+    if (label_input_data_fd == -1)
+      throw;
+
+    // write metadata
+    _u32 metadata[2] = {curr_num_pts, dimension};
+    int  return_value = write(label_input_data_fd, metadata, sizeof(_u32) * 2);
+    if (return_value == -1) {
+      throw;
+    }
+
+    // limits on number of iovec structs per writev means we need to perform
+    // multiple writevs
+    size_t i = 0;
+    while (curr_num_pts > IOV_MAX) {
+      return_value =
+          writev(label_input_data_fd,
+                 (label_to_iovec_map[label] + (IOV_MAX * i)), IOV_MAX);
+      if (return_value == -1) {
+        close(label_input_data_fd);
+        throw;
+      }
+      curr_num_pts -= IOV_MAX;
+      i += 1;
+    }
+    return_value =
+        writev(label_input_data_fd, (label_to_iovec_map[label] + (IOV_MAX * i)),
+               curr_num_pts);
+    if (return_value == -1) {
+      close(label_input_data_fd);
+      throw;
+    }
+
+    free(label_to_iovec_map[label]);
+    close(label_input_data_fd);
+  }
+
+  std::chrono::duration<double> file_writing_time =
+      std::chrono::high_resolution_clock::now() - file_writing_timer;
+  std::cout << "generated " << all_labels.size()
+            << " label-specific vector files for index building in time "
+            << file_writing_time.count() << "\n"
+            << std::endl;
+
+  return label_id_to_orig_id;
+}
+
+// for use on systems without writev (i.e. Windows)
+template<typename T>
+tsl::robin_map<label, std::vector<_u32>>
+generate_label_specific_vector_files_compat(
+    path                        input_data_path,
+    tsl::robin_map<label, _u32> labels_to_number_of_points,
+    std::vector<label_set> point_ids_to_labels, label_set all_labels) {
+  auto          file_writing_timer = std::chrono::high_resolution_clock::now();
+  std::ifstream input_data_stream(input_data_path);
+
+  _u32 number_of_points, dimension;
+  input_data_stream.read((char *) &number_of_points, sizeof(_u32));
+  input_data_stream.read((char *) &dimension, sizeof(_u32));
+  const _u32 VECTOR_SIZE = dimension * sizeof(T);
+  if (number_of_points != point_ids_to_labels.size()) {
+    std::cerr << "Error: number of points in labels file and data file differ."
+              << std::endl;
+    throw;
+  }
+
+  tsl::robin_map<label, char *>            labels_to_vectors;
+  tsl::robin_map<label, _u32>              labels_to_curr_vector;
+  tsl::robin_map<label, std::vector<_u32>> label_id_to_orig_id;
+
+  for (const auto &label : all_labels) {
+    _u32  number_of_label_pts = labels_to_number_of_points[label];
+    char *vectors = (char *) malloc(number_of_label_pts * VECTOR_SIZE);
+    if (vectors == nullptr) {
+      throw;
+    }
+    labels_to_vectors[label] = vectors;
+    labels_to_curr_vector[label] = 0;
+    label_id_to_orig_id[label].reserve(number_of_label_pts);
+  }
+
+  for (_u32 point_id = 0; point_id < number_of_points; point_id++) {
+    char *curr_vector = (char *) malloc(VECTOR_SIZE);
+    input_data_stream.read(curr_vector, VECTOR_SIZE);
+    for (const auto &label : point_ids_to_labels[point_id]) {
+      char *curr_label_vector_ptr =
+          labels_to_vectors[label] +
+          (labels_to_curr_vector[label] * VECTOR_SIZE);
+      memcpy(curr_label_vector_ptr, curr_vector, VECTOR_SIZE);
+      labels_to_curr_vector[label]++;
+      label_id_to_orig_id[label].push_back(point_id);
+    }
+    free(curr_vector);
+  }
+
+  for (const auto &label : all_labels) {
+    path curr_label_input_data_path(input_data_path + "_" +
+                                    std::to_string(label));
+    _u32 number_of_label_pts = labels_to_number_of_points[label];
+
+    std::ofstream label_file_stream;
+    label_file_stream.exceptions(std::ios::badbit | std::ios::failbit);
+    label_file_stream.open(curr_label_input_data_path);
+
+    label_file_stream.write((char *) &number_of_label_pts, sizeof(_u32));
+    label_file_stream.write((char *) &dimension, sizeof(_u32));
+    label_file_stream.write((char *) labels_to_vectors[label],
+                            number_of_label_pts * VECTOR_SIZE);
+
+    free(labels_to_vectors[label]);
+    label_file_stream.close();
+  }
+
+  std::chrono::duration<double> file_writing_time =
+      std::chrono::high_resolution_clock::now() - file_writing_timer;
+  std::cout << "generated " << all_labels.size()
+            << " label-specific vector files for index building in time "
+            << file_writing_time.count() << "\n"
+            << std::endl;
+
+  return label_id_to_orig_id;
+}
+
+/*
+ * Using passed in parameters and files generated from step 3,
+ * builds a vanilla diskANN index for each label.
+ *
+ * Each index is saved under the following path:
+ *  final_index_path_prefix + "_" + label
+ */
+template<typename T>
+void generate_label_indices(path input_data_path, path final_index_path_prefix,
+                            label_set all_labels, unsigned R, unsigned L,
+                            float alpha, unsigned num_threads) {
+  diskann::Parameters label_index_build_parameters;
+  label_index_build_parameters.Set<unsigned>("R", R);
+  label_index_build_parameters.Set<unsigned>("L", L);
+  label_index_build_parameters.Set<unsigned>("C", 750);
+	label_index_build_parameters.Set<unsigned>("Lf", 0);
+  label_index_build_parameters.Set<bool>("saturate_graph", 0);
+  label_index_build_parameters.Set<float>("alpha", alpha);
+  label_index_build_parameters.Set<unsigned>("num_threads", num_threads);
+
+  std::cout << "Generating indices per label..." << std::endl;
+  // for each label, build an index on resp. points
+  double total_indexing_time = 0.0, indexing_percentage = 0.0;
+  std::cout.setstate(std::ios_base::failbit);
+  diskann::cout.setstate(std::ios_base::failbit);
+  for (const auto &label : all_labels) {
+    path curr_label_input_data_path(input_data_path + "_" +
+                                    std::to_string(label));
+    path curr_label_index_path(final_index_path_prefix + "_" +
+                               std::to_string(label));
+
+    size_t number_of_label_points, dimension;
+    diskann::get_bin_metadata(curr_label_input_data_path,
+                              number_of_label_points, dimension);
+    diskann::Index<T> index(diskann::Metric::L2, dimension,
+                            number_of_label_points, false, false);
+
+    auto index_build_timer = std::chrono::high_resolution_clock::now();
+    index.build(curr_label_input_data_path.c_str(), number_of_label_points,
+                label_index_build_parameters);
+    std::chrono::duration<double> current_indexing_time =
+        std::chrono::high_resolution_clock::now() - index_build_timer;
+
+    total_indexing_time += current_indexing_time.count();
+    indexing_percentage += (1 / (double) all_labels.size());
+    print_progress(indexing_percentage);
+
+    index.save(curr_label_index_path.c_str());
+  }
+  std::cout.clear();
+  diskann::cout.clear();
+
+  std::cout << "\nDone. Generated per-label indices in " << total_indexing_time
+            << " seconds\n"
+            << std::endl;
+}
+
+/*
+ * Manually loads a graph index in from a given file.
+ *
+ * Returns both the graph index and the size of the file in bytes.
+ */
+load_label_index_return_values load_label_index(path label_index_path,
+                                                _u32 label_number_of_points) {
+  std::ifstream label_index_stream;
+  label_index_stream.exceptions(std::ios::badbit | std::ios::failbit);
+  label_index_stream.open(label_index_path, std::ios::binary);
+
+  _u64         index_file_size, index_num_frozen_points;
+  _u32         index_max_observed_degree, index_entry_point;
+  const size_t INDEX_METADATA = 2 * sizeof(_u64) + 2 * sizeof(_u32);
+  label_index_stream.read((char *) &index_file_size, sizeof(_u64));
+  label_index_stream.read((char *) &index_max_observed_degree, sizeof(_u32));
+  label_index_stream.read((char *) &index_entry_point, sizeof(_u32));
+  label_index_stream.read((char *) &index_num_frozen_points, sizeof(_u64));
+  size_t bytes_read = INDEX_METADATA;
+
+  std::vector<std::vector<_u32>> label_index(label_number_of_points);
+  _u32                           nodes_read = 0;
+  while (bytes_read != index_file_size) {
+    _u32 current_node_num_neighbors;
+    label_index_stream.read((char *) &current_node_num_neighbors, sizeof(_u32));
+    nodes_read++;
+
+    std::vector<_u32> current_node_neighbors(current_node_num_neighbors);
+    label_index_stream.read((char *) current_node_neighbors.data(),
+                            current_node_num_neighbors * sizeof(_u32));
+    label_index[nodes_read - 1].swap(current_node_neighbors);
+    bytes_read += sizeof(_u32) * (current_node_num_neighbors + 1);
+  }
+
+  return std::make_tuple(label_index, index_file_size);
+}
+
+/*
+ * Custom index save to write the in-memory index to disk.
+ * Also writes required files for diskANN API -
+ *  1. labels_to_medoids
+ *  2. universal_label
+ *  3. data (redundant for static indices)
+ *  4. labels (redundant for static indices)
+ */
+void save_full_index(path final_index_path_prefix, path input_data_path,
+                     _u64                           final_index_size,
+                     std::vector<std::vector<_u32>> stitched_graph,
+                     tsl::robin_map<label, _u32>    entry_points,
+                     label universal_label, path label_data_path) {
+  // aux. file 1
+  auto          saving_index_timer = std::chrono::high_resolution_clock::now();
+  std::ifstream original_label_data_stream;
+  original_label_data_stream.exceptions(std::ios::badbit | std::ios::failbit);
+  original_label_data_stream.open(label_data_path, std::ios::binary);
+  std::ofstream new_label_data_stream;
+  new_label_data_stream.exceptions(std::ios::badbit | std::ios::failbit);
+  new_label_data_stream.open(final_index_path_prefix + "_labels.txt",
+                             std::ios::binary);
+  new_label_data_stream << original_label_data_stream.rdbuf();
+  original_label_data_stream.close();
+  new_label_data_stream.close();
+
+  // aux. file 2
+  std::ifstream original_input_data_stream;
+  original_input_data_stream.exceptions(std::ios::badbit | std::ios::failbit);
+  original_input_data_stream.open(input_data_path, std::ios::binary);
+  std::ofstream new_input_data_stream;
+  new_input_data_stream.exceptions(std::ios::badbit | std::ios::failbit);
+  new_input_data_stream.open(final_index_path_prefix + ".data",
+                             std::ios::binary);
+  new_input_data_stream << original_input_data_stream.rdbuf();
+  original_input_data_stream.close();
+  new_input_data_stream.close();
+
+  // aux. file 3
+  std::ofstream labels_to_medoids_writer;
+  labels_to_medoids_writer.exceptions(std::ios::badbit | std::ios::failbit);
+  labels_to_medoids_writer.open(final_index_path_prefix +
+                                "_labels_to_medoids.txt");
+  for (auto iter : entry_points)
+    labels_to_medoids_writer << iter.first << ", " << iter.second << std::endl;
+  labels_to_medoids_writer.close();
+
+  // aux. file 4 (only if we're using a universal label)
+  if (universal_label != UINT32_MAX) {
+    std::ofstream universal_label_writer;
+    universal_label_writer.exceptions(std::ios::badbit | std::ios::failbit);
+    universal_label_writer.open(final_index_path_prefix +
+                                "_universal_label.txt");
+    universal_label_writer << universal_label << std::endl;
+    universal_label_writer.close();
+  }
+
+  // main index
+  _u64         index_num_frozen_points = 0, index_num_edges = 0;
+  _u32         index_max_observed_degree = 0, index_entry_point = 0;
+  const size_t METADATA = 2 * sizeof(_u64) + 2 * sizeof(_u32);
+  for (auto &point_neighbors : stitched_graph) {
+    index_max_observed_degree =
+        std::max(index_max_observed_degree, (_u32) point_neighbors.size());
+  }
+
+  std::ofstream stitched_graph_writer;
+  stitched_graph_writer.exceptions(std::ios::badbit | std::ios::failbit);
+  stitched_graph_writer.open(final_index_path_prefix);
+
+  stitched_graph_writer.write((char *) &final_index_size, sizeof(_u64));
+  stitched_graph_writer.write((char *) &index_max_observed_degree,
+                              sizeof(_u32));
+  stitched_graph_writer.write((char *) &index_entry_point, sizeof(_u32));
+  stitched_graph_writer.write((char *) &index_num_frozen_points, sizeof(_u64));
+
+  size_t bytes_written = METADATA;
+  for (_u32 node_point = 0; node_point < stitched_graph.size(); node_point++) {
+    _u32 current_node_num_neighbors = stitched_graph[node_point].size();
+    std::vector<_u32> current_node_neighbors = stitched_graph[node_point];
+    stitched_graph_writer.write((char *) &current_node_num_neighbors,
+                                sizeof(_u32));
+    bytes_written += sizeof(_u32);
+    for (const auto &current_node_neighbor : current_node_neighbors) {
+      stitched_graph_writer.write((char *) &current_node_neighbor,
+                                  sizeof(_u32));
+      bytes_written += sizeof(_u32);
+    }
+    index_num_edges += current_node_num_neighbors;
+  }
+
+  if (bytes_written != final_index_size) {
+    std::cerr << "Error: written bytes does not match allocated space"
+              << std::endl;
+    throw;
+  }
+
+  stitched_graph_writer.close();
+
+  std::chrono::duration<double> saving_index_time =
+      std::chrono::high_resolution_clock::now() - saving_index_timer;
+  std::cout << "Stitched graph written in " << saving_index_time.count()
+            << " seconds" << std::endl;
+  std::cout << "Stitched graph average degree: "
+            << ((float) index_num_edges) / ((float) (stitched_graph.size()))
+            << std::endl;
+  std::cout << "Stitched graph max degree: " << index_max_observed_degree
+            << std::endl
+            << std::endl;
+}
+
+/*
+ * Unions the per-label graph indices together via the following policy:
+ *  - any two nodes can only have at most one edge between them -
+ *
+ * Returns the "stitched" graph and its expected file size.
+ */
+template<typename T>
+stitch_indices_return_values stitch_label_indices(
+    path final_index_path_prefix, _u32 total_number_of_points,
+    label_set                                all_labels,
+    tsl::robin_map<label, _u32>              labels_to_number_of_points,
+    tsl::robin_map<label, _u32> &            label_entry_points,
+    tsl::robin_map<label, std::vector<_u32>> label_id_to_orig_id_map) {
+  size_t                         final_index_size = 0;
+  std::vector<std::vector<_u32>> stitched_graph(total_number_of_points);
+
+  auto stitching_index_timer = std::chrono::high_resolution_clock::now();
+  for (const auto &label : all_labels) {
+    path curr_label_index_path(final_index_path_prefix + "_" +
+                               std::to_string(label));
+    std::vector<std::vector<_u32>> curr_label_index;
+    _u64                           curr_label_index_size;
+    _u32                           curr_label_entry_point;
+
+    std::tie(curr_label_index, curr_label_index_size) = load_label_index(
+        curr_label_index_path, labels_to_number_of_points[label]);
+    curr_label_entry_point = random(0, curr_label_index.size());
+    label_entry_points[label] =
+        label_id_to_orig_id_map[label][curr_label_entry_point];
+
+    for (_u32 node_point = 0; node_point < curr_label_index.size();
+         node_point++) {
+      _u32 original_point_id = label_id_to_orig_id_map[label][node_point];
+      for (auto &node_neighbor : curr_label_index[node_point]) {
+        _u32 original_neighbor_id =
+            label_id_to_orig_id_map[label][node_neighbor];
+        std::vector<_u32> curr_point_neighbors =
+            stitched_graph[original_point_id];
+        if (std::find(curr_point_neighbors.begin(), curr_point_neighbors.end(),
+                      original_neighbor_id) == curr_point_neighbors.end()) {
+          stitched_graph[original_point_id].push_back(original_neighbor_id);
+          final_index_size += sizeof(_u32);
+        }
+      }
+    }
+  }
+
+  const size_t METADATA = 2 * sizeof(_u64) + 2 * sizeof(_u32);
+  final_index_size += (total_number_of_points * sizeof(_u32) + METADATA);
+
+  std::chrono::duration<double> stitching_index_time =
+      std::chrono::high_resolution_clock::now() - stitching_index_timer;
+  std::cout << "stitched graph generated in memory in "
+            << stitching_index_time.count() << " seconds" << std::endl;
+
+  return std::make_tuple(stitched_graph, final_index_size);
+}
+
+/*
+ * Applies the prune_neighbors function from src/index.cpp to
+ * every node in the stitched graph.
+ *
+ * This is an optional step, hence the saving of both the full
+ * and pruned graph.
+ */
+template<typename T>
+void prune_and_save(path final_index_path_prefix, path input_data_path,
+                    std::vector<std::vector<_u32>> stitched_graph,
+                    unsigned                       stitched_R,
+                    tsl::robin_map<label, _u32>    label_entry_points,
+                    label universal_label, path label_data_path,
+                    unsigned num_threads) {
+  size_t dimension, number_of_label_points;
+  auto   diskann_cout_buffer = diskann::cout.rdbuf(nullptr);
+  auto   std_cout_buffer = std::cout.rdbuf(nullptr);
+  auto   pruning_index_timer = std::chrono::high_resolution_clock::now();
+
+  diskann::get_bin_metadata(input_data_path, number_of_label_points, dimension);
+  diskann::Index<T> index(diskann::Metric::L2, dimension,
+                          number_of_label_points, false, false);
+
+  // not searching this index, set search_l to 0
+  index.load(final_index_path_prefix.c_str(), num_threads, 1);
+
+  diskann::Parameters paras;
+  paras.Set<unsigned>("R", stitched_R);
+  paras.Set<unsigned>(
+      "C", 750);  // maximum candidate set size during pruning procedure
+  paras.Set<float>("alpha", 1.2);
+  paras.Set<bool>("saturate_graph", 1);
+  std::cout << "parsing labels" << std::endl;
+
+  index.prune_all_nbrs(paras);
+  index.save((final_index_path_prefix + "_pruned").c_str());
+
+  diskann::cout.rdbuf(diskann_cout_buffer);
+  std::cout.rdbuf(std_cout_buffer);
+  std::chrono::duration<double> pruning_index_time =
+      std::chrono::high_resolution_clock::now() - pruning_index_timer;
+  std::cout << "pruning performed in " << pruning_index_time.count()
+            << " seconds\n"
+            << std::endl;
+}
+
+/*
+ * Delete all temporary artifacts.
+ * In the process of creating the stitched index, some temporary artifacts are
+ * created:
+ * 1. the separate bin files for each labels' points
+ * 2. the separate diskANN indices built for each label
+ * 3. the '.data' file created while generating the indices
+ */
+void clean_up_artifacts(path input_data_path, path final_index_path_prefix,
+                        label_set all_labels) {
+  for (const auto &label : all_labels) {
+    path curr_label_input_data_path(input_data_path + "_" +
+                                    std::to_string(label));
+    path curr_label_index_path(final_index_path_prefix + "_" +
+                               std::to_string(label));
+    path curr_label_index_path_data(curr_label_index_path + ".data");
+
+    if (std::remove(curr_label_index_path.c_str()) != 0)
+      throw;
+    if (std::remove(curr_label_input_data_path.c_str()) != 0)
+      throw;
+    if (std::remove(curr_label_index_path_data.c_str()) != 0)
+      throw;
+  }
+}
+
+int main(int argc, char **argv) {
+  // 1. handle cmdline inputs
+  std::string data_type;
+  path        input_data_path, final_index_path_prefix, label_data_path;
+  label       universal_label;
+  unsigned    num_threads, R, L, stitched_R;
+  float       alpha;
+
+  auto index_timer = std::chrono::high_resolution_clock::now();
+  handle_args(argc, argv, data_type, input_data_path, final_index_path_prefix,
+              label_data_path, universal_label, num_threads, R, L, stitched_R,
+              alpha);
+
+  // 2. parse label file and create necessary data structures
+  std::vector<label_set>      point_ids_to_labels;
+  tsl::robin_map<label, _u32> labels_to_number_of_points;
+  label_set                   all_labels;
+
+  std::tie(point_ids_to_labels, labels_to_number_of_points, all_labels) =
+      parse_label_file(label_data_path, universal_label);
+
+  // 3. for each label, make a separate data file
+  tsl::robin_map<label, std::vector<_u32>> label_id_to_orig_id_map;
+  _u32 total_number_of_points = point_ids_to_labels.size();
+
+#ifndef _WINDOWS
+  if (data_type == "uint8")
+    label_id_to_orig_id_map = generate_label_specific_vector_files<uint8_t>(
+        input_data_path, labels_to_number_of_points, point_ids_to_labels,
+        all_labels);
+  else if (data_type == "int8")
+    label_id_to_orig_id_map = generate_label_specific_vector_files<int8_t>(
+        input_data_path, labels_to_number_of_points, point_ids_to_labels,
+        all_labels);
+  else if (data_type == "float")
+    label_id_to_orig_id_map = generate_label_specific_vector_files<float>(
+        input_data_path, labels_to_number_of_points, point_ids_to_labels,
+        all_labels);
+  else
+    throw;
+#else
+  if (data_type == "uint8")
+    label_id_to_orig_id_map =
+        generate_label_specific_vector_files_compat<uint8_t>(
+            input_data_path, labels_to_number_of_points, point_ids_to_labels,
+            all_labels);
+  else if (data_type == "int8")
+    label_id_to_orig_id_map =
+        generate_label_specific_vector_files_compat<int8_t>(
+            input_data_path, labels_to_number_of_points, point_ids_to_labels,
+            all_labels);
+  else if (data_type == "float")
+    label_id_to_orig_id_map =
+        generate_label_specific_vector_files_compat<float>(
+            input_data_path, labels_to_number_of_points, point_ids_to_labels,
+            all_labels);
+  else
+    throw;
+#endif
+
+  // 4. for each created data file, create a vanilla diskANN index
+  if (data_type == "uint8")
+    generate_label_indices<uint8_t>(input_data_path, final_index_path_prefix,
+                                    all_labels, R, L, alpha, num_threads);
+  else if (data_type == "int8")
+    generate_label_indices<int8_t>(input_data_path, final_index_path_prefix,
+                                   all_labels, R, L, alpha, num_threads);
+  else if (data_type == "float")
+    generate_label_indices<float>(input_data_path, final_index_path_prefix,
+                                  all_labels, R, L, alpha, num_threads);
+  else
+    throw;
+
+  // 5. "stitch" the indices together
+  std::vector<std::vector<_u32>> stitched_graph;
+  tsl::robin_map<label, _u32>    label_entry_points;
+  _u64                           stitched_graph_size;
+
+  if (data_type == "uint8")
+    std::tie(stitched_graph, stitched_graph_size) =
+        stitch_label_indices<uint8_t>(
+            final_index_path_prefix, total_number_of_points, all_labels,
+            labels_to_number_of_points, label_entry_points,
+            label_id_to_orig_id_map);
+  else if (data_type == "int8")
+    std::tie(stitched_graph, stitched_graph_size) =
+        stitch_label_indices<int8_t>(
+            final_index_path_prefix, total_number_of_points, all_labels,
+            labels_to_number_of_points, label_entry_points,
+            label_id_to_orig_id_map);
+  else if (data_type == "float")
+    std::tie(stitched_graph, stitched_graph_size) = stitch_label_indices<float>(
+        final_index_path_prefix, total_number_of_points, all_labels,
+        labels_to_number_of_points, label_entry_points,
+        label_id_to_orig_id_map);
+  else
+    throw;
+
+  // 5a. save the stitched graph to disk
+  save_full_index(final_index_path_prefix, input_data_path, stitched_graph_size,
+                  stitched_graph, label_entry_points, universal_label,
+                  label_data_path);
+
+  // 6. run a prune on the stitched index, and save to disk
+  if (data_type == "uint8")
+    prune_and_save<uint8_t>(final_index_path_prefix, input_data_path,
+                            stitched_graph, stitched_R, label_entry_points,
+                            universal_label, label_data_path, num_threads);
+  else if (data_type == "int8")
+    prune_and_save<int8_t>(final_index_path_prefix, input_data_path,
+                           stitched_graph, stitched_R, label_entry_points,
+                           universal_label, label_data_path, num_threads);
+  else if (data_type == "float")
+    prune_and_save<float>(final_index_path_prefix, input_data_path,
+                          stitched_graph, stitched_R, label_entry_points,
+                          universal_label, label_data_path, num_threads);
+  else
+    throw;
+
+  std::chrono::duration<double> index_time =
+      std::chrono::high_resolution_clock::now() - index_timer;
+  std::cout << "pruned/stitched graph generated in " << index_time.count()
+            << " seconds" << std::endl;
+
+  clean_up_artifacts(input_data_path, final_index_path_prefix, all_labels);
+}
diff --git a/tests/search_disk_index.cpp b/tests/search_disk_index.cpp
index a765669ef..2825ac1cc 100644
--- a/tests/search_disk_index.cpp
+++ b/tests/search_disk_index.cpp
@@ -425,4 +425,5 @@ int main(int argc, char** argv) {
     diskann::cerr << "Index search failed." << std::endl;
     return -1;
   }
+
 }
\ No newline at end of file
diff --git a/tests/search_memory_index.cpp b/tests/search_memory_index.cpp
index c0e3c094a..153662988 100644
--- a/tests/search_memory_index.cpp
+++ b/tests/search_memory_index.cpp
@@ -357,4 +357,5 @@ int main(int argc, char** argv) {
     diskann::cerr << "Index search failed." << std::endl;
     return -1;
   }
+
 }
\ No newline at end of file
diff --git a/tests/utils/compute_groundtruth.cpp b/tests/utils/compute_groundtruth.cpp
index d57eb70c3..5de1e7351 100644
--- a/tests/utils/compute_groundtruth.cpp
+++ b/tests/utils/compute_groundtruth.cpp
@@ -327,14 +327,11 @@ inline std::vector<size_t> load_filtered_bin_as_float(
   std::cout << "start and end ids: " << start_id << ", " << end_id << std::endl;
   reader.seekg(start_id * ndims * sizeof(T) + 2 * sizeof(uint32_t),
                std::ios::beg);
-  //    data = new T[nptsuint64_t * ndimsuint64_t];
   T *data_T = new T[nptsuint64_t * ndimsuint64_t];
   reader.read((char *) data_T, sizeof(T) * nptsuint64_t * ndimsuint64_t);
   std::cout << "Finished reading part of the bin file." << std::endl;
   reader.close();
-  //  data =  (nptsuint64_t*ndimsuint64_t, ALIGNMENT);
   data = aligned_malloc<float>(nptsuint64_t * ndimsuint64_t, ALIGNMENT);
-  //#pragma omp parallel for schedule(dynamic, 32768)
   for (int64_t i = 0; i < (int64_t) nptsuint64_t; i++) {
     if (std::find(pts_to_labels[start_id + i].begin(),
                   pts_to_labels[start_id + i].end(),
diff --git a/tests/utils/stats_label_data.cpp b/tests/utils/stats_label_data.cpp
index 3287179c3..b983b02d3 100644
--- a/tests/utils/stats_label_data.cpp
+++ b/tests/utils/stats_label_data.cpp
@@ -1,8 +1,5 @@
-//
-// Created by Neel Karia on 10/01/2021.
-//
-// Generate base file and groundtruth file for points belonging to a particular
-// label
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT license.
 
 #include <iostream>
 #include <fstream>
@@ -33,7 +30,6 @@ void stats_analysis(const std::string labels_file, std::string univeral_label,
                     _u32 density = 10) {
   std::string   token, line;
   std::ifstream labels_stream(labels_file);
-  // int           cnt[51] = {0};
   std::unordered_map<std::string, _u32> label_counts;
   std::string                           label_with_max_points;
   _u32                                  max_points = 0;
@@ -68,20 +64,11 @@ void stats_analysis(const std::string labels_file, std::string univeral_label,
             << (float) dense_pts / (float) labels_per_point.size() << std::endl;
   std::sort(labels_per_point.begin(), labels_per_point.end());
 
-  /*  for (_u32 p = 90; p < 100; p += 1) {
-      std::cout << labels_per_point[(p * labels_per_point.size()) / 100.0]
-                << std::endl;
-    }
-    std::cout << labels_per_point[labels_per_point.size() - 1] << std::endl;
-  */
-
   std::vector<std::pair<std::string, _u32>> label_count_vec;
 
   for (auto it = label_counts.begin(); it != label_counts.end(); it++) {
     auto& lbl = *it;
     label_count_vec.emplace_back(std::make_pair(lbl.first, lbl.second));
-    //    std::cout << "Label " << lbl.first << " has " << lbl.second
-    //              << " data points" << std::endl;
     if (lbl.second > max_points) {
       max_points = lbl.second;
       label_with_max_points = lbl.first;