Mostrando entradas con la etiqueta Python. Mostrar todas las entradas
Mostrando entradas con la etiqueta Python. Mostrar todas las entradas

domingo, 7 de junio de 2020

Optimización en Python con SciPy

Resolver problemas de optimización es una herramienta esencial para la ingeniería y la ciencia y debería ser una prioridad en los programas académicos. La bibliografía para aprender la teoría que podría sugerir es el capítulo 7 del libro gratuito Mathematics for Machine Learning. Por ahora me concentraré a detallar una nota para usar la función minimize() del modulo SciPy. Me estoy basando en gran parte en el ejemplo de APmonitor. No sólo será un traducción sino que añadiré detalles importantes que fueron omitidos en su tutorial.

En la jerga de optimización, se le llama función objetivo a la función que se desea optimizar alguno o algunos de sus argumentos. En este ejemplo la función objetivo a minimizar será:
La solución del problema de minimización consiste en encontrar los valores del vector x para los cuales la salida de la función objetivo sea mínima. En la práctica la función objetivo está sujeta a restricciones por lo que este ejemplo no será la excepción:

Restricción 1. El producto de todas las variables debe ser mayor o igual a 25.
Restricción 2. La suma de los cuadrados de las variables debe ser igual a 40.
Finalmente hay que considerar los intervalos que acotan a cada una de las variables:  
Para escribir el programa de Python usando el submódulo optimize de SciPy necesitamos saber lo siguiente. La función minimize() tomará, para este ejemplo, estos argumentos:

func. Una función que debe definirse con un único argumento vectorial y retornando un único escalar.
method. La lista de métodos de minimización disponibles aparecen en la documentación. Es importante saber que los únicos métodos que permiten optimización con restricciones son ''COBYLA", "SLSQP" y "trust-constr".
bounds. Lista de tuplas con las cotas para cada variable.
constraints. Lista de diccionarios donde cada uno tiene la esctructura:
  • 'type': Tipo de restricción. 'eq' para igualdad o 'ineq' para desigualdad
  • 'func': Nombre de la función que describe la desigualdad (debe crearse) 
Código completo:

lunes, 23 de marzo de 2020

Escribir un archivo binario desde Python y leerlo en C

En algunas ocasiones se requiere preparar un conjuntos de datos desde la comodidad de Python y terminar de procesarlos con la eficiencia de C o Fotran. Aquí un apunte sobre como pasar de un array de numpy de datos float32 a un binario para ser leído desde un programa en C. La escritura de binarios en Python por si misma es muy útil para exportar resultados de una manera eficiente en espacio si es que se maneja un volumen de datos inmanejable para un archivo de texto.

Python:
Para estructuras mas complejas de datos será necesario ingenier un sistema de bytes iniciales que den información de la estructura que contiene el archivo.  

C:

viernes, 14 de febrero de 2020

Atractor de Lorenz en Python usando odeint()

Este programa es un ejemplo de dos temas interesantes: el atractor de Lorenz y la resolución de un sistema de ecuaciones diferenciales ordinarias usando el modulo scipy. El sistema que resolveremos (el sistema de Lorenz) consta de 3 ecuaciones diferenciales acopladas (que funcionan de manera interdependiente). Nótese que he cambiado la notación que aparece en el articulo de Wikipedia a una forma vectorial:
Esta notación facilita el entendimiento de la sintaxis de la función  odeint(). Lo primero que debe hacerse, después de importar los módulos correspondientes, es contruir una función que tome como argumentos el vector u y t y retorne el vector de derivadas de u. En este caso particular las ecuaciones diferenciales no dependen de t (son autónomas), aún así es necesario poner como argumento a t porque lo exige la sintaxis de la función odeint(). Después se establece una condición inicial para el vector u y se construye una lista de valores para t sobre el intervalo en que nos interese tener a la solución. Finalmente se coloca todo lo anterior como argumentos del solver odint() como aparece en el código completo:

Esta será la gráfica resultante (El Atractor Extraño de Lorenz):


viernes, 6 de diciembre de 2019

Juego del gato (tic-tac-toe) vía algoritmo minimax.

