21 feb. 2013

Resolución de sistemas de ecuaciones lineales por descomposión QR usando Numpy

En Diciembre se  publicó un artículo donde se explica como resolver sistemas de ecuaciones.

Este artículo se basa de un artículo en Inglés QR descomposition with numpy .

Si se desea averiguar más sobre la descomposición QR se puede consultar a la página de wikipedia ó de la siguiente página.

Las ecuaciones que se usaron son:
3x+9y-10z  =   24
x-6y+4z      =   -4
10x-2y+8z  =  20

Donde Ax = b.
A = [[3 9 -10][1 -6 4][10 -2 8]] y
B  = [[24][-4[[20]]




El código se muestra a continuación:


#!/usr/bin/env python
# -*- coding: utf-8 -*-

from numpy import *


#Se define los valores de la matriz A
A = array([[3,9,-10],[1,-6,4],[10,-2,8]])

Q,R = linalg.qr(A) # qr decomposition of A

#Se definen los valores de la matriz B
b = array([[24],[-4],[20]])


#resolver Ax=b usando la funcion estandar numpy
x = linalg.solve(A,b)


#resolver Ax = b usando Q y R.
y = dot(Q.T,b)
xQR = linalg.solve(R,y) 

print "\nComparacion de Soluciones:"
print x.T,'Ax=b'
print xQR.T,'Rx=y'

Al ejecutar el script se tiene lo siguiente:
python qr.py  Comparacion de Soluciones: [[ 2.99029126  0.40776699 -1.13592233]] Ax=b [[ 2.99029126  0.40776699 -1.13592233]] Rx=y
Como se puede ver, la solución usando la función estándar de numpy y por la descomposición QR generan el mismo resultado.




20 feb. 2013

Creación de gráfico de burbujas con matplotlib

Los gráficos de burbuja presentan los datos como una serie de burbujas, el tamaño de las cuales es proporcional a la cantidad de datos. 

Un gráfico de este tipo resulta muy efectivo para mostrar el número de productos vendidos en cierta región.

Existe una herramienta que se llama Trendalyzer  desarrollada por la Fundación Gapminder, está herramienta permite convertir series de estadísticas en gráficos interactivos. 

Es Hans Rosling uno de los fundadores de Gapminder quien ha mostrado el potencial de la herramienta.
El siguiente vídeo es una demostración de la herramienta:

Con matplotlib en Python también se pueden crear gráficas de burbuja, el artículo se basa en un artículo en Inglés llamado "How to make Bubble charts with matplotlib". 

Los datos se obtendrán de un archivo csv, estos datos son de estadística del crimen en los Estados Unidos por estado del año 2005. El archivo csv se puede bajar aquí

La información en el archivo csv se tiene de la siguiente forma (se muestra las primeras líneas del archivo):


state,murder,forcible_rape,robbery,aggravated_assault,burglary,larceny_theft,motor_vehicle_theft,population
United States,5.6,31.7,140.7,291.1,726.7,2286.3,416.7,295753151
Alabama,8.2,34.3,141.4,247.8,953.8,2650,288.3,4545049
Alaska,4.8,81.1,80.9,465.1,622.5,2599.1,391,669488
Arizona,7.5,33.8,144.4,327.4,948.4,2965.2,924.4,5974834
Arkansas,6.7,42.9,91.1,386.8,1084.6,2711.2,262.1,2776221
California,6.9,26,176.1,317.3,693.3,1916.5,712.8,35795255



Como se ve la primera línea se tiene el título de cada columna, en la segunda línea se tiene la información total de Estados Unidos y luego se muestra la información por Estado.




Para efectos de programación la primera y segunda línea no son relevantes para la gráfica de burbuja.




A continuación el código del programa:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from pylab import *
from scipy import *

#Leer los datos desde el archivo csv
#Archivo original
#durl = 'http://datasets.flowingdata.com/crimeRatesByState2005.csv'
#Archivo modificado.
durl = 'http://127.0.0.1/descargas/crimeRatesByState2005-es.csv'
#Se captura los datos del archivo csv.
rdata = genfromtxt(durl,dtype='S8,f,f,f,f,f,f,f,i',delimiter=',')

rdata[0] = zeros(8) # se elimina los titulos 
rdata[1] = zeros(8) # se elimina la estadistica total de estados unidos (2da linea)
x = []
y = []
color = []
area = []

#Se captura los datos de cada fila
for data in rdata:
    x.append(data[1]) # asesinatos
    y.append(data[5]) # robo
    color.append(data[6]) # hurtos
    area.append(sqrt(data[8])) # poblacion
    # graficando las primeras 8 letras del nombre del estado
    text(data[1], data[5],data[0],size=11,horizontalalignment='center')
    

# se crea la grafica
sct = scatter(x, y, c=color, s=area, linewidths=2, edgecolor='w')
sct.set_alpha(0.75)
#Se define los ejes
axis([0,11,200,1280])
#Las etiquetas de cada eje
xlabel('Asesinatos por cada 100.000 habitantes')
ylabel('Robos por cada 100.000 habitantes')
show()

La siguiente figura muestra el resultado de la ejecución del programa:



19 feb. 2013

Manejar información de un archivo csv con csvkit

El programa csvkit es una herramienta desarrollada en Python que facilita la manipulación de la información contenida en un archivo con formato csv.

Se utilizará como ejemplo los datos de la página data.gov. Los datos que se utilizará son del Departamento de Asuntos de Veteranos de beneficios educativos de los Estatos Unidos (disculpen la traducción) del año 2009.

Es necesario instalar csvkit, en este caso se usa el comando pip de python:
pip install csvkit

Obtener los datos:
mkdir beneficios

cd beneficios 

Bajar archivo 2009.csv con el comando wget:
wget -O 2009.csv -U "Mozilla/5.0 (X11; U; Linux x86_64; en-US) AppleWebKit/534.16 (JHTML, like Gecko) chrome/10.0.648.205 Safari/534.16" http://www.data.gob/download/4029/csv 

Verificar las primeras 5 líneas del archivo:

head -n 5 2009.csv 

State Name,State Abbreviate,Code,Montgomery GI Bill-Active Duty,Montgomery GI Bill- Selective Reserve,Dependents' Educational Assistance,Reserve Educational Assistance Program,Post-Vietnam Era Veteran's Educational Assistance Program,TOTAL,

ALABAMA,AL,01,"6,718","1,728","2,703","1,269",8,"12,426",

ALASKA,AK,02,776,154,166,60,2,"1,158",

ARIZONA,AZ,04,"26,822","2,005","3,137","2,011",11,"33,986",

ARKANSAS,AR,05,"2,061",988,"1,575",886,3,"5,513",

Se puede usar el mismo comando wget para bajar los archivos de los años 2010.csv, 2011.csv y 2012.csv.

Obtener la información de las columnas con csvcut:

csvcut -n 2009.csv

  1: State Name

  2: State Abbreviate

  3: Code

  4: Montgomery GI Bill-Active Duty

  5: Montgomery GI Bill- Selective Reserve

  6: Dependents' Educational Assistance

  7: Reserve Educational Assistance Program

  8: Post-Vietnam Era Veteran's Educational Assistance Program

  9: TOTAL

 10: 

Como se puede observar, el archivo csv maneja 9 columnas. 

Con el comando csvcut se puede obtener información entre la fila 2 y 3 (State Abbreviate y Code), sólo se desea mostrar las primeras 5 líneas del archivo: 
csvcut -c 2,3 2009.csv | head -n 5
State Abbreviate,Code
AL,01
AK,02
AZ,04
AR,05


Se puede también manejar estadisticas bajo demanda con csvstat. Se genera la estadistica de la información de las columnas 1,4,9 y 10. Para este caso se utiliza csvcut y se pasa la información a csvstat:

csvcut -c 1,4,9,10 2009.csv | csvstat
  1. State Name
<type 'unicode'>
Nulls: True
Unique values: 52
Max length: 17
  2. Montgomery GI Bill-Active Duty
<type 'int'>
Nulls: True
Min: 435
Max: 34942
Sum: 325723
Mean: 6263.90384615
Median: 3548.0
Standard Deviation: 7537.86225373
Unique values: 52
  3. TOTAL
<type 'int'>
Nulls: True
Min: 768
Max: 46897
Sum: 506914
Mean: 9748.34615385
Median: 6520.0
Standard Deviation: 10070.4022127
Unique values: 52
  4. _unnamed
<type 'NoneType'>
Nulls: True
Values: 




Se puede realizar busquedas por filas con csvgrep. En este caso la información total del Estado de Illinois:

csvcut -c 1,"TOTAL" 2009.csv | csvgrep -c 1 -m ILLINOIS

State Name,TOTAL

ILLINOIS,"21,964"

Voltear orden de las columnas con csvcut: 
csvcut -c 9,1 2009.csv | head -n 5
TOTAL,State Name
"12,426",ALABAMA
"1,158",ALASKA
"33,986",ARIZONA
"5,513",ARKANSAS
En este caso se cambia el orden de las columnas al decirle a csvcut el orden de las columnas. 

Ordenar con csvsort: 
csvcut -c 9,1 2009.csv | csvsort -r | head -n 5
TOTAL,State Name
46897,CALIFORNIA
40402,TEXAS
36394,FLORIDA
33986,ARIZONA

Se puede dar un formato para que sea legible la información con csvlook:
csvcut -c 9,1 2009.csv | csvsort -r -l | csvlook
|--------------+-------+--------------------|
|  line_number | TOTAL | State Name         |
|--------------+-------+--------------------|
|  1           | 46897 | CALIFORNIA         |
|  2           | 40402 | TEXAS              |
|  3           | 36394 | FLORIDA            |
|  4           | 33986 | ARIZONA            |
|  5           | 21964 | ILLINOIS           |
|  6           | 20541 | VIRGINIA           |
|  7           | 18236 | GEORGIA            |
|  8           | 15730 | NORTH CAROLINA     |
|  9           | 13967 | NEW YORK           |
|  10          | 13962 | MISSOURI           |
|  11          | 13614 | COLORADO           |
|  12          | 13314 | OHIO               |
|  13          | 13011 | PENNSYLVANIA       |
|  14          | 12426 | ALABAMA            |
|  15          | 11492 | WASHINGTON         |
|  16          | 10085 | MARYLAND           |
|  17          | 9791  | MINNESOTA          |
|  18          | 9344  | MICHIGAN           |
|  19          | 9206  | OKLAHOMA           |
|  20          | 9013  | IOWA               |
|  21          | 8840  | WEST VIRGINIA      |
|  22          | 8757  | TENNESSEE          |
|  23          | 8081  | WISCONSIN          |
|  24          | 7872  | SOUTH CAROLINA     |
|  25          | 7809  | INDIANA            |
|  26          | 6652  | LOUISIANA          |
|  27          | 6388  | KENTUCKY           |
|  28          | 6009  | MASSACHUSETTS      |
|  29          | 5870  | OREGON             |
|  30          | 5513  | ARKANSAS           |
|  31          | 5511  | NEW JERSEY         |
|  32          | 5416  | NEBRASKA           |
|  33          | 5345  | UTAH               |
|  34          | 4947  | KANSAS             |
|  35          | 4551  | NEW MEXICO         |
|  36          | 4424  | PUERTO RICO        |
|  37          | 4299  | MISSISSIPPI        |
|  38          | 3728  | NEVADA             |
|  39          | 2997  | CONNECTICUT        |
|  40          | 2751  | IDAHO              |
|  41          | 2521  | HAWAII             |
|  42          | 1992  | SOUTH DAKOTA       |
|  43          | 1920  | MAINE              |
|  44          | 1795  | MONTANA            |
|  45          | 1778  | NORTH DAKOTA       |
|  46          | 1326  | NEW HAMPSHIRE      |
|  47          | 1175  | RHODE ISLAND       |
|  48          | 1158  | ALASKA             |
|  49          | 1145  | DELAWARE           |
|  50          | 1117  | WYOMING            |
|  51          | 1084  | DIST. OF COLUMBIA  |
|  52          | 768   | VERMONT            |
|  53          |       |                    |
|--------------+-------+--------------------|


Para finalizar se puede salvar el trabajo en un nuevo archivo csv:
csvcut -c 9,1 2009.csv | csvsort -r -l > 2009_ranking.csv

Si se desea aprender más de la herramienta csvkit se puede revisar la página de la documentación

En próximo artículo se mostrará como usar csvkit desde un programa en Python.

Teorema de muestreo explicado con numpy

El teorema demuestra que la reconstrucción exacta de una señal periódica continua en banda base a partir de sus muestras, es matemáticamente posible si la señal está limitada en banda y la tasa de muestreo es superior al doble de su ancho de banda. 
En el artículo se mostrará una tasa de muestreo a diferentes frecuencias, desde el valor doble a la frecuencia base, luego a un valor menor.

Algo más de teoría:
El teorema de muestreo de una señal continua que x (t) limitada en banda a B Hz pueden ser recuperados de sus muestras x [n] = x (n * T), donde n es un número entero, si T es mayor que o igual a 1 / (2B) sin pérdida de ninguna información. Y llamamos 2B la tasa de Nyquist.

El muestreo a una tasa inferior a la tasa de Nyquist se denomina submuestreo, se produce el efecto aliasing.

Si se desea más información sobre el Teorema de Muestreo se puede consultar a wikipedia.

Este artículo se basa en un artículo en Inglés "The sampling theorem explained with numpy" .

El código se muestra a continuación:


#!/usr/bin/env python



#De numpy se importa lo necesario para graficar la

#funcion seno

from numpy import linspace,sin,cos,pi,ceil,floor,arange



#De pylab se importa plot, show y axis. Lo necesario para crear

#la grafica

from pylab import plot,show,axis





#Muestreo de una seganl de ancho de banda 40 hz

# con velocidad de muestreo de 80 Hz

f = 40;  # Hz

#Tiempo minimo y maximo

tmin = -0.3;

tmax = 0.3;

#Se define el tiempo de la segnal.

t = linspace(tmin, tmax, 400);

#Se define la segnal de muestreo

x = cos(2*pi*t) + cos(2*pi*f*t)

#Se grafica el tiempo y la segnal.

plot(t, x)



# sampling the signal with a sampling rate of 80 Hz

# in this case, we are using the Nyquist rate.

#Muestreo de la segnal con una velocidad de muestreo de 80 Hz.

#Periodo de muestreo

T = 1/80.0;

#Tiempo minimo

nmin = ceil(tmin / T);

#Tiempo maximo

nmax = floor(tmax / T);

#Tiempo de la segnal.

n = arange(nmin,nmax);

#Segnal a la velocidad de muestreo

x1 = cos(2*pi*n*T) + cos(2*pi*f*n*T);

#Se grafica la segnal.

plot(n*T, x1, 'bo')



#Muestreo de la segnal con una velocidad de muestreo de 35Hz.

#Note que 35Hz esta por debajo de la velocidad de Nyquist.

T = 1/35.0;

nmin = ceil(tmin / T);

nmax = floor(tmax / T);

n = arange(nmin,nmax);

x2 = cos(2*pi*n*T) + cos(2*pi*f*n*T);

plot(n*T, x2, '-r.',markersize=15)







axis([-0.3, 0.3, -1.5, 2.3])

show()

La gráfica generada es la siguiente:

Con puntos azules se tiene el muestreo a 80Hz, con puntos rojos se tiene el muestreo a 35 Hz, se nota que el muestreo a 80 Hz es suficiente para capturar la oscilación de la señal.

En la siguiente gráfica se tiene un muestreo a 10 Hz que está por debajo de la frecuencia base de la señal (40 Hz).

Ahora se muestra la frecuencia de muestreo a 20 Hz:

Para terminar se muestra la frecuencia de muestreo a 30 Hz:

Para terminar se muestra la gráfica a una frecuencia de muestreo de 320 Hz:
Como puede notarse, mientras menor es la frecuencia de muestreo con respecto a la frecuencia base de la señal no se puede generar la señal original a partir de la muestra, mientras se va a aumentando la señal hasta llegar a la frecuencia base, se nota que se tiene más muestras para dicha recuperación pero sigue sin ser suficiente, es a partir del doble de la frecuencia base que la muestra puede ser generada.


18 feb. 2013

Encontrar el mínimo de una función usando fmin de scipy

Este artículo muestra como se encuentra el valor mínimo de una función con la función fmin de scipy.

El artículo es una versión en Español del artículo en Inglés "How to find the minimum of a function using fmin from scipy".

El código se muestra a continuación:


#!/usr/bin/env python



#Importar numpy,pylab y la función fmin de scipy.optimize.

import numpy

import pylab

from scipy.optimize import fmin



# Se define la funcion a partir de lambda.

rsinc = lambda x: -1 * numpy.sin(x)/x



#Se define un valor x0 de -5.

x0 = -5

#Se calcula el valor minimo de la funcion en el punto x0

xmin0 = fmin(rsinc,x0)



#Se define el punto x1 con valor -4

x1 = -4

#Se calcula el valor minimo de la funcion en el punto x1

xmin1 = fmin(rsinc,x1)



# se grafica la funcion.

x = numpy.linspace(-15,15,100)

y = rsinc(x)

pylab.plot(x,y)



#Se define el punto x0 en la grafica de la funcion

pylab.plot(x0,rsinc(x0),'bd',xmin0,rsinc(xmin0),'bo')

#Se define el punto x1 en la grafica de la funcion

pylab.plot(x1,rsinc(x1),'rd',xmin1,rsinc(xmin1),'ro')

pylab.axis([-15,15,-1.3,0.3])

pylab.show()



La figura a continuación muestra la gráfica de la función:

El punto azul es el mínimo encontrado al inicio desde diamante azul  (x= -5), el punto rojo es el mínimo encontrado iniciando desde diamante rojo. Este punto es el mínimo global encontrado en la función.

Además de la gráfica se genera la siguiente salida que genera la función fmin:


Optimization terminated successfully.
         Current function value: -0.128375
         Iterations: 18
         Function evaluations: 36
Optimization terminated successfully.
         Current function value: -1.000000
         Iterations: 19
         Function evaluations: 38



15 feb. 2013

Graficar la intercepción de 2 funciones

En este artículo se explica como encontrar la intercepción de 2 funciones, en este caso una función senoidal con una cosenoidal.

El artículo se basa en uno en Inglés que se llama "How to find the intersection of two functions".

La gráfica mostrará la intercepción de las funciones sin(x) y cos(x) desde el valor inicial de x igual a cero. Se utilizará la función fsolve la cual devuelve la raíz de una ecuación (no lineal ) definida por func(x) = 0.

A continuación se muestra el código:


#!/usr/bin/env python



from scipy.optimize import fsolve

#Importar pylab

import pylab

#Importar numpy

import numpy



#Se define la funcion que calcula la intercepcion

#de 2 funciones.

def EncontrarIntercepcion(fun1,fun2,x0):
    #Se usa la función lambda de la diferencia de ambas funciones con
    #valor inicial x0
    return fsolve(lambda x : fun1(x) - fun2(x),x0)









if __name__ == '__main__':

    #se calcula el resultado de la intercepcion de 2 funciones

    #sin y cos con valor inicial de 0.

    resultado = EncontrarInterseccion(numpy.sin,numpy.cos,0.0)

    #se genera el rango de -2 a 2.

    x = numpy.linspace(-2,2,50)

    #Se genera la grafica, 

    pylab.plot(x,numpy.sin(x),x,numpy.cos(x),result,numpy.sin(result),'ro')

    #se muestra la grafica.

    pylab.show()

La gráfica se muestra a continuación:
La función seno se encuentra en color azul,la gráfica del coseno en color verde y la intercepción es un punto de color rojo.

Enhanced by Zemanta

14 feb. 2013

Encontrar la raíz de una función con fsolve

Este artículo explica como calcular la raíz de una función utilizando la función fsolve.

El artículo se basa en un artículo en Inglés "How to find the rooots of a function with fsolve".

La función fsolve retorna la raíces de una ecuación no lineal definida por f(x) = 0.
Para este caso se calculará la raíz de la función f(x) = x^3.

A continuación se muestra el código:


#Import fsolve para calcular la raiz de la funcion x^3

from scipy.optimize import fsolve

#Importar pylab

import pylab

#importar numpy

import numpy



#se calcula la potencia 3 de x con la funcion lambda

potencia3 = lambda x : x**3



#Se calcula la raiz de x^3 iniciando con x = 10

resultado = fsolve(potencia3,10) # starting from x = 10

print resultado



#Se define 400 valores de x entre -4 a 4

x = numpy.linspace(-4,4,400)

#Se genera la grafica, pasando el valor de x

#la potencia 3era de x, el valor de resultado, la potencia 3era de resultado

pylab.plot(x,potencia3(x),resultado,potencia3(resultado),'ro')

#Se define el grid

pylab.grid(b=1)

#Se muestra la grafica

pylab.show()

La gráfica muestra el punto donde se encuentra la raíz de la función:




 .

Graficar una función de 2 variables con matplotlib

Este artículo se basa en el artículo en Ingles "How to plot a function of two variables with matplotlib" .

Se tendrá 2 gráficas de una función de 2 variables, la primera será la gráfica de intensidad y la segunda gráfica será una gráfica 3D.

A continuación se muestra el código:

#!/usr/bin/env python



#De numpy se explorta arange y exp

from numpy import exp,arange

#De pylab se importa meshgrid, cm, imshow, contour, clabel, clorbar, axis, title y show

from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show





from mpl_toolkits.mplot3d import Axes3D

from matplotlib import cm

from matplotlib.ticker import LinearLocator, FormatStrFormatter

import matplotlib.pyplot as plt



#Se define la funcion que se va a graficar

def z_func(x,y):

    return (1-(x**2+y**3))*exp(-(x**2+y**2)/2)



def graficaIntencidad(Z):

    #Se dibuja la funcion

    im = imshow(Z,cmap=cm.RdBu)

    

    #Se agrega el contorno de lineas con sus etiquetas

    cset = contour(Z,arange(-1,1.5,0.2),linewidths=2,cmap=cm.Set2)

    clabel(cset,inline=True,fmt='%1.1f',fontsize=10)

    

    #Se agrega la barra de colores a la derecha

    colorbar(im)

    

    #Se crea el titulo de la grafica con estilo latex

    title('$z=(1-x^2+y^3) e^{-(x^2+y^3)/2}$')

    #Se muestra la grafica

    show()



def grafica3D(X,Y,Z):

    fig = plt.figure()

    ax = fig.gca(projection='3d')

    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,

                           cmap=cm.RdBu,linewidth=0, antialiased=False)

    

    ax.zaxis.set_major_locator(LinearLocator(10))

    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink=0.5, aspect=5)

    plt.show()



