Aula 3 - NumPy e SciPy

A parte do numpy tem como principal base o guia oficial do usuário disponibilizado gratuitamente aqui.

Observação: As notas a seguir estão em constante desenvolvimento, aperfeiçoamento e revisões.

Conteúdo

Pré-requisitos

Podemos explorar nossos estudos desse ponto em diante de duas formas:

  • Localmente com Jupyter Notebooks;
  • Remotamente com o Google Colab.

Este último já vem previamente configurado com todas as libs que utilizaremos na parte numérica do curso. Para a primeira opção, para que tudo corra bem, precisamos garantir que todo mundo está utilizando a mesma versão tanto do Python quanto em suas libs. Para isso, vamos utilizar o conda, um gerenciador open source de pacotes e ambientes.

Instalação do conda (Miniconda)

Caso sua máquina não possua o conda instalado, você pode baixar o instalador condizente com seu computador e sistema operacional neste link, basta seguir as instruções lá providenciadas. Com isso, você instalará o Miniconda - a versão com "o mínimo necessário" do conda - em seu computador, que é o que iremos precisar aqui.

Utilizando e instalando ambientes (envs)

Após instalado o conda e reiniciado seu terminal, você deve ser capaz de checar os ambientes (environments, que chamaremos apenas de envs) instalados (execute o comando abaixo sem ! no seu terminal):

In [1]:
!conda env list
# conda environments:
#
base                     /home/volps/miniconda3
brumadinho               /home/volps/miniconda3/envs/brumadinho
daml                     /home/volps/miniconda3/envs/daml
fenics16-1               /home/volps/miniconda3/envs/fenics16-1
fenics16-2               /home/volps/miniconda3/envs/fenics16-2
fenics17-1               /home/volps/miniconda3/envs/fenics17-1
fenics17-2               /home/volps/miniconda3/envs/fenics17-2
fenics18-1               /home/volps/miniconda3/envs/fenics18-1
sci                      /home/volps/miniconda3/envs/sci
sci_env               *  /home/volps/miniconda3/envs/sci_env

Este comando lista os envs instalados. Note que já temos um env que você provavelmente ainda não tem (a menos que vc esteja executando este notebook no Colab). Vamos resolver isso?

Para tanto, siga a sequência de passos abaixo:

  1. Clone o repositório remoto do curso para sua máquina com o comando git clone https://github.com/pedrosiracusa/curso_python_eamc.git e entre no diretório com cd curso_python_eamc

  2. Adicione o canal (conda-forge) que utilizaremos para baixar os pacotes corretos no conda: conda config --add channels conda-forge

  3. Crie o mesmo pacote que temos - o sci_env - por meio do arquivo de dependências sci_env.yml com o seguinte comando: conda env create -f sci_env.yml

  4. Aguarde a instala...zzzzZZZZZZZZZZZZZZZZzzZZZzz

  5. Após a instalação, reinicie seu terminal e acesse o ambiente com o comando source activate sci_env (Linux) ou conda activate sci_env (Windows).

Pronto, com isso temos exatamente o mesmo ambiente, com as mesmas bibliotecas!

NumPy


O que é?

É um pacote fundamental para aplicações de computação científica em Python. Com ele, o desenvolvedor é capaz de:

  • Montar objetos matriciais N-dimensionais e objetos derivados desses;
  • Realizar operações eficientes de Álgebra Linear computacional;
  • Operações básicas de Estatística;
  • Utilizar geração de números pseudo-aleatórios;
  • Operações vetorizadas;
  • Integrar com códigos em C/C++ e Fortran;

E muito mais!

Conceitos básicos

Importando o pacote

In [2]:
import numpy as np

Tipo básico: numpy.ndarray

O tipo básico da numpy é o numpy.ndarray, o quão chamaremos aqui frequentemente de "array".

Inicializando

Um simples vetor:

In [3]:
array1 = np.array([1., 2., 3.])

array1
Out[3]:
array([1., 2., 3.])

Uma matriz:

In [4]:
matrix1 = np.array(
    [[1., 2., 3.],
     [4., 5., 6.],
     [7., 8., 9.]]
)

matrix1
Out[4]:
array([[1., 2., 3.],
       [4., 5., 6.],
       [7., 8., 9.]])

Note que a matriz nada mais é que um array de arrays.

Ao longo do curso, veremos outras formas de inicializar arrays.

Atenção!!!

A entrada para a criação de arrays deve ser listas ou tuplas:

In [5]:
array_tuple = np.array((1, 2, 3))

array_tuple
Out[5]:
array([1, 2, 3])
In [6]:
matrix_mixed = np.array(
    [(1, 2, 3),
    [4, 5, 6]]
)

matrix_mixed
Out[6]:
array([[1, 2, 3],
       [4, 5, 6]])
Inicializando arrays preenchidos com zeros ou um

Em algumas situações, já são conhecidas as demandas e estruturas dos arrays que vão ser utilizados. Neste caso, o numpy oferece funções para criar arrays preenchidos com zeros ou um.

  • Vetor
In [7]:
array_zeros = np.zeros(3)
array_zeros
Out[7]:
array([0., 0., 0.])
  • Matriz
In [8]:
shape = (4, 4)
matrix_zeros = np.zeros(shape)
matrix_zeros
Out[8]:
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

Analogamente, podemos fazer com valor unitário.

  • Vetor
In [9]:
array_ones = np.ones(3)
array_ones
Out[9]:
array([1., 1., 1.])
  • Matriz
In [10]:
shape = (4, 4)
matrix_ones = np.ones(shape)
matrix_ones
Out[10]:
array([[1., 1., 1., 1.],
       [1., 1., 1., 1.],
       [1., 1., 1., 1.],
       [1., 1., 1., 1.]])
Atributos básicos
  • Número de eixos (axes) do array:
In [11]:
array1.ndim
Out[11]:
1
In [12]:
matrix1.ndim
Out[12]:
2

Note que uma matriz, nesse caso, é considerado um objeto bidimensional, já que são arrays dentro de um array.

  • Forma de um array:
In [13]:
array1.shape
Out[13]:
(3,)
In [14]:
matrix1.shape
Out[14]:
(3, 3)