Me dejaron el ejercicio de hacer un juego de gato y escogí el camino de la inteligencia artificial de la vieja escuela. Pensé que era algo trivial pero si me costó varios días. Planeo hacer en vacaciones un video sobre el algoritmo minimax en mi canal porque los que hay tienen códigos que ahora me parecen innecesariamente complejos. La verdad es que, aunque mis amigos creen que lo digo sólo por modestia, soy muy pendejo para entender las cosas a la primera y más aún si no me explican paso a paso. Ya que vi claro todo el proceso recursivo me acordé del final de la película War Games. Les debo el video pero por ahora pueden consultar mi reporte aquí. La implementación que hice hace uso de la siguiente clase para instanciar los estados posibles del juego:

Contiene un método que verifica si el tablero ya esta lleno y otro que verifica si es el turno del jugador X. El pseudocódigo es el siguiente:
[click para agrandar]


La función evaluación() retorna +10 si 'X' gana la partida, -10 si 'O' gana y 0 si hay empate o si aún no hay ganador en ese estado. El código completo en Python lo pueden encontrar aquí. Este programa sólo tiene el fin de ayudar a entender el algoritmo minimax. No es una implementación eficiente ya que resuelve el problema por fuerza bruta repitiendo incluso el paso por varios estados ya visitados. Pueden mejorarlo añadiendo estrategias como poda alfa-beta o memorizando estados ya visitados. 

domingo, 22 de septiembre de 2019

Registro de imágenes con SIFT en OpenCV

Se le llama registro al proceso de trasformación entre varios conjuntos de datos que permiten "encajar" o superponer patrones similares. Quizá el proceso más conocido sea la creación automática de fotos panorámicas a partir de varias fotos independientes. En esta entrada haremos un ejemplo más sencillo. Observen estas dos fotografías de un mismo objeto (uno de mis libros favoritos) tomadas en ángulos diferentes:
Lo que haremos sera colocar la imagen de la izquierda con la misma orientación y escala que la imagen de la derecha. Tomaremos de referencia las características comunes entre ambas imágenes (el libro principalmente). ¿Cómo encontramos estos puntos? En 2004 David Lowe publicó un algoritmo llamado SIFT (Scale-Invariant Feature Transform)  que permite encontrar puntos clave invariantes a las rotaciones y, en buena medida, a los cambios de iluminación. Los puntos clave con mejor correspondencia entre ambas imágenes se ven así:
Este algoritmo está patentado y no aparece en la instalación estándar de OpenCV. Deben seguir las instrucciones que puse en mi entrada anterior para poderlo utilizar. El código que escribí es una modificación para SIFT del código que pueden encontrar acá. Lo he reducido y comentado para mejor claridad:
La imagen final registrada es esta:

sábado, 21 de septiembre de 2019

Instalar OpenCV-Python con SIFT y SURF habilitados

Quiero ahorrarles algunas horas y frustración con esto. Estas notas están hechas para Linux (uso Debian) pero posiblemente funcione también en Windows (no lo garantizo). Voy asumir que ya tienen instalado Anaconda... Primero tienen que crear un nuevo entorno virtual (esto es por seguridad y no dañar su entorno principal de Python por si algo malo ocurre). Para crear un nuevo entorno:

$ conda create --name mi_entorno_OpenCV -y

Pueden nombrar su entrono como quieran pero es buena idea ponerle un nombre que haga referencia a su modulo. Si trabajan con IA seguramente necesitarán un nuevo para Tensorflow. Ahora, si es la primera vez que crean un entorno virtual tendrán que configurar su terminal para permitir que conda administre los entornos virtuales:

$ conda init bash

Si no usan bash debe poner la shell que utilicen. Bien, una vez que hayan hecho esto pueden activar su entorno con:

$ conda activate mi_entorno_OpenCV

Notarán que habra cambiado la etiqueta antes del promp y deberá aparecer el nombre de su entorno. Deberán instalar ahora matplotlib (su nuevo entorno solo tendrá numpy y scipy [si no, deberán instalarlos también]):

$ conda install -c conda-forge matplotlib

Spyder tampoco estará instalado en su nuevo entorno así que deberan instalarlo también:

$ conda install -c anaconda spyder

Recomiendo leer Managing enviroments de la documentación de conda. Ahora si, para instalar una versión de OpenCV que permita utilizar los algoritmos SIFT y SURF deberán instalar específicamente esta:

$ pip install opencv-contrib-python==3.4.2.16 

Listo. Terminada la instalación pueden hacer una prueba con los códigos de ejemplo de la documentación de OpenCV. Recuerden que deben mandar a llamar a Spyder desde la terminal con el entorno virtual que crearon para OpenCV. para salir del entorno utilicen:

$ conda deactivate

 En la próxima entrada pondré un ejemplo utilizando el algoritmo SIFT.

lunes, 15 de julio de 2019

Gráficas de Maplotlib en Tkinter

He estado aprendiendo usar el módulo TKinter para hacer interfaces gráficas en Python y me está encantando. A pesar de no tener un asistente gráfico como pyQT, se me ha hecho bastante fácil de usar. Para incrustar gráficos de Matplotlib en un GUI de Tkinter pongo el siguiente ejemplo en donde cada botón manda a graficar una función diferente:

  

martes, 2 de julio de 2019

Áreas, perímetros y centroides de imágenes binarias con scikit-image

Utilizando las funciones label() y regionprops() del submódulo measure de scikit-image podemos obtener varios parámetros muy útiles de regiones binarias contiguas (blobs binarios). En este ejemplo construiremos una imagen binaria con dos círculos de diferente tamaño de los cuales obtendremos sus propiedades.

Las salidas del programa son:

Región 1:
Area (pixeles): 15361
Perimetro (pixeles): 458.2741699796952
Centroide: (80.0, 100.0)


Región 2:
Area (pixeles): 7825
Perimetro (pixeles): 326.3919189857866
Centroide: (220.0, 220.0)

Es posible  construir un dataframe de pandas a partir del diccionario de propiedades generado con la función regionprops_table(). Si están comenzando con scikit-image pueden revisar mis notas de mi curso introductorio aquí.

lunes, 25 de marzo de 2019

Gráfica de lecturas de un DHT11 en Python

Esta nota corta será para no olvidar como hacer una gráfica en Python de un conjunto de mediciones de temperatura a lo largo del tiempo. El sensor utilizado en un DHT11 (que muy tardé me di cuenta que tenía muy baja resolución y que es mejor el DHT22). Generé un archivo de texto a partir de comunicación serial entre un PIC y una Raspberry Pi dónde cada linea tiene la siguiente estructura:

2019-03-22 23:11:33,22,46
fecha, temperatura, humedad relativa 

Algo muy útilil que aprendí a hacer es trasferir archivos entre un host local y la Raspberry a través de la red con un solo comando de consola. El comando es scp y su sintaxis es la siguiente para copiar del host local a la Raspberry :

$ scp /un/dir/local/archivo usuario@ipremota:/un/dir/remoto/archivo

Y de la Raspberri al host local:

$ scp usuario@ipremota:/un/dir/remoto/archivo /un/dir/local/archivo

Ya con el archivo en mi computadora, utilicé el módulo Pandas para leer el archivo .txt dónde registre cerca de dos dias y medio de datos de temperatura y humedad relativa. Aplique un filtro de mediana con el modulo scipy.signal para quitar los picos abruptos en los datos. Finalmente realice los ajustes necesarios para mostrar únicamente la hora en el eje x. El programa completo es el siguiente:

Gráfica:

lunes, 4 de febrero de 2019

Toolbox de Control para Python

Para aquellos que busquen un sustituto de código abierto para Matlab hay una buena noticia: Existe el módulo Python Control Systems Library. Acabo de probarlo y me parece mucho mejor que otras alternativas como los toolbox de control para Scilab y Octave. Si usan Anaconda, pueden instalarlo desde la terminal usando conda:

$ conda install -c conda-forge control

Quiero mostrarles algunos ejemplos simples utilizando el modelo de un motor NXT de Lego que describo en una entrada anterior. El modelo es el siguiente:
Creamos el objeto del modelo de la siguiente forma:

#Modelo como función de transferencia
num = [193.3]
den = [1.0, 18.735, 0.0]
sys = ctrl.tf(num,den) #Lazo abierto
print(sys)


El print nos mostrara la función de transferencia como una fracción en texto (al igual que en Matlab). Vamos a graficar ahora la respuesta al escalón de la planta en lazo cerrado:

#Modelo en lazo cerrado
sys_cl = ctrl.feedback(sys,1)

#Respuesta al escalón
t,y = ctrl.step_response(sys_cl)


#Graficas
plt.plot(t,y)
plt.grid(True)
plt.title(u"Respuesta al escalón")
plt.xlabel("segundos")


Lo que nos dará el siguiente resultado:
Para obtener el diagrama de bode:

plt.figure(2)
ctrl.bode_plot(sys)


Para visualizar el diagrama de polos y zeros:

from control.pzmap import pzmap 

plt.figure(3)
polos,zeros = pzmap(sys_cl), plt.grid(True)
print(polos[0])
print(polos[1])


¡Trataré de crear notebook's con más ejemplos y subirlos en mi GitHub!

lunes, 21 de enero de 2019

Leer imagenes .BIL del INEGI con Python

Los archivos de Banda Intercalada por Linea (BIL) es un formato de archivo binario en el que las distintas bandas de una imagen se concatenan formando una única matriz de datos. Los disponibles en la página del INEGI suelen tener una sola banda por lo que su extracción resulta aún más sencilla.

Para este ejemplo vamos utilizar fotografías aéreas. Vamos a ir a la pestaña de Datos, luego Mapas y finalmente a Topografía que nos mostrará este mapa. Seleccionamos la escala de 1:20,000 y nos aparecerán cuadriculas sobre el territorio y seleccionamos la que sea de nuestro interés. Al momento en que escribo esta entrada la página tiene un bug y muestra el mensaje de "No existen mapas relacionados con esta carta" cuando se da clic en cualquier zona y se podría pensar que no hay datos disponibles para ella, pero los datos aparecerán más abajo. Voy a utilizar una fotografía que cubre parte del municipio de Jojutla de Juarez, Morelos [archivo E14A69e]. Los archivos vienen en una carpeta comprimida ZIP. La información requerida para extraer los datos (dimensiones y formato) viene en un archivo .txt. En este caso el archivo es e14a69e.txt:

El programa de lectura y visualización es el siguiente:

La visualización de la imagen resultante es esta:

jueves, 7 de septiembre de 2017

Clasificador k-NN supervisado con Scikit-Learn

Scikit-Learn es un módulo para Python que incluye varias rutinas de clasificación, regresión y clustering entre otras herramientas matemáticas utilizadas en machine learning y minería de datos. En esta entrada únicamente trataré un ejemplo sencillo que permitirá introducir a la idea central detrás del clustering o análisis de grupos.

Imaginemos que existe una especie de conejos en la que las hembras tienden a ser grandes y tener pelaje de color gris claro mientras que los machos tienden a ser más pequeños y tener un pelaje más oscuro. Crearemos un programa que implemente un clasificador k-NN (k-nearest neighbors) utilizando el modulo scikit-learn. Este algoritmo de clasificación es muy sencillo. Su trabajo es dividir un espacio de n dimensiones (dónde cada una de ellas en la práctica representa un atributo o característica de un objeto) en N regiones (dónde cada una representa una clase de objetos). El modelo de división del espacio se basa en la razón de clases de los 'k' datos vecinos más cercanos a un punto de dicho espacio. Por ejemplo, para el caso en dónde solo existen n = 2 características y N = 2 clases (triángulos y cuadrados):
Para k = 3 (circulo solido), el elemento desconocido en verde será clasificado dentro de la categoría de los triángulos (2 triángulos vs 1 cuadrado). Para un k = 5 (circulo punteado), el elemento sera clasificado dentro de la categoría de los cuadrados (3 cuadrados vs 2 triángulos). El modelo de clasificación debe poder hacer predicciones para todos los posibles datos que caigan en cualquier punto del espacio de características. El proceso de generación de este modelo se conoce como entrenamiento. Existen dos caminos: supervisado, cuando el conjunto de datos de entrenamiento esta clasificado desde un principio para el algoritmo y no supervisado, cuando no se le dice algoritmo a que clase pertenece cada elemento del conjunto de entrenamiento. En nuestro ejemplo de los conejos utilizaremos un entrenamiento supervisado a partir de un conjunto de datos de entrenamiento que generé con la siguiente distribución (conejos.csv):
El color de los conejos está especificado por valores en escala de grises de 8-bits (0-255). El código se puede dividir en 4 secciones: carga y acondicionado de los datos de entrenmiento, instanciación del clasificador desde el módulo, entrenamiento del clasificador y finalmente la predicción de la clase de un elemento arbitrario. Una vez entrenado el clasificador, debe poder retornar una predicción de la clase a la cual puede pertenecer el vector [color,tamaño]. El código es el siguiente:

Para más detalles recomiendo revisar la documentación de sklearn.neighbors.KNeighborsClassifier.

jueves, 13 de julio de 2017

Ejemplo de programación multithread en Python

Como anexo a la entrada anterior, pongo un ejemplo del mismo programa pero en Python utilizando el módulo threadings. En este programa hago un cambio importante, estoy utilizando un producto punto del módulo numpy, una operación que está optimizada y escrita en C lo que me da tiempos de ejecución menores que los del programa anterior [Ref]:

Ejecutamos y medimos el tiempo de ejecución con:

time (python MatMulThreads.py)

Estos fueron los tiempos de ejecución en una Raspberry Pi 3 para matrices de 2000x2000:

Debe tomarse en cuenta que a pesar de que ambos programas son semanticamente iguales, su implementación no lo es. Es de notarse que en este caso no hay tanta mejoría entre 2 y 4 threads.

martes, 27 de junio de 2017

Segmentación de imágenes por umbral con scikit-image

Segmentar una imagen a través de umbrales es probablemente la forma más sencilla de identificar regiones de interés en una imagen. Esto hace que este método sea vulnerable a ruido, sombras y otros factores más a la hora de utilizarlo en la práctica. Sin embargo es una forma sencilla de introducir a las funciones de segmentación del módulo scikit-image. Primero, para este ejemplo tomaremos la siguiente pintura de Kazimir Malevich:
 Separaremos la imagen en secciones para las cuales obtendremos un array con los índices de cada pixel para cada una de ellas. Para ello cargamos primero la imagen:

#Cargar imagen
im = skimage.io.imread("malevich-map-style.png")
imGS = skimage.color.rgb2gray(im) #Escala de grises


 y creamos una imagen binarizada a partir de un umbral en la imagen en escala de grises:

#Umbral
thresh = 0.7
bw = binary_erosion(closing(imGS < thresh, np.ones((3,3))))


La variable bw es la imagen binarizada corregida y es obtenida aplicando la función closing seguida de binary_erosion (ambas contenidas en skimage.morphology). Closing() realiza una operación de dilatación seguida de una erosión, lo que permite cerrar "agujeros" oscuros en una imagen. El primer argumento que toma en este ejemplo es la imagen binarizada resultante de aplicar el umbral imGS < 0.7 y el segundo es una matriz que determina el tamaño de la dilatación/erosión. Binary_erosion es solo una operación de erosión que elimina las falsas regiones debidas a las zonas brillantes del fondo. Teniendo la imagen binarizada limpia, obtenermos los segmentos:

#Segmentación
labeled_im =
skimage.measure.label(bw)
lab_img = label2rgb(labeled_im, image=imGS,bg_label=0)


dónde labeled_im es una nueva imagen obtenida con la función label()en dónde cada pixel de la imagen en escala de grises original esta clasificado con su número de región . Si esto no es muy claro, pueden visualizar este arreglo de datos con el explorador de variables de Spyder dando doble click en el nombre de la variable:
La imagen lab_img es una visualización gráfica de las etiquetas de región sobre la imagen original que al mostrarla con Matplotlib y utilizando la función plt.text() junto con los datos contenidos en labeled_im, el resultado final es el siguiente:
 Código completo:

# -*- coding: utf-8 -*-
"""
Created on Tue Jun 27 13:27:53 2017

@author: rodolfo
"""

