Multiplicación en punto fijo Pt. 2 [formato fixdt(signed,n,s,b)]

Esta entrada es una continuación del tema de multiplicación punto fijo y tocará el tema de un formato de mapeo lineal que permite ajustar números enteros a un rango lineal arbitrario como puede ser los rangos de voltaje estándar para DAC's y ADC's. El formato es el siguiente:

fixdt(signed,n,s,b) 

• signed: toma el valor de 1 si el entero almacenado es signado y 0 en caso contrario.
• n:  longitud en bits del entero almacenado. 
• s: pendiente de la función de mapeo lineal. 
• b: bias de la función de mapeo lineal.

En esta convención un número con parte entera y decimal se representa de la siguiente manera:

 El rango de este formato queda definido por sus parámetros de la siguiente manera:

 La multiplicación queda definida cómo:

 La multiplicación es mas complicada que en el formato fixed(signed,n,l) pero se simplifica cuando b = 0 siendo en este caso similar a la multiplicación en el otro formato. El siguiente paso en definir el procedimiento para computar los valores de n, s y b según el intervalo continuo en que el que deseemos trabajar. Al igual que el formato anterior ,el parámetro n define los valores máximos para los enteres almacenados. Por tanto, teniendo definidos los límites del intervalo real y del intervalo entero, los valores de s y b se calculan resolviendo el siguiente sistema de ecuaciones:
 

 Realizaremos un ejemplo práctico para demostrar la utilizad de este formato su ventaja con respecto a fixed(signed,n,l). Implementaremos en un FPGA un generador de señal cosenoidal (he escrito una entrada sobre los detalles de estos generadores). Utilizaremos una DAC LTC1661 que tiene una resolución de 10 bits y un rango de salida de 0-Vreff*(1023/1024). Las salidas Vcc de la tarjeta que voy a utilizar son de 3.3 volts por lo que este intervalo será aproximadamente de 0-3.3 V. Ya que esta DAC no puede dar salidas de voltaje negativo el valor para signed sera '0'. Con esta información podemos encontrar los valores para establecer nuestro formato. Podemos obtenerlos mediante este código de Matlab:

minimo = 0; %Volts
maximo = 3.3; % Volts
is_signed = false;
word_length = 10;

[E_min, E_max] = range(fi([], is_signed, word_length, 0));

A = double ([E_min, 1; E_max, 1]);
b = double ([minimo; maximo]);
x = A\b;
format long g
slope = x(1)
bias = x(2)

La ejecución del código nos da una pendiente s = 0.0032258064516129 y un bias b = 0. La señal para este ejemplo tendrá una frecuencia de 0.05 Hz y un periodo de muestreo de 1 segundo (si pueden usar un osciloscopio pueden probar con una frecuencia más alta y un periodo de muestro más corto). Generaremos un código VHDL con HDL Coder a partir de este modelo de Simulink:

Los subsitemas que he llamado Fixed-Point Multiplication contienene la implementación de la multiplicación que hemos descrito aquí. Dado que b = 0 la multiplicación de punto fijo simplemente se define como el producto de los enteros almacenados y la pendiente s:

 Algo que podría resultar confuso y difícil de entender al principio es el por qué el producto con s se realiza en formato fixed(signed,n,l). La respuesta es que ese formato se utiliza para principalmente para la multiplicación de números fraccionales en un entorno que solo permita trabajar con números enteros de forma interna (que es el caso de nuestro sistema de hardware). El formato de mapeo lineal fixdt(signed,n,s,b) sólo se encarga de que las operaciones entre enteros tengan el sentido físico que esperamos de forma externa. Esto tambien explica un segundo hecho que tambien puede resultar confuso: las operaciones entre enteros se realizan en modo signed (para facilitar las restas) aunque se ha definido el modo unsigned. 

Para generar el código nos vamos a Code > HDL Code > Generate HDL. Para la comunicación SPI con la LTC1661 he modificado un código que había publicado anteriormente. Los módulos son los siguiente:

• Divisor
• DAC_SPI
• Generador_Coseno_HDL (Generado por HDL Coder) 

El esquemático del sistema es el siguiente:

Nota: durante la generación del código se indica que el reloj para el modulo del generador sea el doble del establecido. En este caso el divisor tiene una salida de 2 Hz.

Según la simulación, y si implementamos correctamente las cosas, las secuencia de voltaje que deberíamos ver experimentalmente a la salida del DAC es la siguiente:

Aquí muestro un video donde se puede ver que no tengo un buen multímetro y tuve que aumentar a 4 segundos el periodo de muestreo para dar tiempo de estabilizar las lecturas:
Referencias
Compute Slope and Bias, MathWorks
Fixed-Point: Rage and Precision, MathWorks 

Generación de una señal coseno discreta

Podemos generar una función senoidal o cosenoidal de forma discreta mediante una ecuación de diferencias que es obtenida a partir de la transformada Z de una función senoidal o cosenoidal según sea el caso. Esta ecuación de diferencias es fácilmente implementable en un microcontrolador o FPGA. En este ejemplo crearemos un generador de señal cosenoidal. Partimos de la transformada Z de una función coseno discreto considerándola como la respuesta al impulso de la función de transferencia de un sistema Y(z)/X(z):

Dónde W0 es la frecuencia angular en tiempo continuo y Ts es el tiempo de muestreo [ver mi entrada sobre el concepto de frecuencia en tiempo discreto]. Multiplicamos ahora por z^-2/z^-2 para obtener un modelo físicamente sintetizable (debido a que los valores positivos en los exponentes de z implican muestras desde el futuro):

Multiplicamos ambos lados por el denominador de la función de transferencia y despejamos para Y(z):

Aplicamos transformada Z inversa para obtener finalmente nuestra ecuación de diferencias:

Para simular el generador utilizaremos Xcos de Scilab que es una alternativa de código abierto de Simulink [si no están familiarizados con esta herramienta aquí hay un video curso en español]. Antes de colocar los bloques debemos ir a menú superior en Simulation>Set Context y definimos las siguientes variables:

Ts = 0.1
W0 = 2*%pi/10
k0 = cos(W0*Ts)
k1 = 2*k0
 
El diagrama en Xcos es el siguiente:

Recordemos que hemos supuesto desde el principio que x[n] = delta[n]. La función Delta de Kronecker (impulso discreto) es un super bloque formado por el siguiente subsistema:

Dónde el retardo 1/z está inicializado a 1. Obsérvese que en el diagrama principal estamos enviando los resultados de la simulación al workspace. Personalmente, las gráficas nativas Scilab me parecen feas. Para hacer una grafica en Python exportaremos los datos de la simulación a un archivo CSV con el siguiente comando en la consola de Scilab después de ejecutar la simulación:

csvWrite([A.time,A.values],'cos_data.csv',',')

Teniendo el archivo podemos generar una gráfica en Python con el siguiente código:

import numpy as np
import matplotlib.pyplot as plt

A = np.loadtxt('cos_data.csv',delimiter = ',')

plt.stem(A[:,0],A[:,1])
plt.ylim([-1.5,1.5])
plt.grid(True)
plt.xlabel('Ts*n')
plt.ylabel('y[n]')
plt.show()

  Se deja como ejercicio construir un generador para función seno.

Multiplicación en punto fijo [formato fixdt(signed,n,l)]

Aunque al principio pueda chocar con la intuición, es posible realizar de forma digital una multiplicación de dos números con parte entera y decimal utilizando únicamente números enteros. De hecho, para aplicaciones de instrumentación científica o simulaciones numéricas dónde la precisión sea la prioridad es incluso aconsejable evitar el uso del famoso punto flotante. Recomiendo enormemente ver este capítulo del canal PBS Infinite Series para entender desde el punto de vista matemático el problema fundamental de los tipos double y single. En pocas palabras este problema radica en la en las limitaciones en la accesibilidad sobre la recta numérica computable. La resolución del punto variable no es homogénea a lo largo de todo su rango. Una segunda limitación de las operaciones de punto flotante tiene que ver a su costo de implementación en hardware; construir circuitos digitales de punto flotante consumen más recursos en un FPGA.

El formato de punto fijo con el trabajaremos en esta entrada es:

 fixdt(signed,n,l)

• signed : toma el valor de 1 si el entero almacenado es signado y 0 en caso contrario
• n : longitud en bits del entero almacenado 
• l : longitud en bits de la parte fraccional

En esta convención un número con parte entera y decimal se representa de la siguiente manera:

Los parámetros n y signed nos da el rango del formato:

A partir de la relación entre el Valor del Mundo Real (como se le llama en la jerga de punto fijo y denotaremos como Val) y el entero almacenado (que denotaremos por E) podemos definir la multiplicación de punto fijo de la siguiente forma:

Esto significa que la multiplicación se realiza como el producto de dos enteros mas un bitshift a la derecha de l locaciones. Para comprender mejor este formato realizaremos un ejemplo práctico concreto. Supongamos que queremos digitalizar una señal de voltaje sinusoidal de amplitud de 2 V con una frecuencia de 1 Hz mediante un ADC de 12 bits de resolución y un rango de entrada de 0-4 V. Supongamos que una vez digitalizada la señal queremos realizar un escalamiento de 0.321 y enviar esta nueva señal a un DAC de la misma resolución y rango. ¿Cómo resolveríamos el problema sin utilizar operaciones de punto flotante?

Bien, usaremos el formato fixdt(signed,n,l) para este ejemplo. Lo primero que debemos hacer es establecer los parámetros del formato. Ya que nuestro ADC no puede leer voltajes negativos, el parámetro signed será igual a 0. Dado que la resolución es de 12 bits entonces n = 12. Sabiendo que el voltaje máximo de lectura es 4 V encontramos el valor l de la siguiente manera:

Tenemos que el valor de l es 10 dando un rango de 0-3.99 V. Entonces el formato para nuestro ejemplo queda definido así: fixdt(0,12,10). El siguiente paso en convertir la ganancia de 0.321 en formato fixdt(0,12,10):

El equivalente de la ganancia resulta ser: 329. Podemos realizar ahora una simulación de nuestro ejemplo en Simulink y comparar el desempeño de la operación de punto fijo contra la de punto flotante:

Click para agrandar

El opam en el circuito añade el offset necesario para levantar la señal de la fuente y pueda ser adquirida por el ADC. Los bloques de ADC y DAC son simplemente subsistemas que realizan las operaciones de conversión de decimal a fixdt(1,12,10) y de muestro ZOH. Si bien Matlab/Simulink permite realizar las conversiones al formato fixdt(signed,n,l), decidí manejar uint16 para poder generalizar el ejemplo para ser implementado en microcontroladores y no sólo en FPGA's. Nuestro formato fixdt(0,12,10) cabe en un uint16 y eso también permite hacer una intuición más clara de que simplemente trabajamos con formatos enteros convencionales. La gráfica resultante con la comparativa es la siguiente:

Vemos que la correspondencia entre las operaciones de punto fijo y punto flotante son muy buenas. Podemos ahora a aventurarnos a escribir un codigo VHDL para nuestra operación de escalamiento de este ejemplo:

La simulación del módulo mostrando algunos valores de ejemplo para los enteros almacenados de entrada y salida es la siguiente:

 Click para agrandar

Para finalizar quisiera hacer un último comentario. Quizá notaron que tuvimos mucha suerte al calcular el valor del exponente l y nos dio un número que se ajustaba perfectamente al rango del ADC. No fue suerte, yo elegí el rango para que este ejemplo fuera más sencillo. En la practica los DAC's suelen tener rangos de 0-5.2 V o 0-3.3 V. ¿Que se hace en esos casos? Aquí es donde entra un nuevo formato: fixdt(signed,n,s,b). Pero hablaremos de él en otra entrada.

Referencias:

Fixed-Point: Rage and Precision, MathWorks 

Convertir un modelo de Simulink a código VHDL con HDL Coder [Ejemplo 2]

En una entrada anterior exponía el como generar un código VHDL de una función booleana sencilla. Este será también un ejemplo sencillo pero mucho más interesante. Generaremos un controlador P para la posición de un motor DC con caja de engranes y encoder [este modelo en particular] alimentado por un driver Pololu Qik 2s12v10. El esquema del sistema es el siguiente:

 La posición del motor será leída mediante el encoder que viene integrado en el motor. El modelo de Simulink que convertiremos a código VHDL es el siguiente:

Bloques requeridos:
HDL Coder > Commonly Used Blocks > In1
HDL Coder > Commonly Used Blocks > Out1
HDL Coder > Sources > Constant
HDL Coder >  Math Operations > Divide
HDL Coder > Discontinuities > Saturation
HDL Coder > Signal Attributes > Data Type Conversion

Diagrama:

 La operación de división por números que no sean potencias de 2 no está disponible en todos los dispositivos por lo que en este caso usaremos 1/32 como aproximación de 0.03 [es posible implementar en hardware divisiones entre números arbitrarios pero serán tratados en otro ejemplo]. Es necesario establecer en todos los bloques el tipo de dato int16 dando doble click a cada bloque y cambiando el tipo de dato la siguiente forma:

 También sera necesario dar doble click en el bloque Divide, ir a Signal Attributes y activar el modo "Saturate on integer overflow" y verificar que el modo de redondeo sea "zero" o "simplest". Habiendo hecho todo lo anterior generemos el código VHDL yéndonos a la pestaña Code > HDL Code > Generate HDL. Les generará el siguiente código:

Para complementar el controlador usaremos los siguientes módulos VHDL que ya he publicado por acá:

• Módulo UART
• Módulo Qik2s12v10
• Módulo Encoder X4 

Agregando todos los archivos .vhd requeridos al proyecto se puede proceder a conectar los bloques de forma esquemática en ISE de la siguiente manera:

Si no están familiarizados con la creación de proyectos esquemáticos en ISE aquí hay un video-tutorial en español muy bien explicado [utilizan Verilog pero el procedimiento es el mismo]. El bloque selector es simplemente una declaración when que envía una referencia de 0 si el switch esta apagado y de 2400 (90° en cuentas de encoder) cuando esta encendido. En este video muestro el funcionamiento del sistema implementado en una tarjeta Basys2:

Modoulo VHDL para un encoder de cuadratura en modo X4

El código de este módulo es sólo una modificación del código de Dr Dew (Scilab Ninja). La única diferencia es que le he agregado el conteo de ascenso y descenso para un registro de 16 bits. National Instrumentes tiene una nota en español muy completa sobre el manejo de encoders de cuadratura. La ventaja del modo X4 es que permite tener la máxima resolución angular posible para un encoder.

Parálisis

Acabo de tener la pesadilla/parálisis del sueño más extraña y aterradora hasta ahora. Escribiré los detalles mientras los tenga frescos. La fotografía que muestro es relevante y llegaré a ella en un momento. En la primera parte del sueño me encontraba con varias personas en un extraño tren de vapor que no tenía paredes ni techo ni piso, sólo una estructura o esqueleto en el teníamos que movernos con mucho cuidado para no caer. Teníamos que alimentar la caldera con carbón lo cual era muy difícil de hacer por la forma de los vagones. Empecé a cansarme y me que parado recargado en la estructura mirando el suelo y las vías moverse relativamente por debajo de mi. Comencé a hacer un ejercicio de mindfulness que estoy muy acostumbrado a hacer en la vida real. Esto hizo forzar el "renderizado" del sueño y los detalles eran cada vez más y más reales [esto ya me había ocurrido antes]. De pronto el tren se descarrila. Salgo volando y todo se apaga. En ese momento me di cuenta que estaba soñando y que había muerto en el sueño al caer del tren. Lo extraño era que ahora estaba en una especie de limbo en mi mente. De pronto pude ver una imagen. Era justo la de la fotografía, era lo que tenía frente a mi ojo izquierdo (el derecho estaba tapado por una almohada). Sabía ahora que estaba de regreso en el mundo real. Traté de moverme pero estaba paralizado. Estaba tranquilo porque ya me había pasado muchas veces antes. Espere un momento a que la parálisis pasara. Pero pasaron varios minutos sin que lo hiciera. Traté de moverme lo más fuerte que pude hasta que ocurrió algo muy extraño. Empecé a sentir mover mi brazo, pero era como mover un miembro fantasma. Más extraño aún era que podía tocar y sentir el celular y el libro, sentía que podía moverlos pero la imagen que veía o estaba congelada o los movimientos que hacía no eran reales. Estaba seguro que seguía despierto. Esperé varios minutos más a que este estado pasara, pero para mí horror no lo hacía. Traté de calmarme y pensar en qué hacer. Si mi vista se había congelado tal vez podía hacer una llamada y pedir ayuda. Lo intenté, sentía todos los movimientos de mi mano y el celular pero la imagen seguía congelada y no sabía lo que hacía. Me rendí y trate de pensar en otra cosa. Había comenzado a notar algunos cambios en mi vista con cierto retraso (como el "lag" en un vídeo en vivo). Se me ocurrió intentar dar la vuelta y sentí que lo hice. Unos segundos después la imagen cambio abruptamente al techo. Esto ya era demasiado extraño y de pronto me vino la idea de que tal vez algo me había pasado en la realidad, que tal vez mi casa se había venido abajo, tenía un daño cerebral y sólo estaba atrapado en este limbo esperando morir en cualquier momento. Había pasado cerca de media hora. Tenía mucho miedo, comencé a gritar y sacudirme para despertar pero no podía. De pronto volví a ver la imagen de la foto y ya estaba completamente despierto. Me levanté, verifiqué que podía moverme de nuevo. Lo único que se me ocurrió decir fue "the fuck was this?" y comencé a escribir esto para no olvidarlo. Ahora me pregunto si no estoy en una simulación y esto no fue más que un fallo en la matrix que quizá yo mismo induje al hacer mindfulness en un sueño.

 [Originalmente una nota de Facebook]

Módulo de control para un Qik 2s12v10 de Pololu en VHDL

Estoy empezando a aprender el HDL Coder Toolbox de Matlab/Simulink pero antes de empezar hacer pruebas con controladores PID en hardware necesitaba escribir el código de un modulo que me permitiera controlar el driver para motores DC Qik 2s12v10 que había comprado para mi montura altazimutal. El módulo esta diseñado para trabajar junto con este módulo UART del que había hablado en la entrada pasada. El esquema de la entidad del modulo es la siguiente:

Ver código

Esta versión del módulo solo puede enviar los valores de potencia en 7-bits signados y no puede leer la corriente que consumen los motores ni los mensajes de error. Espero en un futuro cree un módulo más completo. Pero la funcionalidad de esta versión es suficiente para realizar muchos proyectos. El siguiente diagrama muestra la implementación esquemática para una tarjeta Basys2 en dónde probé el módulo:

Aquí muestro en un video el funcionamiento del módulo:

Módulo UART en VHDL

Este módulo fue por el profesor Carlos García Lucero de la FCE-BUAP. Aunque debe estar disponible en la página de INTESC considero útil repostearlo por aquí para que tenga mejor accesibilidad para los buscadores. El código esta muy bien comentado y está implementado en un solo archivo vhd.

Ver código

Convertir un modelo de Simulink a código VHDL con HDL Coder [Ejemplo 1]

Para este primer ejemplo, generarémos el código VHDL para una función booleana sencilla de la forma X = AB xor C desde un diagrama en Simulink.

Bloques requeridos
HDL Coder > Commonly Used Blocks > In1
HDL Coder > Commonly Used Blocks > Out1
HDL Coder > Logic and Bit Operations > Logical Operators

Diagrama

Una vez terminado el bloque debemos irnos a la pestaña Code > HDL Code > Generate HDL. Si no hay ningún error deberá aparecer lo siguiente en la ventana de comandos:

### Generating HDL for 'Ejemplo_HDL'.
### Starting HDL check.
### Begin VHDL Code Generation for 'Ejemplo_HDL'.
### Working on Ejemplo_HDL as hdlsrc\Ejemplo_HDL\Ejemplo_HDL.vhd.
### Creating HDL Code Generation Check Report Ejemplo_HDL_report.html
### HDL check for 'Ejemplo_HDL' complete with 0 errors, 3 warnings, and 0 messages.
### HDL code generation complete. 

Por defecto el código generado aparecerá en el directorio  hdlsrc\. El código generado para este ejemplo es el siguiente:

Jojutla

