                            , const std::size_t n_cols
                            , const std::vector<std::size_t> initial_clusters) {
     using Clustering::Tools::min_multiplicator;
+    // data preparation
     std::size_t first_frame_above_threshold;
     double sigma2;
     std::vector<FreeEnergy> fe_sorted;
     std::vector<std::size_t> prev_clustering;
-    std::size_t prev_max_state;
-    // data preparation
+    unsigned int prev_max_state;
     std::tie(prev_clustering
            , first_frame_above_threshold
            , sigma2

    int comp_size = min(bsize, n_rows - offset);
    if (tid < comp_size) {
      for (unsigned int j=0; j < n_cols; ++j) {
        smem[tid*n_cols+j] = sorted_coords[(tid+offset)*n_cols+j];
      }
    }
    __syncthreads();
    if (gid < i_to) {
      unsigned int self_id = gid + 1;
      unsigned int cluster_id = clustering[gid];
      if (cluster_id == 0) {
        cluster_id = self_id;
      }
      unsigned int ref_id = tid+bsize;
      // load reference coordinates for re-use into shared memory
      for (unsigned int j=0; j < n_cols; ++j) {
        smem[ref_id*n_cols+j] = sorted_coords[gid*n_cols+j];
      }
      for (unsigned int i=0; i < comp_size; ++i) {
        float dist2 = 0.0f;
        for (unsigned int j=0; j < n_cols; ++j) {
          float c = smem[ref_id*n_cols+j] - smem[i*n_cols+j];
          dist2 = fma(c, c, dist2);
        }
        if (dist2 < max_dist2) {
          unsigned int min_id = min(cluster_id, i+offset+1);
          if (min_id == self_id) {
            cluster_id = max(cluster_id, i+offset+1);
          } else {
            cluster_id = min_id;
          }
        }
      }
      clustering[gid] = cluster_id;
    }
  }

//  unsigned int
//  follow(const std::vector<unsigned int>& clustering
//       , unsigned int id
//       , unsigned int min_id = 0) {
//    if (id == 0
//     || id == min_id) {
//      return id;
//    } else {
//      return follow(clustering
//                  , clustering[id-1]
//                  , (min_id == 0)
//                      ? id
//                      : std::min(id, min_id));
//    }
//  }

  std::vector<unsigned int>
  sanitize_state_names(std::vector<unsigned int> clustering) {

    //TODO: lump frames into microstates.
    //      be careful: cluster ids do not necessarily
    //                  denote frame ids (especially not with previous
    //                  clustering)

    // map names to ascending index
    unsigned int id=0;
    std::map<unsigned int, unsigned int> name_dict;
    for (unsigned int name: unique_names) {
      ++id;
      name_dict[name] = id;
    }
    // relabel states
    for (unsigned int& i: clustering) {
      i = name_dict[i];
    }
    return clustering;
  }

//TODO remove old code
//    for (unsigned int i=0; i < clustering.size(); ++i) {
//      clustering[i] = follow(clustering
//                           , clustering[i]);
//    }
//    return clustering;
//  }


  std::vector<std::size_t>
  initial_density_clustering(const std::vector<float>& free_energy

                                                        , n_rows
                                                        , initial_clusters);
    // write log
    screening_log(sigma2
                , first_frame_above_threshold
                , fe_sorted);
    float max_dist2 = 4*sigma2;
    // prepare CUDA environment
    int n_gpus = get_num_gpus();
    std::vector<float*> d_coords_sorted(n_gpus);
    std::vector<unsigned int*> d_clustering(n_gpus);
    std::vector<unsigned int*> d_ref_states(n_gpus);
    // sort coords (and previous clustering results)
    // according to free energies
    std::vector<float> tmp_coords_sorted(n_rows * n_cols);
    std::vector<unsigned int> clustering_sorted(n_rows);
    for (unsigned int i=0; i < n_rows; ++i) {
      for (unsigned int j=0; j < n_cols; ++j) {
        tmp_coords_sorted[i*n_cols+j] = coords[fe_sorted[i].first*n_cols+j];
      }
      // intialize energy-sorted clustering results
      if (i < first_frame_above_threshold) {
        clustering_sorted[i] = prev_clustering[fe_sorted[i].first];
      }
    }
    // initialize new (unclustered) frames
    unsigned int prev_last_frame =
      std::distance(clustering_sorted.begin()
                  , std::find(clustering_sorted.begin()
                            , clustering_sorted.end()
                            , 0));
    for (unsigned int i=prev_last_frame; i < first_frame_above_threshold; ++i) {
      clustering_sorted[i] = ++prev_max_state;
    }
    // computational range of single gpu
    unsigned int gpu_rng =
      min_multiplicator(first_frame_above_threshold - prev_last_frame
                      , n_gpus);