ndarray.shape fornece as dimensões de cada eixo. No caso de uma matriz $A \in \mathbb{R}^{m\,\times\,n}$, como temos em matrix1, ndarray.shape retorna justamente $(m, n)$.

  • Tamanho de um array:
In [15]:
array1.size
Out[15]:
3
In [16]:
matrix1.size
Out[16]:
9

Retorna o número de entrada de valores no objeto ndarray.

  • Tipo dos elementos em um array:
In [17]:
array1.dtype
Out[17]:
dtype('float64')
In [18]:
matrix1.dtype
Out[18]:
dtype('float64')

Ou seja, todas as entradas foram compostas por números em representação de ponto flutuante de dupla precisão. Caso seja necessário utilizar outro tipo de representação, podemos especificar isso na hora de declaramos o array por meio do argumento dtype:

In [19]:
matrix1_32 = np.array(
    [[1., 2., 3.],
     [4., 5., 6.],
     [7., 8., 9.]],
    dtype=np.float32
)
In [20]:
matrix1_32.dtype
Out[20]:
dtype('float32')

Ou mesmo podemos criar arrays de números complexos:

In [21]:
array_complex = np.array([[2. + 3.j, 1.]], dtype=np.complex)

array_complex
Out[21]:
array([[2.+3.j, 1.+0.j]])
In [22]:
matrix_complex = np.array(
    [[2. + 1.j, 1.],
     [0., 1.j]], 
    dtype=np.complex
)

matrix_complex
Out[22]:
array([[2.+1.j, 1.+0.j],
       [0.+0.j, 0.+1.j]])
  • Tamanhos dos itens em bytes de um array:
In [23]:
array1.itemsize
Out[23]:
8
In [24]:
matrix1.itemsize
Out[24]:
8
In [25]:
matrix1_32.itemsize
Out[25]:
4

Note que um array de float64 tem $\dfrac{64}{8} = 8$ bytes, enquanto um float32 tem $\dfrac{32}{8} = 4$ bytes.

  • Endereço dos dados em memória:
In [26]:
array1.data
Out[26]:
<memory at 0x7f6fa4667288>
In [27]:
matrix1.data
Out[27]:
<memory at 0x7f6fb8067630>
In [28]:
matrix1_32.data
Out[28]:
<memory at 0x7f6fb8067b40>

Em Python, geralmente não utilizamos esse tipo de acesso ou informação. Utilizamos apenas a parte boa do C.

Operações e utilitários básicos

Criação de arrays a partir de intervalos:

In [29]:
sequence_1 = np.arange(start=2, stop=10, step=2)

sequence_1
Out[29]:
array([2, 4, 6, 8])

ou

In [30]:
sequence_2 = np.linspace(start=2, stop=8, num=4)

sequence_2
Out[30]:
array([2., 4., 6., 8.])

Note que embora ambos criem arrays igualmente espaçados, exite uma diferença na forma de uso. O np.arange te possibilita dizer o início e o fim, sendo que o fim não está incluso (assim como ocorre com os slicing de listas). Além disso, o np.arange controla o incremento que será utilizado até alcançar o fim da sequência.

Em contrapartida, o np.linspace controla precisamente o início e o fim (que estará incluso). No entanto, ele tem como entrada no terceiro argumento (aqui explicitado) a quantidade de pontos que terá entre o início e o fim. O incremento será consequência disso.

Para simplificar, geralmente ambos os métodos são usados com o nome dos argumentos implícitos:

In [31]:
sequence_1 = np.arange(2, 10, 2)

sequence_1
Out[31]:
array([2, 4, 6, 8])
In [32]:
sequence_2 = np.linspace(2, 8, 4)

sequence_2
Out[32]:
array([2., 4., 6., 8.])
In [33]:
sequence_3 = np.linspace(0, 100, 200)

sequence_3
Out[33]:
array([  0.        ,   0.50251256,   1.00502513,   1.50753769,
         2.01005025,   2.51256281,   3.01507538,   3.51758794,
         4.0201005 ,   4.52261307,   5.02512563,   5.52763819,
         6.03015075,   6.53266332,   7.03517588,   7.53768844,
         8.04020101,   8.54271357,   9.04522613,   9.54773869,
        10.05025126,  10.55276382,  11.05527638,  11.55778894,
        12.06030151,  12.56281407,  13.06532663,  13.5678392 ,
        14.07035176,  14.57286432,  15.07537688,  15.57788945,
        16.08040201,  16.58291457,  17.08542714,  17.5879397 ,
        18.09045226,  18.59296482,  19.09547739,  19.59798995,
        20.10050251,  20.60301508,  21.10552764,  21.6080402 ,
        22.11055276,  22.61306533,  23.11557789,  23.61809045,
        24.12060302,  24.62311558,  25.12562814,  25.6281407 ,
        26.13065327,  26.63316583,  27.13567839,  27.63819095,
        28.14070352,  28.64321608,  29.14572864,  29.64824121,
        30.15075377,  30.65326633,  31.15577889,  31.65829146,
        32.16080402,  32.66331658,  33.16582915,  33.66834171,
        34.17085427,  34.67336683,  35.1758794 ,  35.67839196,
        36.18090452,  36.68341709,  37.18592965,  37.68844221,
        38.19095477,  38.69346734,  39.1959799 ,  39.69849246,
        40.20100503,  40.70351759,  41.20603015,  41.70854271,
        42.21105528,  42.71356784,  43.2160804 ,  43.71859296,
        44.22110553,  44.72361809,  45.22613065,  45.72864322,
        46.23115578,  46.73366834,  47.2361809 ,  47.73869347,
        48.24120603,  48.74371859,  49.24623116,  49.74874372,
        50.25125628,  50.75376884,  51.25628141,  51.75879397,
        52.26130653,  52.7638191 ,  53.26633166,  53.76884422,
        54.27135678,  54.77386935,  55.27638191,  55.77889447,
        56.28140704,  56.7839196 ,  57.28643216,  57.78894472,
        58.29145729,  58.79396985,  59.29648241,  59.79899497,
        60.30150754,  60.8040201 ,  61.30653266,  61.80904523,
        62.31155779,  62.81407035,  63.31658291,  63.81909548,
        64.32160804,  64.8241206 ,  65.32663317,  65.82914573,
        66.33165829,  66.83417085,  67.33668342,  67.83919598,
        68.34170854,  68.84422111,  69.34673367,  69.84924623,
        70.35175879,  70.85427136,  71.35678392,  71.85929648,
        72.36180905,  72.86432161,  73.36683417,  73.86934673,
        74.3718593 ,  74.87437186,  75.37688442,  75.87939698,
        76.38190955,  76.88442211,  77.38693467,  77.88944724,
        78.3919598 ,  78.89447236,  79.39698492,  79.89949749,
        80.40201005,  80.90452261,  81.40703518,  81.90954774,
        82.4120603 ,  82.91457286,  83.41708543,  83.91959799,
        84.42211055,  84.92462312,  85.42713568,  85.92964824,
        86.4321608 ,  86.93467337,  87.43718593,  87.93969849,
        88.44221106,  88.94472362,  89.44723618,  89.94974874,
        90.45226131,  90.95477387,  91.45728643,  91.95979899,
        92.46231156,  92.96482412,  93.46733668,  93.96984925,
        94.47236181,  94.97487437,  95.47738693,  95.9798995 ,
        96.48241206,  96.98492462,  97.48743719,  97.98994975,
        98.49246231,  98.99497487,  99.49748744, 100.        ])

