Operações Avançadas no Numpy#

Vale lembrar que sempre que iniciamos um novo arquivo jupyter e o kernel estiver zerado (nada tiver sido compilado) necessitamos compilar os módulos que necessitamos importar:

import numpy as np

Concatenando arrays#

As funções vstack e hstack permitem concatenar dois ou mais arrays verticalmente e horizontalmente, respectivamente:

a = np.arange(6).reshape(3, 2)
b = np.random.random((3, 2))

print(np.vstack((a, b)))
print()

print(np.hstack((a, b)))
[[0.         1.        ]
 [2.         3.        ]
 [4.         5.        ]
 [0.77470779 0.84853832]
 [0.68088904 0.62880377]
 [0.81534892 0.45023848]]

[[0.         1.         0.77470779 0.84853832]
 [2.         3.         0.68088904 0.62880377]
 [4.         5.         0.81534892 0.45023848]]

A função column_stack concatena dois ou mais arrays unidimensionais na forma de colunas em um array 2D resultante:

a = np.arange(3)
b = np.random.random(3)

print(np.column_stack((a, b)))
print()

print(np.hstack((a, b)))  # resultado diferente
[[0.         0.47960465]
 [1.         0.36236234]
 [2.         0.31575287]]

[0.         1.         2.         0.47960465 0.36236234 0.31575287]

Essas funções de concatenação são todas casos especiais de uso mais comum da função concatenate, que permite definir o eixo sobre o qual ocorrerá a concatenação.

Copiando Arrays#

Nota

Simplesmente usar «=» não faz uma cópia, mas, assim como com listas, você terá vários nomes apontando para o mesmo objeto ndarray.

Portanto, precisamos entender se dois arrays, A e B, apontam para:

  • o mesmo array, incluindo forma e espaço de dados/memória

  • o mesmo espaço de dados/memória, mas talvez formas diferentes (uma view)

  • uma cópia separada dos dados (ou seja, armazenada completamente separada na memória)

Todos esses casos são possíveis:

  • B = A

    isso é atribuição. Nenhuma cópia é feita. A e B apontam para os mesmos dados na memória e compartilham a mesma forma, etc. Eles são apenas dois rótulos diferentes para o mesmo objeto na memória.

  • B = A[:]

    isso é uma view ou cópia rasa. As informações de forma de A e B são armazenadas independentemente, mas ambos apontam para o mesmo local de memória para os dados.

  • B = A.copy()

    isso é uma cópia profunda. Um objeto completamente separado será criado na memória, com um local completamente separado na memória.

Vamos ver exemplos:

a = np.arange(10)
print(a)
[0 1 2 3 4 5 6 7 8 9]

Aqui está a atribuição — podemos simplesmente usar o operador is para testar a igualdade.

b = a
b is a
True

Como b e a são o mesmo, alterações na forma de um são refletidas no outro — nenhuma cópia é feita.

b.shape = (2, 5)
print(b)
a.shape
[[0 1 2 3 4]
 [5 6 7 8 9]]
(2, 5)
b is a
True
print(a)
[[0 1 2 3 4]
 [5 6 7 8 9]]

Uma cópia rasa cria uma nova view do array — os dados são os mesmos, mas as propriedades do array podem ser diferentes.

a = np.arange(12)
c = a[:]
a.shape = (3,4)

print(a)
print(c)
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
[ 0  1  2  3  4  5  6  7  8  9 10 11]

Como os dados subjacentes estão na mesma memória, alterar um elemento de um é refletido no outro.

