def cluster(data,inputs,verbose=False):
    """Clusters data
    
    Using the new offset model, this method uses a greedy algorithm to cluster
    the data. It starts with all the data points in separate clusters and tests
    whether combining them increases the overall log-likelihood (LL). It then
    iteratively joins pairs of clusters which cause the greatest increase in
    the LL, until no join increases the LL.
    
    arguments:
    inputs -- the 'X's in a list, one item per cluster
    data -- the 'Y's in a list, one item per cluster
    
    returns a list of the clusters.    
    """
    N=len(data)
    
    
    #Define a set of N active cluster
    active = []
    for p in range(0,N):
        active.append([p])

    loglikes = np.zeros(len(active))
    loglikes[:] = None

    pairloglikes = np.zeros([len(active),len(active)])
    pairloglikes[:] = None
    pairoffset = np.zeros([len(active),len(active)])

    it = 0
    while True:
    
        if verbose:
            it +=1
            print("Iteration %d" % it)
        
        #Compute the log-likelihood of each cluster (add them together)
        for clusti in range(len(active)):
            if verbose:
                sys.stdout.write('.')
                sys.stdout.flush()
            if np.isnan(loglikes[clusti]):
                loglikes[clusti], unused_offset = get_log_likelihood_offset(inputs,data,[clusti])

            #try combining with each other cluster...
            for clustj in range(clusti): #count from 0 to clustj-1
                temp = [clusti,clustj]
                if np.isnan(pairloglikes[clusti,clustj]):
                    pairloglikes[clusti,clustj],pairoffset[clusti,clustj] = get_log_likelihood_offset(inputs,data,temp)

        seploglikes = np.repeat(loglikes[:,None].T,len(loglikes),0)+np.repeat(loglikes[:,None],len(loglikes),1)
        loglikeimprovement = pairloglikes - seploglikes #how much likelihood improves with clustering
        top = np.unravel_index(np.nanargmax(pairloglikes-seploglikes), pairloglikes.shape)

        #if loglikeimprovement.shape[0]<3:
        # #no more clustering to do - this shouldn't happen really unless
        # #we've set the threshold to apply clustering to less than 0
        #    break 
        
        #if theres further clustering to be done...
        if loglikeimprovement[top[0],top[1]]>0:
            active[top[0]].extend(active[top[1]])
            offset=pairoffset[top[0],top[1]]
            inputs[top[0]] = np.vstack([inputs[top[0]],inputs[top[1]]-offset])
            data[top[0]] = np.hstack([data[top[0]],data[top[1]]])
            del inputs[top[1]]
            del data[top[1]]
            del active[top[1]]

            #None = message to say we need to recalculate
            pairloglikes[:,top[0]] = None 
            pairloglikes[top[0],:] = None 
            pairloglikes = np.delete(pairloglikes,top[1],0)
            pairloglikes = np.delete(pairloglikes,top[1],1)
            loglikes[top[0]] = None
            loglikes = np.delete(loglikes,top[1])
        else:
            break
            
        #if loglikeimprovement[top[0],top[1]]>0:
        #    print "joined"
        #    print top
        #    print offset
        #    print offsets
        #    print offsets[top[1]]-offsets[top[0]]

    #TODO Add a way to return the offsets applied to all the time series
    return active