De niño crecí en un lugar que, junto con la formación académica me mis papás, moldeó en gran parte mi amor por la ciencia. Vivía en una unidad habitacional que tiene una vista muy amplia hacia el valle del suroeste de Morelos con dos volcanes en el fondo. El paso de las estaciones es muy notorio en esa zona de México gracias a los cambios en el color de la vegetación en los cerros, el ir y venir de los glaciares de los volcanes y las "lluvias" de tisne del final de la zafra azucarera. La cerros también me daban una referencia para darme cuenta que el Sol salía en un lugar diferente a lo largo del año, efecto que noté muy bien durante la formación de las mañanas en la secundaria. Las fumarolas del Popocatepetl eran muy comunes y siempre visibles por el clima. Las tormentas eléctricas también son particularmente intensas en esa región, algo que sólo noté hasta que me mudé a otros estados. Hasta la fecha no he estado en un lugar de México en dónde haya tantas especies de plantas y animales en un mismo patio como en ese lugar. Jojutla siempre me recordó que la Tierra es un lugar vivo y de constante cambio. Este mes lo ha hecho más que nunca y nos ha dado a muchos un recordatorio de lo frágiles que somos los humanos ante el poder de la naturaleza

* * *

El 26 de septiembre de 2007 publiqué la primera estrada en este, mi primer y único blog continuamente activo. Me hubiera gustado que mi entrada de X aniversario tuviera un tema menos sombrío. El de 19 de septiembre de este mes ocurrió uno de los sismos más devastadores en la historia de mi país. Mis primeras entradas hablaban de la región, ahora en buena parte en ruinas, en la que pasé toda mi infancia. Nunca nos pasó por la mente el vivir en Jojutla un desastre tan grande. Los sismos siempre fueron frecuentes pero nunca pasaron un breve susto. De niños crecimos haciendo colectas de vivieres hacia lugares que nos parecían muy lejanos. Este fue un tema recurrente entre mis amigos de la infancia con quienes estuve apoyado como brigadista de remoción de escombros desde el martes pasado hasta el domingo. Caminar por aquella ciudad ciudad tan querida reducida a una zona de guerra me provocó una sensación extraña de irrealidad. Al mismo tiempo el ver el compromiso de tanta gente me animó de forma inesperada. Por ahora las emergencias ya están atendidas y solo queda trabajar en las demoliciones y reconstrucción. Aún hay un largo camino por trabajar. Me hubiera quedado más tiempo pero desafortunadamente tenía acceso a internet sólo por breves momentos y no me dí cuenta que la BUAP reanudaba labores hasta el 2 de octubre. Sólo tenía que regresar por mis tramites de titulación. Posiblemente regrese a apoyar este jueves.  

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.

De finales y comienzos

Nunca pensé que terminar un tratamiento psicológico se sintiera como terminar un libro, juego o una serie que te gustó mucho. Te libraste del cruel villano después tantas batallas, finalmente, está hecho, pero ahora sabes que no volverás a ver a quienes te acompañaron en esta larga campaña. Quizá haya un nuevo enemigo y los antiguos combatientes deban ser reunidos de nuevo en el futuro. Pero por ahora camino sobre las ruinas del régimen vencido con una mezcla de libertad y soledad en el pecho... Pero me si me detengo a pensarlo, ¿no es este el inicio de una historia completamente nueva?

He Has Left Us Alone But Shafts of Light Sometimes Grace the Corners of Our Rooms

Tenía bastante que no disfrutaba tanto un disco de Post-Rock. El disco debut de Thee Silver Mt. Zion Memorial Orchestra. Poniéndome un poco en contexto, acabo de emborracharme solo con una Calavera Witbier mientras lo escuchaba. La lluvia de verano en Puebla aún suena allá afuera. Puedo ahorrarme una larga reseña de por qué este disco en una obra de arte si únicamente les recomiendo salir a buscar unas buenas cervezas para acompañarlo. Acompañarlo en soledad y con el mejor equipo de sonido que tengan. Descargar

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.

Programación multithread en C (pthread.h)

Un thread o hilo de ejecución es el subconjunto mínimo de instrucciones dentro de un programa que pueden ser ejecutadas secuencialmente [existe el caso particular en el que un solo hilo de ejecución representa al programa entero]. Este subconjunto de instrucciones viene asociado con los registros y locaciones de memoria que son usados durante su ejecución. Si un programa está divido en varios hilos de ejecución, cada uno de ellos puede ejecutarse de forma concurrente (o de forma paralela si se cuenta con múltiples procesadores). Es importante tener claros los conceptos de proceso e hilo de ejecución cuando se va hacer programación paralela. Recomiendo ver los videos sobre hilos de ejecución y gestión de procesos de la UCAM dónde resumen muy bien los temas en menos de 10 minutos (español). De manera conceptual, podemos representar la relación entre los hilos de ejecución (threads) y un proceso de la siguiente manera:

Se conoce como Pthreads o POSIX threads a un modelo de ejecución estandar (independiente del lenguaje de programación) para sistemas Unix. Esta definido como un conjunto de rutinas y tipos en C implementados en el header: pthread.h.

Si un programa en C/C++ no crea nuevos hilos explícitamente, será ejecutado como un proceso de hilo único. Para la crear un nuevo thread con la librería pthreads utilizamos la función pthread_create:

int pthread_create(thread_ID,att,rutina_concurrente,argumento) 

thread_ID: apuntador al identificador de hilo (tipo pthread_t)

att: apuntador a la estructura que contiene los atributos del hilo. NULL para atributos por defecto.

rutina_concurrente: puntero a la función que contiene la rutina que va a ejecutar el hilo.

argumento: argumento opcional a pasar hacia rutina_concurrente (tipo void*). NULL si ninguno es requerido.

Pueden encontrar una descripción más detalla del resto de las funciones de la librería aquí.

Ejemplo: Multiplicación matricial 

En esta entrada quiero hacer énfasis en la creación de hilos y ejecución. Como ejemplo de esto, probaremos un programa que realice una multiplicación de matrices cuadradas de forma paralela. Consideremos el algoritmo general para la multiplicación de matrices:

 Una forma de paralelizar el producto matricial es dividir las operaciones entre filas y columnas entre los hilos (no es la forma más eficiente pero será muy útil para aprender a dividir un programa en hilos independientes). El número total de productos punto entre n filas y  m columnas es igual a n*m. Si usamos matrices de 16x16, el total de operaciones entre filas y columnas sería de 256. Si dividimos estas operaciones en partes iguales entre 4 hilos de ejecución, cada hilo llevará a cabo 64 operaciones. Podemos hacer esto dividiendo las filas de la matriz A a lo largo del índice 'i'. Un ejemplo para matrices 8x8 repartido en dos hilos:

Cada hilo en este caso realiza 8 operaciones de dot(fila,columna) del total de 16. Con 4 hilos, cada uno accedería a una única fila de la matriz A realizado 4 operaciones. La rutina que pasaremos como argumento a la función pthread_create() es la siguiente:

void *matmul(void* id_arg){
  int i,j,k;
  long  id = (long) id_arg;
  int rows_per_thr = col/num_of_threads;
  int start_index = id*rows_per_thr;
  int final_index = (id+1)*rows_per_thr;

  for(i=start_index;i < final_index;i++){
   for(j=0;j < col;j++){
    for(k=0;k < row;k++){
      C[i][j] += A[i][k]*B[k][j]; 
    }
   }
  }
}

Todos los hilos ejecutarán exactamente la misma rutina pero los datos a los que accederán (elementos de 'A','B' y 'C') dependieran del identificador de cada hilo. Ya que los hilos tienen acceso a la memoria compartida del proceso principal podrán leer directamente las locaciones de memoria de las variables globales del programa. Sólo un identificador de hilo será pasado como argumento para la rutina paralela. La creación de threads se hace dentro de main() de la siguiente manera:

  pthread_t tid[num_of_threads];
  long rank;

  //Creación de threads
  for (rank = 0; rank < num_of_threads; rank++)
     pthread_create(&tid[rank], NULL,matmul , (void*) rank);

En este punto del código los hilos comienzan a ejecutar independiente la rutina que se les ha asignado. Sin embargo aquí hay un problema: el programa principal continuara con su ejecución sin esperar a que los hilos terminen sus tareas. Para resolver esto es necesario indicarle explícitamente al programa principal esperar a cada uno de los hilos con la función pthread_join():

  //Unión de threads
  for (rank = 0; rank < num_of_threads; rank++)
      pthread_join(tid[rank], NULL);