Operações aritméticas

Com o numpy, as operações aritméticas são realizadas elemento-por-elemento do array.

In [34]:
a = np.arange(5)

a
Out[34]:
array([0, 1, 2, 3, 4])

(note que se np.arange tiver apenas um argumento, ele gerará uma sequência igualmente espaçada até o número inteiro imediatamente antes ao argumento)

In [35]:
b = np.arange(0, 50, 10)

b
Out[35]:
array([ 0, 10, 20, 30, 40])
In [36]:
soma = a + b
subtracao = b - a
potenciacao = a ** 2
divisao = b[1:] / a[1:]
produto = a * b

print(soma)
print(subtracao)
print(potenciacao)
print(divisao)
print(produto)
[ 0 11 22 33 44]
[ 0  9 18 27 36]
[ 0  1  4  9 16]
[10. 10. 10. 10.]
[  0  10  40  90 160]

Note que a operação de divisão de vetores não é matematicamente definida, mas aqui é possível pois se trata de operações elemento-por-elemento. Além disso, o produto entre vetores matematicamente definido como produto interno não ocorre por meio do operador * (iremos ver isso em breve). Podemos fazer com matrizes de maneira similar:

In [37]:
A = np.array(
    [[1., 3.], 
     [2., 4.]]
)
B = np.array(
    [[10., 0.], 
     [10., 20.]]
)

print(A, '\n')
print(B)
[[1. 3.]
 [2. 4.]] 

[[10.  0.]
 [10. 20.]]
In [38]:
soma = A + B
subtracao = B - A
potenciacao = A ** 2
divisao = B / A
produto = A * B

print(soma, '\n')
print(subtracao, '\n')
print(potenciacao, '\n')
print(divisao, '\n')
print(produto)
[[11.  3.]
 [12. 24.]] 

[[ 9. -3.]
 [ 8. 16.]] 

[[ 1.  9.]
 [ 4. 16.]] 

[[10.  0.]
 [ 5.  5.]] 

[[10.  0.]
 [20. 80.]]

Novamente, note que a operação de divisão em matrizes matematicamente não existe e o produto entre matrizes também não é o usual. As operações que ocorrem acima foram realizadas elemento-por-elemento. Caso desejemos realizar o produto de matrizes, o Python oferece um operador específico para isso:

In [39]:
produto_matrizes = A @ B

produto_matrizes
Out[39]:
array([[40., 60.],
       [60., 80.]])

ou

In [40]:
produto_matrizes_alternativo = A.dot(B)

produto_matrizes_alternativo
Out[40]:
array([[40., 60.],
       [60., 80.]])

sendo esse último uma notação mais verbosa do que está ocorrendo. Podemos, dessa forma, fazer a clássica operação entre matrizes e vetores do tipo $Ax$:

In [41]:
mat = np.array([[2., 1.], [1, 2.]])
vec = np.array([1., 0.])
prod_matvec = mat @ vec

prod_matvec
Out[41]:
array([2., 1.])

Comentários sobre operações com tipos de dados distintos

Quando se opera objetos ndarray com tipos diferentes, o que será que acontece?

In [42]:
a = np.ones(3, dtype=np.int32)

a
Out[42]:
array([1, 1, 1], dtype=int32)
In [43]:
b = np.linspace(0, 5, 3)

b
Out[43]:
array([0. , 2.5, 5. ])
In [44]:
a + b
Out[44]:
array([1. , 3.5, 6. ])

Note que ocorre uma coerção de tipo para o que apresentar maior grau de precisão (ou mais genérico).

Funções universais da NumPy

A NumPy fornece as principais funções matemáticas capazes de operar de forma vetorial. Em outras palavras, capaz de retornar arrays, ao contrário da biblioteca-padrão math. Vamos ver um exemplo simples:

In [45]:
import math


def sin_function(x):
    return math.sin(x)
In [46]:
vec = np.array([1 / 3, 1 / 2])
sin_vec = sin_function(vec)

sin_vec
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-46-f2b7933f87d0> in <module>
      1 vec = np.array([1 / 3, 1 / 2])
----> 2 sin_vec = sin_function(vec)
      3 
      4 sin_vec

<ipython-input-45-f0da795297e3> in sin_function(x)
      3 
      4 def sin_function(x):
----> 5     return math.sin(x)

TypeError: only size-1 arrays can be converted to Python scalars

Pela numpy:

In [47]:
def sin_function_np(x):
    return np.sin(x)
In [48]:
sin_vec = sin_function_np(vec)

sin_vec
Out[48]:
array([0.3271947 , 0.47942554])

Basicamente, as funções fornecidas são as mesmas que as da lib math, porém com a capacidade exibida acima. Lembrando que essas funções também são executadas elemento-por-elemento. A lista completa de funções disponíveis pode ser acessada aqui.

Alguns métodos de numpy.ndarray úteis

No dia-a-dia científico, muitas vezes necessitamos de algumas operações com certa frequência, como somas em vetores, valores máximos ou mínimos e etc. Os objetos ndarray são munidos de métodos que executam esses cálculos. Vejamos alguns deles a seguir.