if __name__ == '__main__':

    #rango de valores para x y y.

    x = arange(-3.0,3.0,0.1)

    y = arange(-3.0,3.0,0.1)

    

    #Se define la grilla de puntos

    X,Y = meshgrid(x, y)

    #Se evalua la funcion segun los valores de X y Y

    Z = z_func(X, Y)

    

    graficaIntencidad(Z)

    grafica3D(X,Y,Z)

Al ejecutar el script se mostrará la primera gráfica, al cerrarla aparecerá la segunda gráfica.

La gráfica de intensidad se muestra a continuación:

La gráfica 3D se muestra a continuación:

13 feb. 2013

Creación de grafos con networkx. Parte 3.

En el artículo anterior se muestra como crear grafos a partir del contenido de un archivo, en el primero se explica lo básico en la creación de grafos, en este se tocará el tema de analizar los grafos, de extraer información de los mismos.

El ejemplo a continuación se crea un grafo con varios nodos (de la "a" a la "z"), se enlazan todos los nodos, luego se muestra la cantidad de nodos, la cantidad de relaciones entre nodos, se le asigna un atributo a un nodo, un atributo entre 2 nodos, se averigua los vecinos de un nodo, el valor del atributo de un nodo, entre 2 nodos, se muestra también la ruta más corta para llegar del nodo m al b; cual es el promedio de la ruta más corta, se muestra la relación de los nodos con respecto a m, cual es la ruta más corta utilizando el Algoritmo de Dijkstra entre m y b. Para terminar se genera el grafo con varios métodos (espectral, aleatoria, circular y normal).

