k-means|| a.k.a. Scalable K-Means++ Bahman Bahmani et al. introduced k-means||, which is a faster version of k-means++.

This algorithm is taken from page 4 of their paper, Bahmani, B., Moseley, B., Vattani, A., Kumar, R., & Vassilvitskii, S. (2012). Scalable k-means++. Proceedings of the VLDB Endowment, 5(7), 622-633.
Unfortunately I don't understand those fancy Greek letters, so I need some help to understand how this works. As far as I understand this algorithm is an improved version of k-means++, and it uses oversampling, to reduce the number of iterations: k-means++ has to iterate $k$ times, where $k$ is the number of desired clusters.
I got a very good explanation through a concrete example of how k-means++ works, so I will use the same example again.
Example
I have the following dataset:
(7,1), (3,4), (1,5), (5,8), (1,3), (7,8), (8,2), (5,9), (8,0)
$k = 3$ (number of desired clusters)
$\ell = 2$ (oversampling factor)

I started to calculate it, but I'm not sure if I got it right, and have no idea about Steps 2, 4, or 5.


*

*Step 1: $\mathcal{C} \leftarrow$ sample a point uniformly at random from $X$ 
Let's say the first centroid is $(8,0)$ (same as in k-means++) 

*Step 2: $\psi \leftarrow \phi_X(\mathcal{C})$
no idea

*Step 3: 


*

*$d^2(x, \mathcal{C}) = [2, 41, 74, 73, 58, 65, 4, 90]$ 
We calculate the squared distances to the closest center to each point. In this case we have only one center so far, $(8,0)$.

*$\ell \cdot d^2(x, \mathcal{C}) = [4, 81, 148, 146, 116, 130, 8, 180]$
(Because $\ell = 2$ in this case.)

*$\text{cumulative } \ell \cdot d^2(x, \mathcal{C}) = [4, 85, 233, 379, 495, 625, 633, 813]$
Pick $\ell = 2$ random numbers in the interval $[0, 813)$. Say you pick $246.90$ and $659.42$. They fall in the ranges $[379, 495)$ and $[633, 813)$ which correspond to the 4th and 8th items respectively.

*Repeat it $\mathcal{O}(\log \psi)$ times, but what is $\psi$ (calculated in Step 2) in this case?  


*Step 4: For $x \in \mathcal{C}$, set $w_x$ to be the number of points in $X$ closer
to $x$ than any other point in $\mathcal{C}$.

*Step 5: Recluster the weighted points in $\mathcal{C}$ into $k$ clusters.
Any help in general or in this particular example would be great. 
 A: 
data points: (7,1),(3,4),(1,5),(5,8),(1,3),(7,8),(8,2),(5,9),(8,0) 
l = 2 // oversampling factor
k = 3 // no. of desired clusters

Step 1:

Suppose the first centroid is $\mathcal{C}$ is $\{ c_1\} = \{ (8,0) \}$. $X = \{x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8\}=\{(7,1),(3,4),(1,5),(5,8),(1,3),(7,8),(8,2),(5,9)\}$

Step 2:
$\phi_X(\mathcal{C})$ is the sum of all smallest 2-norm distances (euclidean distance) from all points from the set $X$ to all points from $\mathcal{C}$. In other words, for each point in $X$ find the distance to the closest point in $\mathcal{C}$, in the end compute the sum of all those minimal distances, one for each point in $X$.
Denote with $d^2_{\mathcal{C}}(x_i)$ as the distance from $x_i$ to the closest point in $\mathcal{C}$. We then have $\psi = \sum_{i=1}^{n}d^2_{\mathcal{C}}(x_i)$.
At step 2, $\mathcal{C}$ contains a single element (see step 1), and $X$ is the set of all elements. Thus in this step the $d^2_{\mathcal{C}}(x_i)$ is simply the distance between the point in $\mathcal{C}$ and $x_i$. Thus $\phi = \sum_{i=1}^{n}{||x_i-c||^2}$.

$\psi = \sum_{i=1}^nd^2(x_i,c_1) = 1.41+6.4+8.6+8.54+7.61+8.06+2+9.4 = 52.128$
  $log(\psi) = log(52.128) = 3.95 = 4 (rounded)$

