Física Computacional - FSC-5705

só um divisor

Python para quem sabe programar

Porque Python?

Pergunta obvia, porque aprender essa linguagem se eu já sei outra?. A reposta bem simples e rápida, por que, em ciências a gente precisas analisar os dados obtidos de um programa, experimento, teste, etc, e aqui onde o python pode ser muito útil para todos. O python tem uma comunidade gigantesca de usuários, muitos deles da área científica e eles tem desenvolvido módulos para o python para realizar a suas tarefas de forma mais eficiente, assim a ideia é pegar carona nos módulos já existentes. Em esse sentido o python é muito parecido ao java (de fato, é orientado a objeto igual que o java) ja que possui uma enorme quantidade de classes já definidas, a diferencia é que o java pertencia a SUN (que foi comprada pela Oracle e está largando o java) e o python é mantido pela comunidade. Mas para chegar a isso vamos ver as nuances dessa linguagem.

Primeiro Script

Python é uma linguagem de script, portanto ela interpretada e no compilada, ou seja, quando ela está sendo executa é que acontece a conversão para código de maquina, isso faz dela que seja muito lenta assim não espere fazer dela uma linguagem para tratar cálculos complexos, mas dependendo da área da física pode ser mais útil para o tratamento de dados do que C ou Fortran.

primeiro script
Figura 01: executando um script python.

Na figura 01 se mostra um tipo script, ele inicia no cabeçalho informando qual é a linguagem, no caso #!/usr/bin/env python. O primeiro que você deve saber é que ele é um arquivo de texto normal, assim poderia ser escrito em qualquer editor simples (usuários windows: eso significa usar o bloco de notas). Note que o # é o mesmo simbolo utilizado para comentário. Para executar o script existem diversas formas, a mais simples é digitar python nome_do_script; a outra forma é converter o script em um executável: chmod u+x nome_do_script e depois executar ele utilizando o comando ./ (execute): ./nome_do_script . O outro jeito é entrar no ipython e digitar run nome_script (note que no devem ser usados acentos!!). Esse último método é o que será utilizado com windows quando se instala o python (x,y) (o problema é que por enquanto só tem para 32 bits)

primeiro script desde o ipython
Figura 02: O ipython.

Para sair do ipython basta digital CTRL+d (a tecla control junto da tecla d). O ipython é um shell interativo em python, você pode rodar comando python diretamente nele. Para rodar o ipython basta digitar no terminal ipython e se você deseja utilizar a matplotlib é melhor digitar no terminal ipython -pylab

O python é extremadamente chato com a tabulação, por favor cuide isso, tabular em python significa entrar num bloco tipo if ou for (Do em Fortran), while, etc.

Variáveis

O python tem a peculiaridade de não ser tipificado, ou seja, se você trabalha em fortran o C (e derivados), você sempre necessita definir as variáveis, se são inteiras, reais (float) ou caráteres (strings), no python isso não é necessário, as variáveis (constantes) são objetos e na execução o python saberá que se trata de uma variável de um tipo determinado.


          #!/usr/bin/env python

          a = 2
          b = a + 1
          a = '2'
          c = 'a vale ' + a

          print b
          print c

          

Similar a outras linguagens as variáveis básicas do python (na verdade os objetos) podem ser strings (palavras - cadeia de caráteres), inteiros, reais (float), complexos (intrínseco a fortran 90 e superiores só), lógicos; mas como as variáveis armazenam objetos elas podem armazenar coisas ainda mais complexas como matrizes, figuras, gráficos, músicas, etc. Um ponto importante (principalmente para programador em C já que os "fortraneros" em geral nem sabem desse problema) é o referente à precisão, recomendo dar uma olhada no módulo decimal do python e em o documento relativo a representação de números de ponto flutuante - reais. Você pode dar a uma variável o nome que desejar, menos aqueles nomes associados às palavras reservadas: and, as, assert, break, class, continue, def, del, elif, else, except, False, finally, for, from, global, if, import, in, is, lambda, None, nonlocal, not, or, pass, raise, return, True, try, with, while

Podemos converter uma variável de um tipo em outro utilizando int (nome_variavel), float (nome_variavel) e str (nome_variavel). Veremos adiante que podemos até realizar conversões para outros tipos de objetos mais complexos.

Operações lógicas

O pytho permite trabalhar diretamente com objetos lógicas: True (verdadeiro) e False (falso) que podem ser atribuídos à uma variável (note que é em maiúscula),

Sobre esse objetos se admitem a seguintes operações lógicas> and, or e not, como exemplo temos

            #!/usr/bin/env python

            a = 1.0
            b = 2.0
            c = 3
            d = True
            e = (a>b)

            print ( e and (b<c) )
            print ( ( not e ) and d )

            t1 = 'mae'
            t2 = 'pae'
            t  = (t1==t2)

            print t
          

note que os parenteses ajudam na compressão da expressão.

Operações matemáticas básicas

Por padrão o python trás definidas as operações básica: soma (+), resta (-), divisão (/) e multiplicação(*), resto da divisão (%) (mod em fortran) e potenciação (**) as quais funciona para inteiros:

            #!/usr/bin/env python

            a = 2
            b = 3

            c = a + b; print '%d + %d = %d' % (a, b, c)
            c = a - b; print '%d - %d = %d' % (a, b, c)
            c = a * b; print '%d * %d = %d' % (a, b, c)
            c = a / b; print '%d / %d = %d' % (a, b, c)
            c = a % b; print '%d %% %d = %d' % (a, b, c)
            c = a**b; print '%d**%d  = %d' % (a, b, c)
          

números reais (floats)

            #!/usr/bin/env python

            a = 2.0
            b = 3.0

            c = a + b; print '%f + %f = %f' % (a, b, c)
            c = a - b; print '%f - %f = %f' % (a, b, c)
            c = a * b; print '%f * %f = %f' % (a, b, c)
            c = a / b; print '%f / %f = %f' % (a, b, c)
            c = a % b; print '%f %% %f = %f' % (a, b, c)
            c = a**b; print '%f**%f  = %f' % (a, b, c)
          

números complexos (note que o print não tem como representar um complexo, assim primeiro foi convertido para string - usando str() - e depois impresso)

            #!/usr/bin/env python
            #!/usr/bin/env python

            a = 2.3 + 4.0j
            b = 3.1 + 6.5j

            c = a + b; print '%s + %s = %s' % (str(a), str(b), str(c))
            c = a - b; print '%s - %s = %s' % (str(a), str(b), str(c))
            c = a * b; print '%s * %s = %s' % (str(a), str(b), str(c))
            c = a / b; print '%s / %s = %s' % (str(a), str(b), str(c))
            c = a % b; print '%s %% %s = %s' % (str(a), str(b), str(c))
            c = a**b; print '%s**%s = %s' % (str(a), str(b), str(c))
          

e a soma com strings (adiante veremos muitas mais operações com string, que é uma classe especial)

            #!/usr/bin/env python

            a = 'ola'
            b = ' mundo'
            c = 2012

            c = a + b
            print '%s + %s = %s' % (a, b, c)
            d = a + b + ' em ' + str(c)
            print c
          

O ponto e vírgula nos exemplos prévios significa fim da linha. Além dessas operações o python também disponibiliza recursos para manipular dados em nível de bits: deslocamento à esquerda (<<), deslocamento à direita (>>), e (&), ou (|), xor (^) e complemento de 1 (~).

Importando um módulo

Para importar um módulo o procedimento é como se mostra no exemplo embaixo

            #!/usr/bin/env python

            import math

            a = 4.0
            b = math.sqrt(a)
            print b
          

As vezes o nome do módulo é comprido resultado enfadonho digitar ele, para isso podemos importar o modulo com outro nome, como mostrado embaixo

            #!/usr/bin/env python

            import math as m

            a = 4.0
            b = m.sqrt(a)
            print b
          

Para evitar a chatice de ter que escrever o nome do modulo, mais um ponto, a fim de acessar à função (atributo) desejada podemos importar todo o módulo, para isso fazemos

            #!/usr/bin/env python

            from math import *

            a = 4.0
            b = sqrt(a)
            print b
          

finalmente, podemos importar só uma função (atributo) que desejemos de um dado módulo

            #!/usr/bin/env python

            from math import sqrt

            a = 4.0
            b = sqrt(a)
            print b
          

Principais Funções matemática

Para acessar às funções matemática devemos importa a classe matemática do python, math, caso seja necessário trabalhar com números complexos devemos utilizar cmath. O módulo math permite acesso a funções como: sin, cos, tan, sinh, cosh, exp, log, pi, e, etc

            #!/usr/bin/env python

            from math import sinh, exp, e, pi

            x  = 2*pi
            r1 = sinh(x)
            r2 = 0.5 * (exp(x) - exp(-x))
            r3 = 0.5 * (e**x - e**(-x))
            print '%.16f %.16f %.16f' % (r1, r2, r3)
          

Contudo, resulta mais interessante importar o módulo numérico do python que é muito mais completo e que além das operações matemáticas do math nos permite trabalhar com matrizes. Esse módulo é o numpy (numerical python) que é parte do esforço das ferramentas do python para matemática, ciências e engenharia conhecido como scipy. Uma lista completa das funções acessíveis do módulo numpy (com exemplos) pode ser encontrada aqui. Igualmente resulta interessante dar uma olhada no livro de receitas do scipy para ter acesso a funções não tradicionais.

            #!/usr/bin/env python

            import numpy as np
            from fractions import Fraction  #aritmetica racional
            from sympy.mpmath import mp     #matematica simbolica

            theta       = 30.0
            thetaRad    = np.deg2rad(theta)
            thetaRad_pi = np.deg2rad(theta) / np.pi

            f = Fraction(thetaRad_pi).limit_denominator()  #cria uma fracao

            print '%.1f graus =  %f radianos = %s pi radianos' % (theta, thetaRad, f)

            senoAng = np.sin(thetaRad)
            print 'o seno de %.1f graus = %f' % (theta, senoAng)

            cosenoAng2 = np.sqrt(3.0) / 2.0
            theta2Rad  = np.arccos(cosenoAng2)  #arccos nao acos
            theta2     = np.rad2deg(theta2Rad)

            print '%.1f graus tem um coseno igua a = %.15f' % (theta2, cosenoAng2)

            mp.dps = 100                                  #numero de digitos
            piLongo = mp.pi
            print 'pi com 100 decimais = %s' % piLongo    #piLongo como uma string
            print 'pi com 100 decimais = %.100f' %piLongo #piLongo como um decimal
          

Vale destacar que enquanto numpy define as funções trigonométricas "inversas" como: arcsin, arccos, arctan, o módulo math utiliza a mesma definição utilizada no C: asin, acos, atan, isso também é válido para as funções hiperbólicas.

Imprimindo no terminal

O print é o comando de escrita ao terminal que já foi utilizado em vários pontos deste texto. O que pode merecer um pouco mais de atenção é a formatação que pode ser dado ao dados. Ao igual que no C, o print converte para string aquilo que vai ser impresso no terminal, assim essa formatação pode ser aplicada às string. Para que a conversão seja feita de forma direta é necessário informar ao print qual é a formatação que será dada ao dados segundo:

print "formatação" % (variável1, variável2, ...)

onde a formatação tem a sintaxe

%[flags][largura][.precisão]código

em que flags é dada por uma das opções da seguinte tabela

flag significado
* argumento especifica a largura ou precisão
- justificado á esquerda
+ mostra o sinal positivo
<spc> deixa um espaço em branco antes de um número positivo
# adiciona zero (0) ao octal ou (0x ou 0X) ao hexadecimal
0 completa com 0 à esquerda a largura do campo ([largura]).

largura diz a quantidade de campos que serão utilizados para imprimir e precisão diz a máxima quantidade de dígitos apos o ponto decimal.

onde código diz respeito a um dos mencionado na coluna 1 da próxima tabela

código significado
s String
c Carácter
d, i Número decimal.
o Número octal
x, X Número inteiro hexadecimal. o x imprime as letras em minuscula e o X em maiúscula.
f imprime um número real em notação decimal.
e, E Real em notação científica. O e imprime em minúscula e o E maiúscula.
g, G Dependendo de qual tiver o menor comprimento na impressão, utiliza f (ou F) ou e (ou E).
% Imprime o caráter porcentagem (%).

No exemplo embaixo se exemplifica o uso da impressão formatada. Aproveitamos e mostramos o módulo de matemática simbólica do python (tipo maple, matemática, maxima, etc)

            #!/usr/bin/env python

            import sympy                     #modulo de comp. symbolica

            x = sympy.Symbol('x')                  #x eh um simbolo
            y = sympy.integrate(x**2 + x + 1, x)   #avalia a integral

            print y                    #imprime em tela o resultado, eh uma string
            sympy.pprint(y)            #pretty print, funcao de sympy
            y2 = y.subs(x, -3)         #subtitui x por -3 e avalia

            print "%+07.2f" % (y2.evalf())  #converte a real, se positivo imprime +
                                            #7 campos, 2 digitos apos zero, prenche
                                            #com zero a esquerda o espaços em branco
          

a título de esclarecimento, na linha 6 foi criado um objeto do tipo sympy de forma que herda todos os métodos da classe sympy, por isso na linha 10 podemos acessar o método subs e na linha 12 o método evalf().

Lendo do terminal

Para ler do terminal o python possui 2 funções (na realidade mais, podemos fazer um open do stdio com opção r): input() e raw_input(). A diferença entre elas é que raw_input() sempre devolve uma string, enquanto que input() avalia a operação realizada entre os objetos de entrada (ou seja realiza uma conversão de tipo), vejamos o exemplo

            #!/usr/bin/env python

            a = raw_input("Por favor, digite sua idade ")
            b = input("Por favor, digite sua idade novamente ")

            print "primeiro informastes que tens %s anos" % a
            print "depois informastes que tens %s anos" % b
          

Se você digitar algo como 21 + 2 a variável a vai armazenar literalmente o que você digitou: 21 + 2, enquanto que a variável b vai armazenar 23, já que o input avalia antes de salvar na variável, mas o resultado em b é uma string. Como a leitura do terminal sempre da uma string o melhor é converter ela para o tipo de dado que se quer, no exemplo anterior é obvio que se quer um inteiro por tanto o script muda para

            #!/usr/bin/env python

            a = int( raw_input("Por favor, digite sua idade ") )
            b = int( input("Por favor, digite sua idade novamente ") )

            print "primeiro informastes que tens %03d anos" % a
            print "depois informastes que tens %03d anos" % b
          

de forma que se você digitar de novo 21 + 2 vamos ter um erro de conversão já na atribuição em a; caso voce digite 23 e depois 21 + 2, não dará ero. Veja que agora podemos formatar a saída já que a e b são de tipo inteiro (imprime um inteiro com no mínimo 3 campo, preenchendo com zeros à esquerda as lacunas).

Strings e modulo de expressões regulares

Para escrever uma string basta colocar o que se deseja entre aspas simples ou duplas, aspas triplas são utilizadas para definir um bloco (um paragrafo). Além das string temos as sequências: \\ (barra invertida), \' (apostrofo), \" (aspas), \a (bip), \b (retrocesso), \f (avanço de formulário), \n (nova linha), r (cariage return), \t (tabulação horizontal), \v (tabulação vertical) e \u que permite escrever em unicode (exemplo \u'\xe4' que é ä -- ver a tabela). Em relação à conversão de unicode (por exemplo UTF-8) para latin o python disponibiliza o comando unicode() que converte de "normal" para unicode e str(), que converte de unicode para normal.

Como em todas as linguagens, strings são vetores chamados de imutáveis, consequentemente cada letra (tudo entre aspas é uma letra), é um elemento desse vetor. O vetor iniciam em 0 e vão até $n-1$. É possível pesquisar o tamanho de uma string utilizando a função len()

            #!/usr/bin/env python

            a = ( raw_input("Por favor, digite algo:\n") )

            print "voce digitou: %0s" % a
            print "que tem um comprimento de %d" % len(a)
            print "o primeiro carater eh: %c" % a[0]
            print "o ultimo carater eh: %c" % a[-1]
            b = len(a) / 2
            print "o carater do meio eh: %c" % a[b]

            parte = a.split()
            print 'tem %d palavras' % len(parte)
            print "o primeiro palavra eh: %s" % parte[0]
            print "a ultimo palavra eh: %s" % parte[-1]
          

O python disponibiliza uma serie de funções que atuam sobre as string como o count do exemplo embaixo. É importante conhecer alguns deles pois devemos lembrar que tudo o que o python lé como arquivo é uma string (via raw_input, input, open ou ate mesmo utilizando o módulo sys), assim por exemplo, se você tem um arquivo de varias linhas e varias colunas você pode ler a linha, mas as colunas você separa usando split(). As funções disponíveis para tratar strings são muitas mas em geral as que a gente necessita para tratar nossos dados são poucas.

            #!/usr/bin/env python


            a = 'Para escrever uma string basta colocar o\
            que se deseja entre aspas simples ou duplas \
            aspas triplas sao utilizadas para definir um bloco'

            b = ( raw_input("Por favor, digite uma letra:\n") )
            c = a.count( b )
            print 'no texto:\n%s \na letra %s apareceu %d vezes' %\
            (a, b, c)
          

Mencionamos quando falamos sobre a formatação do print que podemos utilizar ela para formatar uma string como se mostra no próximo exemplo

            #!/usr/bin/env python

            import os    #modulo de acesso ao sistema operacional

            Rho     = 3.0
            Temp    = 0.543
            Semente = 1

            sRho     = '%06.2f' % Rho
            sTemp    = '%06.3f' % Temp
            sSemente = '%02d' % Semente

            file = 'dado_rho_' + sRho + '_temp_' + sTemp + '_seed_' + sSemente + '_.dat'
            print file, os.path.isfile(file)   #os.path.isfile checa se o arquivo
                                               #'file' existe
          

Expressões Regulares

Uma dos objetivos de se aprender python para auxiliar na pesquisa é utilizar só ele e esquecer o resto das funções existentes no Bash tipo grep, sed, awk, tail, head, etc., funções muito usadas na análise dos resultados. O módulo intrínseco de manipulação de strings faz muito desse trabalho contudo ele se complementa ainda mais com o uso de expressões regulares (existentes em grep, awk, sed, etc.). Não pretendo dar toda uma descrição de todas as funções regulares, mas só exemplificar alguns casos (e simultaneamente introduzir outros métodos do python) a fim de motivar à procurar por mais informação ao respeito

Uma expressão regular é um método formal de especificar um padrão de texto. Ou seja, são utilizados algumas para procurar por padrões em algum texto (string). Em python o módulo que da acesso às funções regulares é o re

Metacaratere significado
. casa qualquer caráter
[ ... ] lista de caráteres a serem casados (os pontos não fazem parte)
[ ^... ] lista de caráteres a não serem casados (os pontos não fazem parte)
? casa zero ou um caráter
* casa zero, um ou mais caráteres
+ casa um ou mais caráteres
{n,m} casa de n até m coisas, exemplo 5{1,3} casa: 5, 55, 555
^ casa no inicio, ex. ^a: abra, ala, ana, arpa, ...
$ casa no fim da linha, ex. $r: radar, camihnar, ...
\b casa no inicio ou fim da palavra
\c escape, transforma em litera o c, ex. \^: literal o ^, \.: literal o .
| ou, exemplo ^a|^b: ana, bola, aras, barca, ...
( ... ) casa o seguinte grupo, ex. (ba): aba, ambar, bah, ..
\1..\9 casa para atrás, ex.

Na tabela embaixo apresento as classes utilizadas para nos referirmos ao elementos sobre as quais podemos aplicar as expressões regulares, primeiramente as classes do tipo caráter

classe caráter significado
\c Control
\s Espaço em branco
\S Outro caráter diferente de espaço em branco
\d Dígitos, ex. 1, 2, 3, 4, ...
\D Outro caráter diferente de dígito
\w Caráter tipo letra - alfanumérico (letra e dígitos) e _, ex.
\W Qualquer caráter diferente do tipo letra
\x Carácter do tipo hexadecimal
\O Carácter do tipo octal

Tabela anterior na formatação POSIX

classe POSIX significado
[:upper:] letras maiúscula
[:lower:] letra minúsculas
[:alpha:] Todas as letras
[:alnum:] Dígitos e letras
[:digit:] Dígitos
[:xdigit:] Dígito hexadecimal
[:punct:] Carateres de pontuação
[:blank:] Espaço e tabulações
[:space:] Espaços
[:cntrl:] Caráter control
[:graph:] caráter imprimível
[:print:] Caráter imprimível e espaços
[:word:] dígitos, letras, "_" (underline)

Como exemplo vejamos o seguinte script python

            #!/usr/bin/env python

            import re    #modulo de expressões regulares

            file  = 'dado_rho_003.00_temp_00.543_seed_01_.dat'

            Parte =  re.split('\D\D+',file) #procura mais de um não dígitos 
                                            #contíguos
            print Parte
          

Blocos de decisão ou blocos if

Similar a outras linguagens o python possui o if para tratar da toma de decisões. A principal diferença é que o python utiliza a indentação para definir o bloco lógico ou seja, o bloco que será executado caso a condição do if seja verdadeira, vejamos um exemplo para entender melhor

            #!/usr/bin/env python

            from math import *
            from cmath import sqrt as csqrt

            coef = raw_input("digite os coeficiente a, b, c \
            da eq de II grau, separados por virgula ")

            Parte = coef.split(',')
            a = float(Parte[0])
            b = float(Parte[1])
            c = float(Parte[2])

            radicando = b**2 - 4.0 * a * c

            if (radicando > 0.0):     #se radicando > 0 faz
              raiz = sqrt(radicando)
            elif (radicando < 0.0):   #caso radicando > 0 faz
              raiz = csqrt(radicando)
            else:                     #se nao faz
              raiz = 0.0

            raiz1 = (-b + raiz) / (2.0 * a)
            raiz2 = (-b - raiz) / (2.0 * a)

            print 'as raizes de ax**2 + bx + c = 0 sao:\n \
            raiz 1 = %s\n raiz 2 = %s' % ( str(raiz1), str(raiz2) )
          

Os parenteses à volta da expressão de teste não são obrigatórios. Note que terminada a opção você coloca : para informar o inicio do bloco, a próxima linha (linha 17, 19, 21) é indentada pois se refere ao bloco que deve ser computado caso a linha acima seja verdadeira.

A linha 4 apresenta uma forma de importação de módulo no antes vista, como tanto o módulo math como o cmath possuem a função (método) sqrt foi importada a correspondente ao método cmath como (as) csqrt. Dessa forma temos ambas funções a disposição.

Blocos de repetição

A estrutura de repetição padrão do python é o while (condição): bloco. A condição é alguma condição lógica válida a qual pode ser até mesmo só um True transformando esse bloco em um "loop" infinito, nesse caso devemos colocar um if que quando seja verdeiro execute um break do loop. A fim de entender melhor vejamos o seguinte exemplo que calcula a sequencia de fibonacci até que um número o próximo número da sequência seja maior do que 10000

            #!/usr/bin/env python

            f1 = 1
            f2 = 1

            proximo = f1 + f2

            while (proximo <= 10000):
              print proximo
              f1      = f2
              f2      = proximo
              proximo = f1 + f2  
          

Observe que podemos ter feito essa mesma sequência utilizando loop infinito

            #!/usr/bin/env python

            f1 = 1
            f2 = 1

            proximo = f1 + f2

            while True:
              print proximo
              f1      = f2
              f2      = proximo
              proximo = f1 + f2
              if (proximo > 10000):
                break
          

Além do while o python disponibiliza o for que é ótimo para percorrer elementos de um objeto complexo, tipo arquivos, vetores, etc. A função range() é uma função ótima para ser utilizada junto com for já que ela gera uma sequência de números que vão de zero até o valor colocado entre as chaves, vejamos uns exemplos

            #!/usr/bin/env python

            import numpy as np

            sigma = 1.0  #desvio standard
            ro    = 0.0  #media


            y = np.random.normal(ro, sigma, 1000) #sera criada um matriz de 1
                                                  # dimensao com 1000 elementos

            print type(y)  #type devolve o tipo do objeto, no caso numpy.ndarray
            print y.shape  #shape devolve a forma da matriz, no  caso(1000, )
            print y.ndim   #ndim devolve a quantidade de dimensoes da matriz
                           #no caso 1
            
            
            for i in range(y.shape[0]):     #shape[0] diz a quantidade de elemento
                                            #na dimensao 0 do objeto y
              print y[i]    #imprime o i esimo elemento do objeto
          

No exemplo acima são gerados 1000 números aleatórios e armazenados num objeto numpy.ndarray. Utilizando a função range obtemos a quantidade de dados na dimensão 0 da matriz y y logo é impresso cada dado do vetor. Outro uso comum do for é para ler as linhas de um arquivo aberto.

Arquivos

O comando de abertura de arquivo é o open('nome_arq', 'modo', buf ), onde nome_arq é o caminho até o arquivo partindo do ponto onde está o script - caso o arquivo esteja no mesmo local é só o mesmo nome do arquivo - e modo pode ser uma das seguintes opções: r - leitura -, w - escritura -, a - anexar -, b - (ex. rb leitura de um binário) binário -, + - (ex. r+) leitura + escritura ou vice versa. Se buf = 1 então python faz buffer do arquivo copiando em memoria e não no disco e só copiando ao disco quando for feito um flush() ou um close(). Si buf = 0 então o arquivo é copiado diretamente ao disco.

O python te permite escrever usando a função write() como mostra o exemplo embaixo

            #!/usr/bin/env python

            arq = open('texto.txt', 'w')
            arq.write('Ola, ')
            arq.write('mundo!')
            arq.close()
          

Podemos utilizar a função print para escrever no arquivo

            #!/usr/bin/env python

            arq = open('fibonacci.dat', 'w')

            f1 = 1
            f2 = 1

            proximo = f1 + f2

            while (proximo <= 10000):
              print >>arq, proximo  #note que o >> envia a saida do print para arq
              f1      = f2
              f2      = proximo
              proximo = f1 + f2

            arq.close()
          

Para ler podemos utilizar o comando read()

            arq = open('texto.txt', 'w')
            arq.read(2)
            'Ol'           #<-- saida no terminal
            arq.read()
            'ola, mundo!'  #<-- saida no terminal
          

Se nosso arquivo tem varias linhas podemos utilizar a função (método) readline() - ler linha:

            #!/usr/bin/env python

            arq = open('textoLongo.txt', 'w')
            While True:
              linha = arq.readline()
              if not linha: break
              temp    = linha.rstrip()  #remove o '\n' final
              colunas = temp.split()    #separa em colunas
            arq.close()
          

o, utilizar um laço for para ler a linha

            #!/usr/bin/env python

            arq = open('dados.dat', 'r')

            for linha in arq:
              temp    = linha.rstrip()
              colunas = temp.split()

            arq.close()
          

Na linha 3 é aberto o arquivo dados.dat para leitura (r de read). O for da linha 5 percorre todas as linhas. Como cada linha acaba em um ENTER (retorno de carro \n) é necessário eliminar esse caráter final isso é feito por rstrip() e a nova linha modificada é armazenada em temp (de temporário). Depois linha é separadas em colunas usando como delimitador o espaço em branco. Finalmente fechamos o arquivo.

Caso queiramos ler uma opção da entrada padrão, ao estilo opção para o programa, utilizamos o modulo sys

            #!/usr/bin/env python

            import sys

            texto    = sys.stdin.read()
            palavras = texto.split()
            quanPal  = len(palavras)
            print 'O total de palavras no texto eh = %d' %quanPal
          

quando executado no terminal como por exemplo python script.py texto.txt ele vai imprimir a quantidade de palavras no texto.

Vetores e Matrizes: versão padrão do python

O python disponibiliza 3 padrões diferentes de tratar vetores/matrizes, essas formas padrões podem não ser as mais adequadas para tratar problemas matemáticos, para isso é melhor utilizar as matrizes do modulo númerico contudo, esse métodos resultam muito adequados e bem vidos na análise dos dados, por isso vamos dar uma olhada nesses vetores intrinsecos.

Listas

Você pode construir uma lista em python de um forma muito simples, como por exemplo

              #!/usr/bin/env python
              
              lista1 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
              lista2 = ['Galileo', 'Newton', 'Descartes', 'Leibniz', 'Huygens', 'Young', 'Faraday', 'Maxwell', 'Boltzmann', 'Langevine', 'Einstein']


              print 'usando range e len para percorrer com a lista'                
              for i in range(len(lista1)):
                print lista1[i]

              print '\n'
              print 'quarto elemento da lista 2 = ', lista2[3]

              print '\n'
              print 'usando for para percorrer a lista'
              for fisico in lista2:
                print fisico

              print '\n'
              print 'numero de elementos na lista 1 = %d' % len(lista1)
              print 'menor elemento da lista 1 = %d'      % min(lista1)
              print 'maior elemento da lista 1 = %d'      % max(lista1)

              print '\n'
              print 'numero de elementos na lista 2 = %d' % len(lista2)
              print 'menor elemento da lista 2 = %s'      % min(lista2)
              print 'maior elemento da lista 2 = %s'      % max(lista2)

              print '\n'
              lista3 = [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
              lista4 = lista3 + lista1
              print 'lista3 + lista1 = ', lista4

              print '\n'
              ordenaLista4  = sorted(lista4)
              print 'lista4 ordenada acendente - sorted = ', ordenaLista4

              print '\n'
              ordenaLista4 = sorted(lista4, reverse=True)
              print 'lista4 ordenada descente - sorted  = ', ordenaLista4

              print '\n'
              lista4.sort()
              print 'lista4 ordenada acendente - sort = ', lista4

              print '\n'
              lista4.sort(reverse=True)
              print 'lista4 ordenada descente - sort  = ', lista4

              print '\n'
              lista2.sort()
              print 'lista2 ordenada acendente - sort = ', lista2

              print '\n'
              lista2[3] = 'Michael Faraday'
              print 'lista2 = ', lista2

              a = lista2[-1]

              del lista2[-1]  #o elemento -1 eh o ultimo, -2 eh o antes do ultimo, etc
              print 'removido o ultimo elemento de lista2 = ', lista2

              print '\n'
              lista2.append('Young')
              print 'anexado um novo elemento na ultima posicao de lista2 = ', lista2

              print '\n'
              print 'posicao na lista 2 do Michael Faraday = ', lista2.index('Michael Faraday')

              print '\n'
              lista2.insert(2, 'Planck')
              print 'anexado um elemento na posicao 2 = ', lista2
            

Nesse exemplo vemos a funcionalidade da lista. Até foi mostrado como ordenar uma lista utilizando a função sort (ordena e salva na mesma variável), e sorted (ordena e salva em outra variável). Na linha 25 juntei duas listas, a forma lembra a aquilo que se faz com as strings. Além dessa função (método) existem muitos outras, por exemplo a função append que permite adicionar mais elementos a uma lista

Tuplas

São iguais às listas, a principal diferença é que são imutáveis. Outra diferencia é que quando são criadas se utiliza parenteses "()" e não colchetes "[]"

              #!/usr/bin/env python
                tupla1 = (1, 'florianopolis', 17, 'Brasil', 5)

                tupla2 = tuple( sorted(tupla1) )

                for item in tupla2:
                  print item
            

O métodos aplicáveis às tuplas são quase os mesmo da lista, você só deve ter em mente que a tupla não se modifica.

Dicionários

São estruturas onde o índice da matriz não necessariamente é um número, e sim um palavra. Dessa forma na verdade o que se tem é um conjunto de dados chave-valor onde o valor pode ser acessado a través da chave, vejamos um exemplo

              #!/usr/bin/env python

              fisicos = [ 'Galileo', 'Newton', 'Descartes', 'Leibniz', 'Huygens', 'Young', 'Faraday', 'Maxwell', 'Boltzmann', 'Langevine', 'Einstein' ]
              descobertas = [ 'relatividade clasica', 'Gravitacao Universal', 'conservacao do momentum', 'vis viva', 'luz como onda', 'interferencia',\
              'inducao eletromagnetica', 'ondas eletromagneticas', 'estatstica: equiparticao, entropia', 'eq. dif. para o mov. Browniano', 'relatividade especial e geral, etc.' ]

              dic1 = dict( zip(fisicos,descobertas) )

              print '\n'
              print ' ------------------- dicionario 1:  -------------------'
              print dic1
              print '\n'


              print ' ------------------- items no dicionario 1:  -------------------'
              print dic1.items()
              print '\n'

              print ' ---------- numero de elementos no dicionario 1:  ---------'
              print len( dic1 )
              print '\n'

              lista1 = dic1.iteritems()

              print ' ------ items no dicionario 1 transformado em lista:  ------'
              for item in lista1:
                print item
              print '\n'

              print ' ----------- valores no dicionario 1:  -----------'
              for dado in dic1.values():
                print dado
              print '\n'

              print ' ------------------- chaves do dicionario 1:  -------------------'
              for dado in dic1.keys():
                print dado
              print '\n'

              dic2 = {'Galileo': 'relatividade clasica', 'Newton': 'Gravitacao Universal', 'Descartes': 'conservacao do momentum', 'Leibniz': 'vis viva', 'Huygens': 'luz como onda', 'Young': 'interferencia', 'Faraday': 'inducao eletromagnetica', 'Maxwell': 'ondas eletromagneticas',  'Boltzmann': 'estatstica: equiparticao, entropia', 'Langevine': 'eq. dif. para o mov. Browniano', 'Einstein': 'relatividade especial e geral, etc.'}

              print ' ------------------- dicionario 2:  -------------------'
              print dic2
              print '\n'

              print ' ------------- remove um item de dicionario 2:  -------------'
              dic2.pop('Langevine')
              print dic2
              print '\n'

              print ' ------ acesso a elementos especificos do dicionario 2:  ------'
              print 'Faraday no dic1 = %s, e no dic2 = %s' % (dic1['Faraday'] , dic2['Faraday'])
              print '\n'

              dic3 = {}
              dic3['Galileo'] = 'relatividade clasica'
              dic3['Newton'] = 'Gravitacao Universal'

              print ' ------------------- dicionario 3:  -------------------'
              print dic3
              print '\n'

              print ' ------------------- apagando dicionario 3:  -------------------'
              dic3.clear()
              print dic3
              print '\n'

              dic4 = {}
              dic4 ['Newton'] = ['Gravitacao Universal', 'F=GMm/r^2']
              dic4['Einstein'] = ['relatividade especial e geral, etc.', 'E = mc^2']

              print ' -------- dicionario 4: dicionario de lista ---------'
              for dado in dic4['Newton']:
                print dado
              print '\n'

              print ' ----- checando a existencia de uma chave no dicionario 4: -----'
              print dic4.has_key('Lagrange')
              print '\n'
            

Um exemplo usando listas e dicionários na analise de dados

Uma aplicação direta na análise de dados consiste em criar uma lista de dicionários. É comum termos situações onde nosso programa gera uma infinidade de dados, pastas, etc. Uma prática comum é colocar no nome do arquivo/pasta informações relevantes a respeito dos parâmetros que caraterizam aquele resultado, dessa forma no pós processamento precisamos associar essa informação contida no nome do arquivo/pasta à informação dentro do arquivo/pasta. A fim de me explicar melhor vejamos o seguinte exemplo, na figura embaixo observamos varias pastas e uns poucos arquivos. Cada pasta tem dentro de se própria os mesmos arquivos mostrados na figura: ecin, epot, est_fin.dat, est_ini.dat, etot, facv.dat.

Analise de dados
Figura 03: Pastas com dados a serem analisados.

Note que o nome da pasta informa a temperatura e densidade, que foram os parâmetros mantidos fixos durante o cálculo das grandezas dentro da pasta. Nosso objetivo é calcular a energia média em função temperatura, mantendo a densidade fixa: para isso temos que temos que:

  1. Escolher somente as pastas, ler seu nome e dividir este de forma a saber qual é a temperatura e densidade.
  2. Entrar dentro da pasta, ler a pasta etot.dat e calcular o valor médio.
  3. Criar um dicionario onde colocamos {'arq': pasta, 'rho': valRho, 'temp': valTemp, 'etot': valEne} ou seja, nome do arquivo, valor de $\rho$ para aquele arquivo, valor da temperatura e o valor da energia média.
  4. Anexar o dicionario anterior a uma lista, isso forma uma matriz de 4 colunas e muitas linhas = número de densidades totais vezes o numero de temperaturas totais.
  5. Ordenar a linha tomando como base a densidade, assim teremos varias linhas contíguas da matriz com a mesma densidade, mas as temperaturas bagunçadas.
  6. Ordenar em temperatura, agora além de termos as linhas contíguas da matriz com a mesma densidade, as temperaturas estão em ordem.
  7. Criar uma pasta com o nome e_vs_rho_fixT_valoDaTemp.dat onde valoDaTemp é a temperatura à qual vão corresponder os dados dentro da pasta. Dentro da pasta serão colocados as densidades (ordenadas) e os valores de energia associados.

O código para isso é listado embaixo

              #!/usr/bin/env python

              import numpy as np
              import re, os, os.path, shutil, glob
              from operator import itemgetter

              dirData            = os.getcwd()
              padraoTeste        = 'den_*'
              listaArqCasaPadrao = glob.glob(padraoTeste)

              matrizDados = []
              chave       = ['arq', 'rho', 'temp', 'etot']

              for pasta in listaArqCasaPadrao:
                if ( os.path.isdir(pasta) ):

                  quebraNome = re.split("_", pasta)
                  sRho  = "0.%s" % (quebraNome[1])
                  sTemp = "%.1f" % ( int( quebraNome[3])*0.1 )

                  etotArq     = os.path.join(dirData, pasta, 'etot.dat')
                  tempo, etot = np.loadtxt(etotArq, unpack=True)
                  ene         = np.mean(etot)

                  val        = [pasta, sRho, sTemp, ene ]
                  linhaDados = dict(zip(chave, val))
                  #linhaDados = {'arq': pasta, 'rho': sRho, 'temp': sTemp, 'etot': ene}
                  matrizDados.append(linhaDados)

              OrdenaRho  = sorted(matrizDados, key=itemgetter('rho'))
              OrdenaTemp = sorted(OrdenaRho, key=itemgetter('temp'))

              tempAtual = '-100.0'
              for i in range( len(OrdenaTemp) ):
                sRho  = OrdenaTemp[i]['rho']
                sTemp = OrdenaTemp[i]['temp']
                etot  = OrdenaTemp[i]['etot']
                if (sTemp != tempAtual):
                  tempAtual = sTemp
                  print sRho, etot
                  fileName = 'e_vs_rho_fixT_' + sTemp + '_.dat'
                  try:
                    fileT    = open(fileName, 'w')
                  except:
                    fileT.close()
                    fileT    = open(fileName, 'w')
                print >>fileT, sRho, etot

            

Nesse código são utilizados alguns módulos novos:

  1. os: utilizado para acessar funções do sistema operativo (operational system) como chown, chmod, setui, setgui, chdir (cd), getcwd (igual a pwd do bash), link, symlink, mkdir, remove, rename, rmdir, kill, system (permite executar um comando do qualquer do sistema).
  2. os.path: modulo que permite manipular o caminho (que é parte do nome) de um arquivo. Possui como função: abspath, basename, exists, expanduser, getatime, getmtime, isabs, isfile, isdir, islink, ismount, join (cria um caminho), samefile, split (contrario a join), supports_unicode_filenames.
  3. shutil: ferramentas do shell. As funções que possui são: copyfile, copymode, copy, move, disk_usage, chown,
  4. glob: procura padrões utilizando caráteres * como curinga. É equivalente a fazer algo como ls *.dat no terminal.
  5. itemgetter: permite extrair um item de um determinado objeto (tipo uma dicionario, etc.)

Analisando o código:

Funções definidas pelo usuário

Funções, como sabemos, são pequenos programas dentro de nosso programa, assim constituem um bloco com caraterísticas próprias desenhado para realizar uma tarefa específica. Para definir uma função utilizamos a estrutura seguinte: def nomeDaFuncao (variaveisDeEntrada): Note os dosi pontos no fim. Por ser um bloco, segue a mesma ideia dos outros blocos, fica definido pela endentação do código embaixo da definição da função. As funções podem ser definidas em qualquer parte do programa, desde que seja antes da primeira vez que for chamada, mas é recomendável sempre definir elas logo no inicio do script. Você pode definir ela em outro arquivo: suponhamos que o arquivo tenha por nome minhasFuncoes.py, e uma das funções que você definiu dentro dele seja a superFuncao, para que seu script possa acessa à superFuncao você deve importar a função como temos realizado ate agora com os módulos, exemplo: from minhasFuncoes import superFuncao. As funções só devolvem um objeto, mas lembre que um objeto em Python pode ser uma estrutura extremadamente complexa como um numero, uma lista, tupla, um array ou mesmo pode não devolver nada (por exemplo uma função que grava dados no disco), o retorno de dados se faz mediante a utilização do comando return estrutura. Vejamos uns exemplos

            #!/usr/bin/env python

            import numpy as np

            def ForcaEletrica(q, r1, r2):
              """Esta funcao recebe como dados de entrada
              um vetor de duas dimensoes com as cargas das
              particulas, a posicao da particula 1 e a posicao
              de particula 2 """

              k = 9.0e9

              Dr  = r2 - r1
              mDr = np.sqrt(np.dot(Dr,Dr))

              F = k * q[0]*q[1] / mDr**3 * Dr

              return F

            #---------------------------------------------

            #Programa Principal
            q  = np.array([10.0, 2.0])
            r1 = np.array([0,0,0])
            r2 = np.array([10,10,10])

            F = ForcaEletrica(q, r1, r2)

            print F
          

Funções Lambda ou anônimas

São funções definidas na linha. Elas são legais porque podem entrar em lugares que uma função comum não consegue, por exemplo, na hora de definir uma lista, você pode, utilizando lambda, definir uma função para cada elemento da lista. A sintaxe é simples:

lambda argumento 1, argumento 2, ..., argumento N : expressão usando argumentos

para esclarecer vejamos um exemplo de uma função definida de forma normal e depois utilizando lambda:

            #!/usr/bin/env python

            #define factorial utilizando uma funcao
            def fatorial_def(x):
              if (x <= 1):
                return 1
              else:
                return x * fatorial_def(x-1)

            #define factorial utilizando lambda
            fatorial_lab = lambda x: 1 if (x <= 1) else x*fatorial_lab(x-1)

            print fatorial_def(5), fatorial_lab(5)
          
Blocos if dentro de uma função lambda debemos em geral se vale de uma propriedade de operadores lógicos do python: Quando se avalia um teste lógico composto (utilizando operadores) o python sempre devolve o objeto do lado esquerdo do operador caso seja certo o teste lógico e se é falso devolve o objeto do lado esquerdo, como se mostra no seguinte exemplo
            #!/usr/bin/env python

            a = 12
            b = 7
            c = 5

            s = ((a > 13) and b) or c
            print s

            s = ((a > 6) and b) or c
            print s
          
Note que equivale a fazer:
            #!/usr/bin/env python

            if (a > 13):
              s = b
            else:
              s = c
            print s

            if (a > c):
              s = b
            else:
              s = c
            print s
          
A Função lambda pode ter comportamentos complexos, como quando está aninhada dentro de um outra função, é como se tivéssemos uma função de função
            #!/usr/bin/env python

            def funcao (x):
              return lambda y: x+ y

            f = funcao(99)
            print f

            print f(10)
          
o resultado desses print são <function <lambda> at 0xdaa5f0> e depois 109.

Vetores e matrizes, o jeito numpy

A principal diferencia com as listas é que as matrizes no numpy são de um tipo só ou seja, nas listas da para misturar string, inteiros, reais, complexos, etc, o que você bem entender; já nos vetores do numpy você só pode usar objetos de um mesmo tipo ou inteiros, ou reais ou complexos (sendo que cada um deles pode ser de diversas precisões, como se mostra na tabela embaixo)

Tipo Descrição
bool verdadeiro ou Falso, 1 bit
inti inteiro caraterístico do sistema (ou int32 ou int64)
int8 Byte ($-128$ a $127$)
int16 Inteiro ($-32768$ a $32767$)
int32 Inteiro $\left(2\times10^{31}\right.$ a $\left.2\times10^{31}-1\right)$
int64 Inteiro $\left(2\times10^{63}\right.$ a $\left.2\times10^{63}-1\right)$
uint8 Byte ($0$ a $255$)
uint16 Inteiro ($0$ a $65535$)
uint32 Inteiro $\left(0\right.$ a $\left.2\times10^{32}-1\right)$
uint64 Inteiro $\left(0\right.$ a $\left.2\times10^{64}-1\right)$
float16 real de precisão média: bit de sinal, expoente de 5 bits e 10 bits de mantisa
float32 real de precisão simples: bit de sinal, expoente de 8 bits e 23 bits de mantisa
float64 ou float real de precisão dupla: bit de sinal, expoente de 11 bits e 52 bits de mantisa
complex64 número complexo formado por 2 float32 - parte real e imaginaria
complex128 ou complex número complexo formado por 2 float64 - parte real e imaginaria

Definindo um array

Para definirmos uma matriz ou um vetor existem diversas forma

          #!/usr/bin/env python

          import numpy as np

          print '\n---- De uma lista python para um vetor numpy -----'
          lista = []
          lista.append(1.0+1j)
          lista.append(2.0+5j)
          lista.append(3.0-3j)
          a = np.array(lista, dtype=np.complex)
          print lista, a, a.dtype, a.shape

          print '\n-------- Usando array --------'
          a = np.array([23, 5, 14])
          print a, a.dtype, a.shape

          a = np.array([True, False, True])
          print a, a.dtype, a.shape

          a = np.array([[1.0, 3.0, 5.0], [7.0, 9.0, 11.0], [13.0, 15.0, 17.0]], dtype=np.float32)
          print a, a.dtype, a.shape

          print '\n-------- Usando arange ----------'
          a = np.arange(5,dtype=np.float16)
          print a, a.dtype, a.shape

          a = np.arange(15,dtype=np.float16).reshape(3,5)
          print a, a.dtype, a.shape

          a = np.arange(1.0, 10.0, 2.0, dtype=np.float)
          print a, a.dtype, a.shape

          print '\n----------- Usando zeros, ones e empty --------------'
          a = np.zeros( 5, dtype=np.int32 )
          print a, a.dtype, a.shape

          a = np.ones( 5, dtype=np.complex )
          print a, a.dtype, a.shape

          a = np.empty( 16, dtype=np.float32 ).reshape(4,4)
          print a, a.dtype, a.shape

          print '\n------ Modificando um elemento da matriz -------'
          np.set_printoptions(precision=1)          
          a[3,3] = 1.0
          a[2,1] = 9.0
          print a
        
[usuario@python:exemplo ]# python exemplo24.py
[(1+1j), (2+5j), (3-3j)] [ 1.+1.j 2.+5.j 3.-3.j] complex128 (3,)

Usando Array
[23 5 14] int64 (3,)
[ True False True] bool (3,)
[[ 1. 3. 5.]
[ 7. 9. 11.]
[ 13. 15. 17.]] float32 (3, 3)

Usando Arange
[ 0. 1. 2. 3. 4.] float16 (5,)
[[ 0. 1. 2. 3. 4.]
[ 5. 6. 7. 8. 9.]
[ 10. 11. 12. 13. 14.]] float16 (3, 5)
[ 1. 3. 5. 7. 9.] float64 (5,)

Usando linspace: gera dados desde inicio ate fim,
igualmente espaciados, incluindo extremos
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
[ 0. 0.25 0.5 0.75 1. 1.25 1.5 1.75 2. 2.25
2.5 2.75 3. 3.25 3.5 3.75 4. 4.25 4.5 4.75 5.
5.25 5.5 5.75 6. 6.25 6.5 6.75 7. 7.25 7.5
7.75 8. 8.25 8.5 8.75 9. 9.25 9.5 9.75 10. ]
Usando Zeros
[0 0 0 0 0] int32 (5,)

Usando Ones
[ 1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j] complex128 (5,)

Usando Empty
[[ -1.57829049e-07 4.58392754e-41 2.81180559e-37 0.00000000e+00]
[ 3.17444459e-37 0.00000000e+00 3.17448405e-37 0.00000000e+00]
[ 3.17452351e-37 0.00000000e+00 3.17456297e-37 0.00000000e+00]
[ 3.17460243e-37 0.00000000e+00 3.17464190e-37 0.00000000e+00]] float32 (4, 4)

modificando um elemento da matriz
[[ -1.6e-07 4.6e-41 2.8e-37 0.0e+00]
[ 3.2e-37 0.0e+00 3.2e-37 0.0e+00]
[ 3.2e-37 9.0e+00 3.2e-37 0.0e+00]
[ 3.2e-37 0.0e+00 3.2e-37 1.0e+00]]

algumas operações, elemento a elemento, que podem ser feitas com matrizes do numpy

          #!/usr/bin/env python

          import numpy as np

          #De uma lista python para um vetor numpy
          print 'Converte de uma lista para array'
          lista = []
          lista.append(1.0+1j)
          lista.append(2.0+5j)
          lista.append(3.0-3j)
          a = np.array(lista, dtype=np.complex)
          print lista, a, a.dtype, a.shape

          #Usando array
          print '\nUsando Array'
          a = np.array([23, 5, 14])
          print a, a.dtype, a.shape

          a = np.array([True, False, True])
          print a, a.dtype, a.shape

          a = np.array([[1.0, 3.0, 5.0], [7.0, 9.0, 11.0], [13.0, 15.0, 17.0]], dtype=np.float32)
          print a, a.dtype, a.shape

          #Usando arange
          print '\nUsando Arange'
          a = np.arange(5,dtype=np.float16)
          print a, a.dtype, a.shape

          a = np.arange(15,dtype=np.float16).reshape(3,5)
          print a, a.dtype, a.shape

          a = np.arange(1.0, 10.0, 2.0, dtype=np.float)
          print a, a.dtype, a.shape

          #Usando linspace
          print '\nUsando linspace: gera dados desde inicio ate fim,'
          print 'igualmente espaciados, incluindo extremos'
          a = np.linspace(0,10,11)
          print a
          a = np.linspace(0,10,41)
          print a

          #usando zeros, ones e empty
          print '\nUsando Zeros'
          a = np.zeros( 5, dtype=np.int32 )
          print a, a.dtype, a.shape

          print '\nUsando Ones'
          a = np.ones( 5, dtype=np.complex )
          print a, a.dtype, a.shape

          print '\nUsando Empty'
          a = np.empty( 16, dtype=np.float32 ).reshape(4,4)
          print a, a.dtype, a.shape

          np.set_printoptions(precision=1)
          print '\nmodificando um elemento da matriz'
          a[3,3] = 1.0
          a[2,1] = 9.0
          print a
        

Lendo e escrevendo no disco

O numpy possui uma função para ler desde um arquivo uma tabela de dado, a função loadtxt() e seu uso é muito simples, como se mostra no exemplo embaixo (use os dados: dados.dat):

          #!/usr/bin/env python

          from numpy import *

          print '\nLendo todos os dados e armazenando em uma matriz'
          DataIn = loadtxt('dados1.dat')
          print DataIn

          print '\nLendo dado a dado da matriz'
          for line in DataIn:
            print line[0], line[1], line[2]


          print '\nUsando unpack'
          x, y, yerr = loadtxt('dados1.dat', unpack=True)

          for i in range(len(x)):
            print x[i], y[i], yerr[i]

          print '\nUsando unpack para ler colunas especificas'
          x, y = loadtxt('dados1.dat', unpack=True, usecols=[0,1])

          for i in range(len(x)):
            print x[i], y[i]
        

No caso da escritura a gente dispõe da função savetxt(). É claro que seu uso é simples, como o da função loadtxt() más ela possui mais opções, como é de se esperar, já que nos podemos querer escrever formatado nossos dados, nesse caso devemos utilizar

%[flags][largura][.precisão]código

onde o código está dado pela tabela acima. Embaixo mostramos um exemplos (criando um gráfico e una animação)

          #!/usr/bin/env python

          import numpy as np
          import matplotlib.pyplot as plt
          from visual import *
          from PIL import Image
          import ImageGrab

          xo    = 0.0            #m
          yo    = 10.0            #m
          g     = -9.8           #m/s^2
          theta = np.radians(30.0)  #converte a radianos
          vo    = 10.0           #m/s

          vox = vo * np.cos(theta)
          voy = vo * np.sin(theta)

          #define os coeficiente da equacao de II grau
          a = 0.5 * g
          b = voy
          c = yo

          #define o polinomio de II grau
          polinomioY = np.poly1d([a, b, c])

          #imprime o polinomio
          print '\nequacao de II grau que define o mov. :'
          print polinomioY

          #tem duas raizes: zero (quando eh lanzado) e outra
          raiz = polinomioY.roots
          tv = raiz[np.where (raiz != 0.0) ]
          print '\ntempo de vo = ', tv[0]

          #calcula a derivada e tira a raiz (assim obtem o maximo)
          derivada = np.poly1d.deriv(polinomioY)
          ts       = derivada.roots
          print '\no tempo de subida eh = ', ts[0]

          #calculamos a posicao da particula em 100 pontos
          t = np.linspace(0.0, tv[0], 100)
          y = polinomioY(t)
          x = xo + vox * t

          ymax = polinomioY(ts[0])
          xmax = xo + vox * tv[0]

          #construindo uma matriz de 100 linhas e 2 colunas
          matrizDados = np.column_stack( (x, y) )

          #salvando num arquivo a matriz
          np.savetxt('projetil.dat', matrizDados, fmt=('%.2f', '%.2f'))

          #plotando os dados
          plt.plot(x, y, '*')
          plt.show()

          #animacao do disparo
          scene         = display(title='Trajetoria do projetil', width=800, height=800, center=(0.5*xmax,0.5*ymax,1.0), background=(1,1,1))
          eixoX         = arrow(pos=(0,0,0), axis=(xmax,0,0), shaftwidth=0.1)
          eixoY         = arrow(pos=(0,0,0), axis=(0,ymax,0), shaftwidth=0.1)
          eixoZ         = arrow(pos=(0,0,0), axis=(0,0,10), shaftwidth=0.1)

          objeto        = sphere()
          objeto.pos    = vector(xo, yo, 0)
          objeto.radius = 1.0
          objeto.color  = (0.5,0.5,0.5)
          objeto.traj   = curve(color=color.black, radius = 0.1)

          i = 0
          while i<100:
            rate(10)

            objeto.pos  = vector(x[i], y[i], 0)
            objeto.traj.append( pos=objeto.pos )
            i += 1
        

O gráfico resultante pode ser visto na figura embaixo:

primeiro plot com matplotLib
Figura 03: Matplotlib.

enquanto que a animação é mostrada a continuação:

trabalhando com matrizes

            #!/usr/bin/env python

            from pylab import *
            import numpy as np

            n = np.array([[3.0, -0.1, -0.2], [0.1, 7.0, -0.3], [0.3, -0.2, 10.0]])

            print '\nMatriz n:'
            print n

            print '\nA II coluna de eh'
            print n[:,1] #lembara que vai de 0 ate 2

            print '\nA II lina de eh'
            print n[1,:]

            print '\nTranposta da matriz '
            print np.transpose(n)

            print '\nvalor minimo contido na matriz'
            print np.amin(n)

            print '\nvalor maximo contido na matriz '
            print np.amax(n)

            print '\ndeterminante da matriz '
            print np.linalg.det(n)  #linalg pacote de algebra linear

            print '\ninversa da matriz'
            print np.linalg.inv(n)

            print '\ndecomposicao de valor singular da matriz n'
            print np.linalg.svd(n)

            A = np.array([[11, 2, 3, 1, 4], [2, 9, 3, 5, 2], [3, 3, 15, 4, 3], [1, 5, 4, 12, 4], [4, 2, 3, 4, 17]])
            print A

            print '\nAuto valores da matriz'
            print np.linalg.eigvals(A)

            #Os autovetores sao as colunas da matriz 1
            print '\nAuto vetores da matriz'
            eve1 = np.linalg.eig(A)[1][:,0]
            eve2 = np.linalg.eig(A)[1][:,1]
            eve3 = np.linalg.eig(A)[1][:,2]
            eve4 = np.linalg.eig(A)[1][:,3]
            eve5 = np.linalg.eig(A)[1][:,4]

            eva1 = np.linalg.eig(A)[0][0]
            eva2 = np.linalg.eig(A)[0][1]
            eva3 = np.linalg.eig(A)[0][2]
            eva4 = np.linalg.eig(A)[0][3]
            eva5 = np.linalg.eig(A)[0][4]

            print '\nAuto valor e auto vetor da matriz 1'
            print eva1, eve1
            print 'Auto valor e auto vetor da matriz 2'
            print eva2, eve2
            print 'Auto valor e auto vetor da matriz 3'
            print eva3, eve3
            print 'Auto valor e auto vetor da matriz 4'
            print eva4, eve4
            print 'Auto valor e auto vetor da matriz 5'
            print eva5, eve5

            prod  = np.dot(A,eve1) #produto matricial
            prod2 = eva1 * eve1  #produto de um escalar com matriz
            print '\nTestando Ax = lx'
            print prod, "=", prod2
          

Sistema linear de equações

Vamos supor o seguinte problema: Três paraquedistas estão ligados uns aos outros através de cordas sem peso enquanto caim, em queda livre, a uma velocidade de 5 m/s. Calcule a tensão em cada seção da corda e a aceleração do time sabendo os seguintes dados: $m_1 = 70$, $m_2 = 60$, $m_3 = 40$, $c_1 = 10$, $c_2 = 14$, $c_3 = 17$ Analisando o diagrama de corpo livre chegamos o seguinte sistema de equações lineares \begin{eqnarray} m_1g - T - c_1v & = & m_1a\\ m_2g - T - c_2v -R & = & m_1a\\ m_3g - c_3v +R & = & m_1a\\ \nonumber \end{eqnarray} onde $c_i$ são os coeficientes de arrastro. Substituindo os valores obtemos o seguinte sistema de equações \[ \nonumber \begin{bmatrix} 70 & 1 & 0\\ 60 & -1 & 1\\ 40 & 0 & -1\\ \end{bmatrix}\, \begin{bmatrix} a \\ T \\ R\\ \end{bmatrix} = \begin{bmatrix} 636 \\ 518 \\ 307\\ \end{bmatrix} \] Resolvendo com pynum usamos o método linalg.solve:

            #!/usr/bin/env python

            import numpy as np

            A = np.array([[70, 1, 0], [60, -1, 1], [40, 0, -1]])
            b = np.array([636, 518, 307])

            print np.linalg.solve(A,b)
          

Interpolação

O numpy e o scipy tem diversos módulos que permitem interpolar dados em um intervalo dado. No próximo exemplo são dados alguns valores de x e seguidamente se calcula o logaritmo. Com essa informação se calcular uma função, resultante da aplicação da interpolação de Lagrange nos dados e com isso se calcula o valor interpolado da função no ponto $x=2$, posteriormente se compara com o valor exato. No fim é mandada a plotar

            #!/usr/bin/env python

            import numpy as np
            import scipy.interpolate
            import matplotlib as mpl
            import matplotlib.pyplot as plt

            x   = np.array([1, 3, 4, 5, 6])
            lnx = np.log(x)

            #polinomio resultante da interpolacao de lagrange
            L = scipy.interpolate.lagrange(x,lnx)

            print 'O polinomio interpolado eh:'
            print L

            print '\nInterpolando em x = 2 obtemos por Lagrange:'
            print L(2.0)

            print '\nErro na interpolacao'
            print np.log(2.0) - L(2.0)

            #Usando o interp1d
            Il = scipy.interpolate.interp1d(x, lnx, kind='linear')
            print '\nInterpolando com interp1d linear em x = 2 obtemos:'
            print Il(2.0)

            Ic = scipy.interpolate.interp1d(x, lnx, kind='cubic')
            print '\nInterpolando com interp1d cubica em x = 2 obtemos:'
            print Ic(2.0)

            x    = np.linspace(1, 6, 25)
            y    = np.log(x)
            lagY = L(x)
            Ilin = Il(x)
            Icub = Ic(x)

            #plotando os dados
            mpl.rcParams['text.usetex'] =  True
            graf = plt.figure()
            fig1 = graf.add_subplot(111)
            fig1.plot(x, y, 'o', markersize=5, color= 'orange', label=r'$f(x) = \log(x)$')
            fig1.plot(x, lagY, '-',  color = 'yellow', label=r'y = Lagrange')
            fig1.plot(x, Ilin, '--', color = 'blue',   label=r'y = Linear'  )
            fig1.plot(x, Icub, '-.', color = 'black',  label=r'y = Cubica'  )
            fig1.set_ylabel('eixo Y')
            fig1.set_xlabel('eixo X')
            plt.legend(loc=2);
            graf.savefig( 'figPlt01.png', format='png' )
            plt.show()
          

O resultado se apresenta na próxima figura

Interpolação
Figura 05:Interpolação.

Outra forma de interpolação muito utilizada para a suavização das curvas é o cubic spltine. Vejamos um exemplo

            #!/usr/bin/env python

            import numpy as np
            import matplotlib as mpl
            import matplotlib.pyplot as plt
            import scipy.interpolate

            x = np.linspace(0, 2*np.pi, 10)
            y = np.sin(x)

            xnew = np.linspace(0, 2*np.pi, 20)

            Ic = scipy.interpolate.interp1d(x, y, kind='cubic')
            yc = Ic(xnew)

            spline = scipy.interpolate.splrep(x, y)
            ys     = scipy.interpolate.splev(xnew, spline)

            mpl.rcParams['text.usetex'] =  True
            plt.plot(x, y, '^', markersize=7, color= '#AAAAAA')
            plt.plot(xnew, yc, '-',  color = 'yellow')
            plt.plot(xnew, ys, '--',  color = 'black')
            plt.legend(['$y = \sin(x)$', 'cubic', 'spline'])
            graf.savefig( 'figPlt01.png', format='png' )
            plt.show()
          
Interpolação Spline
Figura 06:Interpolação por spline.

Ajuste de curvas polinomiais

Primeiramente vejamos como ajustar uma curva à um polinômio. Para exemplificar isso vamos gerar uns dados e adicionaremos a eles um pouco de ruido, a fim de simular dados experimentais, depois vamos ajustar com numpy.polyfit, o resultado dessa função são os coeficiente do polinômio de ordem n que você pressupõe ajuste os dados, para criar o polinômio em si utilizamos numpy.poly1d

              #!/usr/bin/env python

              import numpy as np
              import matplotlib as mpl
              import matplotlib.pyplot as plt
              import scipy.interpolate

              #Gera 41 dados igualmente
              x = np.linspace(0, 20, 41);
              y = x**3 - 20 * x**2 + 10 * x + 2

              #gera 41 dados aleatorios no intervalo
              #(0.0, 1.0]
              xr = np.random.random_sample((41,))
              yr = np.random.random_sample((41,))

              #adiciona ruido aos dados
              xnew = x + 0.2 * x * (xr - 0.5)
              ynew = y + 0.2 * y * (yr - 0.5)

              #calcula os coeficientes do polinomio
              print 'coeficientes do polinomio'
              coef = np.polyfit(x, y, 3)
              print coef

              #cria o polinomio
              print '\nPolinomio ajustado'
              pol = np.poly1d(coef)
              print pol

              #calcula os novos y com o juste
              yajus = pol(x)

              #cria o grafico
              plt.plot(xnew, ynew, 'o')
              plt.plot(x, y, '-')
              plt.plot(x, yajus, '--')
              plt.legend(['dados com ruido', 'curva limpa', 'ajuste'], loc=2)
              plt.savefig( 'figPlt03.png', format='png' )
              plt.show()
            
Que da como resultado:
coeficientes do polinomio
[ 1. -20. 10. 2.]

Polinomio ajustado
   3      2
1 x - 20 x + 10 x + 2
ajuste polinomial
Figura 06:Ajuste polinomial.

Ajuste no linear de curvas

Para se ajustar curvas não lineares utilizamos a função curve_fit do módulo scipy.optimize, de todas as funções até agora vistas está é a mais complexa de todas

    curve_fit(f, dadox, dadoy, p0=None, sigma=None, **kw)
onde Como resultado, a função curve_fit retorna um par $(popt,pcov)$ onde Como exemplo vamos ajustar os dados.dat (da tabela 1 do livo Métodos computacionais da Física de Claudio Scherrer). No exemplo 1.3 do livro citado os dados da I coluna representam as frequências, $\omega$ e os da II coluna a potencia e supõe-se sejam ajustados por uma "Lorentziana" dada por \[ P(\omega) = \frac{A}{(\omega - \omega_0)^2 + \Gamma^2} \nonumber \] onde $A$ é a amplitude, $\omega_0$ é a frequência de ressonância e $\Gamma$ é a largura, portanto nosso parâmetros a serem calculados da forma mais ótima possível. O código em python para esse ajuste é o seguinte
              #!/usr/bin/env python

              import numpy as np
              import scipy.optimize
              import matplotlib.pyplot as plt

              def lorentz(x, a, wo, c):
                p = a / ((x-wo)**2 + c**2)
                return p

              def residuo(parametros, y, x):
                err = y - lorentz(parametros, x)
                return err


              chute = np.array([16.0, 5.0, 1.0])
              x, y  = np.loadtxt('dados_claudio.dat', unpack=True)

              resAjuste, mcov = scipy.optimize.curve_fit(lorentz, x, y, p0=chute)

              #resAjuste = scipy.optimize.leastsq(func=residuo, x0=chute, args=(y, x), maxfev=2000, ftol=1.0e-8)

              print "Resultado do ajuste: ", resAjuste
              print "\nMatriz de covarianca (diagonal eh varianca dos resultados):"
              print mcov

              xnew = np.linspace(0.0, 12.0, 100)
              ynew = np.zeros(100)
              for i in range(len(xnew)):
                ynew[i] = lorentz(xnew[i], resAjuste[0], resAjuste[1], resAjuste[2])

              #mpl.rcParams['text.usetex'] =  True
              plt.plot(x, y, 'o', xnew, ynew, '-')
              plt.savefig( 'ajuste.png', format='png' )
              plt.show()
            
ajuste no linear
Figura 07:Ajuste não linear.
Que da como resultado:
[usuario@pclabfis: ]# python ajusteNaoLinear2.py
Resultado do ajuste: [ 21.62608754 4.68830763 1.14494062]

Matriz de covarianca (diagonal eh varianca dos resultados):
[[ 5.32976167e-01 -5.79257123e-04 1.66331862e-02]
[ -5.79257123e-04 1.97323981e-04 -1.74538139e-05]
[ 1.66331862e-02 -1.74538139e-05 5.54313034e-04]]
É possível ainda fazer ajustes mais complexos utilizando a função optimize.leastsq do modulo scipy permite mudar parâmetros do ajuste como tolerância (ftol, xtol, gtol), numero de interações (maxfev), etc e além disso tem um retorno mais complexo: infodict, mesg, ier, cov_x, vejamos o exemplo
              #!/usr/bin/env python

              import numpy as np
              import scipy.optimize
              import matplotlib.pyplot as plt

              def lorentz(parametros, x):
                a  = parametros[0]
                wo = parametros[1]
                c  = parametros[2]
                p = a / ((x-wo)**2 + c**2)
                return p

              def residuo(parametros, y, x):
                err = y - lorentz(parametros, x)
                return err


              chute = np.array([16.0, 5.0, 1.0])
              x, y  = np.loadtxt('dados_claudio.dat', unpack=True)

              resAjuste, cov_x, infodict, mesg, ier = scipy.optimize.leastsq(func=residuo, x0=chute, args=(y, x), maxfev=2000, ftol=1.0e-8, full_output=1)

              print "Resultado do ajuste: ", resAjuste

              print "\nMatriz de covarianca (diagonal eh varianca dos resultados):"
              print cov_x

              print '\nNumero total de interacoes:'
              print infodict['nfev']

              print '\nSolucao foi ontida se for 1,2,3 ou 4:'
              print ier

              print '\nCausa da falha ():'
              print mesg

              xnew = np.linspace(0.0, 12.0, 100)
              ynew = np.zeros(100)
              for i in range(len(xnew)):
                ynew[i] = lorentz(resAjuste, xnew[i])

              #mpl.rcParams['text.usetex'] =  True
              plt.plot(x, y, 'o', xnew, ynew, '-')
              plt.savefig( 'ajuste.png', format='png' )
              plt.show()
            
Que da como resultado:
[usuario@pclabfis: ]#python ajusteNaoLinear.py
Resultado do ajuste: [ 21.62608754 4.68830763 1.14494062]

Matriz de covarianca (diagonal eh varianca dos resultados):
[[ 8.17762583e+00 -8.88772952e-03 2.55208359e-01]
[ -8.88772952e-03 3.02760571e-03 -2.67799515e-04]
[ 2.55208359e-01 -2.67799515e-04 8.50500429e-03]]

Numero total de interacoes:
21

Solucao foi ontida se for 1,2,3 ou 4:
1

Causa da falha ():
Both actual and predicted relative reductions in the sum of squares
are at most 0.000000
Também imprime um gráfico igual ao 07 mostrado acima.