Valores máximos e mínimos em um array

Vamos antes gerar um vetor de valores aleatórios com o pacote random do numpy:

In [49]:
from numpy import random as rnd

array1 = rnd.random(10)

array1
Out[49]:
array([0.34336861, 0.27182526, 0.13845444, 0.90895815, 0.02422253,
       0.42396994, 0.74019314, 0.9223146 , 0.50825702, 0.16044042])

Criamos um vetor com 10 entradas aleatórias entre 0 e 1. Agora vamos retirar algumas informações sobre ele.

  • Valor máximo
In [50]:
array1.max()
Out[50]:
0.9223146022241675
  • Índice da entrada de valor máximo
In [51]:
imax = array1.argmax()
imax
Out[51]:
7
In [52]:
array1[imax]
Out[52]:
0.9223146022241675
  • Valor mínimo
In [53]:
array1.min()
Out[53]:
0.024222528975309188
  • Índice da entrada de valor mínimo
In [54]:
array1.argmin()
Out[54]:
4

Somatório e produtórios dos valores em um array

  • Soma dos valores
In [55]:
array1.sum()
Out[55]:
4.442004105570682
  • Soma cumulativa
In [56]:
array1.cumsum()
Out[56]:
array([0.34336861, 0.61519386, 0.7536483 , 1.66260645, 1.68682898,
       2.11079892, 2.85099206, 3.77330666, 4.28156368, 4.44200411])

Note que:

In [57]:
np.allclose(array1.sum(), array1.cumsum()[-1])
Out[57]:
True
  • Produto dos valores
In [58]:
array1.prod()
Out[58]:
6.715479688671335e-06
  • Produto cumulativo dos valores
In [59]:
array1.cumprod()
Out[59]:
array([3.43368606e-01, 9.33362595e-02, 1.29228197e-02, 1.17463022e-02,
       2.84525147e-04, 1.20630108e-04, 8.92895791e-05, 8.23530826e-05,
       4.18565325e-05, 6.71547969e-06])

Também é constatado que:

In [60]:
array1.prod() == array1.cumprod()[-1]
Out[60]:
True

E o que acontece em arrays multidimensionais?

Neste caso, o formato não é levado em conta. O método é aplicado ao longo do array como se os valores de entradas fossem uma lista.

In [61]:
shape = (3, 3)
vec1 = rnd.random(shape)

vec1
Out[61]:
array([[0.52757469, 0.44499104, 0.24150775],
       [0.17061932, 0.23133227, 0.48625264],
       [0.79868698, 0.93480275, 0.93662729]])
In [62]:
vec1.min()
Out[62]:
0.1706193224360818

Porém, podemos especificar em qual axis os métodos serão aplicados:

In [63]:
vec1.min(axis=1)
Out[63]:
array([0.24150775, 0.17061932, 0.79868698])

A operação acima retornou o mínimo de cada linha. Analogamente:

In [64]:
vec1.min(axis=0)
Out[64]:
array([0.17061932, 0.23133227, 0.24150775])

retornou o mínimo em cada coluna. O procedimento pode ser aplicado da mesma maneira para as funções expostas anteriormente.

Indexação, fatiamento e iterações em arrays

Indexação

In [65]:
array1 = np.linspace(5, 20, 10)
array1
Out[65]:
array([ 5.        ,  6.66666667,  8.33333333, 10.        , 11.66666667,
       13.33333333, 15.        , 16.66666667, 18.33333333, 20.        ])

Semelhantemente às listas, os índices começam a contar do 0. Portanto, a entrada na posição 3 é acessado como

In [66]:
array1[2]
Out[66]:
8.333333333333334

Fatiamento

Array unidimensional

Em algumas situações, se deseja acessar ou operar em apenas uma parte de um array. No array unidimensional, basta fazer array[n:m]. Isso retornará as entradas do índice n até o m-1.

In [67]:
array1[0:3]
Out[67]:
array([5.        , 6.66666667, 8.33333333])

Neste caso, como o índice inicial é 0, podemos omitir e simplesmente fazer

In [68]:
array1[:3]
Out[68]:
array([5.        , 6.66666667, 8.33333333])

De forma semelhante, se fizermos array[n:] a fatia ocorrerá do índice n até o fim do array:

In [69]:
array1[7:]
Out[69]:
array([16.66666667, 18.33333333, 20.        ])

Outro exemplo:

In [70]:
sliced_array = array1[2:8]
sliced_array
Out[70]:
array([ 8.33333333, 10.        , 11.66666667, 13.33333333, 15.        ,
       16.66666667])

Também podemos fazer operações em apenas trechos dos arrays:

In [71]:
array1[:3] + 10
Out[71]:
array([15.        , 16.66666667, 18.33333333])

Ou operar partes de arrays com outras partes de outras arrays (sem fazer laços!), desde que sejam compatíveis em tamanho e dimensões:

In [72]:
array2 = np.linspace(0, 5, 10)
array3 = array1[:5] + array2[5:]

array3
Out[72]:
array([ 7.77777778, 10.        , 12.22222222, 14.44444444, 16.66666667])

Esse tipo de operação entre arrays, sem laços, é conhecida como uma das coisas mais poderosas da NumPy, e é comumente chamado de vetorização. Vejamos que alguns exemplos de performance:

In [73]:
array1, array2 = rnd.random(1000), rnd.random(1000)
n = 800
In [74]:
%%timeit
array3 = np.zeros(n)
for i in range(n):
    array3[i] = array1[i] + array2[i]
    
array3
661 µs ± 11.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [75]:
%%timeit
array3 = array1[:n] + array2[:n]

array3
1.95 µs ± 393 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Então podemos verificar que há um ganho de performance/complexidade temporal na ordem de 300 com a operação vetorizada.

Array multidimensional

Neste caso, o fatiamento é aplicado em cada eixo (axis) separadamente.

In [76]:
shape = (4, 4)
matrix1 = rnd.randint(5, size=shape)

print(matrix1)
[[1 2 0 4]
 [3 1 4 4]
 [1 1 1 0]
 [2 1 2 2]]
  • Todas as linhas (:) da coluna 2:
In [77]:
matrix1[:, 1]
Out[77]:
array([2, 1, 1, 1])
  • Todas as colunas da linha 4:
In [78]:
matrix1[3, :]
Out[78]:
array([2, 1, 2, 2])
  • Uma matriz formada das linhas 1 até 3 e com as colunas de 2 e 3 matrix1:
In [79]:
matrix2 = matrix1[:3, 1:3]

matrix2
Out[79]:
array([[2, 0],
       [1, 4],
       [1, 1]])

Agora vamos fazer também um comparativo de performance.

In [80]:
shape = (200, 200)
matrix1 = rnd.randint(5, size=shape)
matrix2 = rnd.randint(5, size=shape)
matrix3, matrix3_loop = np.zeros(shape), np.zeros(shape)
In [81]:
%%timeit
for j in range(matrix1.shape[0]):
    for i in range(matrix1.shape[1]):
        matrix3_loop[i, j] = matrix1[i, j] + matrix2[i, j]
23.2 ms ± 477 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [82]:
%%timeit
matrix3[:, :] = matrix1[:, :] + matrix2[:, :]
43.5 µs ± 586 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Apenas fazendo um simples teste de sanidade:

In [83]:
np.array_equal(matrix3, matrix3_loop)
Out[83]:
True

Disso, vemos que as operações foram equivalentes, porém a forma vetorizada teve um ganho em performance/complexidade temporal praticamente de $\mathcal{O}(10^3)$.

Iterações

De forma análoga às listas em Python, podemos iterar sobre objetos numpy.ndarray.

In [84]:
matrix1 = rnd.randint(10, size=(3, 3))
matrix1
Out[84]:
array([[1, 3, 0],
       [9, 6, 3],
       [4, 9, 9]])
In [85]:
for row in matrix1:
    print(row)
[1 3 0]
[9 6 3]
[4 9 9]

Também é válido para percorrer elementos em vetores (arrays unidimensionais):

In [86]:
array1 = rnd.randint(5, size=5)
array1
Out[86]:
array([0, 2, 3, 0, 4])
In [87]:
for element in array1:
    print(element)
0
2
3
0
4

Manipulação de shapes

Em situações práticas, pode ser que seja necessário um rearranjo da configuração de um array. Aqui veremos brevemente como podemos mudar os shapes dos arrays.

In [88]:
shape = (3, 4)
matrix1 = rnd.randint(10, size=shape)

matrix1
Out[88]:
array([[3, 7, 8, 0],
       [1, 2, 7, 7],
       [7, 7, 5, 4]])
In [89]:
matrix1.shape
Out[89]:
(3, 4)

Como já vimos, o atributo shape nos informa quantos elementos há em cada eixo do array multidimensional. Inicialmente, vamos transformar a matrix1 em um array unidimensional:

In [90]:
matrix1.ravel()
Out[90]:
array([3, 7, 8, 0, 1, 2, 7, 7, 7, 7, 5, 4])

Podemos também rearranjar o conteúdo de um array para um outro shape:

In [91]:
new_shape = (2, 6)
matrix1.reshape(new_shape)
Out[91]:
array([[3, 7, 8, 0, 1, 2],
       [7, 7, 7, 7, 5, 4]])

Note que é necessário que shape[0] * shape[1] = new_shape[0] * new_shape[1]. Mais precisamente, o tamanho do array rearranjado deve ser igual ao do array original. Veja o erro abaixo:

In [92]:
new_shape = (3, 3)
matrix1.reshape(new_shape)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-92-e395430efd40> in <module>
      1 new_shape = (3, 3)
----> 2 matrix1.reshape(new_shape)

ValueError: cannot reshape array of size 12 into shape (3,3)

É importante observar que o reshape não altera o array original, apenas devolve um array rearranjado em relação ao original:

In [93]:
matrix1
Out[93]:
array([[3, 7, 8, 0],
       [1, 2, 7, 7],
       [7, 7, 5, 4]])

Portanto, matrix1 se mantém inalterada. Isso não ocorre quando se utiliza o resize, que modifica também o array original:

In [94]:
new_shape = (2, 6)
matrix1.resize(new_shape)

matrix1
Out[94]:
array([[3, 7, 8, 0, 1, 2],
       [7, 7, 7, 7, 5, 4]])

Empilhamento de arrays

Digamos que precisamos unir os dados presentes em dois arrays diferentes. Aqui veremos sucintamente algumas maneiras de fazer isso.

In [95]:
shape = (2, 2)
matrix1, matrix2 = rnd.randint(5, size=shape), rnd.randint(4, size=shape)

print(matrix1, '\n')
print(matrix2)
[[1 3]
 [2 1]] 

[[1 3]
 [2 3]]
  • Empilhamento horizontal

Os arrays são concatenados horizontalmente:

In [96]:
np.hstack([matrix1, matrix2])
Out[96]:
array([[1, 3, 1, 3],
       [2, 1, 2, 3]])
  • Empilhamento vertical

Analogamente:

In [97]:
np.vstack([matrix1, matrix2])
Out[97]:
array([[1, 3],
       [2, 1],
       [1, 3],
       [2, 3]])
  • Empilhamento de um vetor como coluna em uma matriz
In [98]:
array1 = np.array([10., 10.])

np.hstack([matrix1, array1])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-98-9bfd43ac53dc> in <module>
      1 array1 = np.array([10., 10.])
      2 
----> 3 np.hstack([matrix1, array1])

~/miniconda3/envs/sci_env/lib/python3.7/site-packages/numpy/core/shape_base.py in hstack(tup)
    338         return _nx.concatenate(arrs, 0)
    339     else:
--> 340         return _nx.concatenate(arrs, 1)
    341 
    342 

ValueError: all the input arrays must have same number of dimensions

Observe que isso não é possível com o hstack. Porém, temos o column_stack para nos salvar!

In [99]:
np.column_stack([matrix1, array1])
Out[99]:
array([[ 1.,  3., 10.],
       [ 2.,  1., 10.]])

Também podemos utilizar o column_stack para montar um array multidimensional a partir de arrays que irão compor as colunas:

In [100]:
array2 = np.array([15, 30])

np.column_stack([array1, array2])
Out[100]:
array([[10., 15.],
       [10., 30.]])

