Attempt to compute a 'null model' of the matrix given a model
    to base itself on.

def null_model(matrix, positions=None, lengths=None,
               model="uniform", noisy=False, circ=False,
               sparsity=False):
    """Attempt to compute a 'null model' of the matrix given a model
    to base itself on.
    """

    n, m = matrix.shape
    positions_supplied = True
    if positions is None:
        positions = range(n)
        positions_supplied = False
    if lengths is None:
        lengths = np.diff(positions)

    N = np.copy(matrix)

    contigs = np.array(positions_to_contigs(positions))

    def is_inter(i, j):
        return contigs[i] != contigs[j]

    diagonal = np.diag(matrix)

    if model == "uniform":
        if positions_supplied:
            trans_contacts = np.array([matrix[i, j]
                                      for i, j in itertools.product(range(n),
                                      range(m)) if is_inter(i, j)])
            mean_trans_contacts = np.average(trans_contacts)
        else:
            mean_trans_contacts = np.average(matrix) - diagonal / len(diagonal)

        N = np.random.poisson(lam=mean_trans_contacts, size=(n, m))
        np.fill_diagonal(N, diagonal)

    elif model == "distance":
        distances = distance_diagonal_law(matrix, positions)
        N = np.array([[distances[min(abs(i - j), n)]
                       for i in range(n)] for j in range(n)])

    elif model == "rippe":

        trans_contacts = np.array([matrix[i, j] for i, j in itertools.product(
            range(n), range(m)) if is_inter(i, j)])
        mean_trans_contacts = np.average(trans_contacts)
        kuhn, lm, slope, d, A = rippe_parameters(matrix, positions, circ=circ)

        def jc(s, frag):
            dist = s - circ * (s**2) / lengths[frag]
            computed_contacts = (0.53 * A * (kuhn**(-3.)) *
                                 (dist ** slope) *
                                 np.exp((d - 2) / (dist + d)))
            return np.maximum(computed_contacts, mean_trans_contacts)

        for i in range(n):
            for j in range(n):
                if not is_inter(i, j) and i != j:
                    posi, posj = positions[i], positions[j]
                    N[i, j] = jc(np.abs(posi - posj) * lm / kuhn, frag=j)
                else:
                    N[i, j] = mean_trans_contacts

    if sparsity:
        contact_sum = matrix.sum(axis=0)
        n = len(contact_sum)
        try:
            from Bio.Statistics import lowess
            trend = lowess(np.array(range(n), dtype=np.float64),
                           contact_sum, f=0.03)
        except ImportError:
            expected_size = int(np.amax(contact_sum) / np.average(contact_sum))
            w = min(max(expected_size, 20), 100)
            trend = np.array([np.average(contact_sum[i:min(i + w, n)]) for
                              i in range(n)])

        cov_score = np.sqrt((trend - np.average(trend)) / np.std(trend))

        N = ((N * cov_score).T) * cov_score

    if noisy:
        if callable(noisy):
            noise = noisy
        return noise(N)
    else:
        return N