A continuación el código:

#!/usr/bin/env python



import networkx as net



#Se importa la libreria pyplot de matplotlib como plt



import matplotlib.pyplot as plt

#Se crea la instancia vacia del grafo

g=net.Graph()

#Se crean los nodos con sus relaciones

g.add_edge('a','b')

g.add_edge('b','c')

g.add_edge('c','a')

g.add_edge('a','d')

g.add_edge('f','d')

g.add_edge('g','f')

g.add_edge('h','b')

g.add_edge('i','h')

g.add_edge('i','g')

g.add_edge('j','c')

g.add_edge('k','j')

g.add_edge('l','k')

g.add_edge('m','l')

g.add_edge('k','h')

g.add_edge('i','d')

g.add_edge('f','k')

g.add_edge('m','g')

g.add_edge('n','m')

g.add_edge('o','m')

g.add_edge('p','o')

g.add_edge('q','h')

g.add_edge('r','q')

g.add_edge('s','r')

g.add_edge('t','a')

g.add_edge('u','t')

g.add_edge('v','u')

g.add_edge('j','b')

g.add_edge('w','f')

g.add_edge('x','w')

g.add_edge('y','i')

g.add_edge('z','y')

g.add_edge('n','p')

g.add_edge('z','x')

g.add_edge('z','v')

