# Path: BinaryMatrix.py
# def approximate_from_margins_weights(r, c, w, T = None,
#                                      sort_by_wopt_var = True):
#     """Return approximate samples from row/column-conditional binary matrices.
#                                      
# Return a binary matrix (or a list of binary matrices) sampled
# approximately according to the specified Bernoulli weights,
# conditioned on having the specified margins.
# Inputs:
#   r: row margins, length m
#   c: column margins, length n
#   w: weight matrix, (m x n) matrix with values in (0, +infty)
#   T: number of matrices to sample
#   sort_by_wopt_var: when enabled, column ordering depends on w
# Output:
#   B_sample_sparse: (T default) sparse representation of (m x n) binary matrix
#                    (T >= 1) list of (sparse binary matrices, logQ, logP)
# 
# More explicitly, consider independent Bernoulli random variables
# B(i,j) arranged as an m x n matrix B given the m-vector of row sums
# r and the n-vector of column sums c of the sample, i.e., given that
# sum(B_sample, 1) = r and sum(B_sample, 0) = c.
# 
# An error is generated if no binary matrix agrees with r and c.
# 
# B(i,j) is Bernoulli(p(i,j)) where p(i,j) = w(i,j)/(1+w(i,j)), i.e.,
# w(i,j) = p(i,j)/(1-p(i,j)).  [The case p(i,j) = 1 must be handled by
# the user in a preprocessing step, by converting to p(i,j) = 0 and
# decrementing the row and column sums appropriately.]
# 
# The sparse representation used for output is a matrix giving the
# locations of the ones in the sample. If d = sum(r) = sum(c), then
# B_sample_sparse has dimensions (d x 2). If something goes wrong (due
# to undetected improper input), some of the rows of B_sample_sparse
# may [-1,-1], indicating no entry of B_sample.
# 
# B_sample can be recovered from B_sample_sparse via:
# 
#     B_sample = np.zeros((m,n), dtype=np.bool)
#     for i, j in B_sample_sparse:
#         if i == -1: break 
#         B_sample[i,j] = 1
# """
#     r_prune, c_prune, arrays_prune, unprune = _prune(r, c, w)
#     w_prune = arrays_prune[0]
# 
#     _check_margins(r_prune, c_prune)
# 
#     ### Preprocessing
# 
#     # Sizing (making copies of m and n, as they are mutated during sampling)
#     r_init = r_prune.copy()
#     m, n = len(r_prune), len(c_prune)
#     if (m == 0) or (n == 0):
#         if T:
#             return [unprune([np.empty((0,2)), 0, 0]) for t in xrange(T)]
#         else:
#             return np.empty((0,0))
#     m_init, n_init = m, n
#     assert((m,n) == w_prune.shape)
# 
#     # Sort the row margins (descending)
#     rndx_init = np.argsort(-r_prune)
#     rsort = r_prune[rndx_init]
# 
#     # Balance the weights
#     a_scale, b_scale = canonical_scalings(w_prune, r_prune, c_prune)
#     wopt = apply_scale(w_prune, a_scale, b_scale)
# 
#     # Reorder the columns
#     if sort_by_wopt_var:
#         cndx = np.lexsort((-wopt.var(0), c_prune))
#     else:
#         cndx = np.argsort(c_prune)
#     csort = c_prune[cndx]
#     wopt = wopt[:,cndx]
# 
#     # Precompute log weights
#     logw = np.log(w_prune)
# 
#     # Compute G
#     G = _compute_G(r_prune, m, n, wopt)
# 
#     # Generate the inverse index for the row orders to facilitate fast
#     # sorting during the updating
#     irndx_init = np.argsort(rndx_init)
# 
#     # Compute the conjugate of c
#     cconj_init = conjugate(csort, m)
# 
#     # Get the running total of number of ones to assign
#     count_init = np.sum(rsort)
# 
#     def do_sample():
#         sample_prune = _compute_sample(logw,
#                                        count_init, m_init, n_init,
#                                        r_init, rndx_init, irndx_init,
#                                        csort, cndx, cconj_init,
#                                        G)
#         return unprune(sample_prune)
#     
#     if T:
#         return [do_sample() for t in xrange(T)]
#     else:
#         return do_sample()[0]


import numpy as np
from BinaryMatrix import approximate_from_margins_weights as samp
#!/usr/bin/env python



def fill(m, n, sparse):
    Y = np.zeros((m,n), dtype = np.bool)
    for i, j in sparse:
        if i == -1: break
        Y[i,j] = 1
    return Y

m = 8
n = 6
s = 10

t = 5

while True:
    r = np.random.multinomial(s, np.ones(m) / m)
    c = np.random.multinomial(s, np.ones(n) / n)
    if ((r.min() == 0) or (r.max() >= n) or (c.min() == 0) or (c.max() >= m)):
        continue
    try:
        p = np.random.random((m,n))
        w = p / (1 - p)

