0

I'm relatively new to Python and would appreciate your help!

Suppose I have two square matrices - one large M and one smaller K - and an integer array of indices ind, not necessarily sorted. The length of ind matches the dimensions of K. In Octave/MATLAB, I can easily do this:

M(ind, ind) = K

This will distribute all the components of K to those positions of M that correspond to indices ind. This is often used in Finite Element computations.

Is there a way to do the same thing just as elegantly in Python? You may assume my M and K are NumPy arrays that were constructed via the operations:

M = np.zeros((12, 12))
K = np.zeros((6, 6))

I did some work on these matrices and filled them with data. My ind array is a NumPy array as well. However, when I do something like

M[ind, ind] = K

I get shape mismatch as an error. Plugging ind.tolist() instead of ind into M won't change anything.

Thanks for any advice!

MajinSaha
  • 188
  • 1
  • 9
  • 1
    In `numpy` `M[ind, ind]` accesses a 'diagonal', while in MATLAB it defines a block. `M[ind[:,None],ind]` gets the block, broadcasting a column vector against a row (1d array). In effect what's easy in MATLAB is a bit harder in `numpy`. But the reverse is true; MATLAB requires an extra function to get the 'flat' indices for a diagonal. – hpaulj Nov 18 '20 at 21:02

2 Answers2

1

You are looking for this:

M[np.ix_(ind,ind)] = K

example:

M = np.zeros((12, 12))
K = np.ones((6, 6))
ind = np.arange(6)

output:

[[1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
Ehsan
  • 12,072
  • 2
  • 20
  • 33
0

If you are using numpy arrays, It can be donde directly using this syntax:

import numpy as np

M = np.zeros((12, 12))
K = np.ones((6, 6))

M[:6, :6] = K
# This is equivalent to:
# M[0:6, 0:6] = K

M will then look like this:

array([[1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

Here you have more information about slice indexing in python

https://www.pythoncentral.io/how-to-slice-listsarrays-and-tuples-in-python/

Alejo Dev
  • 2,290
  • 4
  • 29
  • 45