g.add_edge('s','x')

g.add_edge('p','v')

g.add_edge('r','u')



print "Nodos: ",g.nodes()

print "Relaciones: ",g.edges()



#Se define un tamagno al nodo a y un peso a la relacion entre el nodo a y b

g.node['a']['size']=1

g['a']['b']['weight']=1



print "Vecinos de m: ", g.neighbors('m')

print "Valor del nodo a: ",g['a']

print "Peso de la relacion entre a y b: ", g['a']['b']



print "Ruta mas corta entre m y b: ",net.algorithms.shortest_path(g,'m','b')

print "Promedio de la ruta mas corta ",net.algorithms.average_shortest_path_length(g)

print "Muestra la relacion de ruta mas corta entre pares de nodos relacionado con m: ", net.algorithms.all_pairs_shortest_path(g)['m']

print "Ruta mas corta usando el algoritmo de Dijkstra entre m y b:",net.algorithms.dijkstra_path(g, 'm', 'b')





#Se dibuja el grafo

net.draw(g)

#net.draw_spectral(g)

#net.draw_circular(g)

#net.draw_random(g)



#Se muestra el grafico

plt.show()


Para generar el grafo de forma espectral, circular o aleatoria se comenta el anterior y se descomenta el que se desea generar.

El grafo normal se muestra en la siguiente figura:
A continuación se muestra el grafo espectral:
A continuación se muestra el grafo circular:
A continuación se muestra el grafo aleatorio:
A parte de las gráficas se imprime lo siguiente en pantalla:


Nodos:  ['a', 'c', 'b', 'd', 'g', 'f', 'i', 'h', 'k', 'j', 'm', 'l', 'o', 'n', 'q', 'p', 's', 'r', 'u', 't', 'w', 'v', 'y', 'x', 'z']

Relaciones:  [('a', 'c'), ('a', 'b'), ('a', 'd'), ('a', 't'), ('c', 'b'), ('c', 'j'), ('b', 'h'), ('b', 'j'), ('d', 'i'), ('d', 'f'), ('g', 'i'), ('g', 'm'), ('g', 'f'), ('f', 'w'), ('f', 'k'), ('i', 'y'), ('i', 'h'), ('h', 'q'), ('h', 'k'), ('k', 'j'), ('k', 'l'), ('m', 'o'), ('m', 'l'), ('m', 'n'), ('o', 'p'), ('n', 'p'), ('q', 'r'), ('p', 'v'), ('s', 'x'), ('s', 'r'), ('r', 'u'), ('u', 't'), ('u', 'v'), ('w', 'x'), ('v', 'z'), ('y', 'z'), ('x', 'z')]

Vecinos de m:  ['o', 'l', 'g', 'n']

Valor del nodo a:  {'c': {}, 'b': {'weight': 1}, 'd': {}, 't': {}}

