     using Clustering::Tools::boxlimits;
     ASSUME_ALIGNED(coords);
     std::sort(radii.begin(), radii.end(), std::greater<float>());
-    std::size_t n_radii = radii.size();
-    std::vector<float> rad2(n_radii);
-    for (std::size_t i=0; i < n_radii; ++i) {
-      rad2[i] = radii[i]*radii[i];
-    }
-    // sort coordinates on first dimension for neighbor pruning
-    std::vector<float> sorted_coords = dim1_sorted_coords(coords
-                                                        , n_rows
-                                                        , n_cols);
-    // box limits for pruning
-    std::vector<float> blimits = boxlimits(sorted_coords
-                                         , BSIZE_POPS
-                                         , n_rows
-                                         , n_cols);
     int n_gpus;
     cudaGetDeviceCount(&n_gpus);
     if (n_gpus == 0) {

      }
    }
    cudaFree(d_coords);
    cudaFree(d_radii2);
    cudaFree(d_pops);
    return pops;
  }


/*
  Pops
  calculate_populations_partial(const float* coords
                              , const std::vector<float>& sorted_coords
                              , const std::vector<float>& blimits
                              , std::size_t n_rows
                              , std::size_t n_cols
                              , std::vector<float> radii
                              , std::size_t i_from
                              , 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>());
    // setup device & streams
    cudaSetDevice(i_gpu);
    cudaStream_t streams[N_STREAMS];
    for (unsigned int s=0; s < N_STREAMS; ++s) {
      cudaStreamCreate(&streams[s]);
    }
    // copy coords to device
    float* d_coords;
    cudaMalloc((void**) &d_coords
             , sizeof(float) * n_rows * n_cols);
    cudaMemcpy(d_coords
             , coords
             , sizeof(float) * n_rows * n_cols
             , cudaMemcpyHostToDevice);
    float* d_sorted_coords;
    cudaMalloc((void**) &d_sorted_coords
             , sizeof(float) * n_rows * n_cols);
    cudaMemcpy(d_sorted_coords
             , sorted_coords.data()
             , sizeof(float) * n_rows * n_cols
             , cudaMemcpyHostToDevice);
    // copy squared radii to device
    float* d_radii2;
    cudaMalloc((void**) &d_radii2
             , sizeof(float) * n_radii);
    std::vector<float> radii2(radii);
    for (float& r: radii2) {
      r *= r;
    }
    cudaMemcpy(d_radii2
             , radii2.data()
             , sizeof(float) * n_radii
             , cudaMemcpyHostToDevice);
    // tmp buffer for in/out info & reference coords (per stream)
    unsigned int* d_in_radius[N_STREAMS];
    for (unsigned int s=0; s < N_STREAMS; ++s) {
      cudaMalloc((void**) &d_in_radius[s]
               , sizeof(unsigned int) * n_rows * n_radii);
    }
    // result buffer
    unsigned int* d_pops;
    cudaMalloc((void**) &d_pops
             , sizeof(unsigned int) * n_rows * n_radii);
    cudaMemset(d_pops
             , 0
             , sizeof(unsigned int) * n_rows * n_radii);
    // populations per frame
    for (std::size_t i=i_from; i < i_to; ++i) {
      unsigned int i_stream = i % N_STREAMS;
      // 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);
      // +1 for subsequent reduction over otherwise uninitialized range
      in_radius <<< rng+1
                  , BSIZE
                  // shared mem for reference coords
                  , sizeof(float) * n_cols
                  , streams[i_stream] >>> (offset
                                         , d_sorted_coords
                                         , d_coords
                                         , i
                                         , n_rows
                                         , n_cols
                                         , d_radii2
                                         , n_radii
                                         , d_in_radius[i_stream]);
      // compute pops per radius
      if (rng % 2 == 0) {
        rng /= 2;
      } else {
        rng = rng/2 + 1;
      }
      for (unsigned int r=0; r < n_radii; ++r) {
        // pops stored col-wise -> just set an offset ...
        offset = r*n_rows;
        reduce_sum_uint<BSIZE> <<< rng
                                 , BSIZE
                                 , 0
                                 , streams[i_stream] >>> (offset
                                                        , d_in_radius[i_stream]
                                                        , n_rows
                                                        , d_pops
                                                        , r*n_rows + i);
      }
    }
    cudaThreadSynchronize();
    // retrieve pops
    std::vector<unsigned int> tmp_pops(n_rows*n_radii);
    cudaMemcpy(tmp_pops.data()
             , d_pops
             , sizeof(unsigned int) * n_rows * n_radii
             , cudaMemcpyDeviceToHost);
    // sort tmp_pops into pops
    Pops pops;
    for (unsigned int r=0; r < n_radii; ++r) {
      pops[radii[r]].resize(n_rows, 0);
      for (unsigned int i=i_from; i < i_to; ++i) {
        pops[radii[r]][i] = tmp_pops[r*n_rows+i];
      }
    }
    cudaFree(d_sorted_coords);
    cudaFree(d_radii2);
    cudaFree(d_pops);
    for (unsigned int s=0; s < N_STREAMS; ++s) {
      cudaFree(d_in_radius[s]);
    }
    return pops;
  }
*/

  Pops
  calculate_populations(const float* coords

                                         , d_coords
                                         , n_rows
                                         , n_cols
                                         , d_rad2
                                         , n_radii
                                         , d_pops
                                         , i_from
