import copy

def ward_cluster(data, distances = None):
  '''Another go at ward clustering. Much of this depends on the fact that we
  merge downwards.'''

  n_data = len(data)
  clusters = [[j] for j in range(n_data)]
  indices = [j for j in range(len(clusters))]

  history = []

  # if we have not already been given computed distances, compute them,
  # based on common sense...

  if not distances:
    distances = { }
    for i in range(n_data - 1):
      for j in range(i + 1, n_data):
        try:
          distance = sum([(di - dj) * (di - dj) for \
                          di, dj in zip(data[i], data[j])])
        except TypeError, e:
          distance = (data[i] - data[j]) * (data[i] - data[j])

        distances[(i, j)] = distance / 2

  while len(clusters) > 1:

    # compute nearest neighbours

    neighbour = { }

    for i in range(len(clusters) - 1):

      dmin = 1.0e100

      for j in range(i + 1, len(clusters)):
        if distances[(indices[i], indices[j])] >= dmin:
          continue
        dmin = distances[(indices[i], indices[j])]
        neighbour[i] = dmin, j

    # find nearest of near neighbours, i.e. the target to merge the source
    # into

    dmin = 1.0e100
    target = None
    source = None

    for i in range(len(clusters) - 1):
      if neighbour[i][0] >= dmin:
        continue

      dmin = neighbour[i][0]
      target = i
      source = neighbour[i][1]

    # perform the merge

    source_cluster = clusters.pop(source)
    source_index = indices.pop(source)
    history.append((copy.deepcopy(clusters[target]), source_cluster, dmin))

    # update distances

    for i in range(len(clusters)):
      if i == target:
        continue

      n_tot = len(clusters[target]) + len(source_cluster) + \
          len(clusters[i])

      if indices[i] < indices[target]:
        d_target = distances[(indices[i], indices[target])]
      else:
        d_target = distances[(indices[target], indices[i])]

      if indices[i] < source_index:
        d_source = distances[(indices[i], source_index)]
      else:
        d_source = distances[(source_index, indices[i])]

      index_low, index_high = min(indices[i], indices[target]), \
          max(indices[i], indices[target])

      distances[(index_low, index_high)] = (
          (len(clusters[i]) + len(clusters[target])) * d_target +
          (len(clusters[i]) + len(source_cluster)) * d_source -
          len(clusters[i]) * dmin) / float(n_tot)

    # delete the source related distances

    for index in sorted(distances):
      if source_index in index:
        del(distances[index])

    clusters[target] += source_cluster

  return history

if __name__ == '__main__':

  data = [2.0, 6.1, 5.1, 6.2, 2.2, 2.3, 2.4, 0.1, 0.2, 0.3]

  history = ward_cluster(data)

  for i, j, d in history:
    print '%.2f' % d, [data[_i] for _i in i], [data[_j] for _j in j]