import numpy as np
from matplotlib import pyplot as plt
import skimage
from skimage.morphology import closing, binary_erosion
from skimage.color import label2rgb

#Cargar imagen
im =
skimage.io.imread("malevich-map-style.png")
imR = im[:,:,0]
imG = im[:,:,1]
imB = im[:,:,2]
imGS = skimage.color.rgb2gray(im)

#Umbral
thresh = 0.7
bw = binary_erosion(closing(imGS < thresh, np.ones((3,3))))

#Segmentación
labeled_im =
skimage.measure.label(bw)
lab_img = label2rgb(labeled_im, image=imGS,bg_label=0)

#Array de arrays de indices por segmento
reg_index = []
for i in range(0, np.max(labeled_im)):
    u = np.where(labeled_im == i)
    reg_index.append(zip(u[0],u[1]))


#Plots
R = []
cr = []
for i in range(0,np.max(labeled_im)):
    R.append(np.where(labeled_im==i+1))
    cr.append((np.mean(R[i][0]),np.mean(R[i][1])))

plt.figure(1)
plt.imshow(im)

plt.figure(2)
plt.imshow(lab_img)
for i in range(0,np.max(labeled_im)):
    plt.text(cr[i][1],cr[i][0],str(i+1),color='w')


sábado, 20 de mayo de 2017

Experimentando con autómatas celulares

Encontré por accidente el libro The New Kind of Science de Stephen Wolfram mientras buscaba libros de estadística en la biblioteca de la FCFM de la BUAP. Después de leer los primeros capítulos no pude evitar probar en Python algunos autómatas celulares. Es increíble el como con reglas tan simples se puede generar tanta complejidad. Sólo observen, tomemos una una linea de n casillas de largo. Cada casilla puede contener un '1' o un '0'. Supongamos que inicialmente la linea contiene solo '0' exceptuando la casilla n/2 que inicializamos con '1'. Nuestro pequeño autómata se encargará de generar una nueva linea de n casillas debajo de la inicial rellenando las casillas de la nueva linea de acuerdo a una función booleana cuyas entradas son casillas de la linea anterior. El autómata tiene un tamaño de 3 casillas y va rellenado las casillas de la nueva linea moviéndose una casilla a la vez hasta terminar todo el array. Un ejemplo para un autómata sería:

 dónde las casillas blancas representan '0' y las negras '1'. Consideremos un autómata un poco más complicado definido por el siguiente conjunto de reglas:

 Podemos expresar la regla anterior como una función booleana de la siguiente forma:
la cual podemos reducir mediante factorización a la forma:
El código en Python para este autómata celular queda:

import numpy as np
import matplotlib.pyplot as plt

n = 500
iline = np.zeros(n)
iline[n/2] = 1;
lines= [iline]

def rule(A,B,C):
    X = A and not(C) or not(A) and C
    return X

for i in range(0,n):
    line = lines[i]
    nline = np.zeros(n)
    for j in range(0,n-2):
        if rule(line[j],line[j+1],line[j+2]) :
            nline[j+1] = 1
    lines.append(nline)
   
lines = np.array(lines)

plt.imshow(lines, cmap = 'gray')


Y el espectacular resultado graficado es:
La figura resultante es un fractal conocido como Triángulo de Sierpinski y el autómata que lo genera se conoce como Regla 90.

domingo, 12 de marzo de 2017

Convertir coordenadas ecuatoriales a altazimutales (horizontales) con astropy

Algunos meses después de entender la geometría y hacer un programa para convertir entre estas coordenadas para mi tesis de licenciatura, descubrí que las conversiones podían hacerse con el módulo de astropy. No fue tiempo perdido pues entendí bastante bien el proceso. Para una conversión de muy buena precisión recomendaría usar astropy aunque sus funciones tienen el inconveniente de ser relativamente lentas (promedian una ejecución de 3 ms).

# -*- coding: utf-8 -*-
"""
Created on Sun Mar 12 20:28:24 2017

@author: rodolfo
"""
from datetime import datetime
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle

#Locación
Puebla = EarthLocation(lat=Angle('19d02m36s'), 
                       lon=Angle('-98d21m6.94s'), 
                       height=2152*u.m)
#Fecha y hora
utc_time = Time(datetimeutcnow(), scale='utc')
time = utc_time 

#Ecuatoriales
Sirio = SkyCoord(Angle('101d17m18.7469s'), 
                 Angle('-16d42m47.3141s'), 
                 frame='icrs')
#Horizontales
Sirio = Sirio.transform_to(AltAz(obstime=time,location=Puebla))

print "Altura: ", Sirio.alt.degree
print "Azimut; ", Sirio.az.degree

martes, 12 de julio de 2016

¿Qué tan lejos están las estrellas que podemos ver a simple vista?


Hace unos años se había abierto una discusión en reddit (también en Quora)a causa de la tira de xkcd de arriba. Hoy me encontré con una discusión similar ahora por un post en Facebook. Como ni las respuestas en reddit ni Quora me convencieron, decidí encontrar la respuesta yo mismo analizando algunos datos.

Primero, descargué la base de datos HYG (Hipparcos, Yale Bright Star, and Gliese) desde aquí, que nos da la información combinada de los 3 catálogos estelares mencionados, una muestra de cerca de 120,000 estrellas. Teniendo el archivo, escribí un programa en Python usando la librería Pandas [pueden bajarlo aquí] para analizar los datos. Y esto fue lo que encontré para la todas las estrellas resueltas [que pueden distinguirse como fuentes puntuales individuales] con una magnitud menor a 6:

De acuerdo al catálogo, el 98% de las estrellas que podemos ver a simple vista en la Tierra (muy lejos de las ciudades y sin luna) se encuentran a menos de 1000 años luz. Aún más interesante, las distancias de estas estrellas tienen una media de sólo 145 años luz con una desviación estándar de 136 años luz . Considerando los porcentajes en una distribución normal:
estamos hablando de que el 95% de esas estrellas se encuentran a menos de 417 años luz. Debe notarse también de que el 2% de las estrellas que se encuentran a más de 1000 años luz en este catalogo tienen una incertidumbre considerable y todas ellas se redondean a 10,000 años luz. En conclusión, cuando uno mira una estrella en la noche, lo más probable es que esté a menos de 400 años luz de distancia.

martes, 26 de abril de 2016

Máquina de Estados Finitos en Python

Una máquina de estado o autómata finito (Finite State Machine en inglés) es un modelo computacional bastante útil en robótica y automatización. Para este ejemplo haremos una implementación por software de una máquina de estados en Python que puede ser fácilmente adaptada para algún proyecto en una Raspberry Pi. Esta implementación es bastante sencilla, para FSM's más robustas basadas en clases recomiendo checar aquí y aquí. El diagrama de estados para FSM de esté ejemplo será el siguiente:
Este tipo particular de FSM es llamada Máquina de Moore ya que la transición de estado sólo depende del estado actual y del valor de la entrada. Esta implementación en Python aprovechará una característica peculiar de los diccionarios que permite asociar una función con una llave [ref]. De esta manera los estados quedan definidos como funciones dónde podemos incluir todas las acciones que nos interese realizar durante dicho estado (o durante las transiciones). El programa completo:

# -*- coding: utf-8 -*-
"""
Created on Thu Feb  4 23:36:28 2016

@author: Rodolfo
"""

from time import sleep
from random import randint

#Variable global
estado = 'i'

#Estados
def EDOi(entrada):
    global estado
    print('Estado Inicial')
    #Transiciónes
    sleep(2)
    if entrada == 0:
        estado = 'i'
    if entrada == 1:
        estado = 0
        print(u'Transición hacia 0...')

def EDO0(entrada):
    global estado
    print('Estado 0')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 1
        print(u'Transición hacia 1...')
    if entrada == 1:
        estado = 2
        print(u'Transición hacia 2...')
        
def EDO1(entrada):
    global estado
    print('Estado 1')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 0
        print(u'Transición hacia 0...')
    if entrada == 1:
        estado = 2
        print(u'Transición hacia 2...')

def EDO2(entrada):
    global estado
    print('Estado 2')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 1
        print(u'Transición hacia 1...')
    if entrada == 1:
        estado = 0
        print(u'Transición hacia 0...')