Peso de la relacion entre a y b:  {'weight': 1}

Ruta mas corta entre m y b:  ['m', 'l', 'k', 'h', 'b']

Promedio de la ruta mas corta  2.99

Muestra la relacion de ruta mas corta entre pares de nodos relacionado con m:  {'a': ['m', 'g', 'i', 'd', 'a'], 'c': ['m', 'l', 'k', 'j', 'c'], 'b': ['m', 'g', 'i', 'h', 'b'], 'd': ['m', 'g', 'i', 'd'], 'g': ['m', 'g'], 'f': ['m', 'g', 'f'], 'i': ['m', 'g', 'i'], 'h': ['m', 'g', 'i', 'h'], 'k': ['m', 'l', 'k'], 'j': ['m', 'l', 'k', 'j'], 'm': ['m'], 'l': ['m', 'l'], 'o': ['m', 'o'], 'n': ['m', 'n'], 'q': ['m', 'g', 'i', 'h', 'q'], 'p': ['m', 'o', 'p'], 's': ['m', 'g', 'f', 'w', 'x', 's'], 'r': ['m', 'g', 'i', 'h', 'q', 'r'], 'u': ['m', 'o', 'p', 'v', 'u'], 't': ['m', 'g', 'i', 'd', 'a', 't'], 'w': ['m', 'g', 'f', 'w'], 'v': ['m', 'o', 'p', 'v'], 'y': ['m', 'g', 'i', 'y'], 'x': ['m', 'g', 'f', 'w', 'x'], 'z': ['m', 'o', 'p', 'v', 'z']}

Ruta mas corta usando el algoritmo de Dijkstra entre m y b: ['m', 'g', 'i', 'h', 'b']

De esta forma se puede extraer información del grafo como la ruta más corta entre m y b (m,g,i,h y b) que se tiene el mismo resultado usando el Algoritmo de Dijkstra.

6 feb. 2013

Graficar el espectro de frecuencia con Numpy

En este artículo se explicará como generar una gráfica de una señal que se tiene en el dominio del tiempo a generar la gráfica en el dominio de la frecuencia.

Esté artículo se basa de un artículo en Inglés  llamado "How to plot the frequency spectrum with scipy" .

El análisis de frecuencia es el proceso de determinar la representación en el dominio de una señal de dominio del tiempo y comunmente empleando la Transformada de Fourier.

La Transformada Discreta de Fourier (DFT) se usa para determinar el contenido de frecuencias de las señales  y la Transformada Rápida de Fourier (FFT) es un método eficiente para calcular la Transformada Discreta de Fourier. Scipy implemente  FFT y se usará para analizar el espectro de frecuencias.

A continuación se muestra el código que gráfica la amplitud en el dominio de frecuencia de la función coseno del dominio del tiempo:


#!/usr/bin/env python

# -*- coding: utf-8 -*-

#Importa  coseno, linspace y pi

from numpy import cos, linspace, pi

#Importa plot, show, title, xlabel, ylabel y subplot para graficar

from pylab import plot, show, title, xlabel, ylabel, subplot

#Importa fft y arange

from scipy import fft, arange





#

def plotSpectrum(y,Fs):

 """

 grafica la amplitud del espectro de y(t)

 """

 n = len(y) # longitud de la señal

 k = arange(n)

 T = n/Fs

 frq = k/T # 2 lados del rango de frecuancia

 frq = frq[range(n/2)] # Un lado del rango de frecuencia



 Y = fft(y)/n # fft calcula la normalizacion

 Y = Y[range(n/2)]



 plot(frq,abs(Y),'r') # grafica el espectro de frecuencia

 xlabel('Frecuencia (Hz)')

 ylabel('|Y(f)|')





if __name__ == '__main__':

    Fs = 150.0;  # rata de muestreo

    Ts = 1.0/Fs; # intevalo de muestreo

    t = arange(0,1,Ts) # vector tiempo

    ff = 5;   # frecuencia de la señal

    y = cos(5*pi*ff*t)

    

    #Proceso de graficar la señal

    subplot(2,1,1)

    plot(t,y)

    xlabel('Tiempo')

    ylabel('Amplitud')

    subplot(2,1,2)

    #Se llama a la funcion con la señal y la rata de muestreo

    plotSpectrum(y,Fs)

    show()

Al ejecutar el script se genera la siguiente gráfica:

Obtener número de serial y número de versión de Raspberry Pi con Python

El Raspberry Pi tiene información de número del número de revisión del PCB, este número indica que se tiene pequeños cambios en el PCB. Estos cambios pueden afectar en el funcionamiento de los programas Python e incluso como está distribuido los pines (alimentación, entrada/salida, etc) en la placa.

Está información se puede obtener desde el archivo /proc/cpuinfo como se indica a continuación:


ernesto@raspberrypi ~ $ cat /proc/cpuinfo 

Processor	: ARMv6-compatible processor rev 7 (v6l)

BogoMIPS	: 697.95

Features	: swp half thumb fastmult vfp edsp java tls 

CPU implementer	: 0x41

CPU architecture: 7

CPU variant	: 0x0

CPU part	: 0xb76

CPU revision	: 7



Hardware	: BCM2708

Revision	: 000f

Serial		: 00000000bcd34f5e

El número de Revisión es 00f y el serial es 00000000bcd34f5e.
Las variantes del Raspberry Pi según número versión se tiene a continuación (tomado de acá):

  • Modelo B Revisión 1.0: El valor de cpuinfo es 0002.
  • Modelo B Revisión 1.0+ECN001(sin fusible, D14 se removió): El valor de cpuinfo es 0003.
  • Modelo B Revisión 2.0: El valor de cpuinfo es 0004,0005, 0006...
El sitio Raspberrypi spy tiene un par de scripts para obtener el número de revisión y el serial de la placa.

A continuación se muestra el script (con una pequeña corrección):
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Script que busca la informacion del serial y numero de revision en el Raspberry Pi.
Este script se obtiene de los siguientes enlaces:
Numero de Revision:
http://www.raspberrypi-spy.co.uk/2012/09/getting-your-raspberry-pi-revision-number-using-python/#more-574
Numero de Serial:
http://www.raspberrypi-spy.co.uk/2012/09/getting-your-raspberry-pi-serial-number-using-python/#more-570