# Path: BinaryMatrix.py
# def approximate_from_margins_weights(r, c, w, T = None,
#                                      sort_by_wopt_var = True):
#     """Return approximate samples from row/column-conditional binary matrices.
#                                      
# Return a binary matrix (or a list of binary matrices) sampled
# approximately according to the specified Bernoulli weights,
# conditioned on having the specified margins.
# Inputs:
#   r: row margins, length m
#   c: column margins, length n
#   w: weight matrix, (m x n) matrix with values in (0, +infty)
#   T: number of matrices to sample
#   sort_by_wopt_var: when enabled, column ordering depends on w
# Output:
#   B_sample_sparse: (T default) sparse representation of (m x n) binary matrix
#                    (T >= 1) list of (sparse binary matrices, logQ, logP)
# 
# More explicitly, consider independent Bernoulli random variables
# B(i,j) arranged as an m x n matrix B given the m-vector of row sums
# r and the n-vector of column sums c of the sample, i.e., given that
# sum(B_sample, 1) = r and sum(B_sample, 0) = c.
# 
# An error is generated if no binary matrix agrees with r and c.
# 
# B(i,j) is Bernoulli(p(i,j)) where p(i,j) = w(i,j)/(1+w(i,j)), i.e.,
# w(i,j) = p(i,j)/(1-p(i,j)).  [The case p(i,j) = 1 must be handled by
# the user in a preprocessing step, by converting to p(i,j) = 0 and
# decrementing the row and column sums appropriately.]
# 
# The sparse representation used for output is a matrix giving the
# locations of the ones in the sample. If d = sum(r) = sum(c), then
# B_sample_sparse has dimensions (d x 2). If something goes wrong (due
# to undetected improper input), some of the rows of B_sample_sparse
# may [-1,-1], indicating no entry of B_sample.
# 
# B_sample can be recovered from B_sample_sparse via:
# 
#     B_sample = np.zeros((m,n), dtype=np.bool)
#     for i, j in B_sample_sparse:
#         if i == -1: break 
#         B_sample[i,j] = 1
# """
#     r_prune, c_prune, arrays_prune, unprune = _prune(r, c, w)
#     w_prune = arrays_prune[0]
# 
#     _check_margins(r_prune, c_prune)
# 
#     ### Preprocessing
# 
#     # Sizing (making copies of m and n, as they are mutated during sampling)
#     r_init = r_prune.copy()
#     m, n = len(r_prune), len(c_prune)
#     if (m == 0) or (n == 0):
#         if T:
#             return [unprune([np.empty((0,2)), 0, 0]) for t in xrange(T)]
#         else:
#             return np.empty((0,0))
#     m_init, n_init = m, n
#     assert((m,n) == w_prune.shape)
# 
#     # Sort the row margins (descending)
#     rndx_init = np.argsort(-r_prune)
#     rsort = r_prune[rndx_init]
# 
#     # Balance the weights
#     a_scale, b_scale = canonical_scalings(w_prune, r_prune, c_prune)
#     wopt = apply_scale(w_prune, a_scale, b_scale)
# 
#     # Reorder the columns
#     if sort_by_wopt_var:
#         cndx = np.lexsort((-wopt.var(0), c_prune))
#     else:
#         cndx = np.argsort(c_prune)
#     csort = c_prune[cndx]
#     wopt = wopt[:,cndx]
# 
#     # Precompute log weights
#     logw = np.log(w_prune)
# 
#     # Compute G
#     G = _compute_G(r_prune, m, n, wopt)
# 
#     # Generate the inverse index for the row orders to facilitate fast
#     # sorting during the updating
#     irndx_init = np.argsort(rndx_init)
# 
#     # Compute the conjugate of c
#     cconj_init = conjugate(csort, m)
# 
#     # Get the running total of number of ones to assign
#     count_init = np.sum(rsort)
# 
#     def do_sample():
#         sample_prune = _compute_sample(logw,
#                                        count_init, m_init, n_init,
#                                        r_init, rndx_init, irndx_init,
#                                        csort, cndx, cconj_init,
#                                        G)
#         return unprune(sample_prune)
#     
#     if T:
#         return [do_sample() for t in xrange(T)]
#     else:
#         return do_sample()[0]


# Path: minitest_5.py
import numpy as np
from BinaryMatrix import approximate_from_margins_weights as samp

#!/usr/bin/env python



def fill(m, n, sparse):
    Y = np.zeros((m,n), dtype = np.bool)
    for i, j in sparse:
        if i == -1: break
        Y[i,j] = 1
    return Y

m = 8
n = 6
s = 10

t = 5

while True:
    r = np.random.multinomial(s, np.ones(m) / m)
    c = np.random.multinomial(s, np.ones(n) / n)
    if ((r.min() == 0) or (r.max() >= n) or (c.min() == 0) or (c.max() >= m)):
        continue
    try:
        p = np.random.random((m,n))
        w = p / (1 - p)

        Z = samp(r, c, w, T = (t + 1))