pyqgis - Measure distances between specific features in QGIS for Desktop

I'm completely new using PyQGIS and developing plugins.

Currently I have a layer with the X and Y coordinates of electric distribution poles, where each pole is a feature. What I need to do is measure the distance between all the features in the map, but it must be done pole by pole. By example, take the first pole and get all the distances between the first pole and all the others, then, take the second pole and repeat until the end of the dataset.

Any ideas or references for doing this?

As a reference, I was doing the same thing in Matlab with a dataset of 20 features with something like this:

for i=1:20

  for j=1:20
      if i~=j
        distance = euclidean(P1(i), P1(j));

Where P1(i) and P1(j) are the X-Y coordinates of two different points (Each one Written as a complex number) and "euclidean" is a function that calculates the distance between these two points.

Answer

By using PyQGIS this code works:

from math import sqrt

layer = iface.activeLayer()


features = layer.getFeatures()

points = []

for feature in features:
    geom = feature.geometry().asPoint()
    points.append(geom)

n = len(points)


for i in range(n-1):
    for j in range(n):
        if i < j:
            print i, j, sqrt(points[i].sqrDist(points[j]))

I tried out the code at the Python Console with the layer point (projected in meters) of next image:

and the result was (for upper triangular distance matrix):

i j  distance(m)

0 1 28415.4381005
0 2 48380.9487158
0 3 63608.9197655
0 4 44111.1282702
0 5 19609.1738225
1 2 27252.6459011
1 3 47850.7529946
1 4 37821.6827176
1 5 27147.7980092
2 3 21588.551469

2 4 23663.3870275
2 5 34764.3290912
3 4 23048.9959669
3 5 45857.2099628
4 5 24822.4756566

Editing Note:

The code with gene suggestions:

from math import sqrt
import itertools


layer = iface.activeLayer()

features = layer.getFeatures()

points = [feature.geometry().asPoint() for feature in features]

n = len(points)

list = range(n)


for i,j in itertools.combinations(list, 2):
    print i, j, sqrt(points[i].sqrDist(points[j]))

and the results in a QGIS plugin with the same suggestions: