code: def pairdists(alignment, subs_model, alpha=None, ncat=4, tolerance=1e-6, verbose=False):
      
  
      # 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
docs: Load an alignment, calculate all pairwise distances and variances
          model parameter must be a Substitution model type from phylo_utils