#Finite State Machine (FSM)   
def FSM(entrada):
    global estado
    switch = {
       'i':EDOi,
        0 :EDO0,
        1 :EDO1,
        2 :EDO2,
    }
    func = switch.get(estado, lambda: None)
    return func(entrada)

#Programa Principal
while True:    
 FSM(randint(0,1))
 sleep(2)


Salida del programa en la consola en Raspberry:
En caso de que te interese una implementación en lenguaje C, aquí tengo otra entrada dónde explico como hacerlo.

domingo, 10 de abril de 2016

Streaming de datos por TCP en Raspberry Pi + LabVIEW

 Si se desea hacer un streaming de datos desde un sistema de instrumentación remota como una estación meteorológica, sismológica o un telescopio robótico, puede resultar muy útil y económico tener un servidor TCP en un Raspberry Pi y visualizar los datos en tiempo real con LabVIEW desde cualquier lugar. Esto puede hacerse muy fácilmente utilizando el modulo socket para Python.

Un socket es una abstracción de programación para la representación de conexiones y permiten realizar una comunicación bidireccional. En Python se crea de la siguiente manera:

import socket
s = socket.socket()    #Instaciamiento para el objeto 's'
host = ''                      #Host local por defecto
port = 8006                #Puerto
s.bind((host,port))     #Asignación de dirección y puerto a la instancia
s.listen(5)                  #Número máximo de conexiones entrantes 

Para este ejemplo simularemos 2 señales que usaremos para enviar por TCP. En una aplicación real estas señales serían adquiridas mediate ADC's u otros dispositivos conectados a los pines GPIO de la tarjeta:

t = np.linspace(0,2*np.pi,100)
y0 = np.float16(np.sin(t))
y1 = np.float16(0.5*np.sin(t)*np.sin(2*t))

Se puede observar que estoy usando la función float16() para reducir la precisión de los datos para ahorrarle un poco de trabajo al cliente en LabVIEW pero pueden omitirla si requieren más dígitos significativos. Ahora, para parte principal del programa tenemos:

print u"Esperando conexión"
s, addr = s.accept()
print u"Conexión desde: ", str(addr)
i = 0
try:
 while True:
     s.send('CH0'+str(y0[i])+'\n')
     s.send('CH1'+str(y1[i])+'\n')
     i += 1
     if i>=99:
         i = 0
     time.sleep(0.01)
except KeyboardInterrupt:
s.close()

 Observen que en el argumento del método send() estoy concatenando las etiquetas de 'CH0' y 'CH1' para poder separar las señales después con LabVIEW. (recuerden que '\n' es una secuencia de escape que indica nueva linea).

El diagrama a bloques del cliente para LabVIEW es bastante sencillo:
(Click para agrandar)
En el bloque TCP Open Conection colocamos la dirección IP* de la Raspberry y el puerto que establecimos para el socket en el programa de Python. Después de hacer la lectura utilizamos los bloques Match Patern para identificar las etiquetas de los canales y separar sus valores. Finalemnte se convierten las cadenas aisladas a numeros de tipo double para enviarlos al graficador. El sistema funcionando se ve así:


*Aquí estoy usando una ip local. En la práctica se requerirá utilizar una IP pública. Hay varias formas de hacerlo: #1, #2 , #3

Referencias:
Python Advanced Tutorial 6, DrapsTV
Raspberry Pi + Arduino + LabVIEW, Alfredo Cruz

sábado, 20 de febrero de 2016

OpenCV en Python

Instalación de Anaconda en Linux:

Descargar instalador aquí y ejecutar con:

$ bash nombredelinstallador.sh

Instalación de OpenCV usando conda (en la terminal):

$ conda install opencv 

Abrir Spyder:

$ spyder

Ejemplo. Cargar una imagen y mostrarla con matplotlib:

import cv2 
from matplotlib import pyplot as plt 

img = cv2.imread('eco.jpg',0) #Cargar imagen 
plt.imshow(img,cmap = 'gray') # Mostrar imagen 
plt.xticks([]), plt.yticks([]) #Ocultar numeración en ejes 
plt.title('Umberto Eco (1932-2016)'