       cudaFree(d_coords_sorted[i_gpu]);
       cudaFree(d_clustering[i_gpu]);
     }
-
-//TODO
-    std::cout << "#  " << clustering_sorted[0] << std::endl;
-
-//    return clustering;
     return normalized_cluster_names(first_frame_above_threshold
                                   , clustering
                                   , fe_sorted);

    }
    __syncthreads();
    if (gid < i_to) {
      unsigned int cluster_id = clustering[gid];
      if (cluster_id == 0) {
        cluster_id = gid + 1;
      }
      unsigned int ref_id = tid+bsize;
      // load reference coordinates for re-use into shared memory

        if (i_gpu == 0) {
          clustering_sorted[i] = tmp_clust[i];
        } else {
          clustering_sorted[i] = std::max(clustering_sorted[i]
                                        , tmp_clust[i]);
        }
      }
    }
    // reduce clustering to min number of ids
    clustering_sorted = sanitize_state_names(clustering_sorted);
    //TODO: refined clustering (by rerunning with alternative kernel)




    // convert state trajectory from
    // FE-sorted order to original order
    std::vector<std::size_t> clustering(n_rows, 0);
    for (unsigned int i=0; i < n_rows; ++i) {
      clustering[fe_sorted[i].first] = clustering_sorted[i];
    }
