       cudaMemset(d_in_radius[i_stream]
                , 0
                , sizeof(float) * n_rows * n_radii);
+//TODO
       // prune range for faster computation
       // (using largest radius in first dimension)
-      auto min_max_box = Clustering::Tools::min_max_box(blimits
-                                                      , coords[i*n_cols]
-                                                      , radii[0]);
-      unsigned int offset = min_max_box.first * BSIZE;
-      unsigned int rng = (min_max_box.second-min_max_box.first+1);
+//      auto min_max_box = Clustering::Tools::min_max_box(blimits
+//                                                      , coords[i*n_cols]
+//                                                      , radii[0]);
+//      unsigned int offset = min_max_box.first * BSIZE;
+//      unsigned int rng = (min_max_box.second-min_max_box.first+1);
+      unsigned int rng = Tools::min_multiplicator(n_rows, BSIZE);
+      unsigned int offset = 0;
       in_radius <<< rng
                   , BSIZE
                   , 0

      // compute squared dist
      for (j=0; j < n_cols; ++j) {
        c = coords[i_ref*n_cols+j] - sorted_coords[(gid + offset)*n_cols+j];
//        dist2 = fma(c, c, dist2);
        dist2 += c*c;
      }
      // write results: 1.0 if in radius, 0.0 if not
      for (r=0; r < n_radii; ++r) {
        if (dist2 <= radii2[j]) {
          in_radius[r*n_rows + gid] = dist2;
//          in_radius[r*n_rows + gid] = 1.0f;
        } else {
          in_radius[r*n_rows + gid] = -dist2;
//          in_radius[r*n_rows + gid] = 0.0f;
        }
      }
    }

                              , std::size_t i_to
                              , int i_gpu) {
    ASSUME_ALIGNED(coords);
    unsigned int n_radii = radii.size();
    // make sure radii are in descending order
    std::sort(radii.begin(), radii.end(), std::greater<float>());