Analogamente, podemos concatenar como linhas:

In [101]:
np.row_stack([array1, array2])
Out[101]:
array([[10., 10.],
       [15., 30.]])

Em tempo, podemos adicionar dois arrays com a função append, que difere um pouco do seu primo nos objetos list:

In [102]:
np.append(array1, array2)
Out[102]:
array([10., 10., 15., 30.])

Isso também pode ser usado no numpy para criar arrays "dinamicamente" (que não existe no numpy, exatamente):

In [103]:
an_array = np.array([])  # cria um array vazio
for i in range(10):
    an_array = np.append(an_array, i)
    
an_array
Out[103]:
array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])

A grosso modo, a função append possui a seguinte sintaxe: np.append(array, value), onde array será adicionado por value, que pode ser também um outro ndarray ou tipos compatíveis, como list ou tuple. O append não modifica o array, retorna um novo ndarray com operado como requerido na chamada da função.

Observação: Criar arrays com append pode ser custoso! Afinal, usa laços em Python.

Cópias e visões de arrays

Durante criação e manipulações de arrays, às vezes cópias são realizadas ou às vezes temos apenas visões dos arrays. Com numpy, temos alguns tipos de cópias e isso pode ser crucial quando desejamos manipular dados apropriadamente ou mantê-los intactos. Basicamente, temos apenas três casos e os veremos a seguir.

Sem cópias

Em geral, declarações não fazem cópias dos arrays.

In [104]:
array1 = np.arange(0, 10, 2)

array1
Out[104]:
array([0, 2, 4, 6, 8])
In [105]:
array2 = array1

array1 is array2
Out[105]:
True

Então os arrays são iguais. Em Python, isso quer dizer que um mesmo objeto possui dois nomes (ou mais). O que acontece se um dos arrays for modificado?

In [106]:
array2[-1] = 10

array1
Out[106]:
array([ 0,  2,  4,  6, 10])

É importante ter muita atenção nisso, pois é uma fonte de bugs frequente. Funções em Python também modificam o valor dos argumentos:

In [107]:
def f_bug(x):
    x[0] = 20.0
    print(x)
In [108]:
f_bug(array1)
[20  2  4  6 10]
In [109]:
array1
Out[109]:
array([20,  2,  4,  6, 10])

Isso é conhecido como passagem por referência de um argumento de uma função.

Visão ou cópia raza (shallow copy)

Neste caso, diferentes objetos compartilham os mesmos dados. O método view cria um novo array com dados iguais aos de outro array.

In [110]:
array2 = array1.view()

array2
Out[110]:
array([20,  2,  4,  6, 10])
In [111]:
array2 is array1
Out[111]:
False
In [112]:
array2.base is array1
Out[112]:
True

O que acontece com o array1 se modificarmos o array2?

In [113]:
array2 = array2 + 5

array1
Out[113]:
array([20,  2,  4,  6, 10])

Aparentemente tudo OK! Esse tipo de cópia ocorre nos fatiamentos, onde são geradas as "visões" dos dados para um novo array:

In [114]:
array3 = array1[1:4]

array3
Out[114]:
array([2, 4, 6])
In [115]:
array3[1] = 40

array1
Out[115]:
array([20,  2, 40,  6, 10])

Então também devemos tomar cuidado com isso!

Cópia profunda (deep copy)

É uma perfeita cópia, totalmente independente dos dados originais (e segura!). É feita por meio do método copy().

In [116]:
array3 = array1.copy()
In [117]:
array3 is array1
Out[117]:
False
In [118]:
array3.base is array1
Out[118]:
False
In [119]:
array3[:3] = 0.0

array1
Out[119]:
array([20,  2, 40,  6, 10])

E assim estamos totalmente a salvo!

Aplicações básicas em Álgebra Linear computacional

Nesta seção, iremos ver algumas das funcionalidades voltadas para Álgebra Linear computacional em alguns problemas frequentes. Para consultas e detalhes das funcionalidades disponibilizadas abaixo, a documentação oficial.

A NumPy fornece apenas algumas funcionalidades básicas. Métodos mais complexos e eficientes podem ser utilizados no SciPy, como veremos mais na frente do curso. Na NumPy, as funções de Álgebra Linear estão no sub-pacote numpy.linalg.

In [120]:
import numpy.linalg as la

Operações básicas

Vetores e matrizes para testes: $u, v \in \mathbb{R}^n$ e $A, B \in \mathbb{R}^{n\,\times\,n}$, $n = 3$.

In [121]:
n = 3
shape = (n, n)
u, v = rnd.random(n), rnd.random(n)
A, B = rnd.random(size=shape), rnd.random(size=shape)

print('u:\n', u, '\n')
print('v:\n', v, '\n')
print('A:\n', A, '\n')
print('B:\n', B)
u:
 [0.57012495 0.27262238 0.71554926] 

v:
 [0.14246231 0.76020427 0.00672093] 

A:
 [[0.36349565 0.62509516 0.40577012]
 [0.71787282 0.91853055 0.71335495]
 [0.55643619 0.10764589 0.7500328 ]] 

B:
 [[0.52474217 0.3235028  0.31363832]
 [0.24603614 0.40909994 0.14241992]
 [0.22472984 0.92455571 0.45881707]]

Transposta: $A^T$

In [122]:
A.T
Out[122]:
array([[0.36349565, 0.71787282, 0.55643619],
       [0.62509516, 0.91853055, 0.10764589],
       [0.40577012, 0.71335495, 0.7500328 ]])

ou

In [123]:
np.transpose(A)
Out[123]:
array([[0.36349565, 0.71787282, 0.55643619],
       [0.62509516, 0.91853055, 0.10764589],
       [0.40577012, 0.71335495, 0.7500328 ]])

Produto interno

$u \cdot v$

In [124]:
np.dot(u, v)
Out[124]:
0.29327917230988837

$A \cdot u \equiv Au$

In [125]:
np.dot(A, u)
Out[125]:
array([0.66800138, 1.17012979, 0.88327025])

$A \cdot B \equiv AB$

In [126]:
np.dot(A, B)
Out[126]:
array([[0.43572615, 0.74847533, 0.38920642],
       [0.76300199, 1.26754104, 0.68326889],
       [0.48702506, 0.9174937 , 0.53397848]])