Nota: en este ejemplo se considera que la mayoría de las implementaciones de Pthreads los hilos tienen el atributo de joinable por defecto. Para mayor portabilidad se recomienda activar este atributo explícitamente. Pueden encontrar un ejemplo de esto aquí.

Para simplificar el código, las matrices serán leídas en nuestro programa desde un archivo de texto. Creamos este archivo con dos matrices de enteros aleatorios (y una matriz de resultado para verificar la salida del programa en C) con el siguiente código en Python:

import numpy as np
n= 16
A = np.random.randint(0,9,size = (n,n))
B = np.random.randint(0,9,size = (n,n))
C = np.matmul(A,B)
namein = "matext"+str(n)+"x"+str(n)+".txt"
nameres = "resultado"+str(n)+"x"+str(n)+".txt"
np.savetxt(namein,np.concatenate((A,B),0),fmt = "%d",delimiter = " ")
np.savetxt(nameres,C,fmt="%d",delimiter= " ")

El código completo en C es el siguiente:

Compilamos y ejecutamos desde terminal con:

 gcc MatMulThreads.c -lpthread ; ./a.out | less 

Para verificar la reducción en el tiempo de ejecución, corrí el programa en una Raspbery Pi 3 que tiene un quad-core ARM Cortex-A53. Generé dos matrices de 1000x1000 y estos fueron los resultados variando el número de threads:

Este programa tiene varias limitaciones, una de ellas es que sólo admite números pares de hilos. Se deja como ejercicio modificarlo para cualquier número de hilos.  

Bien, hasta aquí se han cubierto los aspectos de creación y unión de threads. Existen aún varios temas importantes más como la sincronización de hilos, especialmente en los accesos de memoria (en el ejemplo de la multiplicaciones de matrices no existen conflictos de memoria ya que cada hilo trabaja con una parte diferente de los arreglos). Estos temas son tratados de forma consista en este artículo y en estas diapositivas de la Universidad de Buenos Aires. 

[Escribí el mismo ejemplo de esta entrada en Python utilizando el módulo threading. Pueden revisarlo aquí.]

Referencias
POSIX Threads Programming, Blaise Barney, Lawrence Livermore National Laboratory

Instalación de GNU Scientific Library (Debian)

GNU Scientific Library (GSL) es una librería numérica para C/C++, al parecer, poco conocida fuera del ámbito académico especializado. Tiene un repertorio de rutinas comparable con el popular módulo SciPy para Python. Encuentro muy útil GSL para sistemas con Linux embebido en robótica o IoT, especialmente en aplicaciones que requieran la operación de estos sistemas de tiempo real. Recomendaría echarle un vistazo si trabajan en proyectos serios con Raspberry Pi's o UDO's que requieran cálculos numéricos avanzados. Su documentación está bien estructurada y es muy amigable. El único inconveniente es que no hay muchos ejemplos y ayuda ahí afuera (ni siquiera en Stack Overflow) así que habrá que leer cuidadosamente el manual de referencia.

Para instalar GSL en Debian y derivados:

sudo apt-get install libgsl0ldbl libgsl0-dev 

Probaremos la instalación con el siguiente ejemplo con operaciones de números complejos (revisar documentación):

#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

int main (void){
 
  //Declaración de variables complejas
  gsl_complex z,w,u;
 
  //Asignación de valores
  GSL_SET_COMPLEX(&z,3,4);  //z = 3+4i
  GSL_SET_COMPLEX(&w,4,5);  //w = 4+5i

  //ejemplos de operaciones
  u = gsl_complex_add(z,w);
  printf("z + w = %2.f + i%.2f\n",u.dat[0],u.dat[1]);
  u = gsl_complex_mul(z,w);
  printf("z*w = %.2f + i%.2f\n",u.dat[0],u.dat[1]);
  u = gsl_complex_log(z);
  printf("log(z) = %.2f + i%.2f\n",u.dat[0],u.dat[1]);
  u = gsl_complex_exp(gsl_complex_add(z,w));
  printf("exp(z+w) = %.2f + i%.2f\n",u.dat[0],u.dat[1]);
  return 0;
}

Compilamos y ejecutamos desde terminal con:

gcc GSL_Complex_demo.c -lgsl -lgslcblas ; ./a.out

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')

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.