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
– DYZ, Commented Sep 7, 2021 at 6:40
@DYZ bernoulli.rvs(p) is not equal to bernoulli.rvs(p) ->bernoulli.rvs(p) is not equal to bernoulli.rvs(q) — Mahdi F.
– Mahdi F., Commented Sep 7, 2021 at 6:44

Tzane · Accepted Answer · 2021-09-07 07:10:35Z

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

(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.)
Also, you can start with a bernoulli matrix of 200x200 and then fill in the "holes."
@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.

DYZ · Accepted Answer · 2022-06-05 22:16:07Z

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

You can change the shape of the output for bernouli with the size parameter such as bernoulli.rvs(0.5, size=(5,5))
@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

Collectives™ on Stack Overflow

calculate matrix without using nested if and for loops?

2 Answers 2

3 Comments

3 Comments

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

3 Comments

Linked

Related