def null_model(matrix, positions=None, lengths=None,<EOL>model="<STR_LIT>", noisy=False, circ=False,<EOL>sparsity=False):
n, m = matrix.shape<EOL>positions_supplied = True<EOL>if positions is None:<EOL><INDENT>positions = range(n)<EOL>positions_supplied = False<EOL><DEDENT>if lengths is None:<EOL><INDENT>lengths = np.diff(positions)<EOL><DEDENT>N = np.copy(matrix)<EOL>contigs = np.array(positions_to_contigs(positions))<EOL>def is_inter(i, j):<EOL><INDENT>return contigs[i] != contigs[j]<EOL><DEDENT>diagonal = np.diag(matrix)<EOL>if model == "<STR_LIT>":<EOL><INDENT>if positions_supplied:<EOL><INDENT>trans_contacts = np.array([matrix[i, j]<EOL>for i, j in itertools.product(range(n),<EOL>range(m)) if is_inter(i, j)])<EOL>mean_trans_contacts = np.average(trans_contacts)<EOL><DEDENT>else:<EOL><INDENT>mean_trans_contacts = np.average(matrix) - diagonal / len(diagonal)<EOL><DEDENT>N = np.random.poisson(lam=mean_trans_contacts, size=(n, m))<EOL>np.fill_diagonal(N, diagonal)<EOL><DEDENT>elif model == "<STR_LIT>":<EOL><INDENT>distances = distance_diagonal_law(matrix, positions)<EOL>N = np.array([[distances[min(abs(i - j), n)]<EOL>for i in range(n)] for j in range(n)])<EOL><DEDENT>elif model == "<STR_LIT>":<EOL><INDENT>trans_contacts = np.array([matrix[i, j] for i, j in itertools.product(<EOL>range(n), range(m)) if is_inter(i, j)])<EOL>mean_trans_contacts = np.average(trans_contacts)<EOL>kuhn, lm, slope, d, A = rippe_parameters(matrix, positions, circ=circ)<EOL>def jc(s, frag):<EOL><INDENT>dist = s - circ * (s**<NUM_LIT:2>) / lengths[frag]<EOL>computed_contacts = (<NUM_LIT> * A * (kuhn**(-<NUM_LIT>)) *<EOL>(dist ** slope) *<EOL>np.exp((d - <NUM_LIT:2>) / (dist + d)))<EOL>return np.maximum(computed_contacts, mean_trans_contacts)<EOL><DEDENT>for i in range(n):<EOL><INDENT>for j in range(n):<EOL><INDENT>if not is_inter(i, j) and i != j:<EOL><INDENT>posi, posj = positions[i], positions[j]<EOL>N[i, j] = jc(np.abs(posi - posj) * lm / kuhn, frag=j)<EOL><DEDENT>else:<EOL><INDENT>N[i, j] = mean_trans_contacts<EOL><DEDENT><DEDENT><DEDENT><DEDENT>if sparsity:<EOL><INDENT>contact_sum = matrix.sum(axis=<NUM_LIT:0>)<EOL>n = len(contact_sum)<EOL>try:<EOL><INDENT>from Bio.Statistics import lowess<EOL>trend = lowess(np.array(range(n), dtype=np.float64),<EOL>contact_sum, f=<NUM_LIT>)<EOL><DEDENT>except ImportError:<EOL><INDENT>expected_size = int(np.amax(contact_sum) / np.average(contact_sum))<EOL>w = min(max(expected_size, <NUM_LIT:20>), <NUM_LIT:100>)<EOL>trend = np.array([np.average(contact_sum[i:min(i + w, n)]) for<EOL>i in range(n)])<EOL><DEDENT>cov_score = np.sqrt((trend - np.average(trend)) / np.std(trend))<EOL>N = ((N * cov_score).T) * cov_score<EOL><DEDENT>if noisy:<EOL><INDENT>if callable(noisy):<EOL><INDENT>noise = noisy<EOL><DEDENT>return noise(N)<EOL><DEDENT>else:<EOL><INDENT>return N<EOL><DEDENT>
Attempt to compute a 'null model' of the matrix given a model
    to base itself on.