c[1] = -1
print(a)
[[ 0 -1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

Até mesmo fatiamentos em um array são apenas views, ainda apontando para a mesma memória.

d = c[3:8]
print(d)
[3 4 5 6 7]
d[:] = 0 
print(a)
print(c)
print(d)
[[ 0 -1  2  0]
 [ 0  0  0  0]
 [ 8  9 10 11]]
[ 0 -1  2  0  0  0  0  0  8  9 10 11]
[0 0 0 0 0]

Existem várias maneiras de verificar se dois arrays são iguais, são views, possuem seus próprios dados, etc.

print(c is a)
print(c.base is a)
#print(c.flags.owndata)
#print(a.flags.owndata)
False
True

Para fazer uma cópia dos dados do array com a qual você possa trabalhar independentemente do original, você precisa de uma cópia profunda.

d = a.copy()
d[:,:] = 0.0

print(a)
print(d)
[[ 0 -1  2  0]
 [ 0  0  0  0]
 [ 8  9 10 11]]
[[0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]

Indexação booleana (também conhecida como mascaramento)#

Existem várias maneiras interessantes de indexar arrays para acessar apenas os elementos que atendem a uma determinada condição.

a = np.arange(12).reshape(3,4)
a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

Aqui, definimos todos os elementos do array que são > 4 para zero.

a > 4
array([[False, False, False, False],
       [False,  True,  True,  True],
       [ True,  True,  True,  True]])
a[a > 4]
array([ 5,  6,  7,  8,  9, 10, 11])
a[a > 4] = 0
a
array([[0, 1, 2, 3],
       [4, 0, 0, 0],
       [0, 0, 0, 0]])

E agora, todos os zeros para -1.

a[a == 0] = -1
a
array([[-1,  1,  2,  3],
       [ 4, -1, -1, -1],
       [-1, -1, -1, -1]])
a == -1
array([[ True, False, False, False],
       [False,  True,  True,  True],
       [ True,  True,  True,  True]])

Vamos considerar outro exemplo, envolvendo dois arrays distintos

a = np.arange(12)
b = np.arange(12)*123
b
array([   0,  123,  246,  369,  492,  615,  738,  861,  984, 1107, 1230,
       1353])
a>6
array([False, False, False, False, False, False, False,  True,  True,
        True,  True,  True])
b[a>6]
array([ 861,  984, 1107, 1230, 1353])

Se tivermos 2 testes, precisamos usar logical_and() ou logical_or().

a = np.arange(12).reshape(3,4)
a[np.logical_and(a > 3, a <= 9)] = 0.0
a
array([[ 0,  1,  2,  3],
       [ 0,  0,  0,  0],
       [ 0,  0, 10, 11]])

Evitando loops#

De maneira geral, você quer evitar loops sobre os elementos de um array.

Aqui, vamos criar coordenadas 1D de x e y e depois tentar preencher um array maior.

M = 3200
N = 6400
xmin = ymin = 0.0
xmax = ymax = 1.0

x = np.linspace(xmin, xmax, M, endpoint=False)
y = np.linspace(ymin, ymax, N, endpoint=False)

print(x.shape)
print(y.shape)
(3200,)
(6400,)

Vamos medir o tempo do nosso código.

import time
t0 = time.time()

g = np.zeros((M, N))

# Índices, que ideia terrível
for i in range(M):
    for j in range(N):
        g[i,j] = np.sin(2.0*np.pi*x[i]*y[j])
        
t1 = time.time()
print("tempo decorrido: {} s".format(t1-t0))
tempo decorrido: 26.826146602630615 s

Agora, vamos fazer isso usando toda a sintaxe de arrays. Primeiro, vamos estender nossos arrays de coordenadas 1D para 2D. O Numpy tem uma função para isso (meshgrid()).

Vamos ver como o meshgrid() funciona primeiro.

x2d, y2d = np.meshgrid([0,1,2,3], [10,20,30], indexing="ij")
x2d
array([[0, 0, 0],
       [1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])
y2d
array([[10, 20, 30],
       [10, 20, 30],
       [10, 20, 30],
       [10, 20, 30]])

Vemos que isso cria 2 arrays bidimensionais, um com os valores de x variando ao longo das linhas e outro com os valores de y variando ao longo das colunas. Isso significa que podemos indexar todos os pontos (x[i], y[j]) através do par x2d, y2d.

Agora, vamos fazer o mesmo exemplo usando este método.

t0 = time.time()
x2d, y2d = np.meshgrid(x, y, indexing="ij")
g2 = np.sin(2.0*np.pi*x2d*y2d)
t1 = time.time()
print("tempo decorrido: {} s".format(t1-t0))
tempo decorrido: 0.2992236614227295 s

No meu laptop, isso é cerca de:Uma verificação final para garantir que eles forneçam a mesma resposta.

print(np.max(np.abs(g2-g)))
0.0

Exemplo: Diferenciação Numérica#

Agora, queremos construir uma derivada, $\( \frac{d f}{dx} \)$

x = np.linspace(0, 2*np.pi, 25)
f = np.sin(x)

Queremos fazer isso sem loops — usaremos views em arrays deslocados entre si. Lembre-se de cálculo que uma derivada é aproximadamente: $\( \frac{df}{dx} = \frac{f(x+h) - f(x)}{h} \)\( Aqui, tomaremos \)h$ como um único elemento adjacente.

f
array([ 0.00000000e+00,  2.58819045e-01,  5.00000000e-01,  7.07106781e-01,
        8.66025404e-01,  9.65925826e-01,  1.00000000e+00,  9.65925826e-01,
        8.66025404e-01,  7.07106781e-01,  5.00000000e-01,  2.58819045e-01,
        1.22464680e-16, -2.58819045e-01, -5.00000000e-01, -7.07106781e-01,
       -8.66025404e-01, -9.65925826e-01, -1.00000000e+00, -9.65925826e-01,
       -8.66025404e-01, -7.07106781e-01, -5.00000000e-01, -2.58819045e-01,
       -2.44929360e-16])
f[1:]
array([ 2.58819045e-01,  5.00000000e-01,  7.07106781e-01,  8.66025404e-01,
        9.65925826e-01,  1.00000000e+00,  9.65925826e-01,  8.66025404e-01,
        7.07106781e-01,  5.00000000e-01,  2.58819045e-01,  1.22464680e-16,
       -2.58819045e-01, -5.00000000e-01, -7.07106781e-01, -8.66025404e-01,
       -9.65925826e-01, -1.00000000e+00, -9.65925826e-01, -8.66025404e-01,
       -7.07106781e-01, -5.00000000e-01, -2.58819045e-01, -2.44929360e-16])
f[:-1]
array([ 0.00000000e+00,  2.58819045e-01,  5.00000000e-01,  7.07106781e-01,
        8.66025404e-01,  9.65925826e-01,  1.00000000e+00,  9.65925826e-01,
        8.66025404e-01,  7.07106781e-01,  5.00000000e-01,  2.58819045e-01,
        1.22464680e-16, -2.58819045e-01, -5.00000000e-01, -7.07106781e-01,
       -8.66025404e-01, -9.65925826e-01, -1.00000000e+00, -9.65925826e-01,
       -8.66025404e-01, -7.07106781e-01, -5.00000000e-01, -2.58819045e-01])
print(len(f[:-1]))
len(f[1:])
24
24
dx = x[1] - x[0]
dfdx = (f[1:] - f[:-1])/dx
dfdx
array([ 0.98861593,  0.92124339,  0.79108963,  0.60702442,  0.38159151,
        0.13015376, -0.13015376, -0.38159151, -0.60702442, -0.79108963,
       -0.92124339, -0.98861593, -0.98861593, -0.92124339, -0.79108963,
       -0.60702442, -0.38159151, -0.13015376,  0.13015376,  0.38159151,
        0.60702442,  0.79108963,  0.92124339,  0.98861593])

Ordenando arrays#

O NumPy fornece diversas funções para ordenar arrays. A função principal é o sort, que retorna uma cópia do array ordenado, com base no array fornecido como argumento:

a = np.array([8, 12, 1, 0, -2, -6, 2, 7, 13])
print(np.sort(a))
[-6 -2  0  1  2  7  8 12 13]

Se a for um array multidimensional, é possível utilizar o parâmetro axis para especificar o eixo ao longo do qual a ordenação será realizada.

b = np.array([[8, 12, 1], [0, -2, -6], [2, 7, 13]])
print(np.sort(b, axis=None))  # se axis for None, o array é aplainado antes de ordenar
print(np.sort(b, axis=0))  # axis=0 irá ordenar através das linhas
print(np.sort(b, axis=1))  # axis=1 irá ordenar através das colunas
[-6 -2  0  1  2  7  8 12 13]
[[ 0 -2 -6]
 [ 2  7  1]
 [ 8 12 13]]
[[ 1  8 12]
 [-6 -2  0]
 [ 2  7 13]]

Também é possível escolher o algoritmo de ordenação usando o parâmetro kind (cujo valor padrão é “quicksort”):

print(np.sort(a, kind='mergesort'))
[-6 -2  0  1  2  7  8 12 13]

Uma outra função relevante é argsort, que, em vez de devolver o array ordenado, retorna os índices que, quando aplicados ao array, o ordenariam corretamente. A função argsort recebe os mesmos parâmetros que a função sort.

print(a)
print(np.argsort(a))
print(a[np.argsort(a)])
print()

b = np.array([
    [8, 12, 1], 
    [0, -2, -6], 
    [2, 7, 13]
])

print(np.argsort(b, axis=0))
print()
print(np.argsort(b, axis=1))
[ 8 12  1  0 -2 -6  2  7 13]
[5 4 3 2 6 7 0 1 8]
[-6 -2  0  1  2  7  8 12 13]

[[1 1 1]
 [2 2 0]
 [0 0 2]]

[[2 0 1]
 [2 1 0]
 [0 1 2]]

Localizando valores#

As funções amax e amin são responsáveis por retornar, respectivamente, o maior e o menor valor do array. É possível, naturalmente, especificar o eixo desejado utilizando o parâmetro axis.

b = np.array([
    [8, 12, 1], 
    [0, -2, -6], 
    [2, 7, 13]
])

print(np.amax(b))
print(np.amin(b))
print()

print(np.amax(b, axis=0))  # maior elemento de cada coluna
print(np.amin(b, axis=1))  # menor elemento de cada linha
13
-6

[ 8 12 13]
[ 1 -6  2]

De forma semelhante à função argsort, as funções argmax e argmin retornam os índices, respectivamente, do maior e do menor elemento do array, com a possibilidade de definir o eixo desejado.

print(np.argmax(b))  # considera o array aplainado (b.ravel())
print(np.argmin(b))
print()

print(np.argmax(b, axis=0))  # índice do maior elemento de cada coluna
print(np.argmin(b, axis=1))  # índice do menor elemento de cada linha
8
5

[0 0 2]
[2 2 0]

Uma função adicional relevante para localizar valores em um array é a função where, que retorna uma tupla com os índices onde uma condição é verdadeira.

print(b)
print(b < 7)
print(np.where(b < 7))
[[ 8 12  1]
 [ 0 -2 -6]
 [ 2  7 13]]
[[False False  True]
 [ True  True  True]
 [ True False False]]
(array([0, 1, 1, 1, 2]), array([2, 0, 1, 2, 0]))

Geração de números aleatórios#

O NumPy oferece várias funções para gerar números aleatórios através do pacote random. Por exemplo, para gerar números uniformemente distribuídos no intervalo [0.0, 1.0), pode-se utilizar a função random:

print(np.random.random())  # apenas um número
print(np.random.random(3))  # vetor com 3 números aleatórios
print(np.random.random((3, 2)))  # matriz aleatória 3x2
0.6301598847448017
[0.14638481 0.91042623 0.49656389]
[[0.60030749 0.6365096 ]
 [0.61274727 0.95872926]
 [0.54761529 0.73669038]]

Para gerar números inteiros uniformemente distribuídos no intervalo [a, b), pode-se utilizar a função randint:

print(np.random.randint(2, 5))  # inteiro em [2, 5)
print(np.random.randint(5))  # inteiro em [0, 5)
print(np.random.randint(2, 5, (3, 2)))  # matriz aleatória 3x2 em [2, 5)
3
3
[[4 3]
 [2 4]
 [3 4]]

Para extrair uma amostra aleatória de um array, pode-se utilizar a função choice:

a = np.arange(8)
print(a)
print()
print(np.random.choice(a))  # um elemento aleatório de a
print()
print(np.random.choice(a, 3))  # três elementos aleatórios de a
print()
print(np.random.choice(a, (3, 2)))  # seis elementos aleatórios de a na forma 3x2
print()
print(np.random.choice(a, 3, replace=False))  # três elementos aleatórios de a sem reposição
print()
print(
    np.random.choice(
        a, 
        3, 
        replace=False,
        p=a/np.sum(a)
    )
)  # três elementos aleatórios de a sem reposição e com diferentes probabilidades
print()
print(np.random.choice(5))  # um elemento aleatório no range(5)
[0 1 2 3 4 5 6 7]

1

[6 6 4]

[[3 1]
 [5 0]
 [7 4]]

[5 1 2]

[2 4 7]

1

Para reorganizar aleatoriamente os elementos de um array, pode-se usar a função shuffle, que altera o próprio array, ou seja, não retorna um novo array como resultado. Exemplo:

a = np.arange(15)

np.random.shuffle(a)
print(a)
[10  4  8  1 12  5  6 11 14 13  9  7  0  2  3]

Os arrays multidimensionais são reorganizados aleatoriamente apenas no primeiro eixo:

a = np.arange(12).reshape(3, 4)
print(a)
np.random.shuffle(a)
print(a)
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

Apenas por curiosidade, para reorganizar aleatoriamente as colunas de um array 2D, pode-se utilizar a função de programação funcional apply_along_axis, que aplica uma função ao longo de um eixo do array.

np.apply_along_axis(np.random.shuffle, 0, a)
print(a)
print()
np.apply_along_axis(np.random.shuffle, 1, a)
print(a)
[[ 0  1  6 11]
 [ 8  5  2  3]
 [ 4  9 10  7]]

[[ 0  6  1 11]
 [ 8  2  5  3]
 [ 7  9  4 10]]

De maneira similar à função shuffle, a função permutation também reorganiza aleatoriamente os elementos de um array, mas retorna um novo array como resultado. Quando o array é multidimensional, ela também afeta apenas o primeiro eixo. Uma diferença em relação à função shuffle é que permutation pode receber um número inteiro x como parâmetro, ao invés de um array. Nesse caso, ela retorna uma permutação do range(x).

a = np.arange(12)
print(a)
print(np.random.permutation(a))
print(a)
print()
print(np.random.permutation(10))
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[10  0  6  4  7  9  5  2 11  1  8  3]
[ 0  1  2  3  4  5  6  7  8  9 10 11]

[7 0 4 2 8 1 9 6 3 5]

Em simulações, é essencial garantir que os resultados possam ser reproduzidos. Para isso, é necessário que os valores aleatórios gerados durante a execução do script sejam consistentes a cada execução. Isso pode ser alcançado definindo a «semente» do gerador de números pseudoaleatórios do NumPy, por meio da função random.seed. Por exemplo, o número gerado pelo código abaixo será sempre o mesmo.

np.random.seed(37)
print(np.random.random())
0.9444966028573069

Álgebra Linear#

A maioria das funções de álgebra linear do NumPy estão acessíveis por meio do módulo np.linalg. Entre as funções importantes, destaca-se a de calcular a inversa de uma matriz:

a = np.array([[1.0, 2.0], [3.0, 4.0]])
print(np.linalg.inv(a))
[[-2.   1. ]
 [ 1.5 -0.5]]

Podemos calcular o determinante e o traço de uma matriz:

print(np.linalg.det(a))
print(np.trace(a))  # o traço está disponível no próprio np
-2.0000000000000004
5.0

Obter a diagonal de uma matriz como um vetor:

np.diag(a)
array([1., 4.])

Obter uma matriz diagonal a partir de um vetor:

np.diag(np.diag(a))
array([[1., 0.],
       [0., 4.]])

Resolve um sistema de equações lineares:

y = np.array([[5.0], [7.0]])
np.linalg.solve(a, y)
array([[-3.],
       [ 4.]])

Encontrar os autovalores e os autovetores de uma matriz:

np.linalg.eig(a)  # retorna uma tupla (autovalores, autovetores)
EigResult(eigenvalues=array([-0.37228132,  5.37228132]), eigenvectors=array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]]))