"""



def getrevision():
    # Extrae la informacion del numero de revision del raspberry pi
    #Se asigna a revision un string de 4 ceros.
    revision = "0000"
    #Se abre el archivo cpuinfo
    #Se 
    try:
        f = open('/proc/cpuinfo','r')
        #Se recorre el archivo
        for linea in f:
            #Si existe el contenido Revision en una linea
            if linea[0:8]=='Revision':
                #Se toma la longitud de la linea
                longitud =len(linea)
                #Se agrega la informacion completa de la revision
                revision = linea[11:longitud-1]
        #Se cierra el archivo
        f.close()
#Si no abre el archivo se asigna 4 ceros a revision
    except IOError:
        revision = "0000"
    #Se retorna el valor de revision
    return revision
def getserial():
    # Extrae la informacion del serial desde cpuinfo
    #Se asigna un string con ceros a cpuserial
    cpuserial = "00000000"
    #Se intenta a capturar la informacion del archivo cpuinfo,
    #si no se tiene un mensaje de error.
    try:
        f = open('/proc/cpuinfo','r')
        #se recorre el archivo linea a linea
        for linea in f:
            #Si existe el contenido Serial en la linea
            if linea[0:6]=='Serial':
                #Se agrega la informacion completa del serial del cpu
                cpuserial = linea[10:-1]
        #Se cierra el archivo
        f.close()
    #Si no abre el archivo se asigna el valor de error al serial
    except IOError:
        cpuserial = "ERROR000000000"
    #retorna el valor del serial del cpu
    return cpuserial

if __name__ == '__main__':
    print u"El número de Serial del Raspberry Pi es: %s " %getserial()
    print u"El número de Revision del Raspberry Pi es: %s" %getrevision()


Al ejecutar el script se obtiene lo siguiente:

ernesto@raspberrypi ~ $ ./versionSerial.py 

El número de Serial del Raspberry Pi es: 00000000bcd34f5e 

El número de Revision del Raspberry Pi es: 000f

Ya con esta información se tiene que el Raspberry Pi es de la revisón 2.0.

A continuación se muestra la organización de los pines de la revisión 1.0:

Y la organización de los pines de la revisión 2.0:

Con esta información ya se conoce como realizar las conexiones de los circuitos con el Rasberry Pi.

5 feb. 2013

Instalación de Raspbian en un Raspberry Pi.

 Raspberry Pi es una placa de computadora de bajo costo (Más información en wikipedia) desarrollada por el Reino Unido.

En la siguiente figura se muestra un esquema de la placa:


A continuación 2 fotos del unboxing  del Raspberry Pi que me llegó hace una semana:



Para poner a funcionar la placa es necesario lo siguiente:

  • Un cable de monitor VGA a HDMI.
  • Una tarjeta de Memoria SD (mínimo 8 GB).
  • Un teclado y ratón USB.
  • Cable de red.
  • Cargador de Celular microusb.
  • Una versión de Linux (para este caso Raspbian una versión de Debian Wheezy para Arquitectura ARM).


Para bajar  Rasbian lo pueden hacer desde el siguiente enlace, el número SHA-1 para verificar la integridad del archivo es el siguiente: 514974a5fcbbbea02151d79a715741c2159d4b0a. 

El usuario del sistema es pi y la clave raspberry .


Para instalar la imagen del Raspbian en la memoria SD es necesario colocar la memoria en un equipo, descomprimir el archivo descargado y desde la línea de comandos ejecutar un dd (para este caso la memoria se encuentra en /dev/sdb, es necesario verificar la identificación de la memoria antes de ejecutar el comando):

dd dd bs=1M if=./2012-12-16-wheezy-raspbian.img of=/dev/sdb

1850+0 registros leídos
1850+0 registros escritos
1939865600 bytes (1,9 GB) copiados, 1610,06 s, 1,2 MB/s

El siguiente vídeo muestra el proceso de arranque de raspbian. 


 Vídeo donde se muestra la configuración del raspbian:


Para finalizar se muestra una foto del escritorio LXDE del Raspbian:

4 feb. 2013

Contador incremental con un display 7 segmentos y Arduino

Luego de explicar como se enciende y apaga de manera constante un LED en un Arduino, ahora se mostrará como conectar un Display 7 segmento de cátodo común.

El display incrementará de 1 a 9, luego pasa a valor 0 y vuelve a contar desde 1 a 9, a intervalos de 1 segundo.

Requerimientos:

  1. Arduino Uno.
  2. Protoboard.
  3. Cables.
  4. Display 7 segmentos de cátodo común.
  5. Resistencia de 220 Ohms.

En la siguiente figura se muestra la identificación de los LEDs del display 7 segmentos:


El esquema de conexión es el siguiente:


Conexiones:

  • Pin a - Digital 7
  • Pin b - Digital 8
  • Pin c - Digital 9
  • Pin d - Digital 10
  • Pin e - Digital 11
  • Pin f - Digital 12
  • Pin g - Digital 13


La siguiente figura se muestra el montaje del circuito:
El código del programa se muestra a continuación:


int pausa=1000; // define la pausa entre cada digito



#Se definen los Pines de salida al display 7 segmentos 

void setup()

{

  pinMode(7, OUTPUT);  

  pinMode(8, OUTPUT);

  pinMode(9, OUTPUT);

  pinMode(10, OUTPUT);

  pinMode(11, OUTPUT);

  pinMode(12, OUTPUT);

  pinMode(13, OUTPUT);

}



//Se define la funcion display que recibe 7

//variables y se asignan a cada una de sus salidas

void display (int a, int b, int c, int d, int e, int f, int g)

// Funcion del display

{

  digitalWrite (7,a);   

  digitalWrite (8,b);   

  digitalWrite (9,c);

  digitalWrite (10,d);

  digitalWrite (11,e);

  digitalWrite (12,f);

  digitalWrite (13,g);

}



//Funcion principal que genera un lazo continuo

void loop() 

// Dependiendo de cada dígito, se envía a la función display

// los estados (0 y 1) a cada uno de los segmentos

{

  display (1,1,1,1,1,1,0); //escribe 0

  delay(pausa);

  display (0,1,1,0,0,0,0); //escribe 1

  delay(pausa);

  display (1,1,0,1,1,0,1); //escribe 2

  delay(pausa);

  display (1,1,1,1,0,0,1); //escribe 3

  delay(pausa);

  display (0,1,1,0,0,1,1); //escribe 4

  delay(pausa);

  display (1,0,1,1,0,1,1); //escribe 5

  delay(pausa);

  display (1,0,1,1,1,1,1); //escribe 6

  delay(pausa);

  display (1,1,1,0,0,0,0); //escribe 7

  delay(pausa);

  display (1,1,1,1,1,1,1); //escribe 8

  delay(pausa);

  display (1,1,1,0,0,1,1); //escribe 9

  delay(pausa);

}

En la siguiente figura se muestra la conexión:


Para terminar, se muestra un vídeo con el contador funcionando: