Load an alignment, calculate all pairwise distances and variances
        model parameter must be a Substitution model type from phylo_utils

def pairdists(alignment, subs_model, alpha=None, ncat=4, tolerance=1e-6, verbose=False):
    """ Load an alignment, calculate all pairwise distances and variances
        model parameter must be a Substitution model type from phylo_utils """

    # Check
    if not isinstance(subs_model, phylo_utils.models.Model):
        raise ValueError("Can't handle this model: {}".format(model))

    if alpha is None:
        alpha = 1.0
        ncat = 1

    # Set up markov model
    tm = TransitionMatrix(subs_model)

    gamma_rates = discrete_gamma(alpha, ncat)
    partials = alignment_to_partials(alignment)
    seqnames = alignment.get_names()
    nseq = len(seqnames)
    distances = np.zeros((nseq, nseq))
    variances = np.zeros((nseq, nseq))

    # Check the model has the appropriate size
    if not subs_model.size == partials[seqnames[0]].shape[1]:
        raise ValueError("Model {} expects {} states, but the alignment has {}".format(model.name,
                                                                                       model.size,
                                                                                       partials[seqnames[0]].shape[1]))

    nodes = [phylo_utils.likelihood.LnlModel(tm) for seq in range(nseq)]
    for node, header in zip(nodes, seqnames):
        node.set_partials(partials[header])  # retrieve partial likelihoods from partials dictionary

    for i, j in itertools.combinations(range(nseq), 2):
        brlen, var = brent_optimise(nodes[i], nodes[j], verbose=verbose)
        distances[i, j] = distances[j, i] = brlen
        variances[i, j] = variances[j, i] = var
    dm = DistanceMatrix.from_array(distances, names=seqnames)
    vm = DistanceMatrix.from_array(variances, names=seqnames)
    return dm, vm