calculate matrix without using nested if and for loops?

Question

I'm trying to make a matrix A (200, 200) where A[i][j] = A[j][i] follows Bernoulli distribution of probability p if the element is in a same cluster of index(i.e. 0-49, 50-99, 100-149, 150-199) and probability q if not.

So far I tried to implement it straightforwardly like below,

import numpy as np
from scipy.stats import bernoulli
def link_weight(p=0.05, q=0.04):
    A = np.zeros((200, 200))
    for i in range(len(A)):
        for j in range(len(A)):
            if i == j:
                A[i][j] = 0
            elif i < j:
                if i < 50:
                    if j < 50:
                        A[i][j] = bernoulli.rvs(p)
                        A[j][i] = bernoulli.rvs(p)
                    else:
                        A[i][j] = bernoulli.rvs(q)
                        A[j][i] = bernoulli.rvs(q)
                if 50 <= i <= 99:
                    if 50 <= j <= 99:
                        A[i][j] = bernoulli.rvs(p)
                        A[j][i] = bernoulli.rvs(p)
                    else:
                        A[i][j] = bernoulli.rvs(q)
                        A[j][i] = bernoulli.rvs(q)
                if 100 <= i <= 149:
                    if 100 <= j <= 149:
                        A[i][j] = bernoulli.rvs(p)
                        A[j][i] = bernoulli.rvs(p)
                    else:
                        A[i][j] = bernoulli.rvs(q)
                        A[j][i] = bernoulli.rvs(q)
                if 150 <= i <= 199:
                    if 150 <= j <= 199:
                        A[i][j] = bernoulli.rvs(p)
                        A[j][i] = bernoulli.rvs(p)
                    else:
                        A[i][j] = bernoulli.rvs(q)
                        A[j][i] = bernoulli.rvs(q)
    return A

But I wanna make it more readable and faster. How can I prove this code by avoiding nested if/for loops?

For your reference, `bernoulli.rvs(p)` is __not__ equal to `bernoulli.rvs(p)`. — DYZ, Sep 07 '21 at 06:40
@DYZ bernoulli.rvs(p) is not equal to bernoulli.rvs(p) ->bernoulli.rvs(p) is not equal to bernoulli.rvs(q) — I'mahdi, Sep 07 '21 at 06:44

Tzane · Accepted Answer · 2021-09-07T07:10:35.160

1

Numpy indexing is really powerful so you could use that:

import numpy as np
from scipy.stats import bernoulli

def link_weight(p=0.05, q=0.04):
    # Create bernouli array
    A = bernoulli.rvs(q, size=(200, 200))
    # Randomize values p=p
    A[0:50, 0:50] = bernoulli.rvs(p, size=(50,50))
    A[50:100, 50:100] = bernoulli.rvs(p, size=(50,50))
    A[100:150, 100:150] = bernoulli.rvs(p, size=(50,50))
    A[150:200, 150:200] = bernoulli.rvs(p, size=(50,50))
    # Make symmetric
    A = A + A.T - np.diag(A.diagonal())
    return A

print(link_weight())

Symmertic matrix: https://stackoverflow.com/a/10806947/14536215
Indexing: https://numpy.org/doc/stable/reference/arrays.indexing.html

Edit: Change the "symmetrization" to use this: https://stackoverflow.com/a/2573982/14536215

edited Sep 07 '21 at 07:10

answered Sep 07 '21 at 06:49

Tzane

2,752
1
10
21

1

`(A + A.T) / 2` is not good: it is not bernoulli anymore. (See what happens if one element is 0 and the other is 1.) – DYZ Sep 07 '21 at 07:06
1

Also, you can start with a bernoulli matrix of 200x200 and then fill in the "holes." – DYZ Sep 07 '21 at 07:08
@DYZ I just was editing my post for the `(A + A.T) / 2` error as I saw your comments. Also, its a great idea to start with a bernouli matrix so I implemented as well. – Tzane Sep 07 '21 at 07:15

DYZ · Answer 2 · 2022-06-05T22:16:07.693

0

Create 10 square arrays for the respective groups of indexes: 4 for p and 6 for q.

a00 = bernoulli.rvs(p, size=(50, 50))
a11 = bernoulli.rvs(p, size=(50, 50))
a22 = bernoulli.rvs(p, size=(50, 50))
a33 = bernoulli.rvs(p, size=(50, 50))
a01 = bernoulli.rvs(q, size=(50, 50))
a02 = bernoulli.rvs(q, size=(50, 50))
a03 = bernoulli.rvs(q, size=(50, 50))
a12 = bernoulli.rvs(q, size=(50, 50))
a13 = bernoulli.rvs(q, size=(50, 50))
a23 = bernoulli.rvs(q, size=(50, 50))

(Naturally, this can be automated.)

Now, combine the blocks into a big matrix:

A = np.vstack([np.hstack([a00, a01, a02, a03]),
               np.hstack([a01, a11, a12, a13]),
               np.hstack([a02, a12, a22, a23]),
               np.hstack([a03, a13, a23, a33])])

Finally, make that array symmetric by combining its lower triangular part and the transpose of it, less the main diagonal:

np.tril(A) + np.tril(A, k=-1).T

edited Jun 05 '22 at 22:16

answered Sep 07 '21 at 07:01

DYZ

55,249
10
64
93

You can change the shape of the output for `bernouli` with the `size` parameter such as `bernoulli.rvs(0.5, size=(5,5))` – Tzane Sep 07 '21 at 07:04
@Tzane, yes, I just changed it. – DYZ Sep 07 '21 at 07:05
@DYZ `So far I tried to implement it straightforwardly like below, ... But I wanna make it more readable and faster. How can I prove this code by avoiding nested if/for loops?` why angry? I think you can not understand and translate above text, please use translate and understand what want OP – I'mahdi Sep 07 '21 at 07:10

calculate matrix without using nested if and for loops?

2 Answers2