ou

In [127]:
A @ B
Out[127]:
array([[0.43572615, 0.74847533, 0.38920642],
       [0.76300199, 1.26754104, 0.68326889],
       [0.48702506, 0.9174937 , 0.53397848]])

Produto diádico

$u \otimes v$

In [128]:
np.outer(u, v)
Out[128]:
array([[0.08122132, 0.43341142, 0.00383177],
       [0.03883841, 0.2072487 , 0.00183228],
       [0.1019388 , 0.54396361, 0.00480915]])

$A \otimes B$

In [129]:
np.outer(A, B)
Out[129]:
array([[0.1907415 , 0.11759186, 0.11400616, 0.08943307, 0.14870605,
        0.05176902, 0.08168832, 0.33607197, 0.16677801],
       [0.32801379, 0.20222003, 0.1960538 , 0.153796  , 0.25572639,
        0.089026  , 0.14047753, 0.5779353 , 0.28680433],
       [0.2129247 , 0.13126777, 0.12726506, 0.09983411, 0.16600053,
        0.05778975, 0.09118865, 0.37515708, 0.18617426],
       [0.37669814, 0.23223386, 0.22515242, 0.17662266, 0.29368172,
        0.10223939, 0.16132744, 0.66371341, 0.3293723 ],
       [0.48199172, 0.2971472 , 0.28808638, 0.22599171, 0.37577079,
        0.13081704, 0.20642122, 0.84923266, 0.42143749],
       [0.37432743, 0.23077232, 0.22373545, 0.1755111 , 0.29183346,
        0.10159595, 0.16031214, 0.65953639, 0.32729943],
       [0.29198553, 0.18000866, 0.17451971, 0.13690341, 0.22763801,
        0.07924759, 0.12504781, 0.51445625, 0.25530242],
       [0.05648634, 0.03482375, 0.03376188, 0.02648478, 0.04403793,
        0.01533092, 0.02419124, 0.09952462, 0.04938977],
       [0.39357384, 0.24263771, 0.23523903, 0.18453518, 0.30683837,
        0.10681961, 0.16855475, 0.69344711, 0.34412785]])

Esse é o produto diádico generalizado. Resulta em matriz de blocos $3 \times 3$.

Produto cruzado

$u \times v$

In [130]:
np.cross(u, v)
Out[130]:
array([-0.54213133,  0.09810703,  0.39457301])

Produto de contração tensorial

$A:B$

In [131]:
np.tensordot(A, B)
Out[131]:
array(1.74291627)

Característica de Matrizes e condicionamento

Traço

In [132]:
np.trace(A)
Out[132]:
2.032058999170806

Determinante

In [133]:
la.det(A)
Out[133]:
-0.041970404338626255

O determinante ser diferente de zero é uma característica muito desejada quando se tem que resolver um sistema. Ele garante que o sistema é invertível, e portanto tenha solução única.

Número de condicionamento

In [134]:
la.cond(A)
Out[134]:
39.59718654379296

O número de condicionamento nos informa a sensibilidade numérica do problema/função/sistema. Um sistema com número de condicionamento baixo não é sensível a erros numéricos que podem se propagar nas manipulações e iterações. Ou seja, é um sistema fácil de se resolver. O problema mal-condicionado é justamente o inverso. São problemas sensíveis e que podem convergir para soluções inadequadas, com muito erro, ou mesmo nem convergir (em métodos iterativos de sistemas). Em geral, dizemos que um sistema com número de condicionamento $\kappa(A) \le 1$ (ou próximo de um) é bem condicionado.

Posto de uma matriz

In [135]:
la.matrix_rank(A)
Out[135]:
3

O posto de uma matriz também é um número muito importante, pois informa quantas linhas linearmente independentes a matriz possui. Podemos usar essa informação para comparar com o número de linhas da matriz. Se for igual, então isso é um indicativo que o sistema é determinado.

Norma

In [136]:
la.norm(A)
Out[136]:
1.8545068183606386

A função la.norm calcula a norma de uma matriz ou vetor. Pode ser usada genericamente como la.norm(A, ord=n), sendo n a ordem da norma. No caso da norma Euclidiana, n = 2, por exemplo. A norma do máximo é dada por:

In [137]:
la.norm(A, np.inf)
Out[137]:
2.3497583103569863

ou, em um vetor:

In [138]:
la.norm(u, np.inf)
Out[138]:
0.7155492635807369

que é equivalente a

In [139]:
u.max()
Out[139]:
0.7155492635807369

Autovalores e autovetores

São fornecidos por uma única função:

In [140]:
e_vals, e_vecs = la.eig(A)

e_vals, e_vecs
Out[140]:
(array([ 1.69680194, -0.06222866,  0.39748572]),
 array([[ 0.48460418,  0.81912549, -0.38439744],
        [ 0.79042533, -0.21184004, -0.51725859],
        [ 0.37468199, -0.53306401,  0.76464512]]))

sendo cada autovalor dado por e_vals[i] e seu correspondente autovetor dado por e_vecs[:, i].

Solução de sistemas lineares

As duas funções principais aqui são o solucionador de sistemas e o cálculo da inversa.

Inversa de uma matriz

In [141]:
la.inv(A)
Out[141]:
array([[-14.58499921,  10.13005231,  -1.74414509],
       [  3.37122449,  -1.11622659,  -0.76220188],
       [ 10.33651171,  -7.35510558,   2.73661809]])

lembrando que isso só é possível já que

In [142]:
la.det(A) != 0
Out[142]:
True

Também podemos facilmente verificar que

In [143]:
np.set_printoptions(5)
la.inv(A) @ A
Out[143]:
array([[ 1.00000e+00, -1.09228e-15,  3.18047e-16],
       [ 2.63890e-16,  1.00000e+00,  9.63232e-17],
       [ 9.89491e-16,  6.48609e-16,  1.00000e+00]])
In [144]:
A @ la.inv(A)
Out[144]:
array([[ 1.00000e+00,  7.18267e-16, -2.22578e-16],
       [-6.17238e-17,  1.00000e+00, -4.98952e-17],
       [ 1.46013e-16, -1.39849e-16,  1.00000e+00]])