Note however that in step 3, the general formula is applied since $\mathcal{C}$ will contain more than one point.
Step 3:
The for loop is executed for $log(\psi)$ previously computed. 
The drawings are not like you understood. The drawings are independent, which means you will execute a draw for each point in $X$. So, for each point in $X$, denoted as $x_i$, compute a probability from $p_x = l d^2(x,\mathcal{C})/\phi_X(\mathcal{C})$. Here you have $l$ a factor given as parameter, $d^2(x,\mathcal{C})$ is the distance to the closest center, and $\phi_X(\mathcal{C})$ is explained at step 2. 
The algorithm is simply:


*

*iterate in $X$ to find all $x_i$

*for each $x_i$ compute $p_{x_i}$

*generate an uniform number in $[0, 1]$, if is smaller than $p_{x_i}$ select it to form $\mathcal{C'}$

*after you done all draws include selected points from $\mathcal{C'}$ into $\mathcal{C}$


Note that at each step 3 executed in iteration (line 3 of the original algorithm) you expect to select $l$ points from $X$ (this is easily shown writing directly the formula for expectation).
for(int i=0; i<4; i++) {

  // compute d2 for each x_i
  int[] psi = new int[X.size()];
  for(int i=0; i<X.size(); i++) {
    double min = Double.POSITIVE_INFINITY;
    for(int j=0; j<C.size(); j++) {
      if(min>d2(x[i],c[j])) min = norm2(x[i],c[j]);
    }
    psi[i]=min;
  }

  // compute psi
  double phi_c = 0;
  for(int i=0; i<X.size(); i++) phi_c += psi[i];

  // do the drawings
  for(int i=0; i<X.size(); i++) {
    double p_x = l*psi[i]/phi;
    if(p_x >= Random.nextDouble()) {
      C.add(x[i]);
      X.remove(x[i]);
    }
  }
}
// in the end we have C with all centroid candidates
return C;

Step 4:
A simple algorithm for that is to create a vector $w$ of size equals to the number of elements in $\mathcal{C}$, and initialize all its values with $0$. Now iterate in $X$ (elements not selected in as centroids), and for each $x_i \in X$, find the index $j$ of the closest centroid (element from $\mathcal{C}$) and increment $w[j]$ with $1$. In the end you will have the vector $w$ computed properly.
double[] w = new double[C.size()]; // by default all are zero
for(int i=0; i<X.size(); i++) {
  double min = norm2(X[i], C[0]);
  double index = 0;
  for(int j=1; j<C.size(); j++) {
    if(min>norm2(X[i],C[j])) {
      min = norm2(X[i],C[j]);
      index = j;
    }
  }
  // we found the minimum index, so we increment corresp. weight
  w[index]++;
}

Step 5:
Considering the weights $w$ computed at the previous step, you follow kmeans++ algorithm to select only $k$ points as starting centroids. Thus, you will execute $k$ for loops, at each loop selecting a single element, drawn randomly with probability for each element being $p(i) = w(i)/\sum_{j=1}^m{w_j}$. At each step you select one element, and remove it from candidates, also removing its corresponding weight.
for(int k=0; k<K; k++) {
  // select one centroid from candidates, randomly, 
  // weighted by w
  // see kmeans++ and you first idea (which is wrong for step 3)
  ... 
}  

All the previous steps continues, as in the case of kmeans++, with the normal flow of the clustering algorithm
I hope is clearer now.
[Later, later edit]
I found also a presentation made by authors, where you can not clearly that at each iteration multiple points might be selected. The presentation is here.
[Later edit @pera's issue]
It is obvious that $log(\psi)$ depends on data and the issue you raised would be a real problem if the algorithm would be executed on a single host/machine/computer. However you have to note that this variant of kmeans clustering is dedicated to large problems, and for running on distributed systems. Even more, the authors, in the following paragraphs above the algorithm description state the following:

Notice that the size of $C$ is significantly smaller than the input  size;  the  reclustering  can  therefore  be  done  quickly. For  instance,  in  MapReduce,  since  the  number  of  centers is  small  they  can  all  be  assigned  to  a  single  machine  and any provable approximation algorithm (such as k-means++) can  be  used  to  cluster  the  points  to  obtain
  k centers.
  A MapReduce implementation of Algorithm 2 is discussed in Section 3.5.
  While our algorithm is very simple and lends itself to a natural parallel implementation (in $log(\psi)$ rounds ), the challenging part is to show that it has provable guarantees.

Another thing to note is the following note on the same page which states:

In practice, our experimental results in Section 5 show that
  only a few rounds are enough to reach a good solution.

Which means you could run the algorithm not for $log(\psi)$ times, but for a given constant time.