Solução de sistemas lineares

Vamos criar um vetor do lado direito com valores fáceis de checarmos a qualidade da solução:

In [145]:
b = np.ones(n)
b
Out[145]:
array([1., 1., 1.])

Agora vamos resolver o sistema $A u = b$:

In [146]:
x = la.solve(A, b)
x
Out[146]:
array([-6.19909,  1.4928 ,  5.71802])

Podemos verificar que:

In [147]:
np.dot(A, x)
Out[147]:
array([1., 1., 1.])

Por baixo dos panos, o NumPy utiliza a _gesv do LAPACK, que nada mais é do que uma fatoração LU. Pró: é um método direto, portanto fornece uma solução "exata", sujeita apenas à precisão da representação numérica empregada. Contra: Pode ser bastante custoso em sistemas grandes.

Aplicações básicas em Estatística

A NumPy tem uma grande variedades de funções para Estatística. Não teremos tempo para detalhar aqui neste curso. Esse material demandaria um curso próprio de Análise de Dados com NumPy ou com pandas. Quem tiver interesse, pode verificar as funções disponíveis na documentação oficial online ou também na página de rotinas estatísticas básicas.

Aqui, veremos só o básico do básico inerente apenas a entrada de valores em vetores e matrizes.

Média

In [148]:
u.mean()
Out[148]:
0.51943219773322

Desvio-padrão

In [149]:
u.std()
Out[149]:
0.184342746224072

Variância

In [150]:
u.var()
Out[150]:
0.03398224808543262

que naturalmente é equivalente a

In [151]:
u.std() ** 2
Out[151]:
0.03398224808543262

Escrita e Leitura

O NumPy oferece rotinas convenientes para leitura e escrita de dados numéricos. Basicamente podemos sobreviver com loadtxt() e seu dual savetxt().

Escrita

Vamos começar com a escrita, assim poderemos usar o resultado para testar a leitura depois.

Escrita de um único vetor

In [152]:
np.savetxt('u_out.dat', u)

Vamos agora verificar como ficou o arquivo:

In [153]:
!cat u_out.dat
5.701249488061702619e-01
2.726223808127529180e-01
7.155492635807368940e-01

Podemos ainda mudar o formato da escrita:

In [154]:
np.savetxt('u_out_2.dat', u, fmt='%1.3e')
In [155]:
!cat u_out_2.dat
5.701e-01
2.726e-01
7.155e-01

Escrita de mais de um vetor

In [156]:
np.savetxt('uv_out.dat', (u, v), fmt='%1.4e')
In [157]:
!cat uv_out.dat
5.7012e-01 2.7262e-01 7.1555e-01
1.4246e-01 7.6020e-01 6.7209e-03

Note que a escrita é feita por linhas. Em cada linha do arquivo de saída, temos um dos vetores. Caso queira salvar cada vetor em uma coluna, basta fazer:

In [158]:
np.savetxt('uv_out_2.dat', np.column_stack((u, v)), fmt='%1.4e')
In [159]:
!cat uv_out_2.dat
5.7012e-01 1.4246e-01
2.7262e-01 7.6020e-01
7.1555e-01 6.7209e-03

Escrita de uma matriz

Uma matriz nada mais é do que um aglomerado de vetores colunas, então é meio intuitivo como o savetxt() irá agir.

In [160]:
np.savetxt('A_out.dat', A, fmt='%1.3e')
In [161]:
!cat A_out.dat
3.635e-01 6.251e-01 4.058e-01
7.179e-01 9.185e-01 7.134e-01
5.564e-01 1.076e-01 7.500e-01

Só para verificação:

In [162]:
A
Out[162]:
array([[0.3635 , 0.6251 , 0.40577],
       [0.71787, 0.91853, 0.71335],
       [0.55644, 0.10765, 0.75003]])

Podemos mudar o delimitador simplesmente com

In [163]:
np.savetxt('A_out_2.dat', A, fmt='%1.3e', delimiter=',')
In [164]:
!cat A_out_2.dat
3.635e-01,6.251e-01,4.058e-01
7.179e-01,9.185e-01,7.134e-01
5.564e-01,1.076e-01,7.500e-01

que pode vir a calhar para determinados leitores csv-like.

Leitura

Leitura simples de um vetor em arquivo

In [165]:
x = np.loadtxt('u_out.dat')
x
Out[165]:
array([0.57012, 0.27262, 0.71555])

Leitura de 2 ou mais vetores por linhas

Aqui apenas com 2 para exemplificar. A extensão para mais é trivial.

In [166]:
x, y = np.loadtxt('uv_out.dat')

x, y
Out[166]:
(array([0.57012, 0.27262, 0.71555]), array([0.14246, 0.7602 , 0.00672]))

A leitura aqui foi feita com um array por linha do arquivo.

Leitura de 2 ou mais vetores por coluna

In [167]:
x, y = np.loadtxt('uv_out_2.dat', unpack=True)

x, y
Out[167]:
(array([0.57012, 0.27262, 0.71555]), array([0.14246, 0.7602 , 0.00672]))

Leitura de partes de uma matriz em arrays

Caso se deseje ler dados em formato matricial de apenas determinadas colunas, basta especificar isso por meio do argumento opcional usecols.

In [168]:
x, y = np.loadtxt('A_out.dat', usecols=(0, 1), unpack=True)

x, y
Out[168]:
(array([0.3635, 0.7179, 0.5564]), array([0.6251, 0.9185, 0.1076]))

Verificando:

In [169]:
A
Out[169]:
array([[0.3635 , 0.6251 , 0.40577],
       [0.71787, 0.91853, 0.71335],
       [0.55644, 0.10765, 0.75003]])

Leitura de uma matriz

Basta simplesmente fazer

In [170]:
M = np.loadtxt('A_out.dat')

M
Out[170]:
array([[0.3635, 0.6251, 0.4058],
       [0.7179, 0.9185, 0.7134],
       [0.5564, 0.1076, 0.75  ]])

Próximos passos

Acima, vimos apenas as funcionalidades básicas da NumPy e que considero as mais frequentes. Há ainda um mundo todo para explorar! Quer conhecer mais? Sugiro olhar a documentação sempre que precisar realizar alguma tarefa específica.