El laberinto de Falken

Prompt Engineering, o cómo hablar con una IA.

noreply@blogger.com (Alberto García Serrano) — Fri, 24 May 2024 07:17:00 +0000

En la era de la inteligencia artificial, el prompt engineering se ha convertido en una habilidad crucial para interactuar eficazmente con modelos de lenguaje como GPT-4. Se trata de diseñar y formular instrucciones claras y efectivas (prompts) para obtener respuestas útiles y precisas de estos modelos generativos.

¿Qué es el Prompt Engineering?

El prompt engineering persigue crear instrucciones detalladas y específicas para modelos de lenguaje con el fin de guiar sus respuestas de manera más efectiva. A medida que los modelos de IA se vuelven más sofisticados, la necesidad de proporcionar instrucciones precisas se vuelve esencial para:

Obtener Respuestas Precisas: Prompts bien diseñados ayudan a obtener respuestas más precisas y relevantes de los modelos de IA.
Mejorar la eficiencia en la Comunicación: Minimiza la necesidad de múltiples iteraciones para obtener la información deseada.
Personalización: Permite adaptar las respuestas de la IA a contextos y necesidades específicos.

Principios básicos del Prompt Engineering

Claridad y especificidad

Para obtener la mejor respuesta posible, es esencial que el prompt sea claro y específico. Un prompt ambiguo puede llevar a respuestas confusas o irrelevantes.

Ejemplo:

Prompt vago: "Háblame sobre los gatos."
Prompt claro: "Describe las características físicas y el comportamiento típico de los gatos domésticos."

Contexto relevante

Proveer contexto suficiente permite que la IA comprenda mejor la solicitud y genere respuestas más pertinentes. Por ejemplo:

Sin contexto: "Explícame la teoría de la relatividad."
Con contexto: "Explica la teoría de la relatividad de Einstein de manera simple, como si se la estuvieras contando a un estudiante de secundaria."

Una técnica útil consiste en indicar al LLM cómo debe comportarse:

Indicar rol: "Compórtate como un profesor de secundaria."
Prompt: "Explícame cómo se resuelve una ecuación paso a paso. Pon un ejemplo."

Instrucciones paso a paso

Dividir las solicitudes complejas en pasos más pequeños y manejables puede ayudar a la IA a abordar cada parte de manera ordenada, ya que le ayudará a organizar la información de la forma que más nos convenga.

Ejemplo:

Instrucción compleja: "Escribe un ensayo sobre la revolución industrial."
Instrucción paso a paso:
1. "Define qué fue la revolución industrial."
2. "Describe las principales innovaciones tecnológicas de la revolución industrial."
3. "Analiza el impacto social y económico de la revolución industrial."

Uso de ejemplos

Si le damos ejemplos dentro del prompt podremos guiar a la IA en la dirección correcta y establecer el tono o formato esperado de la respuesta.

Ejemplo:

Sin ejemplo: "Escribe un poema sobre el otoño."
Con ejemplo: "Escribe un poema sobre el otoño similar a este: 'Las hojas caen suavemente, doradas bajo el sol, el viento canta susurrante, anunciando el frío en unión.'"

Herramientas y técnicas avanzadas en Prompt Engineering

Además de seguir los principios básicos, existen técnicas avanzadas y herramientas que pueden mejorar aún más el proceso de prompt engineering.

Uso de plantillas

Las plantillas predefinidas son extremadamente útiles para tareas repetitivas o formatos específicos de respuestas. Utilizar plantillas asegura que las respuestas de la IA sean consistentes y estructuradas. Muy útil para usar los LLMs en entornos de trabajo.

Ejemplo: plantillas para recetas de cocina

Una plantilla puede ayudar a organizar la información de manera coherente y en el formato que más nos interese. Veamos un ejemplo clásico con recetas de cocina.

Plantilla:

Usa la siguiente plantilla para describir recetas de cocina:
"""
### Receta [Nombre de la receta]

## Ingredientes
1 - [ingrediente]
2 - [ingrediente]

## Pasos
1 - [Paso 1]
2 - [Paso 2]

## Calorías aproximadas
[calorías totales de la receta]
"""

Le hemos dicho que lo que sigue es una plantilla, y seguidamente la describimos. Es interesante observar el uso de ciertos marcadores o placeholders que hemos usado para organizar la plantilla. Las triples comillas dobles (""") delimitan la plantilla. Es una forma estructurar prompts complejos, proporcionar contexto adicional y definir claramente las partes de la solicitud. Por ejemplo: Haz un resumen de """El fin de la Eternidad""". Si no ponemos las triples comillas, el LLM podría no entender que El fin de Eternidad es un libro.

Entre corchetes hemos puesto una descripción con la información que debe rellenar, y con el carácter # formateamos el texto. Es decir, que podemos usar el formato markdown.

Ejemplo de prompt usando plantilla:

Prompt: "Se hierven los macarrones durante 12 minutos hasta que estén al dente. Se escurre el agua. Seguidamente ponemos una lata de atún y tomate frito. Una vez emplatado, ponemos queso rallado."

Respuesta de la IA usando plantilla:

### Receta Macarrones con Atún y Tomate

## Ingredientes
1 - 200 gramos de macarrones
2 - 1 lata de atún
3 - 200 ml de tomate frito
4 - Queso rallado al gusto
5 - Agua y sal para hervir los macarrones

## Pasos
1 - Hervir los macarrones en agua con sal durante 12 minutos hasta que estén al dente.
2 - Escurrir el agua de los macarrones.
3 - Mezclar los macarrones con una lata de atún y tomate frito.
4 - Emplatar y añadir queso rallado al gusto por encima.

## Calorías aproximadas
450 calorías (por ración)

El LLM no sólo ha descrito los ingredientes y los pasos. También ha añadido un ingrediente que faltaba (la sal), le ha puesto un nombre a la receta, ha añadido las cantidades de cada ingrediente y ha calculado las calorías a partir de los ingredientes.

Cadena de pensamiento (Chain of Thought)

Esta técnica implica guiar al modelo a través de una serie de pasos lógicos para resolver un problema o generar una respuesta compleja. Al descomponer la tarea en pasos más pequeños, se mejora la coherencia y precisión de la respuesta.

Ejemplo:

Prompt: "Explica cómo se produce la fotosíntesis."

Cadena de pensamiento:

"Describe la función de la clorofila en las plantas."
"Explica cómo la luz solar es utilizada en el proceso."
"Detalla el papel del dióxido de carbono y el agua."
"Concluye describiendo el resultado final del proceso de fotosíntesis."

Uso de etiquetas XML

El uso de etiquetas XML puede ser útil para estructurar prompts complejos, proporcionar contexto adicional y definir claramente las partes de la solicitud.

Etiquetas XML

Las etiquetas XML se utilizan para etiquetar y estructurar diferentes partes del prompt, facilitando la interpretación por parte de la IA.

Ejemplo:

Prompt:

<prompt>
    <contexto>
        Explica cómo se produce la fotosíntesis en las plantas.
    </contexto>
    <paso>
        1. Describe la función de la clorofila.
    </paso>
    <paso>
        2. Explica el uso de la luz solar en el proceso.
    </paso>
    <paso>
        3. Detalla el papel del dióxido de carbono y el agua.
    </paso>
    <paso>
        4. Concluye con el resultado final de la fotosíntesis.
    </paso>
</prompt>

Modificación y expansión de prompts

A veces, un prompt necesita ser modificado o expandido para proporcionar respuestas más completas. Añadir contexto adicional, ejemplos y desglosar instrucciones en pasos más pequeños puede mejorar significativamente la calidad de las respuestas.

Ejemplo:

Prompt inicial: "Explica la teoría de la relatividad."

Prompt modificado: "Explica la teoría de la relatividad de Einstein de manera simple, como si se la estuvieras contando a un estudiante de secundaria. Asegúrate de incluir ejemplos concretos y analogías que faciliten la comprensión."

Ejemplos prácticos de Prompt Engineering

A continuación, hay varios ejemplos prácticos para diferentes escenarios, demostrando cómo el prompt engineering puede mejorar la interacción con modelos de IA.

Ejemplo 1: Resumen de texto

Prompt: "Resume el siguiente artículo en 100 palabras: [inserta texto del artículo]."

Respuesta mejorada: Para obtener un resumen preciso, proporciona un poco más de contexto.

Prompt mejorado: "Resume el siguiente artículo de noticias en 100 palabras, destacando los puntos clave sobre los impactos económicos y sociales: [inserta texto del artículo]."

Ejemplo 2: Generación de contenido creativo

Prompt: "Escribe una historia corta sobre una aventura en el espacio."

Respuesta mejorada: Añadir detalles específicos sobre los personajes y el escenario.

Prompt mejorado: "Escribe una historia corta de 500 palabras sobre una aventura en el espacio, protagonizada por una astronauta llamada Laura y su robot asistente, quienes descubren un planeta desconocido habitado por seres misteriosos."

Ejemplo 3: Resolución de problemas matemáticos

Prompt: "Resuelve este problema matemático: 3x + 5 = 20."

Respuesta mejorada: Guiar a través de los pasos para resolver el problema.

Prompt mejorado: "Resuelve el siguiente problema matemático mostrando todos los pasos: 3x + 5 = 20. Primero, despeja el término constante y luego divide por el coeficiente de x."

Ejemplo 4: Consultoría técnica

Prompt: "Dame consejos para mejorar la seguridad informática."

Respuesta mejorada: Especificar el tipo de consejos y el contexto.

Prompt mejorado: "Proporciona cinco consejos avanzados para mejorar la seguridad informática en una empresa de tecnología, enfocándote en la protección de datos sensibles y la prevención de ataques de phishing."

Ejemplo 5: Análisis de datos

Prompt: "Analiza los datos de ventas del último trimestre."

Respuesta mejorada: Incluir detalles sobre qué aspectos analizar y el formato esperado de la respuesta.

Prompt mejorado: "Analiza los datos de ventas del último trimestre, enfocándote en identificar las tendencias de ventas mensuales, los productos más vendidos y las regiones con mejor desempeño. Presenta los resultados en un formato de informe ejecutivo."

Conclusión

El prompt engineering es una habilidad importante para interactuar de manera efectiva con modelos de lenguaje avanzados. Al seguir principios básicos como la claridad, especificidad, y el uso de ejemplos, junto con técnicas avanzadas como las que hemos visto más arriba, se pueden obtener respuestas más precisas y relevantes.

Scripts en NMAP

noreply@blogger.com (Alberto García Serrano) — Mon, 29 May 2023 17:11:00 +0000

Cuando pensamos en NMAP, pensamos en el escaneo de puertos de un host objetivo al que estamos relizando una prueba de intrusión, pero gracias a las posibilidades que nos ofrecen su Scripting Engine, NMAP es mucho más que eso.

Antes de continuar, un aviso: algunas de posibilidades que nos ofrecen los scripts de NMAP son bastante intrusivas, por lo que recomiendo hacerlas contra hosts propios, máquinas virtuales como las de Metasploitable, o contrato de pentesting mediante. Para este artículo voy a usar las máquinas de Metasploitable3.

No voy a entrar en los detalles sobre el uso básico de NMAP, ya que hay miles de tutoriales en Internet que hablan sobre ello. Lo cierto es que NMAP tiene algunas opciones que permiten obtener información extra, además de qué puertos están abiertos y cuales no. Por ejemplo, la opción -sV trata de obtener el servicio concreto, e incluso la versión del servicio que está corriendo en cada puerto. Otro ejemplo es la opción -O, que intenta averiguar el sistema operativo del host.

sudo nmap -sV -p80 -O 192.168.56.102

En este caso, le hemos pedido que nos dé información sobre el servicio que está corriendo en el puerto 80, y el sistema operativo que tiene el host. Vemos que ha hecho un buen trabajo averiguando la versión del servidor web (Microsoft IIS httpd 7.5), pero con el sistema operativo no llega a ser muy concreto en cuanto la versión exacta que hay instalada.

Si miramos el directorio /usr/share/nmap/scripts, observaremos que hay una gran cantidad de archivos con extensión .nse. Son scripts con diversas funcionalidades. Desde la recogida de información hasta la explotación de vulneravilidades. Pero vamos por pasos. Para poder usar uno de estos scripts, utilizamos la opción --script de NMAP. En el siguiente ejemplo vamos a usar un script llamado smb-os-discovery.nse, cuya misión es intentar obtener información del sistema operativo del host objetivo a través del protocolo SMB. Como dicho protocolo usa los puertos 139 y 445, vamos a usar el script usando la siguiente instrucción.

sudo nmap -p139,445 --script=smb-os-discovery 192.168.56.102

Gracias al script, hemos obtenido el sistema operativo exacto que usa la máquina, además de alguna información adicional, como el nombre del equipo o el grupo de trabajo de NetBIOS. No está mal ¿no?

Categorías

NMAP viene cargado con más de 600 scripts, y para organizarlos, están divididos en categorías, es decir, que cada script pretenece a una o más categorías. Vamos a hacer un rápido repaso.

auth. Scripts relacionados con la autenticación de usuarios.

broadcast. Descubrimiento de hosts mediante mensajes de broadcast en la red local.

brute. Uso de técnicas de fuerza bruta para descubrir crendenciales sobre diferentes protocolos y aplicaciones.

default. Scripts por defecto, que pueden ser lanzados mediante la opción -sC de NMAP. No tengo claro el criterio, pero supongo que son los scripts más usuales.

discovery. Trata de hacer descubrimiento de hosts e información mediante consultas a servicios SMB, SNMP, directorios LDAP, etc.

dos. Relacionados con ataques tipo DOS (denegación de servicio).

exploit. Script pensados para explotar alguna vulnerabilidad.

external. Son scripts que realizan conexiones con servicios externos, como por ejemplo whois.

fuzzer. Los scripts de esta categoría utilizan técnicas de fuzzing.

intrusive. Son scripts extremadamente intrusivos, que podrían incluso bloquear el host objetivo.

malware. Tratan de identificar si el host ha sido comprometido por algún tipo de malware o tiene alguna puerta trasera.

safe. Son script poco intrusivos que, en principio, no deben dañar el host objetivo.

version. Estos scripts son invocados cuando se usa la opción -sV. No se pueden invocar explícitamente el uso de esta categoría.

vuln. Hacen una búsqueda de vulnerabilidades en el host.

Podemos lanzar los scripts individualmente, como ya hemos visto, o por categorías completas. Los scripts en la categoría default los podemos lanzar con la instrucción nmap -sC. Si queremos lanzar los scripts de otra categoría usamos el modificador --script, pero indicando la categoría en lugar de un script concreto.

nmap --script=vuln

Los nombres de los scripts son bastante representativos, pero si queremos obtener detalles de alguno de ellos podemos usar la opción --script-help.

nmap --script-help smb-enum-users

Uso de los scripts de NMAP para pentesting

Ahora que tenemos una visión general, vamos a ponerlos en práctica contra un host Windows del Metaesploitable3 que os indicaba al principio.

Uno de los servicios que muchas veces no está bien configurado en los servidores Windows es el SNMP, así que suele ser un buen lugar donde lanzar la caña. Veamos si podemos recoger algo usando el script snmp-win32-users. En este caso, como SNMP va por UDP, usaremos el modificador -sU para usar este protocolo, y ya que el puerto de SNMP es el 161, usaremos también el modificador -p161.

$sudo nmap -sU -p161 --script=snmp-win32-users 192.168.56.102
Starting Nmap 7.93 ( https://nmap.org ) at 2023-05-29 17:28 CEST
Nmap scan report for 192.168.56.102
Host is up (0.00069s latency).

PORT    STATE SERVICE
161/udp open  snmp
| snmp-win32-users: 
|   Administrator
|   Guest
|   anakin_skywalker
|   artoo_detoo
|   ben_kenobi
|   boba_fett
|   c_three_pio
|   chewbacca
|   darth_vader
|   greedo
|   han_solo
|   jabba_hutt
|   jarjar_binks
|   kylo_ren
|   lando_calrissian
|   leia_organa
|   luke_skywalker
|   sshd
|   sshd_server
|_  vagrant
MAC Address: 08:00:27:DC:53:82 (Oracle VirtualBox virtual NIC)

Nmap done: 1 IP address (1 host up) scanned in 13.50 seconds

Estupendo, ahí tenemos a una serie de usuarios del sistema. Un autentico fan de Starwars. Veamos si podemos averiguar qué software está instalado en la máquina.

$sudo nmap -sU -p161 --script=snmp-win32-software 192.168.56.102
Starting Nmap 7.93 ( https://nmap.org ) at 2023-05-29 17:32 CEST
Nmap scan report for 192.168.56.102
Host is up (0.00067s latency).

PORT    STATE SERVICE
161/udp open  snmp
| snmp-win32-software: 
|   7-Zip 22.01 (x64); 2023-03-19T02:22:16
|   Java 8 Update 251 (64-bit); 2023-03-19T02:23:56
|   Java SE Development Kit 8 Update 211 (64-bit); 2023-03-19T02:25:46
|   Microsoft .NET Framework 4.5.2; 2023-03-19T02:12:48
|   Microsoft .NET Framework 4.5.2; 2023-03-19T02:12:44
|   Microsoft Visual C++ 2008 Redistributable - x64 9.0.30729.6161; 2023-03-19T02:27:48
|   OpenSSH for Windows 7.1p1-1 (remove only); 2023-03-19T02:06:06
|_  Oracle VM VirtualBox Guest Additions 6.0.8; 2023-03-19T02:45:06
MAC Address: 08:00:27:DC:53:82 (Oracle VirtualBox virtual NIC)

Nmap done: 1 IP address (1 host up) scanned in 13.45 seconds

Y de camino, los servicios que están corriendo en la máquina.

$sudo nmap -sU -p161 --script=snmp-win32-services 192.168.56.102
Starting Nmap 7.93 ( https://nmap.org ) at 2023-05-29 17:33 CEST
Nmap scan report for 192.168.56.102
Host is up (0.00084s latency).

PORT    STATE SERVICE
161/udp open  snmp
| snmp-win32-services: 
|   Apache Tomcat 8.0 Tomcat8
|   Application Host Helper Service
|   Background Intelligent Transfer Service
|   Base Filtering Engine
|   COM+ Event System
|   Certificate Propagation
|   Cryptographic Services
|   DCOM Server Process Launcher
|   DHCP Client
|   DNS Client
|   Desktop Window Manager Session Manager
|   Diagnostic Policy Service
|   Diagnostic System Host
|   Distributed Link Tracking Client
|   Distributed Transaction Coordinator
|   Elasticsearch 1.1.1 (elasticsearch-service-x64)
|   Group Policy Client
|   IKE and AuthIP IPsec Keying Modules
|   IP Helper
|   IPsec Policy Agent
|   MEDC Server Component - Apache
|   MEDC Server Component - Notification Server
|   ManageEngine Desktop Central Server
|   Microsoft FTP Service
|   Network Connections
|   Network List Service
|   Network Location Awareness
|   Network Store Interface Service
|   OpenSSH Server
|   Plug and Play
|   Power
|   Print Spooler
|   RPC Endpoint Mapper
|   Remote Desktop Configuration
|   Remote Desktop Services
|   Remote Desktop Services UserMode Port Redirector
|   Remote Procedure Call (RPC)
|   Remote Registry
|   SNMP Service
|   Security Accounts Manager
|   Server
|   Shell Hardware Detection
|   Software Protection
|   System Event Notification Service
|   TCP/IP NetBIOS Helper
|   Task Scheduler
|   User Profile Service
|   VirtualBox Guest Additions Service
|   Windows Event Log
|   Windows Firewall
|   Windows Font Cache Service
|   Windows Licensing Monitoring Service
|   Windows Management Instrumentation
|   Windows Process Activation Service
|   Windows Remote Management (WS-Management)
|   Windows Update
|   Workstation
|   World Wide Web Publishing Service
|   domain1 GlassFish Server
|   jenkins
|   jmx
|   wampapache
|_  wampmysqld
MAC Address: 08:00:27:DC:53:82 (Oracle VirtualBox virtual NIC)

Nmap done: 1 IP address (1 host up) scanned in 13.55 seconds

Con este listado, ya podríamos consultar la base de datos CVE para ver si hay alguna vulnerabilidad en alguno de los programas o servicios, pero ya que estamos, vamos a pedirle a NMAP que nos adelante algo de trabajo pidiendole que lance los scripts de la categoría vuln contra el host.

sudo nmap --script=vuln 192.168.56.102
Starting Nmap 7.93 ( https://nmap.org ) at 2023-05-29 17:35 CEST
Nmap scan report for 192.168.56.102
Host is up (0.00055s latency).
Not shown: 981 closed tcp ports (reset)
PORT      STATE SERVICE
21/tcp    open  ftp
22/tcp    open  ssh
80/tcp    open  http
|_http-csrf: Couldn't find any CSRF vulnerabilities.
|_http-dombased-xss: Couldn't find any DOM based XSS.
|_http-stored-xss: Couldn't find any stored XSS vulnerabilities.
135/tcp   open  msrpc
139/tcp   open  netbios-ssn
445/tcp   open  microsoft-ds
3306/tcp  open  mysql
3389/tcp  open  ms-wbt-server
4848/tcp  open  appserv-http
| ssl-dh-params: 
|   VULNERABLE:
|   Diffie-Hellman Key Exchange Insufficient Group Strength
|     State: VULNERABLE
|       Transport Layer Security (TLS) services that use Diffie-Hellman groups
|       of insufficient strength, especially those using one of a few commonly
|       shared groups, may be susceptible to passive eavesdropping attacks.
|     Check results:
|       WEAK DH GROUP 1
|             Cipher Suite: TLS_DHE_RSA_WITH_AES_128_CBC_SHA
|             Modulus Type: Safe prime
|             Modulus Source: RFC2409/Oakley Group 2
|             Modulus Length: 1024
|             Generator Length: 8
|             Public Key Length: 1024
|     References:
|_      https://weakdh.org
7676/tcp  open  imqbrokerd
8009/tcp  open  ajp13
8080/tcp  open  http-proxy
| http-enum: 
|_  ..%2f..%2f..%2f..%2f..%2f..%2f..%2f..%2f/var/mobile/Library/AddressBook/AddressBook.sqlitedb: Possible iPhone/iPod/iPad generic file sharing app Directory Traversal (iOS)
8181/tcp  open  intermapper
| ssl-dh-params: 
|   VULNERABLE:
|   Diffie-Hellman Key Exchange Insufficient Group Strength
|     State: VULNERABLE
|       Transport Layer Security (TLS) services that use Diffie-Hellman groups
|       of insufficient strength, especially those using one of a few commonly
|       shared groups, may be susceptible to passive eavesdropping attacks.
|     Check results:
|       WEAK DH GROUP 1
|             Cipher Suite: TLS_DHE_RSA_WITH_AES_128_CBC_SHA
|             Modulus Type: Safe prime
|             Modulus Source: RFC2409/Oakley Group 2
|             Modulus Length: 1024
|             Generator Length: 8
|             Public Key Length: 1024
|     References:
|_      https://weakdh.org
8383/tcp  open  m2mservices
| ssl-dh-params: 
|   VULNERABLE:
|   Diffie-Hellman Key Exchange Insufficient Group Strength
|     State: VULNERABLE
|       Transport Layer Security (TLS) services that use Diffie-Hellman groups
|       of insufficient strength, especially those using one of a few commonly
|       shared groups, may be susceptible to passive eavesdropping attacks.
|     Check results:
|       WEAK DH GROUP 1
|             Cipher Suite: TLS_DHE_RSA_WITH_AES_128_CBC_SHA
|             Modulus Type: Safe prime
|             Modulus Source: RFC2409/Oakley Group 2
|             Modulus Length: 1024
|             Generator Length: 8
|             Public Key Length: 1024
|     References:
|_      https://weakdh.org
9200/tcp  open  elasticsearch
49152/tcp open  unknown
49153/tcp open  unknown
49154/tcp open  unknown
49176/tcp open  unknown
MAC Address: 08:00:27:DC:53:82 (Oracle VirtualBox virtual NIC)

Host script results:
|_smb-vuln-ms10-054: false
|_smb-vuln-ms10-061: NT_STATUS_ACCESS_DENIED
| smb-vuln-ms17-010: 
|   VULNERABLE:
|   Remote Code Execution vulnerability in Microsoft SMBv1 servers (ms17-010)
|     State: VULNERABLE
|     IDs:  CVE:CVE-2017-0143
|     Risk factor: HIGH
|       A critical remote code execution vulnerability exists in Microsoft SMBv1
|        servers (ms17-010).
|           
|     Disclosure date: 2017-03-14
|     References:
|       https://blogs.technet.microsoft.com/msrc/2017/05/12/customer-guidance-for-wannacrypt-attacks/
|       https://technet.microsoft.com/en-us/library/security/ms17-010.aspx
|_      https://cve.mitre.org/cgi-bin/cvename.cgi?name=CVE-2017-0143
|_samba-vuln-cve-2012-1182: NT_STATUS_ACCESS_DENIED

Nmap done: 1 IP address (1 host up) scanned in 152.97 seconds

NMAP ha encontrado por nosotros algunos problemas, como posibles debilidades criptográficas, y lo mejor de todo, vulnerabilidades como la CVE-2017-0143 que tiene SMBv1, que permite la ejecución remota de código. Sólo con esta información, ya tendríamos material interesante para pasar a la fase de explotación de nuestro test de intrusión. Aunque de eso ya hablaremos otro día. Por hoy podemos descansar, que nos lo hemos ganado.

Detección de caras con OpenCV

noreply@blogger.com (Alberto García Serrano) — Sat, 08 Jan 2022 00:59:00 +0000

En visión artificial, una de las tareas básicas es el reconocimiento de objetos, ya que es muy interesante disponer de algún modo de discriminar entre clases de objetos. Por ejemplo, un coche autónomo ha de ser capaz de reconocer a otros coches y distinguirlo de un peatón. Otra utilidad más mundana y que seguro conoces, es el reconocimiento de caras que hace tu móvil cuando vas a hacer una foto.

En este artículos vamos a ver cómo reconocer caras usando OpenCV. Ojo, cuando hablo de reconocer caras, me refiero a distinguir una cara de cualquier otra cosa, no de distinguir entre caras para, por ejemplo, reconocer a quién pertenece esa cara. Ese es otro problema diferente del que os hablaré en otro momento. Como quiero abordar el problema desde un punto de vista práctico, vamos a usar la librería OpenCV con Python. Como imagen de ejemplo vamos a tomar la siguiente.

Nuestra misión es detectar las caras que hay en esta imagen, así que fije el rumbo señor Sulu.

import cv2
from matplotlib import pyplot as plt
img = cv2.imread('TOS-crew.jpg')

Tras importar las librerías OpenCV y Matplotlib para visualizar y manipular la imagen, comprobamos que se ha cargado correctamente y verificamos su tamaño.

img.shape

(1115, 2329, 3)

Si queremos visualizar la imagen con Matplotlib primero debemos convertirla a la codificación RGB, ya que OpenCV por defecto trabaja con imagenes BGR.

img_RGB = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
plt.imshow(img_RGB)

Si vamos a trabajar con imágenes grandes, a veces es mejor redimensionarlas para hacer las más pequeñas e incluso convertirlas a escala de grises para aliviar su tamaño. Como ejemplo, las líneas siguientes realizan estas dos operaciones.

img_gray = cv2.cvtColor(img_RGB, cv2.COLOR_RGB2GRAY)
img_small = cv2.resize(img_gray, (800, 400))
plt.imshow(img_small)

Ahora vamos a realizar las busqueda de las caras en la imagen. Con OpenCV es de lo más sencillo. Sólo necesitamos dos líneas de código Python. Sin entrar en detalles, os cuento que el clasificador que vamos a usar se conoce como Haar Cascade o también como algoritmo de Viola-Jones. Vamos a usarlo porque es de los clásicos (que no el mejor). Si queréis conocer más a fondo como funciona el algoritmo podéis consultar el artículo original de 2004 donde Paul Viola y Michael J. Jones presentaron este método en sociedad. Para que este algoritmo funcione hay que entrenarlo mostrandole caras, y otras cosas que no sean caras, para que aprenda a distinguirlos. El proceso puede ser tedioso, ya que son necesaria una relativamente grande cantidad de imágenes. Por suerte ya hay clasificadores preentrenados para reconocer caras (y otras muchas cosas), que están almacenas en formato .xml y que podemos cargar desde OpenCV. El que vamos a usar se llama haarcascade_frontalface_default.xml y se puede descargar de aquí. Está entrenado con caras mirando al frente, cómo las que queremos detectar. En esta página hay más clasificadores preentrenados.

# cargamos el clasificador preentrenado
detector_caras = cv2.CascadeClassifier('haarcascade_frontalface_default.xml')
detecciones = detector_caras.detectMultiScale(img_RGB)
# array de detecciones. una fila por cada cara. Las columnas con [pos_x, pos_y, size_x, size_y]
detecciones

array([[ 936,  211,   58,   58],
       [ 495,   76,  191,  191],
       [ 708,  109,  235,  235],
       [1199,  120,  195,  195],
       [ 221,  124,  212,  212],
       [1810,  225,  216,  216],
       [1464,  325,  215,  215],
       [1540,  326,  215,  215],
       [1435,  564,  352,  352],
       [ 622,  572,  348,  348],
       [1296,  992,   62,   62]], dtype=int32)

El array detecciones contiene una fila por cada cara detectada, con cuatro columnas cada una indicando la posición X e Y, así como la anchura y altura del recuadro imaginario que contiene la cara que ha detectado. En este caso vemos que ha detectado 11 caras, cuando en nuestra imagen sólo hay seis personajes. Algo no parece ir bien, así que vamos a mostrar la imagen sobre la que vamos a dibujar estos cuadrados que nos ha devuelto el clasificador.

for (x, y, w, h) in detecciones:
    cv2.rectangle(img_RGB, (x, y), (x+w, y+h), (0,255,0), 10)

plt.imshow(img_RGB)

La buena noticia es que ha detectado todas las caras (no ha habido falsos negativos), pero además de las autenticas caras, ha detectado como caras otras partes de la imagen (falsos positivos). Vamos a tratar de realizar algunos ajustes a ver si conseguimos mejorar el resultado. La función detectMultiScale busca las caras realizando varios barridos de la imagen con ventanas de diferentes tamaños (ya que puede haber caras grandes cercanas o pequeñas que estén más lejanas). El tamaño en que se incrementa la ventana en cada barrido se puede cambiar con el parámetro scaleFactor. Por ejemplo un valor de 1.1 (valor por defecto) quiere decir que la ventana seingrementa en un 10% en cada barrido, y un valor de 1.05, aumentaría la ventana en un 5%. Aquí hay que jugar con el valor según sea el tamaño de las caras que aparecen en la imagen. Para nuestra imagen, un valor de 1.7 parece funcionar bien

detecciones = detector_caras.detectMultiScale(img_RGB, scaleFactor = 1.7)

Otro parámetro con el que podemos jugar es minNeighbors. Cuando el algoritmo va recorriendo la imagen, acaba reconociendo la misma cara en múltiples ocasiones (aunque sólo nos devuelva una ventana). Por ejemplo, puede reconocer la misma cara con la misma ventana al desplazarla sólo unos pocos pixeles, o en sus múltiples pasadas con diferentes tamaños de ventana. Con minNeighbors podemos indicar el número mínimo de detecciones de una cara que han de darse para que se considere una detección válida. Por ejemplo, la siguiente línea es capaz de detectar correctamente las caras a pesar de no utilizar el parámetro scaleFactor.

detecciones = detector_caras.detectMultiScale(img_RGB, minNeighbors = 20)

Por último tenemos los parámetros minSize y maxSize, que permiten indicar el tamaño mínimo y máximo de la cara a detectar. Podemos ver la utilidad práctica con un ejemplo. Veamos qué detecta el clasificador con los siguientes parámetros.

detecciones = detector_caras.detectMultiScale(img_RGB, minNeighbors = 10)

Vemos que hay un falso positivo en la oreja del capitán Kirk. Una forma de manejar este tipo de problemas, si sabemos que las caras van a ser mayores (o menores) de cierto tamaño, es indicar un tamaño mínimo (o máximo) para la detección.

detecciones = detector_caras.detectMultiScale(img_RGB, minNeighbors = 10, minSize = (100,100))

Como dije al principio, hay otros tipos de clasificadores que tienen mayor rendimiento, pero este es bastante sencillo y eficiente incluso con baja capacidad de proceso, por lo que es muy interesante. En otra ocasión veremos si es posible, no sólo reconocer que hay una cara, sino, además, intentar saber de quién es la cara. Es lo que se denomina reconocimiento facial, o en general reconocimiento de objetos (o de caracteres, o de tipos de naves de la federación, etc.). Hasta entonces... larga vida y prosperidad.

Cómo leer el sensor de temperatura de la Raspberry Pi Pico

noreply@blogger.com (Alberto García Serrano) — Wed, 10 Feb 2021 17:20:50 +0000

Leyendo las especificaciones de la Raspberry Pi Pico me llamó la atención que venía equipada con un sensor de temperatura en el propio chip, así que me puse a buscar información. En este artículo os quiero mostrar cómo acceder a él y de camino cómo utilizar el conversor analógico/digital (ADC).

El siguiente esquema muestra cómo está conectado el ADC a los GPIO y al sensor de temperatura interno.

Es decir, tenemos un multiplexor de cuatro entradas que permite seleccionar la señal de entrada para el ADC. En nuestro caso, la que nos interesa es la 4, que es la fuente de tensión proveniente del sensor. Para seleccionar la entrada en el multiplexor usamos la siguiente función.

adc_select_input(4);

Para hacernos una idea, en la página 574 del datasheet del RP4020 nos dice que esta tensión es de 0.706V a 27 grados centígrados. La fórmula para aproximar la temperatura con respecto a esta tensión es

T = 27 - (ADC_voltage - 0.706)/0.001721

También hay que poner a 1 el registro TS_EN, para conectar el sensor. El SDK nos permite hacer esto fácilmente con la función

adc_set_temp_sensor_enabled(true);

Veamos pues cómo sería el proceso completo con un ejemplo en el que vamos a usar un display LCD 1602 para mostrar la temperatura. El manejo del LCD se hace por I2C y no voy a entrar en mucho detalle aquí, ya que hay bastante documentación al respecto. Además tenéis un ejemplo en la carpeta pico-examples/i2c/lcd_1602_i2c/ que acompaña al SDK.

Ya os conté en otro artículo los pasos para crear un proyecto para la Raspberry Pi Pico, así que ya sabemos que necesitamos como mínimo dos archivos: CMakeLists.txt y el código fuente, que en este caso lo he llamado termometro.c.

cmake_minimum_required(VERSION 3.12)
#enable_language( C ASM )
include(~/tools/pico/pico-sdk/external/pico_sdk_import.cmake)
#include(pico_sdk_import.cmake)
project(termometro C CXX ASM)
set(CMAKE_C_STANDARD 11)
set(CMAKE_CXX_STANDARD 17)
pico_sdk_init()
add_executable(termometro termometro.c)
target_link_libraries(termometro pico_stdlib hardware_i2c hardware_adc)
pico_add_extra_outputs(termometro)

Este es el contenido del archivo CMakeLists.txt, que es muy similar al que ya vimos en el anterior artículo. Se diferencia en la siguiente línea, que indica al compilador que debe añadir las librería para la gestión del bus I2C (para el LCD) y del ADC.

target_link_libraries(termometro pico_stdlib hardware_i2c hardware_adc)

Antes de poner el código fuente completo, voy a centrarme en los pasos para leer el sensor de temperatura. El código mínimo para leer el sensor tendría una forma como ésta.

#include "hardware/adc.h"
int main() {
    // 12 bits. Vref=3.3V
    const float conv_factor = 3.3f / (1<<12); 
	
    uint16_t lectura;
    float lect_conv, temperatura;
    adc_init();
    adc_set_temp_sensor_enabled(true);
    adc_select_input(4); // imput 4 - sensor temperatura
	    
    while(1) {
        // Consultar ADC
        lectura = adc_read();
        lect_conv = lectura * conv_factor;
        temperatura = 27-(lect_conv-0.706)/0.001721;
        // mostrar temperatura
		printf("%2.1f", temperatura);
		    
        sleep_ms(20000);
    }
}

A la lectura del ADC hay que aplicarle un factor de conversión. En este caso, la tensión máxima es de 3.3V y nuestro ADC es de 12 bits, por lo que este factor es 3.3f/(1<<12). Al resultado le aplicamos la fórmula que habíamos visto más arriba para aproximar la temperatura.

El programa completo es algo más complicado, ya que incluye el código para el manejo del LCD, y sería el siguiente.

#include <stdio.h>
#include "pico/stdlib.h"
#include "hardware/i2c.h"
#include "hardware/adc.h"
#include "pico/binary_info.h"
/* panel LCD 16x2 con adaptador I2C PCF8574 
   GPIO 20 (pin 26)-> SDA
   GPIO 21 (pin 27)-> SCL
   VSYS (pin 39) -> VCC
   GND (pin 38)  -> GND
*/
// algunas definiciones para controlar el LCD
const int LCD_CLEARDISPLAY = 0x01;
const int LCD_RETURNHOME = 0x02;
const int LCD_ENTRYMODESET = 0x04;
const int LCD_DISPLAYCONTROL = 0x08;
const int LCD_CURSORSHIFT = 0x10;
const int LCD_FUNCTIONSET = 0x20;
const int LCD_SETCGRAMADDR = 0x40;
const int LCD_SETDDRAMADDR = 0x80;
const int LCD_ENTRYSHIFTINCREMENT = 0x01;
const int LCD_ENTRYLEFT = 0x02;
const int LCD_BLINKON = 0x01;
const int LCD_CURSORON = 0x02;
const int LCD_DISPLAYON = 0x04;
const int LCD_MOVERIGHT = 0x04;
const int LCD_DISPLAYMOVE = 0x08;
const int LCD_5x10DOTS = 0x04;
const int LCD_2LINE = 0x08;
const int LCD_8BITMODE = 0x10;
const int LCD_BACKLIGHT = 0x08;
const int LCD_ENABLE_BIT = 0x04;
// dirección del bus I2C
// podría cambiar según modelos
#define I2C_PORT_LCD i2c0
static int addr_led = 0x27;
#define LCD_CHARACTER  1
#define LCD_COMMAND    0
#define MAX_LINES      2
#define MAX_CHARS      16
// pines I2C en la placa
#define LCD_SDA 20
#define LCD_SCL 21
void i2c_write_byte(uint8_t val) {
    i2c_write_blocking(I2C_PORT_LCD, addr_led, &val, 1, false);
}
void lcd_toggle_enable(uint8_t val) {
#define DELAY_US 600
    sleep_us(DELAY_US);
    i2c_write_byte(val | LCD_ENABLE_BIT);
    sleep_us(DELAY_US);
    i2c_write_byte(val & ~LCD_ENABLE_BIT);
    sleep_us(DELAY_US);
}void lcd_send_byte(uint8_t val, int mode) {
    uint8_t high = mode | (val & 0xF0) | LCD_BACKLIGHT;
    uint8_t low = mode | ((val << 4) & 0xF0) | LCD_BACKLIGHT;
    i2c_write_byte(high);
    lcd_toggle_enable(high);
    i2c_write_byte(low);
    lcd_toggle_enable(low);
}
void lcd_clear(void) {
    lcd_send_byte(LCD_CLEARDISPLAY, LCD_COMMAND);
}
void lcd_set_cursor(int line, int position) {
    int val = (line == 0) ? 0x80 + position : 0xC0 + position;
    lcd_send_byte(val, LCD_COMMAND);
}
static void inline lcd_char(char val) {
    lcd_send_byte(val, LCD_CHARACTER);
}
void lcd_string(const char *s) {
    while (*s) {
        lcd_char(*s++);
    }
}
void lcd_init() {
    lcd_send_byte(0x03, LCD_COMMAND);
    lcd_send_byte(0x03, LCD_COMMAND);
    lcd_send_byte(0x03, LCD_COMMAND);
    lcd_send_byte(0x02, LCD_COMMAND);
    lcd_send_byte(LCD_ENTRYMODESET | LCD_ENTRYLEFT, LCD_COMMAND);
    lcd_send_byte(LCD_FUNCTIONSET | LCD_2LINE, LCD_COMMAND);
    lcd_send_byte(LCD_DISPLAYCONTROL | LCD_DISPLAYON, LCD_COMMAND);
    lcd_clear();
}
int main() {
	uint16_t lectura;
	float lect_conv, temperatura;
	char s[10];
	const float conv_factor = 3.3f / (1<<12); 
    i2c_init(I2C_PORT_LCD, 100 * 1000);
    gpio_set_function(LCD_SDA, GPIO_FUNC_I2C);
    gpio_set_function(LCD_SCL, GPIO_FUNC_I2C);
    gpio_pull_up(LCD_SDA);
    gpio_pull_up(LCD_SCL);
    bi_decl( bi_2pins_with_func(LCD_SDA, LCD_SCL, GPIO_FUNC_I2C));
	
    lcd_init();
    adc_init();
	adc_set_temp_sensor_enabled(true);
	adc_select_input(4); // imput 4 - sensor temperatura
	
	while(1) {
		// Consultar ADC
		lectura = adc_read();
		lect_conv = lectura * conv_factor;
		temperatura = 27-(lect_conv-0.706)/0.001721;
        // mostrar temperatura en el LCD
		sprintf(s, "%2.1f", temperatura);
		lcd_set_cursor(0,0);
		lcd_string("Temperatura:");
		lcd_set_cursor(1,0);
		lcd_string(s);
		sleep_ms(20000);
		lcd_clear();
	}
    return 0;
}

Las conexiones entre el LCD y la Raspberry Pi Pico son las siguientes.

   GPIO 20 (pin 26)-> SDA en LCD
   GPIO 21 (pin 27)-> SCL en LCD
   VSYS (pin 39) -> VCC en LCD
   GND (pin 38)  -> GND en LCD

Ya vimos que para poder compilar el programa tendríamos que invocar a cmake para que se genere el makefile correspondiente.

cmake .

Y finalmente compilamos con make. El resultado en el siguiente video, donde se observa cómo sube la temperatura al tapar el microcontrolador con el dedo.

Primeros pasos con Raspberry Pi Pico

noreply@blogger.com (Alberto García Serrano) — Sat, 30 Jan 2021 11:50:59 +0000

Hace pocos días la fundación Raspberry Pi sacaba al mercado una curiosa plaquita llamada Raspberry Pi Pico. Curiosa porque se separa de lo que venían haciendo.

Es decir, nos encontramos ante una placa al estilo del Arduino nano, pero comandada por el microcontrolador RP4020, basado en el ARM Cortex M0+ de doble núcleo, que ha sido también diseñado por ellos. En la placa de desarrollo nos ofrecen 26 GPIOs con soporte para SPI, I2C, ADC, UART y PWM. Acompañado de una completísima documentación y un entorno de desarrollo basado en Visual Studio Code. Y lo que es mejor, sólo por 4 dólares.

La instalación del entorno de desarrollo oficial es muy sencilla y viene bien explicada en la documentación. Aunque si, como yo, prefieres la línea de comando, en este artículo os explicaré como poner en marcha un proyecto para este microcontrolador y programarlo con vuestro editor favorito. En este caso, para equipos con Debian y derivados. Antes de nada necesitaremos instalar el toolchain para ARM, git para obtener el SDK y los ejemplos del repositorio de Raspberry, y cmake, que es la herramienta en la que está basada el proceso de compilación en el entorno de desarrollo.

sudo apt install cmake gcc-arm-none-eabi build-essential git

Creamos un directorio que contendrá el SDK y los ejemplos, siempre útiles a la hora de afrontar un proyecto propio.

mkdir pico
cd pico
git clone -b master https://github.com/raspberrypi/pico-sdk.git
cd pico-sdk/
git submodule update --init
cd ..
git clone -b master https://github.com/raspberrypi/pico-examples.git

El directorio pico-sdk contiene las librerías necesarias para la programación del microcontrolador, y la carpeta pico-examples, unos cuantos ejemplos. En la documentación oficial y en los miles de tutoriales que aparecerán estos días puedes ver cómo compilar y ejecutar los ejemplos, pero yo voy a ir directo a explicar como crear un proyecto desde cero. He creado una carpeta llamada proyectos dentro del directorio pico, pero esto ya es cuestión de gustos y de la forma que tengas de organizarte.

mkdir proyectos
cd proyectos
mkdir blinky
cd blinky

En este caso, el proyecto que vamos a crear se llama blinky, y como adivinarás por el nombre, lo que vamos a hacer es hacer parpadear un LED, el hola mundo de los microcontroladores. De hecho el código C es el mismo que viene en los ejemplos, así que simplemente crea un archivo llamado blink.c e introduce el siguiente código, que lo que hace parpadear el LED que hay conectado al GPIO25 de la placa.

/**
 * Copyright (c) 2020 Raspberry Pi (Trading) Ltd.
 *
 * SPDX-License-Identifier: BSD-3-Clause
 *
 #include "pico/stdlib.h"
 int main() {
    const uint LED_PIN = 25;
    gpio_init(LED_PIN);
    gpio_set_dir(LED_PIN, GPIO_OUT);
    while (true) {
        gpio_put(LED_PIN, 1);
        sleep_ms(250);
        gpio_put(LED_PIN, 0);
        sleep_ms(250);
    }
}

Como ya se ha dicho, el entorno de desarrollo está basado en cmake, así que necesitamos crear un archivo llamado CMakeLists.txt donde se configura el proyecto. En este caso es el siguiente.

cmake_minimum_required(VERSION 3.12)
include(~/pico/pico-sdk/external/pico_sdk_import.cmake)
project(blink C CXX ASM)
set(CMAKE_C_STANDARD 11)
set(CMAKE_CXX_STANDARD 17)
pico_sdk_init()
add_executable(blink blink.c)
target_link_libraries(blink pico_stdlib)
pico_add_extra_outputs(blink)

En la línea del include debéis poner la ruta donde instalasteis el SDK, que puede no coincidir con la mía. Sin entrar en mucho detalle, porque cmake es un universo en sí mismo, en la línea project, pondremos el nombre de nuestro proyecto.

En add_executable indicamos el nombre del ejecutable que se va a generar y el archivo con el código fuente. También hay que añadir este nombre en larget_link_libraries, para que se enlace con las librerías que vayamos a usar del SDK, en este caso sólo pico_stdlib, pero su fuéramos a usar otras, como, por ejemplo, las de i2c, añadiríamos también hardware_i2c detrás de pico_stdlib.

Con pico_add_extra_outputs, indicamos que queremos obtener el ejecutable, además de en formato elf, en formato .uf2, que es el que usa el microcontrolador.

Ahora le decimos a cmake que configure el proyecto y genere los archivos necesarios para poder compilar con make. Esto lo hacemos con cmake . (importante poner el . después de cmake). Pero antes tendremos que configurar la variable de entorno PICO_SDK_PATH con la ruta del SDK.

export PICO_SDK_PATH=~/pico/pico-sdk
cmake .

Deberíamos ver algo parecido a esto.

Using PICO_SDK_PATH from environment ('/home/alberto/tools/pico/pico-sdk')
Pico SDK is located at /home/alberto/tools/pico/pico-sdk
Defaulting PICO_PLATFORM to rp2040 since not specified.
Defaulting PICO platform compiler to pico_arm_gcc since not specified.
-- Defaulting build type to 'Release' since not specified.
PICO compiler is pico_arm_gcc
PICO_GCC_TRIPLE defaulted to arm-none-eabi
-- The C compiler identification is GNU 9.2.1
-- The CXX compiler identification is GNU 9.2.1
-- The ASM compiler identification is GNU
-- Found assembler: /usr/bin/arm-none-eabi-gcc
Defaulting PICO target board to pico since not specified.
Using board configuration from /home/alberto/tools/pico/pico-sdk/src/boards/include/boards/pico.h
-- Found Python3: /usr/bin/python3.8 (found version "3.8.5") found components: Interpreter 
TinyUSB available at /home/alberto/tools/pico/pico-sdk/lib/tinyusb/src/portable/raspberrypi/rp2040; adding USB support.
-- Could NOT find Doxygen (missing: DOXYGEN_EXECUTABLE) 
ELF2UF2 will need to be built
-- Configuring done
-- Generating done
-- Build files have been written to: /home/alberto/tools/pico/proyectos/blinky

Si nos fijamos ahora en el directorio, veremos que se han creado una serie de archivos nuevos, entre ellos, por ejemplo, un Makefile. Así que nos resta sólo ejecutar make para que nuestro proyecto se compile y genere un ejecutable.

make

Y deberíamos ver algo como lo siguiente.

Scanning dependencies of target ELF2UF2Build
[  1%] Creating directories for 'ELF2UF2Build'
[  3%] No download step for 'ELF2UF2Build'
[  5%] No patch step for 'ELF2UF2Build'
[  6%] No update step for 'ELF2UF2Build'
[  8%] Performing configure step for 'ELF2UF2Build'
-- The C compiler identification is GNU 9.3.0
-- The CXX compiler identification is GNU 9.3.0
-- Check for working C compiler: /usr/bin/cc
-- Check for working C compiler: /usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /usr/bin/c++
-- Check for working CXX compiler: /usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Configuring done
-- Generating done
-- Build files have been written to: /home/alberto/tools/pico/proyectos/blinky/elf2uf2
[ 10%] Performing build step for 'ELF2UF2Build'
Scanning dependencies of target elf2uf2
[ 50%] Building CXX object CMakeFiles/elf2uf2.dir/main.cpp.o
[100%] Linking CXX executable elf2uf2
[100%] Built target elf2uf2
[ 11%] No install step for 'ELF2UF2Build'
[ 13%] Completed 'ELF2UF2Build'
[ 13%] Built target ELF2UF2Build
Scanning dependencies of target bs2_default
[ 15%] Building ASM object pico_sdk/src/rp2_common/boot_stage2/CMakeFiles/bs2_default.dir/boot2_w25q080.S.obj
[ 16%] Linking ASM executable bs2_default.elf
[ 16%] Built target bs2_default
Scanning dependencies of target bs2_default_bin
[ 18%] Generating bs2_default.bin
[ 18%] Built target bs2_default_bin
Scanning dependencies of target bs2_default_padded_checksummed_asm
[ 20%] Generating bs2_default_padded_checksummed.S
[ 20%] Built target bs2_default_padded_checksummed_asm
Scanning dependencies of target blink
[ 22%] Building C object CMakeFiles/blink.dir/blink.c.obj
[ 23%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_stdlib/stdlib.c.obj
[ 25%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_gpio/gpio.c.obj
[ 27%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_claim/claim.c.obj
[ 28%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_sync/sync.c.obj
[ 30%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_platform/platform.c.obj
[ 32%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_uart/uart.c.obj
[ 33%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_divider/divider.S.obj
[ 35%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_time/time.c.obj
[ 37%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_time/timeout_helper.c.obj
[ 38%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_timer/timer.c.obj
[ 40%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_sync/sem.c.obj
[ 42%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_sync/lock_core.c.obj
[ 44%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_sync/mutex.c.obj
[ 45%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_sync/critical_section.c.obj
[ 47%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_util/datetime.c.obj
[ 49%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_util/pheap.c.obj
[ 50%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/common/pico_util/queue.c.obj
[ 52%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_runtime/runtime.c.obj
[ 54%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_clocks/clocks.c.obj
[ 55%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_watchdog/watchdog.c.obj
[ 57%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_xosc/xosc.c.obj
[ 59%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_pll/pll.c.obj
[ 61%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_vreg/vreg.c.obj
[ 62%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_irq/irq.c.obj
[ 64%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/hardware_irq/irq_handler_chain.S.obj
[ 66%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_printf/printf.c.obj
[ 67%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_bit_ops/bit_ops_aeabi.S.obj
[ 69%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_bootrom/bootrom.c.obj
[ 71%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_divider/divider.S.obj
[ 72%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_double/double_aeabi.S.obj
[ 74%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_double/double_init_rom.c.obj
[ 76%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_double/double_math.c.obj
[ 77%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_double/double_v1_rom_shim.S.obj
[ 79%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_int64_ops/pico_int64_ops_aeabi.S.obj
[ 81%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_float/float_aeabi.S.obj
[ 83%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_float/float_init_rom.c.obj
[ 84%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_float/float_math.c.obj
[ 86%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_float/float_v1_rom_shim.S.obj
[ 88%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_malloc/pico_malloc.c.obj
[ 89%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_mem_ops/mem_ops_aeabi.S.obj
[ 91%] Building ASM object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_standard_link/crt0.S.obj
[ 93%] Building CXX object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_standard_link/new_delete.cpp.obj
[ 94%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_standard_link/binary_info.c.obj
[ 96%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_stdio/stdio.c.obj
[ 98%] Building C object CMakeFiles/blink.dir/home/alberto/tools/pico/pico-sdk/src/rp2_common/pico_stdio_uart/stdio_uart.c.obj
[100%] Linking CXX executable blink.elf
[100%] Built target blink

Si al final vemos el mensaje Built target blink es que todo ha ido bien y estamos listos para pasar el ejecutable al microcontrolador. Aquí también han sido originales los chicos de Raspberry. Para mover el ejecutable a la placa hay que seguir este orden. Primero pulsamos el pequeño botón que hay en la placa (hay que mantenerlo pulsado).

Seguidamente se conecta al ordenador mediante el cable USB y, tras esto, se suelta el botón. La tarjeta se monta como si fuera una unidad más en la que encontramos dos archivos.

El archivo INDEX.HTM nos lleva directos a la documentación y el archivo INFO_UF2.TXT contiene información sobre la versión de la tarjeta. Sólo tenemos que arrastrar el archivo blink.uf2 a esta ventana, e inmediatamente se desmontará la unidad y comenzará a ejecutarse. El LED de la placa debería empezar a parpadear como en el siguiente vídeo.

La intención de este artículo era simple, ver cómo poner en marcha un pequeño proyecto y programar la placa. Espero en un próximo artículo profundizar algo más y crear un proyecto algo más completo y realista. Hasta entonces, a disfrutar.

Manejo de grafos con NetworkX en Python

noreply@blogger.com (Alberto García Serrano) — Wed, 12 Feb 2020 11:55:00 +0000

El aprendizaje computacional es un área de investigación que en los últimos años ha tenido un auge importante, sobre todo gracias al aprendizaje profundo (Deep Learning). Pero no todo son redes neuronales. Paralelamente a estas técnicas, más bien basadas en el aprendizaje de patrones, también hay un auge de otras técnicas, digamos, más basadas en el aprendizaje simbólico. Si echamos la vista algunos años atrás, podemos considerar que quizá, la promesa de la web semántica como gran base de conocimiento ha fracasado, pero no es tan así. Ha ido transmutándose y evolucionando hacia bases de conocimiento basadas en ontologías a partir de las cuales es posible obtener nuevo conocimiento. Es lo que llamamos razonamiento automático y empresas como Google ya lo utilizan para ofrecerte información adicional sobre tus búsquedas. Ellos lo llaman Grafos de Conocimiento o Knowledge Graphs.

Gracias a estos grafos de conocimiento, Google puede ofrecerte información adicional sobre tu búsqueda, además de los enlaces a las páginas web. Por ejemplo, en la imagen siguiente vemos que si buscamos Marie Curie, Google nos ofrece datos sobre Marie, y otras personas con las que está relacionada.

Para entender como funcionan estos grafos de conocimiento necesitamos tener algunos conceptos previos bien asentados de los que quiero hablaros hoy. Así que para empezar y, como toda esta información se almacena en grandes grafos, hoy vengo hablaros de eso, de grafos.

¿Qué es un grafo?

Matemáticamente, un grafo es un par ordenado $G=(V,A)$ donde V es un conjunto de vértices (o nodos) y A un conjunto de aristas que relacionan elementos entre sí. Gráficamente se representan como círculos (que son los vértices) unidos por líneas (que son las aristas). Mediante los grafos podemos describir multitud de procesos, situaciones, relaciones, estructuras, etc. Por ejemplo, un grafo nos puede servir para representar los enlaces moleculares de un elemento químico.

Hay situaciones en las que nos interesa especificar que una relación (arista) tiene una dirección concreta. Por ejemplo, en el siguiente grafo se representan los diferentes estados de un semáforo, cuyas aristas tienen definidas una dirección concreta. Por ejemplo, es posible pasar del color amarillo al rojo, pero no al contrario. Este tipo de grafos se llaman digrafos o grafos dirigidos.

Las aristas, además de ser dirigidas, pueden tener un valor o un peso que aporta cierta información acerca de las relaciones. Son los llamados grafos ponderados. En el siguiente grafo, que representa una red de carreteras, vemos que cada arista contiene la distancia en kilómetros que separa cada ciudad. Esto es útil para, por ejemplo, poder calcular la ruta más corta entre dos puntos del grafo.

Os resumo algunas de las métricas básicas de un grafo: el orden del grafo es el número vértices que lo forman, mientras que la medida del grafo es el número de aristas.
El grado de un vértice es el número de aristas que inciden en él. Coincide con el número de vértices adyacentes (por ejemplo, un vértice de grado=0 sería un vértice aislado).
La teoría de grafos es un campo amplio de las matemáticas, así que en vez de aburriros con más términos, iremos definiendo el resto de conceptos que vayamos necesitando sobre la marcha.

Manejo de grafos con NetworkX en Python

Hechas las presentaciones formales, vamos a remangarnos y a ver cómo podemos gestionar grafos en Python, para lo que nos valdremos de la librería NetworkX (hay otras, pero como había que elegir una me he quedado con la que me parecía más sencilla, aunque no necesariamente mejor que las otras). Con el siguiente código vamos a crear un pequeño grafo que relaciona a Kevin Bacon con otros actores con los que ha trabajado (basándonos en la curiosa teoría de los seis grados de separación que cualquier actor tiene con Kevin Bacon).

import networkx as nx

G = nx.Graph() # crear un grafo

#Añadir nodos
G.add_node("Kevin Bacon")
G.add_node("Tom Hanks")
G.add_nodes_from(["Meg Ryan", "Parker Posey", "Lisa Kudrow"])

#Añadir aristas
G.add_edge("Kevin Bacon", "Tom Hanks")
G.add_edge("Kevin Bacon", "Meg Ryan")
G.add_edges_from([("Tom Hanks", "Meg Ryan"), ("Tom Hanks", "Parker Posey")])
G.add_edges_from([("Parker Posey", "Meg Ryan"), ("Parker Posey", "Lisa Kudrow")])

Con nx.Graph() creamos el grafo (podemos usar nx.DiGraph() para crear un grafo dirigido). Para añadir nodos usamos add_node() para añadir un sólo nodo o add_nodes_from() para añadir varios de una vez. Para las aristas usamos add_edge() o add_edges_from() dependiendo de si queremos agregar una o varias aristas.
Si queremos saber el número de nodos y aristas, o qué elementos lo componen, podemos hacerlo fácilmente con:

print(len(G.nodes))
print(len(G.edges))
print(G.nodes)
print(G.edges)

5
6
['Kevin Bacon', 'Tom Hanks', 'Meg Ryan', 'Parker Posey', 'Lisa Kudrow']
[('Kevin Bacon', 'Tom Hanks'), ('Kevin Bacon', 'Meg Ryan'), ('Tom Hanks', 'Meg Ryan'), ('Tom Hanks', 'Parker Posey'), ('Meg Ryan', 'Parker Posey'), ('Parker Posey', 'Lisa Kudrow')]

Es posible asignar atributos tanto a los nodos como a las aristas. Por ejemplo, podemos especificar que Tom Hanks tiene dos Oscars. Respecto a las aristas, como representan que un actor ha trabajado con otro, podemos añadir un atributo que sea la película en la que ambos han coincidido. Un uso muy habitual que tienen los atributos sobre las aristas es indicar el peso si se trata de un grafo ponderado, mediante el atributo weight.

# asignar atributos a nodos y aristas
G.nodes["Tom Hanks"]["oscars"] = 2
G.edges["Kevin Bacon", "Tom Hanks"]["pelicula"] = "Apolo 13"
G.edges["Kevin Bacon", "Meg Ryan"]["pelicula"] = "En carne viva"
G.edges["Parker Posey", "Meg Ryan"]["pelicula"] = "Algo para recordar"
G.edges["Parker Posey", "Tom Hanks"]["pelicula"] = "Tienes un email"
G.edges["Parker Posey", "Lisa Kudrow"]["pelicula"] = "Esperando la hora"
print(G.nodes(data=True))
print(G.edges(data=True))
print(G["Kevin Bacon"])

[('Kevin Bacon', {}), ('Tom Hanks', {'oscars': 2}), ('Meg Ryan', {}), ('Parker Posey', {}), ('Lisa Kudrow', {})]
[('Kevin Bacon', 'Tom Hanks', {'pelicula': 'Apolo 13'}), ('Kevin Bacon', 'Meg Ryan', {'pelicula': 'En carne viva'}), ('Tom Hanks', 'Meg Ryan', {}), ('Tom Hanks', 'Parker Posey', {'pelicula': 'Tienes un email'}), ('Meg Ryan', 'Parker Posey', {'pelicula': 'Algo para recordar'}), ('Parker Posey', 'Lisa Kudrow', {'pelicula': 'Esperando la hora'})]
{'Tom Hanks': {'pelicula': 'Apolo 13'}, 'Meg Ryan': {'pelicula': 'En carne viva'}}

Si quieremos mostrar gráficamente el resultado podemos usar graphviz.

import graphviz

A = nx.nx_agraph.to_agraph(G)
A.layout('dot')
#A.draw('salida.png') # guardar como png

graphviz.Source(A.to_string()) # mostrar en jupyter

Si queremos almacenar nuestro grafo en disco, NetworX tiene una familia de funicones write_*() que facilitan esta operación. Aunque el más habitual es graphml se pueden almacenar en diferentes formatos (como GML, JSON, Pajek, YAML, etc.)

# guardar grafo
nx.write_graphml(G, "actores.graphml")

Para cargar el grafo usamos

# cargar grafo
G2 = nx.read_graphml("actores.graphml")
print(G2.nodes)

Os he contado los rudimentos básicos del manejo de la librería NetworkX. Obviamente esto es sólo una pequeña parte de lo que puede hacer, pero nos vale como introducción. En próximos artículos os seguiré hablando de grafos y os contaré más cosas interesantes.

Regresión polinómica para Machine Learning

noreply@blogger.com (Alberto García Serrano) — Thu, 05 Dec 2019 14:04:00 +0000

La regresión lineal es una herramienta extremadamente potente para realizar análisis de datos y construir modelos de Machine Learning capaces de hacer inferencias sobre datos nuevos. Sin embargo, esto es sólo cierto para datos que tienen un comportamiento lineal, lo que no siempre ocurre. Vamos a trabajar con un dataset sencillo, pero cuyos datos no son lineales: https://www.kaggle.com/akram24/position-salaries
Este dataset contiene los salarios según el nivel del empleado en dentro de una empresa.

Tras descargarlo de Kaggle y descomprimir el fichero tenemos un archivo CSV que podemos cargar con Pandas.

import pandas as pd

datos = pd.read_csv('Position_Salaries.csv')
print (datos)

    Position            Level  Salary
0   Business Analyst      1    45000
1  Junior Consultant      2    50000
2  Senior Consultant      3    60000
3            Manager      4    80000
4    Country Manager      5   110000
5     Region Manager      6   150000
6            Partner      7   200000
7     Senior Partner      8   300000
8            C-level      9   500000
9                CEO     10  1000000

Son muy pocos datos, lo que hace a este conjunto de datos interesante para usarlo como ejemplo. Como siempre, empezamos por el principio. Vamos a visualizar qué forma tienen los datos.

from matplotlib import pyplot as plt

x = datos['Level'].values.reshape(-1, 1) # necesitamos un array de 2D para SkLearn
y = datos['Salary'].values.reshape(-1, 1)
plt.scatter(x,y)

Vamos a ver qué tal se ajusta a un modelo lineal, para lo cuál aplicamos una regresión lineal valiéndonos de la librería SkLearn.

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(x, y)
y_pred = model.predict(x)

plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

No parece que el modelo lineal se ajuste demasiado bien. Una forma de medir la bondad del ajuste es calcular la media de la raíz del error cuadrático (root mean square error) que nos da una medida del error cometido por el modelo. En concreto se calcula la media de la desviación de los valores estimados. Otra métrica interesante es la medida R2 (R cuadrado), cuyo valor está entre 0 y 1, lo que la hace mejor a la hora de interpretar su valor, y es la fracción de la suma total de cuadrados que se 'explica por' la regresión.

import numpy as np
from sklearn.metrics import mean_squared_error, r2_score

rmse = np.sqrt(mean_squared_error(y,y_pred))
r2 = r2_score(y,y_pred)
print ('RMSE: ' + str(rmse))
print ('R2: ' + str(r2))

RMSE: 163388.73519272613
R2: 0.6690412331929895

Si nos fijamos, sobre todo en R2, vemos que el ajuste no es bueno (mientras más cercano a 1, mejor ajuste). Estaremos de acuerdo en que una curva se ajustaría mejor. Si queremos ajustar una curva a los datos tendremos que trabajar con más dimensiones. Es decir, tendremos que intentar ajustar un polinomio (de segundo grado, por ejemplo). Para ello fijémonos en la ecuación de la recta que hemos ajustado con la regresión lineal.
\begin{align}
y=a+bx
\end{align}
Si quisieramos ajustar los datos con una función de segundo grado, de decir, una curva, necesitamos añadir un término cuadrático. Recordemos que x es un atributo de nuestros datos, así que vamos a añadir a la función este mismo atributo pero al cuadrado.
\begin{align}
y=a+bx+cx^2
\end{align}
Esta transformación la podemos realizar fácilmente son SkLearn gracias PolynomialFeatures(). A la que le indicamos el grado de la función que queremos obtener. Haciendo uso del método fit_transform() obtendremos el término cuadrático que andamos buscando para nuestro modelo.

from sklearn.preprocessing import PolynomialFeatures

poly = PolynomialFeatures(degree=2, include_bias=False)
x_poly = poly.fit_transform(x)
print(x)
print(x_poly)

[[ 1]
 [ 2]
 [ 3]
 [ 4]
 [ 5]
 [ 6]
 [ 7]
 [ 8]
 [ 9]
 [10]]

[[  1.   1.]
 [  2.   4.]
 [  3.   9.]
 [  4.  16.]
 [  5.  25.]
 [  6.  36.]
 [  7.  49.]
 [  8.  64.]
 [  9.  81.]
 [ 10. 100.]]

Vemos que lo que ha hecho esta función es transformar una función lineal en otra cuadrática simplemente añadiendo el cuadrado del primer atributo. Vamos a calcular la regresión sobre esta nueva función y a comprobar si se ajusta mejor.

model.fit(x_poly, y)
y_pred = model.predict(x_poly)

plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

rmse = np.sqrt(mean_squared_error(y,y_pred))
r2 = r2_score(y,y_pred)
print ('RMSE: ' + str(rmse))
print ('R2: ' + str(r2))

RMSE: 82212.12400451244
R2: 0.9162082221443943

Definitivamente el ajuste es mucho mejor, y el valor de R2 está mucho más cerca de 1. Si te estás preguntando qué ocurriría si en vez de una función cuadrática usáramos una función cúbica, no te cortes y pruébalo.

poly = PolynomialFeatures(degree=3, include_bias=False)
x_poly = poly.fit_transform(x)

model.fit(x_poly, y)
y_pred = model.predict(x_poly)

plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

rmse = np.sqrt(mean_squared_error(y,y_pred))
r2 = r2_score(y,y_pred)
print ('RMSE: ' + str(rmse))
print ('R2: ' + str(r2))

RMSE: 38931.50401232937
R2: 0.9812097727913365

Mucho mejor. ¿Y si seguimos aumentando el grado?

poly = PolynomialFeatures(degree=4, include_bias=False)
x_poly = poly.fit_transform(x)

model.fit(x_poly, y)
y_pred = model.predict(x_poly)

plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

rmse = np.sqrt(mean_squared_error(y,y_pred))
r2 = r2_score(y,y_pred)
print ('RMSE: ' + str(rmse))
print ('R2: ' + str(r2))

RMSE: 14503.2349096276
R2: 0.9973922891706611

Guau, el ajuste es casi perfecto. Pues ya puestos vamos a usar una función de grado 100 y será mucho mejor ¿no?
En realidad no. Como no tenemos datos de validación, estamos un poco a ciegas en cuanto a comprobar si nuestro modelo está sobreaprendiendo (overfitting), aunque un ajuste tan perfecto nos tiene que hacer sospechar. Si tuviéramos que usar el modelo para hacer predicciones sobre datos nuevos, seguramente las dos medidas que hemos usado para comprobar la bondad del ajuste (RMSE y R2) empezarían a degradarse considerablemente, ya que lo que ha hecho nuestro modelo es "aprenderse" los datos de entrenamiento, lo que significa que generalizará mal con nuevos datos. En un caso como este, es mejor un modelo más genérico (como el cuadrático) que uno que ajuste demasiado bien los datos de entrenamiento.

Como siempre, os dejo un enlace al notebook de Jupyter con el código del artículo: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/RegresionPolinomica2.ipynb

Análisis de sentimiento y procesado del lenguaje natural

noreply@blogger.com (Alberto García Serrano) — Wed, 13 Nov 2019 19:15:00 +0000

El análisis de sentimiento es una técnica de machine learning, basada en el procesado del lenguaje natural, que pretende obtener información subjetiva de una serie de textos o documentos. Un ejemplo clásico de aplicación consiste en dilucidar si un artículo periodístico es favorable o no favorable en relación a un tema determinado. Las aplicaciones en el mundo real son muchas y variadas. No hace mucho leía un artículo sobre la aplicación de estas técnicas para usarlas como un dato más a la hora de predecir índices bursátiles, analizando lo que decía la prensa económica sobre determinadas empresas. El campo del tratamiento del lenguaje natural es una área de estudio bastante activa actualmente, y utiliza técnicas que pueden llegar a ser bastante complejas. Sin embargo, no es necesario complicarse mucho la vida para obtener resultados bastante decentes.

Como ejemplo, vamos a tratar de predecir si un mensaje de Twitter es positivo o negativo. Como datos de entrenamiento usaremos un dataset con 1.600.000 tweets (no los usaremos todos) con su correspondiente etiqueta. Fuente del dataset: https://www.kaggle.com/kazanova/sentiment140. Se puede abordar el problema desde diferentes enfoques, pero en este caso entrenaremos una red neuronal con los datos del dataset. El primer paso será descargar el conjunto de datos y cargarlo en memoria con Pandas.

import pandas as pd

columnas = ['target','ids','date','flag','user','text']
data = pd.read_csv('dataset/training.1600000.processed.noemoticon.csv', encoding = "latin-1", names=columnas)
# nos quedamos con una selección de 50.000 tweets aleatorios
data = data.sample(n=50000, random_state=42)
data = data.reset_index(drop=True) # reconstruimos los índices para que sean consecutivos
data.head()

	ids	date	flag	user	text
0	2200003196	Tue Jun 16 18:18:12 PDT 2009	NO_QUERY	LaLaLindsey0609	@chrishasboobs AHHH I HOPE YOUR OK!!!
1	1467998485	Mon Apr 06 23:11:14 PDT 2009	NO_QUERY	sexygrneyes	@misstoriblack cool , i have no tweet apps fo...
2	2300048954	Tue Jun 23 13:40:11 PDT 2009	NO_QUERY	sammydearr	@TiannaChaos i know just family drama. its la...
3	1993474027	Mon Jun 01 10:26:07 PDT 2009	NO_QUERY	Lamb_Leanne	School email won't open and I have geography ...
4	2256550904	Sat Jun 20 12:56:51 PDT 2009	NO_QUERY	yogicerdito	upper airways problem

Los campos son los siguientes (extraído directamente de la documentación del dataset):

target: the polarity of the tweet (0 = negative, 2 = neutral, 4 = positive)
ids: The id of the tweet ( 2087)
date: the date of the tweet (Sat May 16 23:58:44 UTC 2009)
flag: The query (lyx). If there is no query, then this value is NO_QUERY.
user: the user that tweeted (robotickilldozr)
text: the text of the tweet (Lyx is cool)

Para este ejemplo voy a usar sólo la columna text (contenido del tweet) y target (la etiqueta para entrenar la red neuronal).
Nos hemos quedado sólo con 50.000 tweets aleatorios, ya que ocuparían demasiada memoria para un ordenador de andar por casa.
Vamos a echar un vistazo a cómo se reparten las etiquetas (campo target).

from matplotlib import pyplot as plt
plt.hist(data['target'], bins=5)

Parece que el dataset sólo etiqueta los tweets positivos(=0) y los negativos(=4), ya que no hay ningún tweet etiquetado como neutral. Además, para alimentar la red neuronal, es más conveniente que tengan valores continuos y que no haya valores que no tengan ningún significado (como el 1 o el 3). Vamos a reasignar las etiquetas como 0=negativo y 1=positivo.

# ahora la etiqueta para los tweets positivos es el 1
data.loc[data['target'] == 4, 'target'] = 1
plt.hist(data['target'], bins=3)

Además vemos que hay aproximadamente el mismo número de ejemplos de tweets negativos y positivos, lo que es muy interesante a la hora de entrenar la red.

Tratamiento de los datos

Vamos a diseñar una red neuronal capaz de aprender si un texto es positivo o negativo (sentimiento). Para ello necesitamos codificar el texto de alguna manera para que la red neuronal, que trabaja con valores numéricos preferentemente, sea capaz de entenderlo. Existen varias técnicas enmarcadas en lo que se conoce como tratamiento del lenguaje natural. En este ejemplo voy a usar una técnica más o menos sencilla pero efectiva. La técnica consiste en codificar un conjunto amplio de palabras (corpus) como valores numéricos, es decir, creamos un índice en el que hacemos corresponder un valor a cada palabra. Cada palabra se codificará luego como un vector donde la posición correspondiente al valor de la palabra contendrá un 1 y las demás posiciones un 0. Es decir, usando una codificación one-hot. Existen corpus de palabras ya disponibles para usar, pero vamos a complicarnos un poco más la vida y vamos a crear el nuestro a partir del dataset de tweets del que disponemos. La siguiente función construye un corpus de palabras a partir de los tweets. Tenemos cuidado de no añadir las referencias a los usuarios (empiezan por @), los enlaces (contienen http) y eliminamos las palabras menores de 3 caracteres, que no aportan demasiada información.

from keras.preprocessing.text import Tokenizer

# codificar tweets -> corpus -> one-hot
# con limit podemos limitar el número de palabras del corpus
def construir_corpus(tweets, limit=10000):
    corpus=[]
    i=0
    for t in tweets:
        if i>limit:
            break
        for w in t.lower().split():
            # si la palabra no está en el corpus y no empieza por @, no contiene http y es mayor de 3 caracteres.
            if w[0] not in ('@') and 'http' not in w and len(w)>3 and not w in corpus:
                corpus.append(w)
                i=i+1
    
    # codificar palabras como enteros
    t = Tokenizer()
    t.fit_on_texts(corpus)
    encoded_corpus = t.texts_to_matrix(corpus, mode='count')
    
    return corpus[:limit], encoded_corpus[:limit]

Vamos a probarlo.

# Probamos con tres frases para confirmar que todo va bien
tweets = ['Sólo sé que no sé nada', 'Pienso luego existo', 'viva @socrates y @descartes']

corpus, encoded_corpus = construir_corpus(tweets)
print (corpus)
print (encoded_corpus)

['sólo', 'nada', 'pienso', 'luego', 'existo', 'viva']
[[0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 1.]]

Se ha construido el siguiente corpus a partir de las tres frases: ['sólo', 'nada', 'pienso', 'luego', 'existo', 'viva'].
También se muestra como se codifica cada palabra. por ejemplo, la primera (sólo) se codifica como: [0. 1. 0. 0. 0. 0. 0.].
Obsérvese que si tenemos 10.000 palabras, nuestra matriz de codificación tendrá 10.000 filas. En realidad sólo las he calculado para que veas qué forma tiene la codificación. Como verás más adelante, no vamos a usar esta matriz para nada.
Ahora vamos a construir el corpus real que vamos a utilizar.

# Ahora usamos los tweets del dataset
tweets = data['text']
corpus, encoded_corpus = construir_corpus(tweets)

Para codificar un tweet, por cada palabra de la que esté compuesto, se pone un 1 en la posición correspondiente a dicha palabra en el vector. Así, conforme al ejemplo de arriba (el de Sócrates y Descartes), el tweet con el texto "Sólo pienso" se codificaría como [0. 1. 0. 1. 0. 0. 0.]. Una vez entendida la idea de cómo se codifica un tweet para que sea útil a la hora de alimentar una red neuronal, vamos a crear una función que codifique todo el conjunto de tweets. Y también necesitamos codificar las etiquetas para el entrenamiento, para lo cuál también vamos a usar una codificación one-hot valiéndonos de la función to_categorical() de Keras.

import numpy as np

from keras.utils import to_categorical

def codifica_tweets(tweets, corpus, corpus_size=10000):
    coded = np.zeros((len(tweets), corpus_size))
    for i, tweet_text in enumerate(tweets):
        words = tweet_text.lower().split()
        for w in words:
            if w in corpus:
                coded[i,corpus.index(w)] = 1

    return coded

x_train = codifica_tweets(data['text'], corpus)
y_train = to_categorical(data['target'])

Seguidamente prepararemos los datos con los que vamos a entrenar la red neuronal, separando del conjunto de entrenaiento una porción de los datos para crear el conjunto de validación y de test. También vamos a barajar los datos para añadir un poco de entropía.

# barajamos el dataset
np.random.seed(42)    
permutation = np.random.permutation(x_train.shape[0])
x_train = x_train[permutation]
y_train = y_train[permutation]
     
# obtenemos el conjunto de validación
num_val = 10000
x_val = x_train[:num_val]
x_train = x_train[num_val:]
y_val = y_train[:num_val]
y_train = y_train[num_val:]

# obtenemos el conjunto de test
x_test = x_train[:num_val]
x_train = x_train[num_val:]
y_test = y_train[:num_val]
y_train = y_train[num_val:]

Una vez preparados los datos, vamos a definir la arquitectura de la red neuronal. Como los tweets se codifican en vectores de dimensión 10.000 (el tamaño de nuestro corpus de palabras) y queremos que clasifique los tweets en dos categorías (positivos y negativos), nuestra red tendra 10.000 entradas y dos salidas en la última capa. En esta última capa usaremos una función de activación de tipo softmax, que nos devolverá una distribución de probabilidad entre las dos salidas (una para la positiva y otra para la negativa), indicándonos cuál es la probabilidad de que corresponda a cada una de las dos clases. Usaremos dos capas ocultas de 16 neuronas densamente conectadas (todas con todas).

from keras import models, layers

model = models.Sequential()
model.add(layers.Dense(16, activation='relu', input_shape=(10000,)))
model.add(layers.Dense(16, activation='relu'))
model.add(layers.Dense(2, activation='softmax'))

model.compile(optimizer='rmsprop',
             loss='categorical_crossentropy',
             metrics=['acc'])

train_log = model.fit(x_train, y_train,
                     epochs=5, batch_size=512,
                     validation_data=(x_val, y_val))

Train on 30000 samples, validate on 10000 samples
Epoch 1/5
30000/30000 [==============================] - 1s 45us/step - loss: 0.6503 - acc: 0.6775 - val_loss: 0.5995 - val_acc: 0.7202
Epoch 2/5
30000/30000 [==============================] - 1s 36us/step - loss: 0.5471 - acc: 0.7519 - val_loss: 0.5465 - val_acc: 0.7296
Epoch 3/5
30000/30000 [==============================] - 1s 40us/step - loss: 0.4922 - acc: 0.7733 - val_loss: 0.5416 - val_acc: 0.7338
Epoch 4/5
30000/30000 [==============================] - 1s 35us/step - loss: 0.4658 - acc: 0.7841 - val_loss: 0.5487 - val_acc: 0.7313
Epoch 5/5
30000/30000 [==============================] - 1s 36us/step - loss: 0.4471 - acc: 0.7949 - val_loss: 0.5568 - val_acc: 0.7285

Si miramos la salida durante el entrenamiento, vemos que, con el conjunto de validación, los mejores resultados se producen en la iteración 3/5. Después parece que se produce un sobreentrenamiento (overfitting). Si queremos asegurarnos podemos mostrar una gráfica.

import matplotlib.pyplot as plt

acc = train_log.history['acc']
val_acc = train_log.history['val_acc']
loss = train_log.history['loss']
val_loss = train_log.history['val_loss']
epochs = range(1, len(acc) + 1)
plt.plot(epochs, acc, 'r', label='Training acc')
plt.plot(epochs, val_acc, 'b', label='Validation acc')
plt.title('Training and validation accuracy')
plt.legend()
plt.figure()
plt.plot(epochs, loss, 'r', label='Training loss')
plt.plot(epochs, val_loss, 'b', label='Validation loss')
plt.title('Training and validation loss')
plt.legend()
plt.show()

Vamos a evaluar el modelo contra el conjunto de test.

test_loss, test_accuracy = model.evaluate(x_test, y_test)
print(test_accuracy)

10000/10000 [==============================] - 1s 68us/step
0.7278000116348267

Hemos conseguido un 72% de precisión. Si hubieramos cortado el entrenamiento en la tercera iteración seguramente hubieramos llegado al 73% (te dejo que lo pruebes). No es un mal dato teniendo en cuenta que sólo hemos usado una pequeña fracción de los tweets y el corpus que hemos generado tampoco es demasiado bueno. Con esta técnica, usando todo el dataset y un buen corpus de palabras podríamos esperar superar el 80%, pero como primera aproximación no está mal.
Ya con el modelo entrenado podemos hacer predicciones sobre cualquier tweet.

predictions = model.predict(x_test)
print (np.argmax(predictions[1]))# clase más probable
print (y_test[1])

0
[1. 0.]

Como os decía al principio, el procesado del lenguaje natural es un campo amplio y actualmente se están realizando bastantes esfuerzos de investigación en esta materia. Como ves, con técnicas más o menos simples podemos acomenter problemas de cierta complejidad. Espero haberte despertado la curiosidad para que te atrevas con técnicas más avanzadas y también más interesantes.
Os dejo un enlace al notebook de Jupyter con el código de este artículo: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/SentimentAnalysis.ipynb

Visión artificial con redes convolucionales (CNN)

noreply@blogger.com (Alberto García Serrano) — Mon, 21 Oct 2019 10:46:00 +0000

Las redes neuronales convolucionales o ConvNets son útiles en variadas aplicaciones dentro del machine learning, sin embargo, donde ha conseguido revolucionar el estado del arte es en la visión artificial y el procesado de imágenes. Lo que caracteriza a este tipo de redes es que es capaz de aprender en las primeras capas una serie de características básicas de la imagen, como son líneas o formas más o menos simples. En posteriores fases aprende a discriminar entre elementos más complejos, como ojos u orejas. Finalmente, en capas más profundas son capaces de diferenciar objetos complejos, como personas, perros, gatos o coches. En la siguiente imagen se plasma visualmente lo que quiero decir.

Antes de hablar de las redes convolucionales convendrá explicar qué es eso de la convolución y para qué sirve. En el ámbito del tratamiento de imágenes, la convolución es una operación que involucra dos matrices. Por un lado la propia imagen a la que queremos aplicar la convolución (en forma de matriz, claro), y otra matriz más pequeña llamada kernel. El resultado de la convolución es otra matriz de las mismas dimensiones que la imagen original. Para calcular la matriz de convolución vamos desplazando el kernel por cada uno de los elementos de la matriz principal y hacemos la suma de los productos de cada uno de los elementos de ambas matrices. En realidad suena peor de los que es. Visualmente podemos imaginarlo así.

El centro del kernel, pues pasaría por cada uno de los elementos de la matriz para ir generando la matriz de convolución. Mejor verlo con un ejemplo.

$$
\begin{vmatrix}
P_{1} & P_{2} & P_{3} & P_{4} & P_{5} \\
P_{6} & P_{7} & P_{8} & P_{9} & P_{10} \\
P_{11} & P_{12} & P_{13} & P_{14} & P_{15} \\
P_{16} & P_{17} & P_{18} & P_{19} & P_{20} \\
P_{21} & P_{22} & P_{23} & P_{24} & P_{25}
\end{vmatrix}
\circledast
\begin{vmatrix}
A & B & C \\
D & E & F \\
G & H & I
\end{vmatrix}
=
\begin{vmatrix}
C_{1} & C_{2} & C_{3} & C_{4} & C_{5} \\
C_{6} & C_{7} & C_{8} & C_{9} & C_{10} \\
C_{11} & C_{12} & C_{13} & C_{14} & C_{15} \\
C_{16} & C_{17} & C_{18} & C_{19} & C_{20} \\
C_{21} & C_{22} & C_{23} & C_{24} & C_{25}
\end{vmatrix}
$$

Para calcular, por ejemplo, el elemento $C_{13}$ de la matriz de convolución haríamos el siguiente cálculo (y así con todos los elementos).

$$C_{13}=A \cdot P_{7} + B \cdot P_{8} + C \cdot P_{9} + D \cdot P_{12} + E \cdot P_{13} + F \cdot C_{14} + G \cdot P_{17} + H \cdot P_{18} + I \cdot P_{19}$$

Un caso especial es cuando aplicamos la convolución sobre un elemento que está en el borde de la matriz. En ese caso, lo más común, es imaginarnos que la matriz está rodeada de ceros.
¿Y para qué sirve esto? Podemos imaginar un kernel como una especie de filtro que nos permite hacer cosas como detectar en las imágenes bordes, líneas verticales y horizontales, así como desenfocar y enfocar la imagen. Algunos kernels interesantes son los siguientes.

Identidad

$$
\begin{bmatrix}
0 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 0
\end{bmatrix}
$$

Detección de bordes

$$ \begin{bmatrix} \ \ 1 & 0 & -1 \\ \ \ 0 & 0 & \ \ 0 \\ -1 & 0 & \ \ 1 \end{bmatrix} $$
$$ \begin{bmatrix} 0 & \ \ 1 & 0 \\ 1 & -4 & 1 \\ 0 & \ \ 1 & 0 \end{bmatrix} $$
$$ \begin{bmatrix} -1 & -1 & -1 \\ -1 & \ \ 8 & -1 \\ -1 & -1 & -1 \end{bmatrix} $$

Afilado

$$
\begin{bmatrix}
\ \ 0 & -1 & \ \ 0 \\
-1 & \ \ 5 & -1 \\
\ \ 0 & -1 & \ \ 0
\end{bmatrix}
$$

Difuminado

$$
\frac{1}{9}
\begin{bmatrix}
1 & 1 & 1 \\
1 & 1 & 1 \\
1 & 1 & 1
\end{bmatrix}
$$

Desenfoque gausiano 3x3

$$
\frac{1}{16}
\begin{bmatrix}
1 & 2 & 1 \\
2 & 4 & 2 \\
1 & 2 & 1
\end{bmatrix}
$$

Desenfoque gausiano 5x5

$$
\frac{1}{256}
\begin{bmatrix}
1 & 4 & 6 & 4 & 1 \\
4 & 16 & 24 & 16 & 4 \\
6 & 24 & 36 & 24 & 6 \\
4 & 16 & 24 & 16 & 4 \\
1 & 4 & 6 & 4 & 1
\end{bmatrix}
$$

Enfoque

$$
\frac{-1}{256}
\begin{bmatrix}
1 & 4 & \ \ 6 & 4 & 1 \\
4 & 16 & \ \ 24 & 16 & 4 \\
6 & 24 & -476 & 24 & 6 \\
4 & 16 & \ \ 24 & 16 & 4 \\
1 & 4 & \ \ 6 & 4 & 1
\end{bmatrix}
$$

En python podemos calcular la matriz de convolución a mano (no es complicado), pero la librería SciPy ya lo hace por nosotros.

from scipy import signal

A=[[255, 7, 3],
   [212, 240, 4],
   [218, 216, 230]]

B= [[1, -1],
    [-1,1]]

conv = signal.convolve(A, B)
print(conv)

[[ 255 -248   -4   -3]
 [ -43  276 -232   -1]
 [   6  -30  250 -226]
 [-218    2  -14  230]]

Con la librería PIL podemos aplicar la operación de convolución a una imagen. Primero la cargamos.

from PIL import Image, ImageFilter
import numpy as np

imgpath = '/home/alberto/Imágenes/imagen.png'
img = Image.open(imgpath)
img = img.convert('L')
display(img)

Hemos convertido la imagen a escala de grises, por lo que nuestra imagen se puede operar como si fuera una sola matriz. Si fuera una imagen en color tendríamos tres matrices, una por cada canal de color RGB.

imgmatrix = np.asarray(img, dtype=np.uint8)
print (imgmatrix)

[[156 164 156 ... 164 156 132]
 [164 164 164 ... 172 156 124]
 [164 156 164 ... 164 156 132]
 ...
 [ 44  44  52 ... 108 100 100]
 [ 44  44  60 ... 100 108 108]
 [ 44  44  52 ... 108 100 108]]

Así pues, como tenemos una matriz entre manos, podemos aplicar la convolución tal que así.

kernelValues = [-2,-1,0,-1,1,1,0,1,2] 
kernel = ImageFilter.Kernel((3,3), kernelValues)
 
im2 = img.filter(kernel)

display(im2)

Una vez tenemos claro qué hace la operación de convolución sobre una imagen, ya podemos hablar de las redes neuronales convolucionales o CNN (Convolutional Neural Networs).

Redes neuronales convolucionales (CNN).

La arquitectura típica de una CNN dedicada al reconocimiento y clasificación de imágenes, está dividida en dos partes: La etapa convolucional y la etapa de clasificación (las CNN pueden usarse en otras aplicaciones, pero aquí os mostraré su uso en la clasificación de imágenes). La fase convolucional utiliza dos tipos de capas: de convolución y de max pooling. Las capas de convolución que en Keras se crean con la función Conv2D() del paquete layers de Keras. Veamos un ejemplo.

Conv2D(64, (3, 3), activation='relu')

El primer parámetro es el número de filtros (kernels) que se va a aplicar a la imagen, y el segundo parámetro es el tamaño de los kernels. Como hemos visto, cada kernel genera un conjunto de características en la imagen inicial. Por ejemplo, algunos kernels detectarán las líneas verticales, otros las horizontales, otros las esquinas, etc. Como resultado vamos a obtener una serie de imágenes (tantas como kernels apliquemos) en las que la red será capaz de reconocer formas simples como líneas y otros tipos de características. Si vamos repitiendo el proceso, y pasamos el resultado por otra capa convolucional, esta será capar de ir reconociendo elementos más complejos, como una oreja o un ojo, y así hasta ser capaz de reconocer un gato o un perro. Evidentemente esto es una simplificación, pero es una idea bastante intuitiva de cómo funciona una red convolucional. Es habitual usar la función de activación relu en este tipo de redes.
El otro tipo de capa que usa una red convolucional para reconocimiento de imágenes es MaxPooling2D().

MaxPooling2D(pool_size=(2, 2))

Esta capa aplica un filtro que reduce la dimensionalidad de la imagen, esto es, la hace más pequeña. El objetivo es, por un lado, reducir el coste computacional. Además se minimiza la posibilidad de overfitting y se consigue aumentar la abstracción sobre los datos de entrada. Lo más fácil será verlo con un ejemplo.

La imagen muestra el resultado de aplicar una operación de max pooling de 2x2. Esto es, se van cogiendo grupos de 2x2 y nos quedamos con el máximo valor del grupo de cuatro píxeles. De esta forma reducimos la imagen que originariamente era de 4x4 a una imagen de 2x2.
Una vez aclarados estos conceptos, vamos a ponernos manos a la obra.

Clasificar las imágenes del dataset CIFAR10

Vamos a usar el dataset CIFAR10. Este dataset contiene 60000 imágenes de tamaño 32x32 dividido en 10 clases diferentes (6000 imágenes por clase). El dataset se divide en 50000 imágenes para entrenamiento y 10000 para test.
Las clases son:

0 - airplane
1 - automobile
2 - bird
3 - cat
4 - deer
5 - dog
6 - frog
7 - horse
8 - ship
9 - truck

Obtendremos el dataset importándolo desde Keras, que ya lo incorpora, por lo que no necesitamos descargarlo desde la web.

from keras.datasets import cifar10

(x_train, y_train), (x_test, y_test) = cifar10.load_data()

Examinemos algunas imágenes para ver qué pinta tienen.

import matplotlib.pyplot as plt

for i in range(9):
    plt.subplot(330 + 1 + i)
    plt.imshow(x_train[i])
    
plt.show()

Vamos a transformar las etiquetas de números enteros a un vector binario, que es más conveniente a la hora de entrenar la red. Esto lo hacemos con la función to_categorical de keras. Esta forma de codificar las etiquetas se denomina one-hot.

from keras.utils import to_categorical

y_train = to_categorical(y_train)
y_test = to_categorical(y_test)

Creamos un conjunto de validación que nos va a servir para ir controlando la calidad del entrenamiento para evitar el overfitting.

import numpy as np

np.random.seed(42)

# barajamos el dataset
permutation = np.random.permutation(x_train.shape[0])
x_train = x_train[permutation]
y_train = y_train[permutation]

# obtenemos el conjunto de validación
num_val = 10000
x_val = x_train[:num_val]
x_train = x_train[num_val:]
y_val = y_train[:num_val]
y_train = y_train[num_val:]

Finalmente definimos el modelo de nuestra red neuronal.

from keras import layers
from keras import models

model = models.Sequential()
model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same', input_shape=(32, 32, 3)))
model.add(layers.MaxPooling2D(pool_size=(2, 2)))
model.add(layers.Conv2D(64, (3, 3), activation='relu', padding='same'))
model.add(layers.MaxPooling2D(pool_size=(2, 2)))
model.add(layers.Flatten())
model.add(layers.Dense(512, activation='relu'))
model.add(layers.Dense(10, activation='softmax'))

Esquemáticamente tenemos lo siguiente.

Detengámonos un poco en el modelo que acabamos de definir. A diferencia de la red que vimos en el anterior artículo sobre clasificación con deep learning, la entrada de la red convolucional es una imagen, en este caso de 32x32, y con tres capas, una por cada canal de color. Hay varias cosas que señalar. La primera es que el tamaño de los kernels suele ser de tamaño impar (3x3, 5x5, etc). Por otro lado, es habitual que el número de kernels que se aplica en cada capa convolucional vaya en aumento y tome valores (normalmente potencias de 2) cada vez más grandes según se profundiza en la red. En nuestro ejemplo hay pocas capas para lo que suele ser habitual (para permitir que este código se ejecute en equipo relativamente modestos). En una red más compleja se suele llegar a usar hasta 512, 1024, 2048 y hasta más número de filtros. También es habitual incluir una capa de max pooling cada pocas capas de convolución (en nuestro ejemplo, entre cada capa).
La última fase de la red es la de clasificación. Se utiliza una capa densamente conectada. Previamente se aplanan con flatten() las salidas de la fase convolucional, que recordemos que es una imagen en dos dimensiones, para convertirlo en un vector. En la última fase se usa la función de activación softmax de diez salidas (una por cada tipo de imagen a clasificar) que nos devuelve una distribución de probabilidad indicando cuál es la probabilidad de que la imagen pertenezca a cada una de las clases.
Ahora podemos compilar y entrenar el modelo.

from keras import optimizers

model.compile(loss='categorical_crossentropy', 
              optimizer=optimizers.RMSprop(lr=1e-4),
              metrics=['acc'])

train_log = model.fit(x_train, y_train, 
                      epochs=100, batch_size=128,
                     validation_data=(x_val, y_val))

Es buena práctica crear un par de gráficas para ver cómo se comporta el modelo durante el entrenamiento. Veamos la gráfica de ajuste con la función de coste y la de exactitud, tanto para el grupo de datos de entrenamiento como los de validación.

import matplotlib.pyplot as plt

acc = train_log.history['acc']
val_acc = train_log.history['val_acc']
loss = train_log.history['loss']
val_loss = train_log.history['val_loss']
epochs = range(1, len(acc) + 1)
plt.plot(epochs, acc, 'r', label='Training acc')
plt.plot(epochs, val_acc, 'b', label='Validation acc')
plt.title('Training and validation accuracy')
plt.legend()
plt.figure()
plt.plot(epochs, loss, 'r', label='Training loss')
plt.plot(epochs, val_loss, 'b', label='Validation loss')
plt.title('Training and validation loss')
plt.legend()
plt.show()

test_loss, test_accuracy = model.evaluate(x_test, y_test)
print(test_accuracy)

10000/10000 [==============================] - 0s 44us/step
0.6504999995231628

En la segunda gráfica observamos que hay un problema. La red neuronal comienza a tener problemas con el conjunto de validación y vemos como la función de coste crece en vez de minimizarse. Es un claro problema de overfitting. Vamos a recurrir a una técnica habitual en estos casos, que consiste en eliminar complejidad de la red. Al disminuir esta, conseguimos que se simplifique el modelo y que logre generalizar mejor. Una forma de reducir esta complejidad es eliminar algunas conexiones de forma aleatoria entre dos capas. Esto podemos hacerlo con una capa Dropout(). Como parámetro pasamos el porcentaje de conexiones que queremos anular (por ejemplo un valor de 0.25 eliminará el 25% de las conexiones).

model = models.Sequential()
model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same', input_shape=(32, 32, 3)))
model.add(layers.MaxPooling2D(pool_size=(2, 2)))
model.add(layers.Dropout(0.25))
model.add(layers.Conv2D(64, (3, 3), activation='relu', padding='same'))
model.add(layers.MaxPooling2D(pool_size=(2, 2)))
model.add(layers.Dropout(0.25))
model.add(layers.Flatten())
model.add(layers.Dense(512, activation='relu'))
model.add(layers.Dropout(0.5))
model.add(layers.Dense(10, activation='softmax'))

model.compile(loss='categorical_crossentropy', 
              optimizer=optimizers.RMSprop(lr=1e-4),
              metrics=['acc'])

train_log = model.fit(x_train, y_train, 
                      epochs=100, batch_size=128,
                     validation_data=(x_val, y_val))

Como vemos, ahora la cosa está mucho mejor. Y además llegamos al 71% de exactitud con el conjunto de test frente al 65% que teníamos antes.

Aun así, un 71% no es muy espectacular. ¿Se puede mejorar? Rotundamente sí. El modelo que hemos usado es demasiado sencillo para este tipo de tareas. Una red de este tipo tiene muchas más capas, pero claro, también necesitaremos una máquina y una GPU más potente. Por ejemplo, con el siguiente modelo, que sigue siendo relativamente sencillo, llegamos al 77% (eso sí, no lo intentéis en casa si no vais armados con una buena GPU). Con una red aún más compleja he llegado a superar una exactitud del 95%, pero usando técnicas que ya os contaré otro día.

model = models.Sequential()
model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same', input_shape=(32, 32, 3)))
model.add(layers.Conv2D(32, (3, 3)), activation='relu')
model.add(layers.MaxPooling2D(pool_size=(2, 2)))
model.add(layers.Dropout(0.25))
model.add(layers.Conv2D(64, (3, 3), activation='relu', padding='same')
model.add(layers.Conv2D(64, (3, 3), activation='relu')
model.add(layers.MaxPooling2D(pool_size=(2, 2)))
model.add(layers.Dropout(0.25))
model.add(layers.Flatten())
model.add(layers.Dense(512), activation='relu')
model.add(layers.Dropout(0.5))
model.add(layers.Dense(10), activation='softmax')

model.compile(loss='categorical_crossentropy', 
              optimizer=optimizers.RMSprop(lr=1e-4),
              metrics=['acc'])

train_log = model.fit(x_train, y_train, 
                      epochs=100, batch_size=128,
                     validation_data=(x_val, y_val))

En la práctica se utilizan modelos más complejos creados por investigadores y a las que se le asigna un nombre para diferenciarlos de otras arquitecturas. Algunos ejemplos son VGG16, AlexNet, ResNet o DenseNet, por nombrar algunas de las arquitecturas más usadas. Como ejemplo os muestro la arquitectura de RasNet, para que os hagáis una idea de la complejidad de los modelos con los que se trabaja habitualmente.

Os invito a experimentar añadiendo y quitando capas, modificando sus parámetros, añadiendo o quitando capas dropout, etc. Es un ejercicio interesante para ir viendo cómo mejora o empeora el resultado.

Clasificación multiclase con deep learning y keras

noreply@blogger.com (Alberto García Serrano) — Sun, 06 Oct 2019 00:28:00 +0000

Dentro de los diferentes tipos de problemas que trata de abordar el aprendizaje automático, uno de los más interesantes es la clasificación. Básicamente consiste en, dado un conjunto de elementos, asignarle a una etiqueta, donde cada etiqueta corresponde a una clase diferente. Por ejemplo, una aplicación interesante y bastante actual es clasificar a los visitantes de una web en dos categorías, los que, con alta probabilidad, pulsarán sobre un anuncio concreto y los que no. Este es un ejemplo de clasificación binaria, porque sólo hay dos categorías. Otro ejemplo puede ser el que ya vimos en el artículo sobre reconocimiento de caracteres manuscritos usando K-NN. En este caso hablábamos de un problema de clasificación multiclase en el que asignábamos a cada imagen una sola etiqueta de un conjunto de diez posibles (los dígitos del 0 al 9). Vimos que K-NN es un enfoque sencillo y efectivo para este tipo de problemas, pero hay otras alternativas. Vamos a enfrentarnos a un problema similar pero con un dataset y una técnica diferente.
Vamos a usar una red neuronal para resolver un problema de clasificación multiclase similar a MNIST (si no tienes muy claro qué es eso de una red neuronal puedes leer un artículo que escribí en mi antiguo blog sobre el perceptrón). Se trata del conjunto fashion MNIST, que creó la empresa Zalando con imágenes de prendas de vestir. Este dataset contiene 60.000 imágenes de prendas (28x28 en escala de grises), y 10.000 imágenes para test. Las categorías con su etiqueta son las siguientes:

0 T-shirt/top
1 Trouser
2 Pullover
3 Dress
4 Coat
5 Sandal
6 Shirt
7 Sneaker
8 Bag
9 Ankle boot

Se trata de, dada una imagen, asignarle su etiqueta numérica correspondiente. Para crear y entrenar la red neuronal vamos a recurrir a la librería keras, que nos ofrece una capa de abstracción sobre TensorFlow. Keras, además de una interfaz de alto nivel para el diseño de modelos de aprendizaje basado en redes de neuronas, nos permite utilizar una serie de datasets que están ya integrados. Las siguiente líneas nos permiten descargar el dataset fashion MNIST y ponerlo en memoria.

from keras.datasets import fashion_mnist
(x_train, y_train), (x_test, y_test) = fashion_mnist.load_data()

Si quiere visualizar una de las imágenes para ver qué pinta tienen puede hacerlo con matplotlib.

import matplotlib.pyplot as plt

plt.imshow(x_test[0], cmap="gray")
plt.show()

Son imágenes de 28x28 píxeles en escala de grises (con 256 niveles de gris). No son muy vistosas, pero tienen la ventaja de que son muy manejables a la hora de alimentar nuestra red neuronal. Como las imágenes son de 28x28 necesitamos 784 entradas (28x28=784) ya que a cada pixel le corresponde una entrada en la red. Así que lo primero que tenemos que hacer es adaptar las imágenes a un formato más manejable.
Primero transformamos cada imagen de una matriz de 28x28 a un vector de tamaño 784. Por otra parte, cada pixel es un entero en el rango 0-255. Lo vamos a normalizar para que cada pixel sea un número real que esté en el rango 0-1.
También vamos a transformar las etiquetas de números enteros a un vector binario (one-hot), que es más conveniente a la hora de entrenar la red. Esto lo hacemos con la función to_categorical() de keras, que nos devuelve una matriz en la que cada elemento es un vector de diez posiciones donde todos sus valores están a 0 menos uno de ellos. Por ejemplo, la etiqueta con el valor 2 quedaría como 0010000000 (no se trata del valor binario de la etiqueta, sino de asignar una posición concreta al 1 de cada etiqueta).

from keras.utils import to_categorical

x_train = x_train.reshape((60000, 784))
x_train = x_train.astype('float32') / 255
y_train = to_categorical(y_train)

x_test = x_test.reshape((10000, 784))
x_test = x_test.astype('float32') / 255
y_test = to_categorical(y_test)

Del conjunto de entrenamiento vamos a separar una porción de ellos como conjunto de validación. Este subconjunto lo vamos a utilizar para, durante el entrenamiento, ir viendo si el modelo sigue funcionando bien con datos que no están en el conjunto de entrenamiento. Esto lo hacemos para evitar el sobreentrenamiento u overfitting, que ocurre cuando un modelo termina "aprendiéndose de memoria" los datos de entrenamiento, con lo cual funciona muy bien con estos datos, pero no con otros datos de fuera del conjunto de entrenamiento. Es decir, el modelo no es capaz de generalizar.

import numpy as np

np.random.seed(42)

# barajamos el dataset
permutation = np.random.permutation(x_train.shape[0])
x_train = x_train[permutation]
y_train = y_train[permutation]

# obtenemos el conjunto de validación
num_val = 10000
x_val = x_train[:num_val]
x_train = x_train[num_val:]
y_val = y_train[:num_val]
y_train = y_train[num_val:]

Ahora llega el momento de diseñar una arquitectura para nuestra red neuronal. No hay reglas de oro a la hora de diseñar una arquitectura, hay mucho de arte, de intuición y de experiencia. Vamos a partir de una arquitectura con una capa de entrada, dos capas ocultas y una de salida.

from keras import models, layers

model = models.Sequential()
model.add(layers.Dense(256, activation='relu', input_shape=(784,))) # 28x28=784
model.add(layers.Dense(128, activation='relu'))
model.add(layers.Dense(10, activation='softmax'))

Esquemáticamente nuestra red tiene la siguiente arquitectura.

Sin entrar en mucho detalle, nuestra red tiene 784 entradas (28x28 píxeles), que se conectan a una primera capa oculta de 256 neuronas conectadas densamente (todas las entradas se conectan a todas las neuronas de esta capa). Como función de activación se usa la función rectificador o relu, cuya función es

\[f(x) = \max(0, x)\]

Hay diversas funciones de activación posibles, aunque en los últimos años la función relu ha ido ganando popularidad frente a la más clásica función de activación sigmoide.
La segunda capa oculta, también densamente conectada y usando la función de activación relu, tiene 128 neuronas. Finalmente, la capa de salida tiene 10 neuronas, una por cada clase que queremos clasificar. En esta capa usamos la función de activación softmax, que nos devuelve una distribución de probabilidad entre las salidas, de forma que cada una de estas salidas nos indica la probabilidad de que la imagen de entrada pertenezca a dicha clase. Es decir, que la suma de todas las salidas sumarán 1.

Las siguientes líneas compilan el modelo y lo parametrizan. En este caso, para optimizar la función de costes, que es la que queremos minimizar, usamos el optimizador rmsprop, que es un tipo de optimización basada en el descenso de gradiente que ya os presenté en otro artículo anterior (con algunas diferencias). como función de coste (loss function) usaremos categorical_crossentropy, que funciona especialmente bien en problemas de clasificación donde cada ejemplo pertenece exclusivamente a una sola categoría.

model.compile(optimizer = 'rmsprop',
             loss = 'categorical_crossentropy',
             metrics = ['accuracy'])

Ya sólo nos resta entrenar el modelo. Vamos a entrenarlo usando lotes (batch) de 128 ejemplos durante 30 iteraciones (epochs). Además le indicamos cuál es nuestro conjunto de validación.

train_log = model.fit(x_train, y_train, 
                      epochs=30, batch_size=128,
                     validation_data=(x_val, y_val))

Finalizado el entrenamiento, podemos probar el modelo contra el conjunto de test y obtener la exactitud de nuestro modelo resultante.

test_loss, test_accuracy = model.evaluate(x_test, y_test)
print(test_accuracy)

10000/10000 [==============================] - 0s 12us/step
0.8788999915122986

Es decir, hemos conseguido una exactitud del 87% asignando las imágenes del conjunto de test a su clase correspondiente. No está mal. Sin embargo, para ir ajustando el modelo necesitamos alguna herramienta que nos permita medir qué tan bien o mal se comporta el modelo. Durante el entrenamiento se genera un histórico con los valores que ha ido tomando la función de coste y la precisión tanto del conjunto de entrenamiento como del de validación. Podemos obtener una representación gráfica.

# grafica con la función de coste

loss = train_log.history['loss']
val_loss = train_log.history['val_loss']
epochs = range(1, len(loss) + 1)
plt.plot(epochs, loss, 'b', label='Training loss')
plt.plot(epochs, val_loss, 'r', label='Validation loss')
plt.title('Training and validation loss')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()
plt.show()

También nos puede ser útil ver cómo ha ido evolucionando la precisión durante el entrenamiento.

# grafica de la exactitud

plt.clf()
acc = train_log.history['accuracy']
val_acc = train_log.history['val_accuracy']
epochs = range(1, len(loss) + 1)
plt.plot(epochs, acc, 'b', label='Training acc')
plt.plot(epochs, val_acc, 'r', label='Validation acc')
plt.title('Training and validation accuracy')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()
plt.show()

En la primera gráfica, la de la función de costes, vemos que a partir de la octava iteración (epoch), la función de coste sobre el conjunto de validación comienza a crecer (empeora), mientras que en el caso del conjunto de entrenamiento sigue decreciendo (mejorando). Eso es señal inequívoca de que nuestro modelo esta sobreentrenado (overfitting). Es decir, que predice muy bien en el conjunto de entrenamiento, pero mal en datos externos a este. Y eso es malo. Podemos hacer varias cosas para evitarlo, pero las dos más sencillas son, quitar neuronas para simplificar el modelo y/o reducir el número de iteraciones. Vamos a probar con la segunda opción (hay otras un poco más complejas y efectivas, pero ya os las cuento otro día). Vamos a ver qué pasa si lo dejamos en 8 iteraciones (que es donde empieza a empeorar la función de coste del conjunto de validación).

10000/10000 [==============================] - 0s 13us/step
0.8809999823570251

Tenemos un 88% de exactitud, lo que de por sí ya es una pequeña mejora, pero además nos hemos ahorrado 22 iteraciones, lo que es una buena cantidad de tiempo de cómputo. Si no tienes una buena GPU pueden llegar a ser algunos minutos en un ordenador no muy potente. La gráfica de la evolución de la función de coste ahora es esta.

y la de exactitud es esta.

Ahora que tenemos un modelo entrenado, podemos hacer predicciones. Por ejemplo sobre el conjunto de test.

predictions = model.predict(x_test)
np.argmax(predictions[0]) # clase más probable

Esto nos devuelve el valor 9 (Ankle boot), que es precisamente la categoría a la que pertenece la primera imagen del conjunto de test, que ya os mostré al principio del artículo.

Se puede seguir jugando con el modelo añadiendo o quitando neuronas, o capas completas. También cambiando el tipo de función de activación, etc. Aquí os dejo el notebook de jupyter por si queréis trastear un poco con la red neuronal que acabamos de crear.

Display de 7 segmentos con 74LS247

noreply@blogger.com (Alberto García Serrano) — Wed, 12 Jun 2019 16:18:49 +0000

Ya hemos hablado del control de un display de siete segmentos en un artículo anterior (implementando la lógica de control con una FPGA).

Las bases son las mismas, y si buscas, encontraras por la red miles de tutoriales sobre como gestionar estos display con Arduino.

La verdad es que no es difícil, pero le veo un par de inconvenientes. Por un lado, si te fijas en la figura de arriba, para gestionar un display son necesarias ocho conexiones digitales entre el Arduino y el display. Y eso sólo para un display, si necesitas más es probable que no tengamos suficientes pines para gestionarlos. Otra pega es que tendremos que programar la lógica del display en el propio Arduino. Es decir, que si queremos mostrar un número, el programa tendrá que evaluar qué LEDs han de encenderse y cuales han de apagarse, lo que requiere tiempo de proceso y un consumo de energía relativamente alto para ser una gestión secundaria (querremos que el microcontrolador se ocupe de tareas menos mundanas que decidir qué LED encender para mostrar, por ejemplo, la temperatura de la sala). Para la gestión de estos displays hay circuitos integrados especializados que son bastante más eficientes gestionando los LED de lo que lo pueda ser un microcontrolador. Son los llamados decodificadores de siete segmentos. Por ejemplo, el 74LS247 (en este caso he usado un sn74LS247n). Estos decodificadores aceptan un valor en formato BCD y activan las salidas correspondientes en el display para mostrar dicho valor.

Lo primero a tener en cuenta es el tipo de display que vamos a usar. Hay dos tipos, de ánodo común y de cátodo común. Como puedes observar en el esquema que más abajo, la diferencia es que en el primero todos los LED comparten la conexión a Vcc mientras que en el segundo, comparten la conexión a tierra. Esto es importante porque dependiendo del tipo de display, tendremos que buscar un integrado que sea capaz de gestionarlo. En este caso he usado un display con ánodo común.

El 74LS247 viene en formato DIP con 16 patillas. Además de la patilla GND y Vcc (5V), disponemos de cuatro pines de entrada (A, B, C y D) que usaremos para indicar el valor que queremos mostrar en el display, en formato BCD. Las salidas (desde a hasta g) se conectan a las patillas correspondientes en el display (tienen el mismo nombre como puedes observar en los esquemas de más arriba). Además disponemos de tres entradas de control: LT (lamp test), RBI (RB input) y RBO (RB outout). Entre otras cosas, son útiles para comprobar el buen estado del display: cuando RBO está a nivel alto y LT está a nivel bajo, se activan todas las salidas, encendiendo todos los LED del display, lo que nos permite comprobar que todos los LED funcionan correctamente. Tienen más funcionalidades, pero aquí no las vamos a usar por el momento.

Así pues, hechas las presentaciones, vamos a ver cómo se conectan el display y el 74LS247.

Las entradas A, B, C y D las he conectado al 74LS247 con unos interruptores DIP mediante unas resistencias pull-down (200Ω). Las salidas de la pastilla las he conectado directamente al display mediante resistencias de 200Ω. Sólo tenemos que manipular los cuatro primeros interruptores para poner el valor binario del dígito que queremos mostrar en el display. Os dejo un vídeo para que veáis como funciona.

Juegos con Arduino y LCD 5110

noreply@blogger.com (Alberto García Serrano) — Tue, 09 Apr 2019 20:24:18 +0000

¿Recuerdas el mítico Nokia 5110? El teléfono, obviamente, ya no se vende, pero su pantalla sigue siendo muy valorada por los makers y aficionados a la electrónica. La razón es que, gracias a su interfaz SPI es muy sencilla de manejar desde un microcontrolador. Como hace tiempo que tengo uno dando vueltas por mis cajones me he animado a probarlo, así que os traigo un pequeño montaje realizado con un Arduino Nano, un potenciómetro y, por supuesto, un LCD 5110. Y que mejor que un juego para probar este pequeño LCD.

Las siguientes imágenes muestran la pantalla LCD por delante y también por la parte trasera donde se pueden ver las patillas de conexión. Hay que tener cuidado, ya que he encontrado este mismo LCD con diferente patillaje. El que yo he usado es como éste.

La pantalla es monocroma con una resolución de 84x48 píxeles. Utiliza el controlador PCD8544 y dispone de LEDs para iluminar la pantalla. Se alimenta con una tensión de 3.3V, y sus los pines digitales también funcionan a 3.3V, por lo que hay que adaptar las señales digitales del Arduino, que funcionan a 5V. Aunque no es lo más óptimo, con unas cuantas resistencias de 10KΩ podemos apañarnos. El esquema completo queda así.

Las resistencias son de 10KΩ, excepto la primera por la izquierda que conecta el pin 7 del Arduino, que es de 1KΩ, y la última de la derecha, que conecta el LCD con el pin 9 del Arduino, que es de 330Ω. Finalmente las conexiones entre el 5110 y el Arduino quedan de la siguiente manera.

       1-VCC       ----------------  3.3V
       2-GND       ----------------  GND
       3-SCE       ----------------  7
       4-RST       ----------------  6
       5-D/C       ----------------  5
       6-DN(MOSI)  ----------------  11
       7-SCLK      ----------------  13
       8-LED       ----------------  9

El código completo del juego es el siguiente.

#include <Arduino.h>
#include <SPI.h>

/* Comandos PCD8544 */
#define LCD_COMMAND  0
#define LCD_DATA     1

/* LCD 84x48 */
#define LCD_WIDTH   84 
#define LCD_HEIGHT  48 
#define WHITE       0  
#define BLACK       1  

/* Conexionado de pines */
const int scePin = 7;   // SCE - Chip select, pin 3 en LCD.
const int rstPin = 6;   // RST - Reset, pin 4 en LCD.
const int dcPin = 5;    // DC - Data/Command, pin 5 en LCD.
const int sdinPin = 11;  // DN(MOSI) - Serial data, pin 6 en LCD.
const int sclkPin = 13;  // SCLK - Serial clock, pin 7 en LCD.
const int blPin = 9;    // LED - Backlight LED, pin 8 en LCD.
const int analog_in = A0; // entrada analogica potenciometro

/* frame buffer para el display */
byte displayMap[LCD_WIDTH * LCD_HEIGHT / 8];

/* Envio de datos al LCD
 *  El primer parámetro indica si estamos enviando un comando o un dato
 *  El segundo parámetro es el comando o el dato
 */
void LCDWrite(byte data_or_command, byte data)
{
  digitalWrite(dcPin, data_or_command);
  digitalWrite(scePin, LOW);
  SPI.transfer(data); 
  digitalWrite(scePin, HIGH);
}

/* establece fila columna
 *  x - columna
 *  y - fila
 */
void gotoXY(int x, int y)
{
  LCDWrite(0, 0x80 | x);
  LCDWrite(0, 0x40 | y);
}

/* Cambia el estado de un pixel en el framebuffer
 *  x, y - posicion
 *  bw - BLACK o WHITE
 */
void setPixel(int x, int y, boolean bw)
{
  if ((x >= 0) && (x < LCD_WIDTH) && (y >= 0) && (y < LCD_HEIGHT))
  {
    byte shift = y % 8;
    if (bw) 
      displayMap[x + (y/8)*LCD_WIDTH] |= 1<<shift;
    else   
      displayMap[x + (y/8)*LCD_WIDTH] &= ~(1<<shift);
  }
}

/* Dibuja una línea de x0,y0 a x1, y1
 * bw - BLACK o WHITE
 */
void setLine(int x0, int y0, int x1, int y1, boolean bw)
{
  int dy = y1 - y0; 
  int dx = x1 - x0; 
  int stepx, stepy;
  if (dy < 0)
  {
    dy = -dy;
    stepy = -1;
  }
  else
    stepy = 1;
    if (dx < 0)
  {
    dx = -dx;
    stepx = -1;
  }
  else
    stepx = 1;
    dy <<= 1; // dy = 2*dy
  	dx <<= 1; // dx = 2*dx
  	setPixel(x0, y0, bw); 
  	if (dx > dy)
  		{
    	int fraction = dy - (dx >> 1);
    	while (x0 != x1)
    	{
      		if (fraction >= 0)
      		{
        	y0 += stepy;
        	fraction -= dx;
      		}
      	x0 += stepx;
      	fraction += dy;
      	setPixel(x0, y0, bw);
    	}
  }
  else
  {
    int fraction = dx - (dy >> 1);
    while (y0 != y1)
    {
      if (fraction >= 0)
      {
        x0 += stepx;
        fraction -= dy;
      }
      y0 += stepy;
      fraction += dx;
      setPixel(x0, y0, bw);
    }
  }
}

/* dibuja un rectángulo
 *  x0, y0 - esquina superior izquierda
 *  x1, y1 - esquina inferior derecha
 *  fill - 1: relleno, 0: sólo perímetro
 *  bw - BLACK o WHITE
 */
void setRect(int x0, int y0, int x1, int y1, boolean fill, boolean bw)
{
  if (fill == 1)
  {
    int xDiff;
    if(x0 > x1)
      xDiff = x0 - x1;
    else
      xDiff = x1 - x0;
      while(xDiff > 0)
    	{
      	setLine(x0, y0, x0, y1, bw);ç
      	if(x0 > x1)
        	x0--;
      	else
        	x0++;
            xDiff--;
    }
  }
  else
  {
    setLine(x0, y0, x1, y0, bw);
    setLine(x0, y1, x1, y1, bw);
    setLine(x0, y0, x0, y1, bw);
    setLine(x1, y0, x1, y1, bw);
  }
}

/* 
 *  vuelca el contenido del framebuffer a la memoria del LCD.
 */
void updateDisplay()
{
  gotoXY(0, 0);
  for (int i=0; i < (LCD_WIDTH * LCD_HEIGHT / 8); i++)
  {
    LCDWrite(LCD_DATA, displayMap[i]);
  }
}

/* establece el contraste del LCD
 *  normalmente entre 40 y 60
 */
void setContrast(byte contrast)
{
  LCDWrite(LCD_COMMAND, 0x21); 
  LCDWrite(LCD_COMMAND, 0x80 | contrast); 
  LCDWrite(LCD_COMMAND, 0x20); 
}

//Envía comandos mágicos al PCD8544
/* 
 *  Inicializar el LCD 5110
 */
void lcdBegin(void)
{
  pinMode(scePin, OUTPUT);
  pinMode(rstPin, OUTPUT);
  pinMode(dcPin, OUTPUT);
  pinMode(sdinPin, OUTPUT);
  pinMode(sclkPin, OUTPUT);
  pinMode(blPin, OUTPUT);
  analogWrite(blPin, 255);
  SPI.begin();
  SPI.setDataMode(SPI_MODE0);
  SPI.setBitOrder(MSBFIRST);
  digitalWrite(rstPin, LOW);
  digitalWrite(rstPin, HIGH);
  LCDWrite(LCD_COMMAND, 0x21);
  LCDWrite(LCD_COMMAND, 0xB0);
  LCDWrite(LCD_COMMAND, 0x04);
  LCDWrite(LCD_COMMAND, 0x14);
  LCDWrite(LCD_COMMAND, 0x20);
  LCDWrite(LCD_COMMAND, 0x0C);
}

/* borrar pantalla
 *  bw - BLACK o WHITE
 */
void clearDisplay(boolean bw)
{
  for (int i=0; i<(LCD_WIDTH * LCD_HEIGHT / 8); i++)
  {
    if (bw)
      displayMap[i] = 0xFF;
    else
      displayMap[i] = 0;
  }
}

int val, pos, bloq_size = 5;
int way_lenght=5;
int way_pos = (LCD_WIDTH/2)-(way_lenght*bloq_size); // posición inicial
int way_dir;
int way[10] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}; // tramos de carretera
int i;

void setup() {
  lcdBegin();
  setContrast(40);
}

void loop() {
  val = analogRead(analog_in); // leer potenciometro
  pos = map(val, 0, 780, LCD_WIDTH-bloq_size, 0);
  way_dir = random(-1,2);
  way_pos = way_pos + (way_dir*bloq_size); // nueva posioción de la carretera
  /* evitar salirnos de la pantalla */
  if (way_pos<0) way_pos = 0;
  if (way_pos>LCD_WIDTH-(bloq_size*way_lenght)) way_pos = LCD_WIDTH-(bloq_size*way_lenght); 
  /* scroll de la carretera */
  for (i=9; i>=1; i--) {
    way[i] = way[i-1];
  }
  way[0] = way_pos;
  
  /* dibujar pantalla */
  clearDisplay(WHITE);
  
  /* dibujar coche */
  setRect(pos, bloq_size*8, pos+bloq_size, bloq_size*8-10+bloq_size, 1, BLACK);
  
  /* dibujar carretera */
  for (i=0; i<=9; i++) {
    if (way[i] > 0) {
      setRect(way[i], i*bloq_size, way[i]+bloq_size, (i*bloq_size)+bloq_size, 0, BLACK);
      setRect(way[i]+(way_lenght*bloq_size), i*bloq_size, way[i]+(way_lenght*bloq_size)+bloq_size, (i*bloq_size)+bloq_size, 0, BLACK);
    }
  }
  
  /* comprobar colisión */
  if ((pos<way[7]+bloq_size) || (pos+bloq_size>way[7]+bloq_size*way_lenght) && way[7]>0) {
    // hay colision
    setRect(pos, bloq_size*8, pos+bloq_size, bloq_size*8-10+bloq_size, 1, WHITE);
    setRect(pos, bloq_size*8, pos+bloq_size, bloq_size*8-10+bloq_size, 0, BLACK);  
    updateDisplay();
    delay(5000);
    // reiniciar el juego
    for (i=0; i<=9; i++) {
      way[i] = -1;
    }
  }
  
  updateDisplay();
  delay(500);
}

Para las funciones de control del LCD me he basado en un código de Sparkfun. Las operaciones de dibujado se realizan sobre una zona de memoria del Arduino (frame buffer). En este caso esa zona de memoria se almacena en la variable displayMap, que es un array de representa el contenido de la pantalla. Cuando deseamos volcar el contenido del frame buffer en la pantalla física, usamos la función LCDWrite. Esta función hace uso del bus SPI para enviar comandos y el contenido del frame buffer al 5110. El resto de funciones son bastante auto explicativas y tienen un pequeño encabezado en su cabecera.

Gestión de interrupciones con Arduino

noreply@blogger.com (Alberto García Serrano) — Sun, 31 Mar 2019 20:50:40 +0000

Cuando se diseña un programa, lo habitual es pensarlo como un flujo de instrucciones que se ejecutan una detrás de otra siguiendo una ejecución secuencial. En el contexto de los microprocesadores y los microcontroladores, si estuviéramos limitados sólo a este modelo secuencial, el desarrollo de aplicaciones sería más complicado. Os pongo un ejemplo muy simplificado: si el sistema que estamos diseñando tiene que atender una serie de actuadores (por ejemplo un botón), tendremos que andar sondeándolos periódicamente para ver si se pulsó el botón. ¿Pero qué ocurre si se pulsa un botón mientras el programa anda ocupado en otra tarea y justo en ese momento no se está sondeando la entrada del actuador? En este caso el sistema no será capaz de detectar la pulsación. La solución en este caso es el uso de interrupciones. Cuándo se dispara una interrupción, el procesador detiene la ejecución del programa que está corriendo, almacena el estado de los registros (incluido el contador de programa), y ejecuta un subprograma, llamado rutina de servicio de interrupción, que atiende la interrupción (en este caso la gestión de la pulsación del botón). Una vez finalizada la rutina, el microprocesador restaura los registros y continúa la ejecución del programa principal por el mismo lugar en el que fue interrumpido. Vamos a ver cómo podemos crear una rutina de interrupción en Arduino y hacer que se ejecute tras la pulsación del botón.

La siguiente figura muestra el esquema del circuito que vamos a usar. Como se puede observar, el circuito no puede ser más simple. Se trata de un botón conectado al puerto del pin 2 del Arduino y un LED conectado al puerto del pin 13 de Arduino.

El objetivo es hacer un programa en Arduino que encienda el LED cuando se pulse el botón. El siguiente código realiza esta acción sin usar interrupciones, es decir, a la manera clásica.

const int pin = 2;
const int pout = 13;
int estadoBoton;
 
void setup() {
  pinMode(pin, INPUT);
  pinMode(pout, OUTPUT);
}
void loop() {
  estadoBoton = digitalRead(pin);
  if (estadoBoton == HIGH) {
    digitalWrite(pout, HIGH);
  } else {
    digitalWrite(pout, LOW);
  }
}

El código es bastante simple: se configura el puerto 2 como entrada (aquí se conecta el botón) y el 13 como salida (aquí se conecta el LED). Dentro del bucle loop leemos constantemente el valor del puerto al que está conectado el botón y si está pulsado encendemos el LED. En caso contrario se apaga. Vamos a hacer algunos cambios para que el programa realice la misma función, pero usando interrupciones.

const int pin = 2;
const int pout = 13;
int estadoBoton;
 
void setup() {
  pinMode(pin, INPUT);
  pinMode(pout, OUTPUT);
  attachInterrupt(digitalPinToInterrupt(pin), enciendeLED, CHANGE);
}
 
void loop() {
  // nada que hacer aquí. circulen.
}
 
void enciendeLED() {
  estadoBoton = digitalRead(pin);
  if (estadoBoton == HIGH) {
    digitalWrite(pout, HIGH);
  } else {
    digitalWrite(pout, LOW);
  }
}

Lo primero que nos llama la atención es que en el bucle principal (función loop) no hay código. Centrémonos en la siguiente línea de código.

attachInterrupt(digitalPinToInterrupt(pin), enciendeLED, CHANGE);

La función attachInterrupt() es la encargada de asociar un evento de interrupción con una rutina de servicio de interrupción. El primer parámetro es la interrupción que vamos a tratar. Con digitalPinToInterrupt convertimos una entrada (un pin) en una interrupción. En este caso estamos creando una interrupción asociada a la pulsación del botón que está en el pin 2. El segundo parámetro permite indicar cuál es la rutina de servicio de interrupción, o lo que es lo mismo, la función que se invoca cuando se dispara dicha interrupción. Por último, podemos indicar cuándo se activa la interrupción. Las opciones son:

CHANGE: La interrupción se activa cuando cambia el estado del pin.
LOW: La interrupción se activa cuando el pin está a nivel bajo.
RISING: La interrupción se activa cuando el pin pasa de nivel bajo a alto.
FALLING: La interrupción se activa cuando el pin pasa de nivel alto a bajo.
HIGH: La interrupción se activa cuando el pin está en nivel alto. Sólo funciona en Arduino Due, Zero y MKR1000.

Si pruebas este código observarás que el LED se enciende tal y como esperábamos. Sin embargo, si pruebas a conectar el botón a otro pin, digamos el A0, comprobarás que, aunque modifiques el código para usar este pin, el código no funciona ¿por qué?

Básicamente porque aunque el microcontrolador soporta interrupciones en todos sus pines de entrada y salida, la plataforma Arduino sólo lo implementa en algunos de ellos. En concreto en los siguientes pines según el modelo de Arduino:

Uno, Nano y Mini : pines 2 y 3
Mega (todos los modelos): 2, 3, 18, 19, 20, 21
Micro, Leonardo: 0, 1, 2, 3, 7
Zero: todos los pines excepto el 4
MKR1000: 0, 1, 4, 5, 6, 7, 8, 9, A1, A2
Due: todos

¿Quiere eso decir que sólo podemos gestionar dos botones en un Arduino Uno? Para nada. Tenemos dos opciones. La primera es acceder a los registros del microcontrolador a bajo nivel, lo cuál puede llegar a ser algo enrevesado según lo que queramos hacer. La segunda opción es usar una librería externa. Hay varias, en este caso vamos a usar PinChangeInterrupt, que podemos instalar fácilmente desde el propio entorno de Arduino. Para ello accedemos a Programa > Incluir librería > Gestionar librería.

Buscamos pinchangeinterrupt y pulsamos el botón instalar.

Una vez instalada ya podemos usar interrupciones en cualquier pin usando la librería. El siguiente código hace los mismo que en anterior pero usando el pin A0.

#include <PinChangeInterrupt.h>
const int pin = A0;
const int pout = 13;
int estadoBoton;
 
void setup() {
  pinMode(pin, INPUT);
  pinMode(pout, OUTPUT);
  attachPCINT(digitalPinToPCINT(pin), enciendeLED, CHANGE);
}
 
void loop() {
  // nada que hacer aquí. circulen.
}
 
void enciendeLED() {
  estadoBoton = digitalRead(pin);
  if (estadoBoton == HIGH) {
    digitalWrite(pout, HIGH);
  } else {
    digitalWrite(pout, LOW);
  }
}

El código es casi idéntico salvo por la primera línea, que incluye las cabeceras de la librería. Y también la línea siguiente, que configura la interrupción.

attachPCINT(digitalPinToPCINT(pin), enciendeLED, CHANGE);

El formato de la función attachPCINT es igual a la de la función attachInterrupt que ya hemos visto más arriba. En el vídeo siguiente podemos ver el código anterior funcionando usando el pin A0.

Regresión lineal y descenso de gradiente con Python

noreply@blogger.com (Alberto García Serrano) — Tue, 26 Mar 2019 21:14:00 +0000

En machine learning, el objetivo principal es encontrar un modelo que explique el comportamiento de un sistema (en el amplio sentido de la palabra). A partir de unos datos de entrenamiento, un sistema de aprendizaje automático ha de ser capaz de inferir un modelo capaz de explicar, al menos en su mayoría, los efectos observados. Pero también aplicar ese aprendizaje. Por ejemplo, un sistema de machine learning muy lucrativo para las empresas anunciantes es aquél que dado un perfil de usuario (datos de entrada A), sea capaz de predecir si pinchará o no (salida B) sobre un anuncio publicitario de, por ejemplo, comida para gatos. No es sencillo crear un modelo capaz de predecir el comportamiento del usuario (o sí), pero en todo caso, existen diferentes técnicas que nos permiten abordar el problema. En el caso del ejemplo que acabamos de ver, el modelo debería ser capaz de clasificar a los usuarios en dos clases diferentes, los que pulsarán y los que no pulsarán el anuncio de comida de gatos. Otro tipo de problema es aquél en el que se pretende crear un modelo capaz de predecir un valor o una cantidad. Por ejemplo, si conozco los metros cuadrados, el número de habitaciones y el número de baños de una vivienda en venta, ¿seré capaz de predecir por qué precio se venderá? Uno podría hacerse una idea observando el precio de venta de alguna casa igual o parecida que se haya vendido antes. Si no tenemos ese dato, podremos extrapolar el precio comparando el precio de venta de otras casas y comparando esa casa con la nuestra. Intuitivamente sabemos que a más metros, más precio ¿pero cuanto más? Es probable que entre los metros de la casa y el precio haya una relación lineal, así que aprovechándonos de ello, es relativamente sencillo construir un modelo (lineal, claro) que nos ayude a predecir el precio de una casa. Así que manos a la obra.
Vamos a basarnos en una técnica estadística llamada regresión lineal, que conceptualmente es bastante sencilla, pero es muy usada en este tipo de problemas. Antes de seguir vamos a generar un conjunto de datos ficticios con los que empezar a trabajar.

import random
import matplotlib.pyplot as plt

def gen_data(n, bias, varianza):
    x = []
    y = []
    for i in range(0, n):
        x.append(i)
        y.append((i + bias) + random.uniform(0, 1) * varianza)
    return x, y

x,y = gen_data(100, 25, 50)
plt.scatter(x, y)

Este código genera una serie de datos que sigue básicamente una recta alrededor de la cual se crean los puntos de forma aleatoria usando una función de distribución uniforme con una varianza dada (50 en este caso). Echemos un vistazo a como quedan los datos.

Por seguir con nuestro ejemplo de las casas puede imaginar que en el eje de abscisas se representan los metros cuadrados que tiene la casa y en el las ordenadas el precio en miles de euros por el que se vendió la casa.
La idea es, a partir de todos estos datos, crear un modelo que sea capaz de predecir el precio de cualquier casa dados sus metros cuadrados. En este caso es muy obvio que un buen predictor es una línea recta (por eso lo de lineal) que se ajuste lo más posible a la nube de puntos. Así pues, nuestro problema ahora es encontrar esa recta, la de mejor ajuste.
La ecuación de una recta en el plano es la siguiente.
\[y=a+bx\]
Donde a se corresponde con el punto de corte de la recta con el eje de ordenadas y b es la pendiente de la recta, es decir, lo rápido que crece o decrece. La variable independiente x representaría en este caso los metros cuadrados de la casa y la variable y el precio de venta. Por lo tanto, si somos capaces de buscar unos buenos valores para a y b, tendremos un modelo lineal para predecir el precio.
Podemos hacerlo de diferentes formas. Típicamente, la recta de regresión se ha calculado a partir de métodos estadísticos, pero (siempre hay un pero) si el número de casas que tenemos es muy grande (no lo llaman big data por gusto) o en vez de una sola variable independiente tenemos varias (metros, número de habitaciones, número de baños, metros de jardín, etc.) puede llegar a ser tedioso cuando no inabordable. Así pues vamos a usar un método para ajustar la recta llamado descenso de gradiente (gradient descent en la lengua de Shakespeare). Pero vayamos por partes. Lo primero que necesitamos es una forma de medir lo bien o mal que se ajusta una recta a la nube de puntos.
A cada par de valores que probemos para a y b lo vamos a llamar hipótesis, y seguidamente comprobaremos cómo se ajusta esa hipótesis calculando un factor de error, que es la suma de lo que, de media, se equivoca dicha hipótesis (recta) al predecir cada uno de los puntos de la nube. Este error se calcula con.
\[e=\frac{1}{2m}\sum_{i=1}^{m}(h(x^{i})-y^{i})^{2}\]
Donde $h(x^{i})$ es la hipótesis que queremos probar. Por ejemplo, si asignamos los valores $a=1$ y $b=2$, tenemos que $h(x^{i})=1+2x^{i}$ siendo $x^{i}$ cada uno de los valores de x para los metros cuadrados de las casas e $y^{i}$ los precios de dicha casa. Como se puede observar, lo que estamos haciendo es restar al precio de la casa que predice nuestra hipótesis (recta) el valor de venta real. Esa diferencia es el error, que se va sumando para todas las casas de las que disponemos de datos. El resultado de la resta (error) se eleva al cuadrado para evitar que unos errores se compensen con otros, de forma que el error siempre sea positivo. Resumiendo, lo que buscamos son los valores de a y b (nuestra hipótesis) que minimice el error devuelto por esta función. Ahí es donde entra en juego nuestro algoritmo del descenso de gradiente.
El algoritmo de descenso de gradiente hace uso de las derivadas (sí, eso que en el instituto pensaste que jamás te serviría para nada) y más concretamente de las derivadas parciales (sí, eso que en la universidad pensaste que jamás te serviría para nada). La función de error que acabamos de ver es una función convexa (no voy a demostrarlo aquí), lo que es muy interesante, ya que sólo tiene un valor mínimo que es global a toda la función. Puedes imaginarla con una forma similar a esta.

Imagina que el punto A es nuestra hipótesis inicial, y como queremos encontrar el valor de a y b que minimicen la función de error, lo que se busca es llegar al punto B. Conceptualmente, la derivada nos indica la velocidad a la que crece o decrece una función en uno de sus puntos. Por lo tanto, lo que debemos hacer es ir moviendo el punto A en la dirección de mayor pendiente para llegar a B. Como hay dos variables en juego, a y b, debemos calcular las derivadas parciales de la función de error con respecto a las dos variables.
\[a=a-\alpha \frac{\partial }{\partial a}e(a,b)\]
\[b=b-\alpha \frac{\partial }{\partial b}e(a,b)\]
El parámetro $\alpha$ es lo que se llama learning rate, y permite controlar la velocidad de descenso hacia el punto mínimo de la función. Un valor bajo hará que la convergencia sea lenta, mientras que un valor alto puede hacer que nos pasemos de frenada y acabemos dando santos de un lado a otro de la función sin alcanzar nunca el mínimo. Valores típicos son 0.1, 0.001, etc.
Si resolvemos las derivadas parciales obtenemos lo siguiente.
\[a=a-\alpha \frac{1}{m}\sum_{i=1}^{m}(h(x^{i})-y^{i})\]
\[b=b-\alpha \frac{1}{m}\sum_{i=1}^{m}((h(x^{i})-y^{i})x^{i})\]
En teoría, cada vez que recalculamos estas dos ecuaciones, estaremos más y más cerca del mínimo de la función de error, y por lo tanto, nuestra recta se ajustará mejor a la nube de puntos.
He hecho una implementación en Python del algoritmo tal y como os lo he contado, para que sea más claro. Sin embargo, este tipo de implementación no es muy eficiente. Normalmente se usa cálculo matricial, pero me ha parecido más pedagógico aplicar directamente las ecuaciones que acabamos de ver. Empecemos con la función para calcular el error (también llamada función de coste o loss function).

def coste(x, y, a, b):
    m = len(x)
    error = 0.0
    for i in range(m):
        hipotesis = a+b*x[i]
        error +=  (y[i] - hipotesis) ** 2
    return error / (2*m)

La función obtiene como parámetros una lista de valores x (metros cuadrados de las casas) e y (precio de las casas), así como los valores a y b para la hipótesis. A partir de ahí simplemente se calcula la ecuación de error que hemos visto más arriba y se devuelve el resultado.

def descenso_gradiente(x, y, a, b, alpha, epochs):
    m = len(x)
    hist_coste = []
    for ep in range(epochs):
        b_deriv = 0
        a_deriv = 0
        for i in range(m):
            hipotesis = a+b*x[i]
            a_deriv += hipotesis - y[i]
            b_deriv += (hipotesis - y[i]) * x[i]
            hist_coste.append(coste(x, y, a, b))
        a -= (a_deriv / m) * alpha
        b -= (b_deriv / m) * alpha
        
    return a, b, hist_coste

La función de descenso de gradiente, además de los parámetros de la función de coste, recibe alpha (learning rate) y el número de iteraciones, o lo que es lo mismo, el número de veces que vamos a recalcular y aplicar la pendiente de a y b para ir acercándonos al mínimo de la función de error. La función devuelve los valores óptimos de a y b que mejor ajustan la recta a la nube de puntos (o al menos un valor cercano) y una lista con el histórico de la función de coste, para poder estudiar cómo ha ido decreciendo el error. Vamos a probar nuestro algoritmo.

a=1
b=1
alpha = 0.0001
iters = 100000
a,b, hist_coste = descenso_gradiente(x, y, a, b, alpha, iters)
a,b

(43.41993518666135, 1.08466250771225)

Por lo tanto, nuestra recta será $y=43.41+1.08x$, que como vemos en la siguiente figura, ajusta bastante bien.

Si representamos gráficamente el error que hemos ido cometiendo en cada iteración podemos ver cómo va decreciendo (en este caso he omitido las primeras iteraciones porque al principio decrece muy bruscamente y la escala de la gráfica no permite observar bien la curva).

Así pues, si queremos saber cuanto deberíamos pagar por una casa de, digamos 100 metros, sólo hay que calcularlo sustituyendo b en la ecuación de la recta: $y=43.41+1.08 \cdot 100=151.41$
Es importante elegir bien el número de iteraciones del algoritmo. Si nos quedamos cortos nuestra recta no se ajustará bien, y si nos pasamos llega un momento en el que la mejora no es apreciable. Un método posible es detener las iteraciones cuando se compruebe que la función ya no se minimiza más allá de un umbral pequeño en cada iteración. Para que os hagáis una idea de como afecta este parámetro os dejo una animación en la que se muestra la recta calculada para diferentes valores del número de iteraciones.

Este artículo trata de ser didáctico, por lo que si quieres usar regresión lineal en la vida real, te aconsejo que mejor uses alguna buena librería como Scikit-learn que tiene un estupendo soporte para hacer este tipo de regresión.
Os dejo el enlace a mi GitHub con el notepad de Jupyter que contiene todo el código de este artículo: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/Regresion%20lineal.ipynb

Manipulación de datos con pandas

noreply@blogger.com (Alberto García Serrano) — Mon, 11 Mar 2019 15:36:00 +0000

Cuando uno lee un libro o un artículo sobre machine learning encuentra multitud de explicaciones sobre el algoritmo tal o cual. Sin embargo, no se habla demasiado sobre la manipulación y el limpiado de los datos, que bajo mi punto de vista es tan o incluso más importante que utilizar el algoritmo adecuado. Nuestro aliado en esta tarea es la librería pandas.
En lugar de hacer un recorrido exhaustivo por las funcionalidades de la librería, he preferido hacer uso de ella con un dataset real, para poder ver así, no sólo cuáles son sus funcionalidades, sino también cómo se aplican a datos reales. Así pues, en lugar de usar un dataset de los clásicos, he recurrido a uno real sacado de la web de datos abiertos del Ayuntamiento de Málaga. En concreto vamos a trabajar con el siguiente dataset, que se corresponde con las lecturas energéticas de los cuadros eléctricos durante el mes de marzo de 2017: https://datosabiertos.malaga.eu/dataset/lecturas-cuadros-electricos-marzo-2017.
Como siempre os dejo un enlace al notebook python con todo lo que me dispongo a contaros: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/Manipulacion%20datos%20Pandas.ipynb
El dataset es un archivo CSV y tiene los siguientes campos:

_id: Id Autonumérico.
ID_NODO: Identificador del nodo.
ID_DIA: Identificador del día de la toma (número de día de mes).
FECHA_TOMA_DATO: Fecha y hora de la toma del dato.
ENERGIA_ACTIVA_ACUMULADA: Energía activa acumulada.
LECTURA_ENERGIA_ACTIVA: Energía activa en el momento de la toma.
ENERGIA_REACTIVA_ACUMULADA: Energía reactiva acumulada.
LECTURA_ENERGIA_REACTIVA: Energía reactiva en el momento de la toma.
POTENCIA_TOTAL_ACUMULADA: Acumulado de la potencia total manejada por el cuadro (expresada en W).
LECTURA_POTENCIA_TOTAL: Potencia total manejada por el cuadro (expresada en W).

Comenzamos importando y cargando el archivo CSV. Para realizar la carga se utiliza la función read_csv(). Por defecto, esta función espera los datos separados por comas, pero si echas un ojo a nuestro fichero veras que el separador es el punto y coma. Por eso usamos el parámetro sep=";". En la variable data_mar se almacenan los datos. Pandas devuelve los datos en una estructura llamada DataFrame. Con la función head() podemos ver las primeras filas de los datos que hemos cargado.

import pandas as pd
data_mar = pd.read_csv("201703.csv", sep=";")
data_mar.head()

	ID_NODO	ID_DIA	FECHA_TOMA_DATO	ENERGIA_ACTIVA_ACUMULADA	LECTURA_ENERGIA_ACTIVA	ENERGIA_REACTIVA_ACUMULADA	LECTURA_ENERGIA_REACTIVA	POTENCIA_TOTAL_ACUMULADA
0	81848	28	2017-02-28 23:11:54	66884.601562	0.000000	46464.699219	0.000000	0.000000
1	98372	28	2017-02-28 23:11:54	264010.593750	0.000000	169612.296875	0.000000	0.000000
2	98526	28	2017-02-28 23:11:54	56888.000000	0.398438	61891.300781	0.000000	14.000000
3	98680	28	2017-02-28 23:11:54	35347.699219	0.398438	37424.699219	0.398438	77.006493
4	98834	28	2017-02-28 23:11:54	93080.703125	0.000000	47944.500000	0.000000	0.000000

Podemos empezar por obtener algunos estadísticos básicos con el objetivo de empezar a conocer con qué datos nos estamos enfrentando. Para ello usamos la función describe().

data_mar.describe()

	ID_NODO	ID_DIA	ENERGIA_ACTIVA_ACUMULADA	LECTURA_ENERGIA_ACTIVA	ENERGIA_REACTIVA_ACUMULADA	LECTURA_ENERGIA_REACTIVA	POTENCIA_TOTAL_ACUMULADA	LECTURA_POTENCIA_TOTAL
count	413108.000000	413108.000000	4.131080e+05	4.131080e+05	4.131080e+05	4.131080e+05	413108.000000	413108.000000
mean	71962.043165	15.982503	2.330070e+05	1.268065e+02	1.679989e+05	1.002321e+02	2291.475042	42.950668
std	26718.517878	8.936663	1.130742e+06	2.271224e+04	8.109683e+05	1.819730e+04	12654.772009	1263.035618
min	540.000000	1.000000	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000	0.000000
25%	70605.000000	8.000000	0.000000e+00	0.000000e+00	0.000000e+00	0.000000e+00	0.000000	0.000000
50%	76150.000000	16.000000	2.475400e+03	0.000000e+00	3.182920e+04	0.000000e+00	0.000000	0.000000
75%	81540.000000	24.000000	1.568925e+05	0.000000e+00	9.713650e+04	0.000000e+00	3.000000	0.000000
max	103620.000000	31.000000	1.267291e+07	1.148814e+07	9.325143e+06	9.117862e+06	234922.625488	216123.260318

Por cada variable numérica obtenemos el número de elementos, la media, la desviación típica, el máximo y mínimo valor, así como los percentiles 25%, 50% y 75%. También podemos obtener el número de elementos de la tabla, así como el número de columnas.

data_mar.shape

(413108, 9)

El dataset contiene información sobre lecturas de contadores eléctricos. En concreto mide la energía activa y la reactiva. Muchas veces, cuando uno estudia un dataset, no domina el campo sobre el cuál tratan los datos, por lo toca investigar un poco. Os resumo un poco en qué consiste esto de la energía activa y reactiva:

Energía activa: los receptores eléctricos alimentados por corriente eléctrica transforman la energía eléctrica en trabajo mecánico y en calor. A este efecto útil se le denomina “energía activa” y se mide en kWh.
Energía reactiva: existen numerosos receptores, tales como motores, transformadores, reactancias, etc., que para funcionar necesitan que se formen campos magnéticos. Estos equipos, en general inductivos, absorben energía de la red para crear los campos magnéticos y la devuelven mientras desaparecen. Con este intercambio de energía, se provoca un consumo suplementario que no es aprovechable por los receptores. Se mide en kVArh. La energía reactiva provoca una sobrecarga en líneas, transformadores y generadores, sin llegar a producir un rendimiento útil. Sin embargo, la factura de energía sí la contabiliza, por lo que puede llegar a incrementarla en cantidades importantes

Vamos a centrarnos en el estudio de la energía acumulada, así que por ahora no nos van a hacer falta los campos de lectura instantánea de energía activa y reactiva. Una buena práctica, sobre todo si el dataset es muy grande y andamos cortos de memoria, es eliminar las columnas que no vamos a necesitar, así que vamos a eliminarlas.

# Borramos los campos que no vamos a necesitar: 
# LECTURA_ENERGIA_ACTIVA, LECTURA_ENERGIA_REACTIVA, LECTURA_POTENCIA_TOTAL
campos = data_mar.columns.values.tolist()
eliminar = ["LECTURA_ENERGIA_ACTIVA", "LECTURA_ENERGIA_REACTIVA", \
    "LECTURA_POTENCIA_TOTAL"]
data_mar = data_mar[[e for e in campos if e not in eliminar]]

También nos puede interesar a veces crear nuevos campos en los que la misma información esté en un formato diferente que nos facilite la operación los datos. Por ejemplo, el campo de fecha y hora están en formato string, y nos puede interesar que sea un campo numérico para poder representarlo cómodamente en una gráfica. Como ejemplo, voy a crear un campo nuevo con la fecha en formato numérico.

# vamos a crear un campo nuevo para la fecha en formato timestamp, 
# que nos será más cómodo a la hora de representar las gráficas.
data_mar["FECHA_NUM"] = pd.to_numeric(pd.to_datetime(data_mar["FECHA_TOMA_DATO"]))
data_mar.head()

	ID_NODO	ID_DIA	FECHA_TOMA_DATO	ENERGIA_ACTIVA_ACUMULADA	ENERGIA_REACTIVA_ACUMULADA	POTENCIA_TOTAL_ACUMULADA	FECHA_NUM
0	81848	28	2017-02-28 23:11:54	66884.601562	46464.699219	0.000000	1488323514000000000
1	98372	28	2017-02-28 23:11:54	264010.593750	169612.296875	0.000000	1488323514000000000
2	98526	28	2017-02-28 23:11:54	56888.000000	61891.300781	14.000000	1488323514000000000
3	98680	28	2017-02-28 23:11:54	35347.699219	37424.699219	77.006493	1488323514000000000
4	98834	28	2017-02-28 23:11:54	93080.703125	47944.500000	0.000000	1488323514000000000

Vamos a examinar un nodo para ver cómo se comportan los datos para un sólo individuo de la muestra. Primero necesitamos averiguar cuántos nodos (contadores eléctricos) hay. Con la función unique() podemos filtrar el campo ID_NODO para ver cuántos nodos distintos hay.

# ¿cuantos nodos tenemos?
len(data_mar["ID_NODO"].unique())

139

Vamos a estudiar un nodo en concreto. El ID_NODO=542

# vamos a estudiar un nodo al azar: el ID_NODO = 542
data_mar[data_mar["ID_NODO"] == 542].plot(x="FECHA_NUM", y="ENERGIA_ACTIVA_ACUMULADA")

Observamos en la gráfica que a partir de cierta fecha el acumulado de energía activa es cero. Si probamos con otros nodos vemos que el patrón se repite. Además, si examinamos las fechas de las lecturas vemos que hay datos del mes anterior (febrero). Un posible interpretación que podemos hacer es que el periodo de lecturas no comienza exactamente el día 1 del mes en curso. Podemos consultar para qué fecha y hora la lectura del nodo 542 es cero.

# veamos las fechas para las que las lecturas son cero
data_mar[(data_mar["ENERGIA_ACTIVA_ACUMULADA"]==0) & (data_mar["ID_NODO"]==542)].head()

	ID_NODO	ID_DIA	FECHA_TOMA_DATO	FECHA_NUM
337194	542	26	2017-03-26 06:27:10	1490509630000000000
337342	542	26	2017-03-26 06:42:10	1490510530000000000
337468	542	26	2017-03-26 06:57:10	1490511430000000000
337492	542	26	2017-03-26 07:12:10	1490512330000000000
337759	542	26	2017-03-26 07:27:10	1490513230000000000

Tenemos dos opciones. Podemos borrar los datos anteriores o posteriores a una fecha concreta. Otra opción es eliminar aquellas lecturas cuyo valor para el campo ENERGIA_ACTIVA_ACUMULADA sea 0. Para ver un ejemplo de cada caso, vamos a borrar las filas que pertenecen al mes anterior y seguidamente vamos a eliminar todas las filas cuyo campo ENERGIA_ACTIVA_ACUMULADA sea 0.

# borramos datos del mes anterior
data_mar = data_mar[data_mar["FECHA_TOMA_DATO"]>"2017-03"]
data_mar.head()

	ID_NODO	ID_DIA	FECHA_TOMA_DATO	ENERGIA_ACTIVA_ACUMULADA	ENERGIA_REACTIVA_ACUMULADA	POTENCIA_TOTAL_ACUMULADA	FECHA_NUM
556	542	1	2017-03-01 00:11:54	360591.906250	121981.601562	0.000000	1488327114000000000
557	543	1	2017-03-01 00:11:54	300521.894531	242683.304688	0.000000	1488327114000000000
558	544	1	2017-03-01 00:11:54	2475.300049	51700.800781	80.024994	1488327114000000000
559	545	1	2017-03-01 00:11:54	0.000000	26725.900391	0.000000	1488327114000000000
560	546	1	2017-03-01 00:11:54	0.000000	0.000000	0.000000	1488327114000000000

# borramos las filas donde se cumpla que 
# ENERGIA_ACTIVA_ACUMULADA = 0 para todos los nodos
data_mar = data_mar[data_mar["ENERGIA_ACTIVA_ACUMULADA"]>0]
data_mar.shape

(237946, 7)

Podemos volver a examinar la gráfica que, ahora sí, nos permite comprender mejor cómo es la demanda de consumo para ese nodo.

# creamos un dataframe con los datos del nodo 542
data_mar_542 = data_mar[data_mar["ID_NODO"] == 542]
data_mar_542.plot(x="FECHA_NUM", y="ENERGIA_ACTIVA_ACUMULADA")

Vamos a comprar la energía activa y la reactiva poniéndolas juntas en una gráfica.

data_mar_542.plot(x="FECHA_NUM", y=["ENERGIA_ACTIVA_ACUMULADA", \
    "ENERGIA_REACTIVA_ACUMULADA"], subplots=True )

Aunque los valores de la energía reactiva son muy inferiores a los de la energía activa, parece que hay una correlación bastante acusada entre ambos tipos de energía. De hecho, si usamos la función corr() de Pandas para calcular la correlación entre ambos campos, vemos que es una correlación casi perfecta (0.99).

data_mar_542["ENERGIA_ACTIVA_ACUMULADA"].corr(data_mar_542["ENERGIA_REACTIVA_ACUMULADA"])

0.99990652807701275

La energía reactiva es un tipo de energía eléctrica, que absorben de la red algunos equipos eléctricos pero que luego la devuelven, por lo que no supone un consumo, aunque sí hay que generarla y transportarla hasta los equipos. Por ello las compañías eléctricas penalizan en factura el consumo de energía reactiva a partir de un cierto valor.
Si quisieramos hacer un análisis de en qué nodos hay un consumo anomalo de energía reactiva, podríamos ver si la relación entre la energía activa y reactiva tienen una baja correlación (aquí estoy haciendo muchas suposiciones que podrían no ser ciertas a nivel eléctrico, pero valga como ejemplo).
¿Se observará la misma correlación en todos los nodos? Vamos a crear un nuevo DataFrame con la correlación entre ambas lecturas para cada contador.

# creamos un nuevo dataframe con los campos NODO_ID y la correlación
c = []
for i in data_mar["ID_NODO"].unique().tolist():
    c.append([i, data_mar[data_mar["ID_NODO"]==i]["ENERGIA_ACTIVA_ACUMULADA"].corr(data_mar["ENERGIA_REACTIVA_ACUMULADA"])])

correlaciones = pd.DataFrame(c, columns=["ID_NODO", "CORRELACION"])
correlaciones.head(10)

	ID_NODO	CORRELACION
0	542	0.999907
1	543	0.961227
2	544	0.999266
3	547	0.999963
4	33756	0.998668
5	71530	NaN
6	71838	NaN
7	71992	NaN
8	72454	0.999997
9	73224	NaN

Observamos que hay correlaciones que toman el valor NaN (Not A Number). ¿Qué ha pasado? Si generamos la gráfica de uno de los nodos cuya correlación es NaN vemos que el valor de ambas variables es siempre constante, es decir, no crece en todo el periodo de tiempo que estamos observando. Por defecto, Pandas calcula la correlación de pearson, cuya formula es cor(i,j) = cov(i,j)/[stdev(i)*stdev(j)]. Como la desviación estandad es cero y un divisor no puede ser cero, obtenemos un bonito NaN. Lo cuál, dicho sea de paso, nos viene de perlas para ver cómo tratar con valores nulos en un DataFrame.

data_mar[data_mar["ID_NODO"]==71530].plot(x="FECHA_NUM", \
    y=["ENERGIA_ACTIVA_ACUMULADA","ENERGIA_REACTIVA_ACUMULADA"], subplots=True )

Vamos a ver qué opciones tenemos para tratar con los valores nulos. Podemos:

Eliminar las filas con valor nulo (con el método dropna())
Rellenar las filas con un valor calculado (con el método fillna())

# rellenar filas NaN con un valor concreto (0 en este caso)
cor_cero = correlaciones["CORRELACION"].fillna(0)
cor_cero.head(10)

0    0.999907
1    0.961227
2    0.999266
3    0.999963
4    0.998668
5    0.000000
6    0.000000
7    0.000000
8    0.999997
9    0.000000
Name: CORRELACION, dtype: float64

O podemos rellenar con la media.

# rellenar filas con el valor medio de la columna (o cualquier otra función))
cor_media = correlaciones["CORRELACION"].fillna(correlaciones["CORRELACION"].mean())
cor_media.head(10)

0    0.999907
1    0.961227
2    0.999266
3    0.999963
4    0.998668
5    0.971874
6    0.971874
7    0.971874
8    0.999997
9    0.971874
Name: CORRELACION, dtype: float64

Si preferimos eliminar las filas con NaN podemos usar el método dropna().

# pero en este caso vamos a optar por eliminarlas
# si axis=0 borra por filas. Si axis=1 borra por columnas
correlaciones = correlaciones.dropna(axis=0)
correlaciones.head(10)

	ID_NODO	CORRELACION
0	542	0.999907
1	543	0.961227
2	544	0.999266
3	547	0.999963
4	33756	0.998668
8	72454	0.999997
10	73840	0.999995
11	74918	0.999972
12	75072	0.999946
13	75534	0.999648

Si lo que queremos es estudiar cómo se comporta el consumo en los cuadros eléctricos donde hay baja correlación, podemos ordenar la lista para ver cuáles son más bajas.

correlaciones = correlaciones.sort_values(by="CORRELACION")
correlaciones.head()

	ID_NODO	CORRELACION
90	52061	0.514461
61	75688	0.617365
84	52775	0.649950
67	77382	0.826462
40	100682	0.885768

Echemos un vistazo a las gráficas de energía activa y reactiva para el nodo con la correlación más baja.

data_mar[data_mar["ID_NODO"]==52061].plot(x="FECHA_NUM", \ 
    y=["ENERGIA_ACTIVA_ACUMULADA","ENERGIA_REACTIVA_ACUMULADA"], subplots=True )

A partir de esta información podríamos deducir que a ese cuadro eléctrico se conecta maquinaria que hace uso intensivo de campos eléctricos, como bobinas o motores.
Vamos ahora a ver cómo se comporta el campo POTENCIA_TOTAL_ACUMULADA. Miramos el primer nodo.

data_mar[data_mar["ID_NODO"]==52061].plot(x="FECHA_NUM", \
    y=["POTENCIA_TOTAL_ACUMULADA"], subplots=True )

Está a cero, y si observamos otros nodos vemos que también. Como no tenemos más detalles sobre el dataset sólo podemos hacer suposiciones (como que por error no se ha rellenado este campo). De todas formas nos viene bien para ver cómo podemos combinar valores de dos o más campos para volcarlos en otro. Vamos a suponer que en este campo, el valor que debería haber es la suma de la energía activa y reactiva (en realidad ambos valores están en unidades distintas, pero por simplicidad supondremos que ambos valores son Kw/h y que podemos sumarlos alegremente). Vamos a combinar pues las suma de los dos valores en el campo POTENCIA_TOTAL_ACUMULADA.

data_mar["POTENCIA_TOTAL_ACUMULADA"] = data_mar["ENERGIA_ACTIVA_ACUMULADA"] \
    + data_mar["ENERGIA_REACTIVA_ACUMULADA"]
data_mar.head()

	ID_NODO	ID_DIA	FECHA_TOMA_DATO	ENERGIA_ACTIVA_ACUMULADA	ENERGIA_REACTIVA_ACUMULADA	POTENCIA_TOTAL_ACUMULADA	FECHA_NUM
556	542	1	2017-03-01 00:11:54	3.605919e+05	1.219816e+05	4.825735e+05	1488327114000000000
557	543	1	2017-03-01 00:11:54	3.005219e+05	2.426833e+05	5.432052e+05	1488327114000000000
558	544	1	2017-03-01 00:11:54	2.475300e+03	5.170080e+04	5.417610e+04	1488327114000000000
561	547	1	2017-03-01 00:11:54	4.085408e+06	3.286907e+06	7.372315e+06	1488327114000000000
566	33756	1	2017-03-01 00:11:54	1.144733e+07	8.198532e+06	1.964586e+07	1488327114000000000

En este momento tenemos, por un lado, el DataFrame con las lecturas de los cuadros y por otro un DataFrame con las correlaciones que hemos calculado. ¿Podemos agregar las correlaciones a la tabla principal? Claro. Esta operación se llama unión (join en inglés). Hay cuatro tipos: inner join, left join, right join y outer join. No voy a entrar en detalle ya que son los mismos que los de SQL. Si queréis saber más podéis consultar aquí: https://es.wikipedia.org/wiki/Join

# unimos los DataFrames con left join, ya que el 
# segundo dataframe (correlaciones) tiene menos filas 
# y por lo tanto, no todos los nodos disponen de datos de correlación 
# (recuerda que borramos los NaN)
data_mar = pd.merge(left=data_mar, right=correlaciones, how="left", \
    left_on="ID_NODO", right_on="ID_NODO")
data_mar.head()

	ID_NODO	ID_DIA	FECHA_TOMA_DATO	ENERGIA_ACTIVA_ACUMULADA	ENERGIA_REACTIVA_ACUMULADA	POTENCIA_TOTAL_ACUMULADA	FECHA_NUM	CORRELACION
0	542	1	2017-03-01 00:11:54	3.605919e+05	1.219816e+05	4.825735e+05	1488327114000000000	0.999907
1	543	1	2017-03-01 00:11:54	3.005219e+05	2.426833e+05	5.432052e+05	1488327114000000000	0.961227
2	544	1	2017-03-01 00:11:54	2.475300e+03	5.170080e+04	5.417610e+04	1488327114000000000	0.999266
3	547	1	2017-03-01 00:11:54	4.085408e+06	3.286907e+06	7.372315e+06	1488327114000000000	0.999963
4	33756	1	2017-03-01 00:11:54	1.144733e+07	8.198532e+06	1.964586e+07	1488327114000000000	0.998668

Los parámetros left y right indican los DataFrames que vamos a fundir. El parámetro how nos permite indicar el tipo de join, en este caso left join. Finalmente left_on y right_on nos sirve para indicar cuál es el campo clave para unir ambos DataFrames. En este caso el campo común es ID_NODO.
Como dije al principio del artículo, mi intención no era hacer una descripción exhaustiva de toda la funcionalidad de pandas, ya que creo que podía ser más clarificador ver las operaciones básicas sobre un dataset real.

Reconocimiento de caracteres manuscritos con k-nn

noreply@blogger.com (Alberto García Serrano) — Thu, 28 Feb 2019 22:10:00 +0000

En la entrada anterior introducimos el algoritmo k nearest neighbors o de los vecinos más cercanos. Os conté como funciona el algoritmo y os presenté un ejemplo que clasificaba puntitos en un espacio en dos dimensiones. Como sé que os quedasteis con ganas de más (sí, aún os oigo suplicar que os cuente más) vamos a poner a prueba al algoritmo k-nn para ver si es capaz de aprender a reconocer caracteres manuscritos (problema conocido como OCR o reconocimiento óptico de caracteres). Como en el anterior artículo, os dejo un enlace al notebook de Jupyter que he usado para escribir este post: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/KNN%20Image%20Classifier.ipynb
Para entrenar nuestro reconocedor de caracteres vamos a usar un dataset bastante conocido: MNIST. Desde esta web podéis descargar los cuatro archivos que lo componen:
train-images-idx3-ubyte.gz: imágenes de entrenamiento
train-labels-idx1-ubyte.gz: etiquetas de las imágenes de entrenamiento
t10k-images-idx3-ubyte.gz: conjunto de imágenes para test
t10k-labels-idx1-ubyte.gz: conjunto de etiquetas de las imágenes de test

Este dataset contiene miles de imágenes con dígitos escritos a mano. En la imagen de cabecera de este artículo puedes ver qué pinta tienen los dígitos. Las imágenes están codificadas dentro del archivo con un formato que está descrito en la web. No voy a entrar en detalles porque todo está ahí. Simplemente os pongo el código que he usado para cargarlas e introducirlas en matrices de NumPy.

def loadMNIST( prefix, folder ):
    intType = np.dtype( 'int32' ).newbyteorder( '>' )
    nMetaDataBytes = 4 * intType.itemsize

    data = np.fromfile( folder + "/" + prefix + '-images.idx3-ubyte', dtype = 'ubyte' )
    magicBytes, nImages, width, height = np.frombuffer( data[:nMetaDataBytes].tobytes(), intType )
    data = data[nMetaDataBytes:].astype( dtype = 'float32' ).reshape( [ nImages, width, height ] )

    labels = np.fromfile( folder + "/" + prefix + '-labels.idx1-ubyte',
                          dtype = 'ubyte' )[2 * intType.itemsize:]

    return data, labels

images_tr, labels_tr = loadMNIST( "train", "./mnist" )
images_te, labels_te = loadMNIST( "t10k", "./mnist" )
# imagenes en array de 60000 x 28 x 28 -> 60000 imagenes de 28x28

El conjunto de entrenamiento tiene 60.000 imágenes y el de test, es decir, el que usaremos para poner a prueba el sistema una vez entrenado, dispone de 10.000 imágenes. Cada una de estas imágenes es una matriz de 28x28 píxeles con un dígito del 0 al 9 dibujado a mano. Si recordáis, en el ejemplo del artículo anterior manejábamos datos de dos dimensiones (por ejemplo, el sueldo y la deuda a la hora de conceder o no un préstamo). Ahora trabajamos con datos de 784 dimensiones (28*28 píxeles). Así que por lo visto vamos a tener que sumar, restar, elevar al cuadrado y hacer raíces cuadradas de matrices de 784 elementos. No está mal ¿verdad? veremos como le sienta esto a NumPy.

# coger un dígito aleatorio del grupo de test
i = random.randint(0,images_te.shape[0])
img_test = images_te[i].flatten()
label_test = labels_te[i]

Empezamos por elegir una imagen aleatoria del conjunto de test para ver si somos capaces de clasificarla correctamente y reconocer el dígito. Para operar con la imagen la aplanamos (flatten), es decir, ponemos toda la matriz en una sola fila (matriz de 1x784). También guardamos su etiqueta para compararla con la clasificación que haga el algoritmo.

# buscamos los vecinos más cercanos (KNN)
k = 5 #número de vecinos

distances = []
for i in range(images_tr.shape[0]):
    dist = np.sqrt(np.sum(np.square(images_tr[i].flatten() - img_test)))
    distances.append((dist, labels_tr[i])) # guardamos las etiquetas y la distancia

#ordenamos por distancia y nos quedamos con los k vecinos más cercanos
distances.sort(key=lambda x: x[0])
neighbors = distances[:k]

Ahora buscamos a los k vecinos más cercanos (en este caso k=5). Si recuerdas el ejemplo anterior, el código era algo más complejo debido a que estábamos clasificando varios datos de test. En este caso sólo clasificamos un dato (el dígito que hemos seleccionado aleatoriamente en el código de más arriba). No entro en mucho detalle ya que en el artículo anterior ya os conté cómo funcionaba. Si no entiendes qué hace este fragmento de código dale un repaso a lo que ya os conté sobre k-nn.

# contamos los votos para ver qué etiqueta gana
votes = [0,0,0,0,0,0,0,0,0,0]
for neighbor in neighbors:
    votes[neighbor[1]] = votes[neighbor[1]] + 1
# obtenemos la etiqueta ganadora
pred_label = votes.index(max(votes))

Una vez elegidos los k vecinos más cercanos toca votar (k-nn es una democracia). Por lo tanto, recontamos las etiquetas de los k vecinos más cercanos y la ganadora será nuestra predicción.

print ("Predicted: " + str(pred_label))
print ("Real: " + str(label_test))
img = plt.imshow(img_test.reshape(28,28), cmap="gray")
display(img)

La siguiente imagen muestra una captura de la ejecución del código.

Y ¡parece que funciona!
Pero ¿qué tasa real de acierto consigue k-nn con el dataset MNIST? Si tenéis curiosidad os propongo que probéis con todo el conjunto de test. Eso sí, son 10.000 imágenes de prueba, así que un poco de paciencia cuando lancéis la ejecución. En la prueba que yo he hecho contra todo el conjunto de test, de las diez mil imágenes el algoritmo ha clasificado correctamente 9.688. Esto es, una tasa de acierto de 0.968 (9.688 / 10.000). O lo que es lo mismo, ha acertado un 96,8% de las veces. No está nada mal ¿verdad?

Clasificación con k-nearest neighbors

noreply@blogger.com (Alberto García Serrano) — Fri, 22 Feb 2019 22:10:00 +0000

En un artículo anterior hablamos de NumPy y me hubiera gustado completarlo con algún ejemplo concreto. Así que voy a aprovechar este artículo para presentaros un algoritmo muy usado en machine learning y de camino a ilustrarlo con código Python (y la librería NumPy).
El algoritmo k-nearest neighbors (abreviado como k-nn) o en español, k vecinos más cercanos, es un algoritmo conceptualmente sencillo y también fácil de implementar, pero que para ciertas tareas puede tener una eficacia más que razonable. En lo que sigue os voy a contar cómo funciona y lo vamos a implementar en Python. En un siguiente artículo, aplicaremos el algoritmo a un problema del mundo real.
K-nn pertenece a la familia de algoritmos de clasificación supervisada, es decir, que necesitan la intervención humana para realizar su aprendizaje. Esta familia de algoritmos son capaces de clasificar elementos dentro de una serie de categorías prefijadas. De forma genérica, esta tarea se llama entrenamiento, y consiste en ofrecer al sistema ejemplos de cada una de las categorías posibles (a este conjunto de datos se le suele llamar dataset o conjunto de entrenamiento). Cuando se ingresan nuevos elementos al sistema éste lo asignará una de las categorías basándose en lo "aprendido" durante el entrenamiento. El proceso de asignar una categoría se llama etiquetado. Hay multitud de algoritmos de clasificación siendo K-nn uno de los más sencillos, por lo que es ideal para empezar a dar los primeros pasos en el aprendizaje automático.
Cuando entrenamos un clasificador hay que elegir en qué características del conjunto de entrenamiento ha de fijarse a la hora de clasificarlos. Por ejemplo, si un banco quiere crear un sistema que ayude a decidir si se aprueba la solicitud de préstamo de un cliente, un par de datos interesantes es el sueldo y las deudas que tenga el cliente (en realidad habría que considerar más, pero vamos a simplificar). En machine learning las características que vamos a usar para entrenar al sistema se llaman atributos (features en inglés). Como tenemos dos atributos vamos a trabajar con datos de dimensión dos. Esto quiere decir que, por cada ejemplo introducido durante el entrenamiento (obtenidos de antiguos clientes), informaremos de los dos atributos (sueldo y deuda) además de a la clase a la que pertenece, que en este caso podríamos etiquetar como "Pagó" o "No pagó" para indicar si ese cliente en concreto pagó el préstamo o por el contrario tuvo problemas para hacerlo.
Dado que sólo tenemos dos atributos (sueldo y deuda) podríamos representar fácilmente en una gráfica. Si el eje de ordenadas representa el sueldo y el de las abscisas representa la deuda, podemos dibujar un punto por cada cliente. Siendo el punto de un color diferente dependiendo de la clase a la que pertenece. En la siguiente figura hemos representado con puntos azules los clientes que pagaron y de rojo los que no (los valores de los ejes están en miles de euros). La pregunta es ¿debemos conceder el préstamo al nuevo cliente representado por el punto verde?

He separado visualmente las dos clases para que quede más claro. Como no siempre será obvio cómo se debe etiquetar el nuevo dato, necesitamos un algoritmo que nos diga si el punto verde debe ser categorizado en una clase o en la otra. El algoritmo es muy sencillo:

k = número de vecinos a tener en cuenta
Calculamos la distancia entre el nuevo punto y todos los del conjunto de entrenamiento
Nos quedamos con los k vecinos más cercanos y etiquetamos el nuevo punto según la clase a la que pertenezcan la mayoría de sus k vecinos

La idea es que el nuevo punto pertenecerá a la misma clase que la de sus vecinos más cercanos, así que la primera decisión a tomar es ¿cuántos vecinos son suficientes vecinos? No hay una respuesta fácil y tampoco quiero entrar en mucho detalle aquí para no confundirte demasiado. La respuesta corta es: prueba con distintos valores para ver cuál funciona mejor con tu conjunto de datos. la respuesta larga la dejamos para otro día. Si observas la siguiente figura, la etiqueta que le corresponde al nuevo punto no es la misma si usamos tres vecinos que si usamos cinco.

Otra concepto que no he definido es el de la distancia entre puntos. Hay varias medidas de distancia, cada una con sus propiedades, sus ventajas y desventajas. Tampoco quiero entrar ahora a explicar cuál es más conveniente para cada aplicación, así que nos quedamos con una que seguramente te suena de tus días de instituto: la distancia euclidiana.

Te refresco la memoria. La distancia entre dos puntos P y Q en un plano (si Pitágoras no andaba muy desencaminado) se obtiene con la siguiente ecuación.
\[d_E(P,Q)=\sqrt{(p_1-q_1)^2 + (p_2-q_2)^2}\]
Si te estás preguntando qué ocurre si nuestros atributos tienen una dimensionalidad mayor de dos, tranquilo. Lo bueno del Álgebra es que funciona bien en dos dimensiones, en tres, en cuatro o en las que haga falta, así que podemos generalizar la ecuación para cualquier cantidad de atributos.
\[d_E(P,Q)=\sqrt{\sum_{i=1}^n (p_i-q_i)^2}\]
Lo bueno de esta ecuación es que usando cálculo matricial un ordenador puede resolverla en un plis plas. Y aquí es donde NumPy nos va a mostrar su verdadero poder (bueno, con dos atributos no lo vamos a notar mucho, pero cuando sean miles será de agradecer).
Hechas las presentaciones, vamos al lío. He implementado el algoritmo en un notebook para Jupyter que os he dejado en el siguiente enlace: https://github.com/albgarse/InteligenciaArtificial/blob/master/Machine%20Learning/KNN.ipynb

Empezamos importando las bibliotecas que vamos a necesitar.

import numpy as np
import random
from IPython.display import display
import matplotlib.pyplot as plt

Vamos a utilizar NumPy para las operaciones con matrices y Matplotlib para generar algunas gráficas. Para mostrar las gráficas en el notebook Jupyter se importa también la función display.

n = 50 # 50 muestras aleatorias
# generamos dos centros para las dos clases de puntos para entrenamiento 
c1 = [random.randint(0,1000), random.randint(0,1000)]
c2 = [random.randint(0,1000), random.randint(0,1000)]
# generamos las muestras aleatorias alrededor de los centros
tuplasC1 = []
tuplasC2 = []
labelsC1 = []
labelsC2 = []
for i in range(int(n/2)):
    tuplasC1.append([c1[0] + random.randint(-100,100), c1[1] + random.randint(-100,100)])
    labelsC1.append(0)
    tuplasC2.append([c2[0] + random.randint(-100,100), c2[1] + random.randint(-100,100)])
    labelsC2.append(1)

labels = labelsC1 + labelsC2
puntos = np.matrix(tuplasC1 + tuplasC2)
# dibujamos los puntos
plt.scatter([puntos[:int(n/2),0]], [puntos[:int(n/2),1]], c="b")
plt.scatter([puntos[int(n/2):,0]], [puntos[int(n/2):,1]], c="g")
display(plt)

En machine learning tan importante como el algoritmo que hayamos decidido usar son los datos con los que lo vamos a alimentar. En este caso vamos a fabricarlos de forma aleatoria. El código anterior genera cincuenta muestras (o puntos o ejemplos, como prefieras) pero he hecho un poco de trampa. Los he creado alrededor de dos puntos distantes para que visualmente queden bien diferenciados. Para mantener la generalidad en todo momento he etiquetado con el valor 0 los puntos pertenecientes a la primera clase y con 1 los de la segunda. Una vez tenemos los datos los ponemos en una matriz de NumPy para poder operar con ellos. Finalmente mostramos la gráfica con los puntos, en azul los de la clase 0 y en verde los de la clase 1. Este es el resultado.

También necesitamos datos nuevos, es decir, aquellos que queremos clasificar.

# generamos puntos aleatorios nuevos para clasificarlos
n_test = 10
tuplas = []
for i in range(n_test):
    tuplas.append([random.randint(0,1000), random.randint(0,1000)])
    
puntos_test = np.matrix(tuplas)
# dibujamos los nuevos puntos junto con los anteriores
plt.scatter([puntos[:int(n/2),0]], [puntos[:int(n/2),1]], c="b")
plt.scatter([puntos[int(n/2):,0]], [puntos[int(n/2):,1]], c="g")
plt.scatter([puntos_test[:,0]], [puntos_test[:,1]], c='r', marker='x')
display(plt)

Generamos diez puntos nuevos de forma totalmente aleatoria y mostramos la gráfica para ver por donde han caído (se muestran como X de color rojo).

Ha llegado el momento de la verdad, veamos cómo se comporta k-nn con nuestros datos.

# usamos KNN para clasificar los nuevos puntos
k = 5 # número de vecinos

pred_label = []
# clasificamos cada uno de los puntos nuevos
for i in range(puntos_test.shape[0]):
    distances = []
    # por cada punto calculamos la distancia con los puntos de entrenamiento 
    for j in range(puntos.shape[0]):
        dist = np.sqrt(np.sum(np.square(puntos[j] - puntos_test[i])))
        distances.append((dist, labels[j])) # guardamos las etiquetas y la distancia

    # ordenamos por distancia y nos quedamos con los k vecinos más cercanos
    distances.sort(key=lambda x: x[0])
    neighbors = distances[:k]
    # contamos los votos para ver qué etiqueta gana
    votes = [0,0]
    for neighbor in neighbors:
        votes[neighbor[1]] = votes[neighbor[1]] + 1
    # obtenemos la etiqueta ganadora
    pred_label.append(votes.index(max(votes)))

He preferido hacer una implementación lo más clara posible en detrimento de la eficiencia, ya que con NumPy se puede realizar toda la operación en un sólo paso sin necesidad de usar los dos bucles que he usado en el código, pero creo que así es mucho más intuitivo. Pues bien, el primer bucle i va recorriendo la matriz de los puntos a etiquetar (recuerda que eran 10) y en el bucle interno j los comparamos con cada uno de los ejemplos o puntos del dataset de entrenamiento (recuerda que eran 50). Como ves no hay mucho misterio: calculamos la distancia entre ambos puntos operando matricialmente con NumPy y almacenamos dicha distancia junto a la clase (etiqueta) a la que pertenece el punto con el que estamos comparando. Una vez hemos medido todas las distancias las ordenamos y nos quedamos con las k menores (en este caso k=5). Ahora necesitamos saber a qué clases pertenecen los k vecinos más cercanos. Como hemos almacenado su etiqueta junto con la distancia, sometemos la decisión a un proceso de votación. Cada uno de los k vecinos vota con su etiqueta y la ganadora se asigna al nuevo punto.

# separar los datos clasificados para dibujarlos
g0 = []
g1 = []
for i in range(len(pred_label)):
    if pred_label[i] == 0:
        g0.append([puntos_test[i,0], puntos_test[i,1]])
    else:
        g1.append([puntos_test[i,0], puntos_test[i,1]])
        
grupo0 = np.matrix(g0)
grupo1 = np.matrix(g1)

# mostrar datos ya clasificados
plt.scatter([puntos[:int(n/2),0]], [puntos[:int(n/2),1]], c="b")
plt.scatter([puntos[int(n/2):,0]], [puntos[int(n/2):,1]], c="g")
plt.scatter([grupo0[:,0]], [grupo0[:,1]], c='b', marker='x')
plt.scatter([grupo1[:,0]], [grupo1[:,1]], c='g', marker='x')
display(plt)

Ahora que hemos clasificados los nuevos puntos los metemos en una matriz simplemente para representarlos en una gráfica, junto con el color definitivo que le corresponde según han sido clasificados por k-nn. De esta forma comprobamos visualmente si la clasificación ha ido según lo esperado.

Para ser un algoritmo que conceptualmente no es nada complejo parece que realiza un buen trabajo. Aunque claro, nuestro ejemplo no es muy impresionante ¿verdad? te preguntarás si todo esto que te he contado sirve realmente para algo útil. Una de las aplicaciones más interesantes de este algoritmo son los sistemas de recomendación. Eso que hace Netflix cuando te recomienda una película según tu historial de visionado o Amazon cuando te recomienda productos en base a tus compras anteriores. Pero también tiene otras utilidades, algunas de ellas bastante más espectaculares que la que acabamos de ver, pero te dejo con la intriga hasta el próximo artículo en el que espero impresionarte un poco más usando a nuestro nuevo amigo k-nn.

Breve introducción a las matrices con NumPy

noreply@blogger.com (Alberto García Serrano) — Sat, 16 Feb 2019 23:47:00 +0000

La librería NumPy traslada la potencia de entornos como Matlab u Octave al lenguaje de programación Python. Acometer un proyecto de Machine Learning o IA requiere un manejo intensivo y eficaz de cálculo matricial, y NumPy pone a nuestra disposición toda esa potencia. La libraría dispone de un montón de herramientas y capacidades que no es posible acometer en un artículo como este, sin embargo, me dispongo a contaros lo más importante en cuanto al manejo de matrices con NumPy.

Empecemos por un ejemplo sencillo: vamos a crear una matriz.

import numpy as np

matrix1 = np.matrix([[1, 2, 3, 4], [5, 6, 7, 8]])
print (matrix)
[[1 2 3 4]
 [5 6 7 8]]

La primera línea importa la librería para poder usarla en el programa y seguidamente se crea una matriz de dimensión 2x4. Alternativamente, podemos usar una sintaxis al estilo de Matlab.

matrix1 = np.matrix('1 2 3 4; 5 6 7 8')

Existen otras tres funciones que nos facilitan la creación de matrices: zeros() crea una matriz inicializada con el valor 0, ones() hace lo mismo pero con todos sus valores a 1. Finalmente, eye() crea la matriz identidad del tamaño que especifiquemos como parámetro.

zeros = np.zeros([2,3])
print (zeros)
[[ 0.  0.  0.]
 [ 0.  0.  0.]]

ones = np.ones([3, 2])
print (ones)
[[ 1.  1.]
 [ 1.  1.]
 [ 1.  1.]]

eye = np.eye(4)
print (eye)
[[ 1.  0.  0.  0.]
 [ 0.  1.  0.  0.]
 [ 0.  0.  1.  0.]
 [ 0.  0.  0.  1.]]

Podemos acceder a valores concretos de la matriz o a porciones completas. Vamos a verlo con algunos ejemplos:

matrix1 = np.matrix('1 2 3; 4 5 6; 7 8 9')
print (matrix1)
[[1 2 3]
 [4 5 6]
 [7 8 9]]

print (matrix1[0])
[[1 2 3]]

En este ejemplo, estamos accediendo a la primera fila (fila 0) de la matriz. Si queremos acceder a un sólo elemento podemos hacerlo así:

print (matrix1[0,1])
2

También podemos acceder a una porción de una fila:

print (matrix1[1, 1:3])
[[5 6]]

Y lo mismo es aplicable a las columnas:

print (matrix1[0:3, 0])
[[1]
 [4]
 [7]]

Las funciones size() y shape() nos devuelven el número de elementos y la dimensión de la matriz respectivamente.

print (matrix1.size)
9

print (matrix1.shape)
(3, 3)

Otra información útil que podemos obtener es el valor máximo, el mínimo y la media de los valores de la matriz:

print (matrix1.max())
9

print (matrix1.min())
1

print (matrix1.mean())
5.0

No sólo es posible acceder a los elementos de la matriz para consultarlos, también pueden ser modificados:

matrix1 = np.matrix('1 2; 3 4')
print (matrix1)
[[1 2]
 [3 4]]

matrix1[1, 0] = 10
print(matrix1)
[[ 1  2]
 [10  4]]

matrix1[0] = [6, 7]
print (matrix1)
[[ 6  7]
 [10  4]]

El primer ejemplo modifica un valor de la matriz y el segundo sustituye una fila entera. Como se puede observar, se sigue la misma notación que hemos usado para acceder a los valores de la matriz un poco más arriba. Otra función útil es flatten(), que convierte la matriz en un array de una dimensión con todos los elementos de la matriz.

print (matrix1.flatten())
[[ 6  7 10  4]]

Hasta ahora hemos visto como se crean las matrices, como se accede a ellas y como cambiamos sus valores. Abordamos ahora las operaciones matemáticas. Empezamos por una que se usa bastante, la matriz traspuesta.

print(matrix1.T) 
[[ 6 10]
 [ 7  4]]

Sobre las matrices se pueden aplicar las operaciones matemáticas básicas, tanto entre matrices como entre una matriz y un escalar. Los siguientes ejemplos ilustran la suma de una matriz y un escalar y la suma entre dos matrices.

matrix1 = np.matrix('1 2; 3 4')
matrix2 = np.matrix('5 6; 7 8')
print (matrix1 + 2)
[[3 4]
 [5 6]]

print (matrix1 + matrix2)
[[ 6  8]
 [10 12]]

Nos resta por ver la multiplicación de matrices, la más usada junto con la traspuesta.

matrix1 = np.matrix('1 2 3; 4 5 6; 7 8 9')
matrix2 = np.matrix('1 2; 3 4; 5 6')
print (matrix1 * matrix2)
[[ 22  28]
 [ 49  64]
 [ 76 100]]

NumPy tiene muchas más posibilidades, pero considero que esto es lo básico para empezar a trabajar con ella. En próximos artículos tendremos la oportunidad de hacer uso de ella en algún que otro ejemplo práctico.

Semáforos y más semáforos

noreply@blogger.com (Alberto García Serrano) — Mon, 14 Jan 2019 22:41:36 +0000

Que levante la mano quién no haya hecho el típico circuito de semáforo para probar una plataforma o un microcontrolador nuevo. El semáforo se podría considerar como el HolaMundo de la electrónica. Cualquier alumno de tecnología o de ingeniería podrá corroborar lo que digo.

Sobre este particular, José Antonio Vacas dio una interesante charla en Malakabot 2017 que os dejo aquí debajo.

Está claro que como ejemplo inicial, el semáforo puede ser un diseño útil, pero después hay que ir más allá. Si andáis faltos de imaginación, José Antonio os deja en su GitHub algunas ideas que os pueden servir de inspiración

Pero si pensabais que tras todo este movimiento se habían acabado los semáforos, estáis muy equivocados. Tras esto de los semáforos se encuentra un interesante reto propuesto por Cesar García del podcast La Hora Maker y David Cuartielles. Se trata de un proyecto colectivo para realizar semáforos en la mayor cantidad de modos posible. Dejo que ellos mismos os presenten el proyecto.

Si queréis participar sólo tenéis que acudir al GitHub de David y seguir las instrucciones. Yo por mi parte me he animado a enviar dos diseños: Uno con el TI MSP430 y otro con FreeRTOS.

El diseño del semáforo es muy simple y se muestra en el siguiente esquema, donde se indica a qué pines han de conectarse los LED. En este caso he usado una placa de desarrollo LaunchPad basada en el microcontrolador MSP430. También os dejo el código C del semáforo.

#include <msp430.h>
enum estados { ROJO, AMARILLO, VERDE };
int color;

int main(void)
{
	WDTCTL = WDTPW + WDTHOLD; 	// Detiene el watchdog
	P3DIR |= BIT0;            	// P3.0 (LED rojo) como salida digital
	P3DIR |= BIT1;            	// P3.1 (LED amarillo) como salida digital
	P2DIR |= BIT6;            	// P2.6 (LED verde) como salida digital
	P1REN |= BIT1;		  	// Resistencia en P1.1 (pulsador derecho)
	P3OUT |= BIT1;		  	// Modo pull-up
	P2OUT |= BIT1;
	P1IE  |= BIT1;			// habilita la interrupción de P1.1
	P1IES |= BIT1;			// Interrupción por flanco de bajada
	P1IFG &= ~(BIT1);		// baja flag de interrupción de P1.1 (por si acaso)
    color = ROJO;
	P3OUT |= BIT0; 			// Empezamos encendiendo el rojo
	P3OUT &= ~BIT1;			// apagamos el amarillo
	P2OUT &= ~BIT6;			// apagamos el verde
    // Activar interrupciones
	__bis_SR_register(GIE); 	// también puede usarse: __enable_interrupt()
    // Modo bajo consumo pero atendemos interrupciones
	__low_power_mode_0();		// Entrada LPM0, habilita interrupciones
	__no_operation();         	// Para depuración
    }
    
// Rutina de gestión de interrupciones de P1
#pragma vector=PORT1_VECTOR
__interrupt
void P1_ISR(void) {
	extern int color;
	if (P1IFG & BIT1) {		// ¿fuente de la interrupción es P1.1 (pulsador)?
		switch (color) {
			case ROJO:
				color = AMARILLO;
				P3OUT &= ~BIT0;	// apagamos el rojo
				P3OUT |= BIT1;	// encendemos el amarillo
				break;
			case AMARILLO:
				color = VERDE;
				P3OUT &= ~BIT1;	// apagamos el amarillo
				P2OUT |= BIT6;	// encendemos el verde
				break;
			case VERDE:
				color = ROJO;
				P2OUT &= ~BIT6;	// apagamos el verde
				P3OUT |= BIT0;	// encendemos el rojo
				break;
		}
		P1IFG &= ~BIT1;		// puesta a cero del flag de interrupción P1.1
	}
}

Para el ejemplo de semáforo con FreeRTOS he usado la placa de desarrollo LaunchPad CC3200 que utiliza una MCU ARM Cortex A4. El esquema de conexionado con la placa y el código son los siguientes.

#include <stdint.h>
#include <stdbool.h>
#include <stdlib.h>
#include "hw_types.h"
#include "hw_ints.h"
#include "hw_memmap.h"
#include "interrupt.h"
#include "rom_map.h"
#include "prcm.h"
#include "gpio.h"
#include "pinmux.h"
#include "FreeRTOS.h"
#include "task.h"

// prototipos de función
void BoardInit(void);

// variables globales
extern void (* const g_pfnVectors[])(void);

//*****************************************************************************
// FreeRTOS User Hook Functions enabled in FreeRTOSConfig.h
//*****************************************************************************
void
vAssertCalled( const char *pcFile, unsigned long ulLine )
{
    //Handle Assert here
    while(1)
    {
    }
}

void vApplicationMallocFailedHook()
{
    //Handle Memory Allocation Errors
    while(1)
    {
    }
}

void vApplicationStackOverflowHook( void *pxTask,
                                   signed char *pcTaskName)
{
    //Handle FreeRTOS Stack Overflow
    while(1)
    {
    }
}

// Board Initialization & Configuration
void BoardInit(void)
{
// Set vector table base
    IntVTableBaseSet((unsigned long)&g_pfnVectors[0]);
    // Enable Processor
    MAP_IntMasterEnable();
    MAP_IntEnable(FAULT_SYSTICK);
    PRCMCC3200MCUInit();
}
  
/*-----------------------------------------------------------*/
//Esta tarea enceinde las luces del semáforo
static void SemaforoBlink( void *pvParameters )
{
//La tarea debe tener un bucle infinito...
	for( ;; )
	{
		// Encender LED rojo
		GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_3,GPIO_PIN_3);
		vTaskDelay(0.5*configTICK_RATE_HZ);
		GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_3,0);
        // Encender LED amarillo
		GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_2,GPIO_PIN_2);
		vTaskDelay(0.5*configTICK_RATE_HZ);
		GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_2,0);
        // Encender LED verde
		GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_1,GPIO_PIN_1);
		vTaskDelay(0.5*configTICK_RATE_HZ);
		GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_1,0);
	}
}
/*-----------------------------------------------------------*/

// Función main
void main()
{ 
	BoardInit();	// Configuración de la placa
    PinMuxConfig();	// Pinmux donde se configuras los puertos para los LED como salida digital
    // Apagamos los LED
    GPIOPinWrite(GPIOA1_BASE,GPIO_PIN_3|GPIO_PIN_2|GPIO_PIN_1,0);
    
	/* creamos la tarea para el semaforo */
	if (xTaskCreate( SemaforoBlink, "Semaforo", configMINIMAL_STACK_SIZE, NULL, tskIDLE_PRIORITY+1, NULL )!=pdPASS)
		while (1);
        vTaskStartScheduler(); // Inicia el planificador de tareas
}

Os animo a que participéis, ya que es un proyecto, además de interesante, bastante divertido.

Selección de atributos relevantes usando la entropía de Shannon

noreply@blogger.com (Alberto García Serrano) — Wed, 19 Sep 2018 10:29:00 +0000

Hoy en día, la detección de amenazas está cada vez más ligada al aprendizaje automático. Si queremos que un IDS o un antivirus pueda lidiar con amenazas nuevas, además de métodos heurísticos hay que recurrir técnicas de clasificación, una parte importante de la IA. En los problemas de clasificación, una de las tareas a la que nos enfrentamos es seleccionar aquellos atributos que mejor describen la variable objetivo. Esto es, elegir aquellos atributos que aportan más información y tienen mayor correlación con la variable explicada. Por ejemplo, si queremos establecer la probabilidad de que alguien padezca cáncer de pulmón o no, una variable que va a ser muy significativa es si es fumadora. En el otro extremo, el color de los ojos de la persona, por poner un ejemplo muy claro, no aporta información relevante.

Por desgracia, no siempre es tan obvia la cantidad de información que aporta un atributo, así que vamos a ver un método relativamente sencillo de medirlo.

En 1948, Claude Shannon -la fotografía de arriba es un grafitti con su retrato- publicó un artículo llamado A Mathematical Theory of Communication donde introducía el concepto de entropía en el contexto de la teoría de la información (de la que Shannon es padre). Simplificando mucho, la entropía es una medida del desorden de un conjunto datos, y se define como:

\[H(X) = -\sum_{i}^{n} p_i \times \log_2(p_i)\]

Donde p_i es la probabilidad relativa de aparición de la propiedad i en el conjunto de datos. Por ejemplo, p₁ podría ser la probabilidad de tener cáncer de pulmón y p₂ la probabilidad de no tenerlo.

Para entenderlo mejor vamos a ir viéndolo con un ejemplo práctico en R. Vamos a usar el conjunto de datos Titanic incluido en R. Carguemos el conjunto de datos y veamos cuáles son sus atributos:

> data(Titanic)
> dimnames(Titanic)
$Class
[1] "1st"  "2nd"  "3rd"  "Crew"

$Sex
[1] "Male"   "Female"

$Age
[1] "Child" "Adult"

$Survived
[1] "No"  "Yes"

> margin.table(Titanic)
[1] 2201

Así pues, tenemos un conjunto de datos de 2201 pasajeros del Titanic que dispone de tres atributos o variables independientes que son:
- La clase en la que viajaba cada persona que puede tomar los valores: primera, segunda, tercera y tripulación.
- El sexo: Hombre o mujer.
- La edad: Niño o adulto.

Tenemos una cuarta variable, que en este caso será nuestra variable explicada o dependiente, que nos dice si la persona sobrevivió o no al hundimiento del barco.
Así pues, nuestro objetivo es encontrar alguna relación entre los atributos antes mencionados y la probabilidad de que un pasajero sobreviva o no al naufragio.

Intuitivamente podemos pensar que el índice de mujeres y niños supervivientes debe ser mayor que el de los hombres por aquello de "las mujeres y los niños primero". Vemos si es cierto.

> margin.table(Titanic, c(2,4))
Survived
Sex        No  Yes
Male   1364  367
Female  126  344

Parece que, en efecto, las mujeres tenían más tendencia a sobrevivir que los hombres. Veámoslo gráficamente.

A ver los niños:

> margin.table(Titanic,c(3,4))
       Survived
Age       No  Yes
  Child   52   57
  Adult 1438  654

Parece que aquí la cosa ya no está tan clara. ¿Y qué ocurre con la clase en la que viajaba el pasajero?

> margin.table(Titanic,c(1,4))
      Survived
Class   No Yes
  1st  122 203
  2nd  167 118
  3rd  528 178
  Crew 673 212

Se observa también cierta correlación (a los de primera parece que les fue algo mejor) pero no tan clara como con el sexo.

Vamos a tratar de medir, usando el concepto de entropía de Shannon, qué cantidad de información aporta cada una de las tres variables independientes.

Calculamos la entropía del conjunto de datos respecto a la variable explicada Survived.
p₁ -> No sobrevive.
p₂ -> Sobrevive.

Comenzamos calculando las probabilidades de p₁ y p₂.

> p<-margin.table(Titanic, c(4))
> p.norm<-p/sum(p)
> p.norm
Survived
No      Yes 
0.676965 0.323035

Así pues, tenemos que:
p₁ = 0.676965
p₂ = 0.323035

Con las probabilidades ya calculadas, nos resta hacer el sumatorio según la fórmula de Shannon.

> entropia<-(-sum(log2(p.norm)*p.norm))
> entropia
[1] 0.9076514

Lo que nos dice este valor es la medida de desorden de la variable explicada con un valor que va de 0 a 1. Un valor de 0 indica orden total, o lo que es lo mismo, o todos sobrevivieron o todos murieron, pero todos los valores son iguales. Cuanto más cerca esté de 1, mayor desorden. Un valor 1 querría decir que el 50% de pasajeros murieron y el otro 50% sobrevivió (hay que tener en cuenta que la entropía no es una función lineal).
Usaremos la notación H(survived) para nombrar la entropía producida por la variable survived. También vamos a definir el concepto de entropía condicionada a la entropía generada fijando, a priori, el valor de una segunda variable. Por ejemplo, H(survived | Sex) es la entropía de la variable survived condicionada al valor de la variable sexo. Se calcula con:

\[H(X|Y) = \sum_{i}^{n} p_x \times H(X|Y=y)\]

P_x es la probabilidad relativa de aparición de la propiedad x en el conjunto de los datos en los que Y=y y H(X|Y=y) es la entropía de la propiedad y en el conjunto de los datos en los que Y=y.

Por ejemplo, tomemos el caso del sexo y calculemos H(Survived|Sex=Male) y H(Survived|Sex=Female).

La entropía para el caso de las mujeres la calculamos de la forma que ya conocemos.

> p<-margin.table(Titanic,c(2,4))
> p
Survived
Sex        No  Yes
Male   1364  367
Female  126  344

> p.male<-p[1,]
> p.female<-p[2,]

> p.female.norm<-p.female/sum(p.female)
> p.female.norm
       No       Yes 
0.2680851 0.7319149 
> entropia.female<-(-sum(log2(p.female.norm)*p.female.norm))
> entropia.female
[1] 0.8387034

Ahora calculamos la entropía para el caso de los hombres.

> p.male.norm<-p.male/sum(p.male)
> p.male.norm
       No       Yes 
0.7879838 0.2120162 
> entropia.male<-(-sum(log2(p.male.norm)*p.male.norm))
> entropia.male
[1] 0.745319

Finalmente, calculamos la entropía H(Survived|Sex) aplicando la fórmula de la entropía condicionada que acabamos de ver.

> p.sex<-margin.table(Titanic,2)
> p.norm.sex<-p.sex/sum(p.sex)
> p.norm.sex
Sex
     Male    Female 
0.7864607 0.2135393 
> entropia.sex<-entropia.male*as.numeric(p.norm.sex[1])+
   entropia.female*as.numeric(p.norm.sex[2])
> entropia.sex 
[1] 0.7652602

Así pues, ya podemos calcular la entropía, pero este valor por sí mismo tampoco nos aclara del todo el asunto de qué atributo aporta más información.

Vamos a definir, basándonos en la entropía, un nuevo concepto al que llamaremos ganancia de información como:

\[GI(X|Y) = H(X) - H(X|Y)\]

En este caso H(X) es la entropía de la variable explicada y H(X|Y) la entropía condicionada que acabamos de calcular.

La ganancia de información nos dice cómo se reduce la entropía cuando añadimos la nueva variable. Una mayor reducción implica una correlación mayor con la variable explicada. Calculémosla con los datos que ya tenemos:

> gi.sex=entropia-entropia.sex
> gi.sex
[1] 0.1423912

Analicemos qué ocurre con los otros dos atributos. Hacemos el cálculo para el atributo edad:

> p<-margin.table(Titanic,c(3,4))
> p.child<-p[1,]
> p.adult<-p[2,]
> p.child.norm<-p.child/sum(p.child)
> p.adult.norm<-p.adult/sum(p.adult)
> entropia.child<-(-sum(log2(p.child.norm)*p.child.norm))
> entropia.adult<-(-sum(log2(p.adult.norm)*p.adult.norm))
> p.age<-margin.table(Titanic,3)
> p.norm.age<-p.age/sum(p.age)
> entropia.age<-entropia.child*as.numeric(p.norm.age[1])+
   entropia.adult*as.numeric(p.norm.age[2])
> gi.age<-entropia-entropia.age
> gi.age
[1] 0.006410718

Y finalmente hacemos el cálculo para el atributo clase:

> p<-margin.table(Titanic,c(1,4))
> p.1st<-p[1,]
> p.2nd<-p[2,]
> p.3rd<-p[3,]
> p.crew<-p[4,]
> p.1st.norm<-p.1st/sum(p.1st)
> p.2nd.norm<-p.2nd/sum(p.2nd)
> p.3rd.norm<-p.3rd/sum(p.3rd)
> p.crew.norm<-p.crew/sum(p.crew)
> entropia.1st<-(-sum(log2(p.1st.norm)*p.1st.norm))
> entropia.2nd<-(-sum(log2(p.2nd.norm)*p.2nd.norm))
> entropia.3rd<-(-sum(log2(p.3rd.norm)*p.3rd.norm))
> entropia.crew<-(-sum(log2(p.crew.norm)*p.crew.norm))
> p.class<-margin.table(Titanic,1)
> p.norm.class<-p.class/sum(p.class)
> entropia.class<-entropia.1st*as.numeric(p.norm.class[1])+
   entropia.2nd*as.numeric(p.norm.class[2])+
   entropia.3rd*as.numeric(p.norm.class[3])+
   entropia.crew*as.numeric(p.norm.class[4])
> gi.class<-entropia-entropia.class
> gi.class
[1] 0.05928794

En resumen:
Ganancia de información para el sexo: 0.1423912
Ganancia de información para la edad: 0.006410718
Ganancia de información para la clase: 0.05928794

Esto significa que el atributo que más información aporta, con diferencia, es el sexo. La clase aporta bastante menos información, y la edad apenas nada.
En la práctica, esto significa que al construir nuestro modelo predictivo, podríamos descartar tranquilamente el atributo edad sin miedo a peder capacidad predictiva.

Me ha parecido más didáctico desarrollar manualmente los cálculos, sin embargo R ya dispone de un paquete para el cálculo de la entropía de Shannon que nos facilitará sin duda la vida.

Juegos con matrices de LEDs y Arduino

noreply@blogger.com (Alberto García Serrano) — Sun, 23 Jul 2017 18:07:03 +0000

En una entrada anterior ya utilizamos una matriz de LEDs junto con el integrado max7219. Una de las virtudes de este integrado es que facilita la conexión de varias matrices usando un bus serie SPI, así como otros tipos de displays, como los de siete segmentos. En este artículo vamos a utilizar dos matrices para construir un pequeño juego tipo Pong, Arkanoid o, en fin, ya sabes a lo que me refiero, una raqueta que controla el jugador y una pelotita que hay que evitar que se salga por la parte inferior.

He montado el circuito sobre una protoboard usando un Arduino Nano, dos matrices de LEDs y un potenciómetro conectado a una entrada analógica del Arduino para controlar la raqueta. El esquema del circuito es el siguiente.

Como ves, el circuito no tiene nada de particular. Los pines D8, D9 y D10 de Arduino se conectan a los pines CS, CLK y DIN del bus SPI de la matriz de LEDs. Además, obviamente, de conectar VCC y GND. El integrado max7219 nos provee, en la parte superior de la matriz de LEDS, de las conexiones necesarias para conectar una segunda matriz (y una tercera, cuarta, etc.). También he aprovechado las salidas VCC y GND de la segunda matriz para conectarlos al potenciómetro y así dejar más limpia de cables la protoboard. Finalmente, he conectado el potenciómetro a la entrada analógica A3 de Arduino. La tensión de entrada en el pin se convierte en digital mediante un ADC de 10 bits, así que tenemos un valor entre 0 y 1023, que usando la función map() de la librería de Arduino, adaptaremos a nuestras necesidades.

El código para Arduino es el siguiente.

#include "LedControl.h" 
#define WIDTH 16
#define HEIGHT 8
#define PADLENGTH 3
 
// Parámetros: pines DIN, CLK, CS y 2 (2 displays)
LedControl lc=LedControl(8,10,9,2);
const int analog_in = A3; // entrada analogica potenciometro
int val, pos; 
byte pantalla[8][16]; // pantalla
int x, y; // posición de la bola
int dx, dy; // delta
int gameover = 0;

void plot(int x, int y)
{
  int screen;
  
  if (x>=0 && x==0 && y<HEIGHT) { 
    if (x > ((WIDTH/2)-1)) {
      x = x - HEIGHT;
      screen = 1;
    } else {
      screen = 0;
    }
    lc.setLed(screen, x, y, 1);
  }
}

void draw_screen(int x, int y)
{
  int i;     
  // draw PAD
  for (i=0; i<PADLENGTH; i++) {
    plot(pos+i, 0);
  }
  // draw ball
  plot(x, y);
  
}

void move()
{
  x = x + dx;
  y = y + dy;
  if (x<1 || x>=WIDTH-1) {
    dx = -dx;
  } 
  
  if (y>=HEIGHT-1) {
    dy = -dy;
  }
  
  if (y==1 && x>=pos && x<=pos+PADLENGTH) {
    // rebota con la raqueta
    dy = -dy;
  } 
  
  if (y < 0) {
    // pelota sale por abajo. Game over
    gameover = 1;  
  }
}
 
void setup() {
  lc.shutdown(0,false);    // Iniciamos la matriz led
  lc.setIntensity(0,5);    // Intensidad de los leds
  lc.clearDisplay(0);      // Apagamos los leds
  lc.shutdown(1,false);    // Iniciamos la matriz led
  lc.setIntensity(1,5);    // Intensidad de los leds
  lc.clearDisplay(1);      // Apagamos los leds
  x = y = 0;
  dx = dy = 1;
  pos = 0;
  Serial.begin(9600);
}
 
void loop() {      
  int i;
  
  while(true) {   
    if (!gameover) {
      // leer entrada
      val = analogRead(analog_in);
      pos = map(val, 0, 1023, WIDTH - PADLENGTH, -1);
      Serial.println(pos);
  
      // actualizar matriz de LEDs
      lc.clearDisplay(0);
      lc.clearDisplay(1);
      draw_screen(x,y);
      move();
      delay(200);
    } else {
      // game over
      for (i=0; i<HEIGHT; i++) {
        lc.setColumn(0,i,B11111111);
        lc.setColumn(1,HEIGHT-i-1,B11111111);
        delay(200);  
      }
      x = y = 0;
      dx = dy = 1;
      gameover = 0;
      
    }
  }
}

La función plot() se encarga de activar un pixel en la pantalla "virtual" de 16x8 (las dos matrices juntas). Es decir, calcula en cuál de las dos matrices está el LED que se quiere activar y cuál es el LED físico que le corresponde dentro de esa matriz. La función draw_screen() dibuja la raqueta y la pelota en cada iteración del bucle principal y, finalmente, la función move() calcula la nueva posición de la pelota, verifica si ha rebotado o si ha tocado la base, en cuyo caso, hemos perdido y se muestra una pequeña animación en las matrices.

Aquí os dejo un vídeo con el resultado final.

Máquinas de estado finito en VHDL

noreply@blogger.com (Alberto García Serrano) — Sat, 25 Mar 2017 18:06:54 +0000

Las máquinas de estado finito, más conocidas por su acrónimo en inglés FSM (Finite State Machine), se utilizan ampliamente en el diseño de circuitos digitales (además de en otros ámbitos de la ingeniería, como la programación), para describir el comportamiento de un sistema según el valor de sus entradas y de cómo van cambiando en el tiempo. Ésta es una definición parcial pero que nos permite hacernos una primera idea intuitiva.

Desde el punto de vista de las FSM, un sistema está compuesto de estados por los que va pasando el sistema, de señales de entrada que modifican esos estados y de señales de salida que pueden utilizarse para conocer el estado del sistema y actuar en consecuencia. Un ejemplo muy visual podría ser un semáforo, el cuál dispone de tres estados diferentes, uno para cada color. Las entradas del sistema las podría generar un temporizador que activa una señal cada cierto tiempo, indicando que hay que pasar al siguiente estado. Por último, las salidas del sistema podrían ser tres señales que indiquen qué lámpara, de las tres disponibles, tiene que encenderse. Para representar una máquina de estados se utilizan diagramas de estados como el siguiente.

Cada círculo representa un estado. En este caso tenemos cuatro estados que se llaman S1, S2, S3 y S4. El estado S1 tiene una flecha que indica que es el estado el inicial, al que se entra tras un reset del sistema. El S4 tiene un doble círculo, que indica que es un estado final. Los estados están unidos por una flechas llamadas transiciones, que indican cómo evoluciona el sistema de un estado a otro según se activan las señales de entrada. En este caso concreto hay tres señales de entrada (a, b, c y d). Si estando en el estado S1 se activa la señal a, el sistema evolucionará al estado S2. Del mismo modo, si estando en el estado S2 se activa la señal b, el sistema evolucionará de nuevo al estado S1. En este caso no hay señales de salida. Si las hay pueden ponerse dentro del círculo o, para no sobrecargar el diagrama, se pueden definir fuera de éste, por ejemplo en una tabla. Vamos a crear una FSM muy básica y a implementarla en VHDL, para poder llevarla a una FPGA. La máquina que vamos a implementar es la mostrada en el siguiente diagrama de estados.

Tenemos cuatro estados (S0, S1, S2 y S3) y dos señales de entrada (A y B). No hemos puesto en el diagrama las señales de salida, pero las definimos en la siguiente tabla:

S0 -> Y0

S1 -> Y1

S2 -> Y2

S3 -> Y3

Por lo tanto, tenemos cuatro señales de salida. La señal Y0 se activa cuando el sistema está en el estado S0, La Y1 cuando el sistema está en el estado S1 y así sucesivamente. Como se puede observar, la salida del sistema depende exclusivamente del estado en el que se encuentra el sistema. Cuando ocurre esto decimos que se trata de una máquina de estados de tipo Moore. Si la salida depende del estado y además de las entradas actuales del sistema, se trataría de una máquina de estados de tipo Mealy. Se puede demostrar que ambos modelos son equivalentes, así que aquí vamos a quedarnos con el primer tipo que es más sencillo.

A nivel físico, una de las formas más efectivas de implementar una FSM es según el siguiente esquema.

Obviamente necesitamos algún elemento de memoria para almacenar el estado actual del sistema. Nosotros usaremos un registro de biestables de tipo D (es por lo tanto un bloque secuencial). Los otros dos bloques son puramente combinacionales. El primero se encarga de generar el estado siguiente a partir del estado actual y de las entradas del sistema. Además, necesitamos otro circuito para generar las señales de salida a partir del estado actual (recordemos que es una máquina de tipo Moore).

Voy a llevar a la práctica el circuito usando una FPGA. En concreto usaré la placa BASYS 3 con una FPGA Artix 7 de Xilinx. Voy a utilizar el botón 0 como señal de reset y los botones 4 y 2 como señales A y B respectivamente. Las salidas Y0 a Y3 las voy a asociar a los LEDs 0 a 3. Para evitar utilizar circuitos antirebote para los pulsadores, he diseñado el diagrama de forma que dos pulsaciones seguidas del mismo botón no hagan evolucionar el sistema a otro estado.

Como tenemos cuatro estados, podríamos usar dos bit para codificar los estados. Por ejemplo: S0->00, S1->01, S2->10 y S3->11. Sin embargo, lo más habitual es dejar que el sintetizador lo haga por nosotros definiendo un tipo para los estados tal y como se ve en la siguiente definición. Así, si añadimos nuevos estados, no tenemos que andar redimensionando el registro de estados.

    -- declaraciones modelo FSM
    type STATES is (S0, S1, S2, S3);
    signal state_reg, state_next: STATES;

Veamos como se implementa en VHDL cada uno de estos tres bloques. Empecemos por el registro de estados.

    -- registro de estados
    process(CLK)
    begin
        if CLK'event and CLK='1' then
            if RST='1' then 
                state_reg <= s0;
            else 
                state_reg <= state_next;
            end if;
         end if;
    end process;

El registro de estados tiene una señal de reset síncrona que lleva al sistema al estado S0. Si la lógica de codificación del estado siguiente genera un nuevo estado, este se almacenará en el registro, si no, la salida del registro se mantiene. Veamos cómo es el circuito que genera el estado siguiente.

    -- Lógica de estado siguiente (circuito combinacional)
    process (state_reg, A, B)
    begin
        state_next <= state_reg; 
        case state_reg is 
        when S0 =>
            if A='1' then
                state_next <= S1; 
            end if; 
        when S1 =>
            if B='1' then
                state_next <= S2; 
            end if; 
        when S2 => 
            if A='1' then
                state_next <= S3; 
            end if; 
        when S3 =>
            state_next <= S3;
        end case;
    end process;

El estado siguiente se genera a partir del estado actual y de las entradas A y B, por lo que ponemos estas tres señales en la lista de sensibilidad del proceso. Para evitar que el sintetizador genere algún latch o elemento de memoria, asignamos por defecto el estado actual a la señal de estado siguiente. El código VHDL describe los estados y las transiciones asociadas a loas entradas según el diagrama de estados que definimos más arriba, de forma que, de ser necesario, asigna el valor del estado siguiente a la variable state_next cuando hay un cambio en las entradas o en el estado. Por último echamos un vistazo al circuito de decodificación para las salidas.

    -- salida tipo Moore
    process (state_reg)
    begin
        -- estableciendo la salida por defecto
        -- nos aseguramos de crear un circuito
        -- combinacional sin latches.
        Y0 <= '0';
        Y1 <= '0';
        Y2 <= '0';
        Y3 <= '0'; 
        case state_reg is 
        when S0 => Y0 <= '1'; 
        when S1 => Y1 <= '1'; 
        when S2 => Y2 <= '1'; 
        when S3 => Y3 <= '1';
        end case; 
    end process;

De nuevo, para evitar que el sintetizador cree elementos de memoria, asignamos un valor por defecto a las salidas (en este caso 0), y según el estado actual, se activa la salida correspondiente según definimos en la tabla de arriba.

El código VHDL completo de nuestra FSM es el siguiente.

library IEEE;
use IEEE.STD_LOGIC_1164.ALL;

entity FSM is
    Port ( CLK : in STD_LOGIC;
           BTN : in STD_LOGIC_VECTOR (4 downto 0);
           LED : out STD_LOGIC_VECTOR (15 downto 0));
end FSM;
architecture Behavioral of FSM is
    -- alias
    alias RST : STD_LOGIC is BTN(0);
    alias A : STD_LOGIC is BTN(4);
    alias B : STD_LOGIC is BTN(2);
    alias Y0 : STD_LOGIC is LED(0);
    alias Y1 : STD_LOGIC is LED(1);
    alias Y2 : STD_LOGIC is LED(2);
    alias Y3 : STD_LOGIC is LED(3);
    
    -- declaraciones modelo FSM
    type STATES is (S0, S1, S2, S3);
    signal state_reg, state_next: STATES;
    begin
    -- registro de estados
    process(CLK)
    begin
        if CLK'event and CLK='1' then
            if RST='1' then 
                state_reg <= s0;
            else 
                state_reg <= state_next;
            end if;
         end if;
    end process;
    
    -- Lógica de estado siguiente (circuito combinacional)
    process (state_reg, A, B)
    begin
        state_next <= state_reg; 
        case state_reg is 
        when S0 =>
            if A='1' then
                state_next <= S1; 
            end if; 
        when S1 =>
            if B='1' then
                state_next <= S2; 
            end if; 
        when S2 => 
            if A='1' then
                state_next <= S3; 
            end if; 
        when S3 =>
            state_next <= S3;
        end case;
    end process;
    
    -- salida tipo Moore
    process (state_reg)
    begin
        -- estableciendo la salida por defecto
        -- nos aseguramos de crear un circuito
        -- combinacional sin latches.
        Y0 <= '0';
        Y1 <= '0';
        Y2 <= '0';
        Y3 <= '0'; 
        case state_reg is 
        when S0 => Y0 <= '1'; 
        when S1 => Y1 <= '1'; 
        when S2 => Y2 <= '1'; 
        when S3 => Y3 <= '1';
        end case; 
    end process;    
    end Behavioral;

Os dejo el archivo de restricciones para la placa BASYS 3 que he usado para este ejemplo, donde se definen los puertos de los pulsadores, los LEDs y el reloj del sistema.

## LEDs
set_property PACKAGE_PIN U16 [get_ports {LED[0]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[0]}]
set_property PACKAGE_PIN E19 [get_ports {LED[1]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[1]}]
set_property PACKAGE_PIN U19 [get_ports {LED[2]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[2]}]
set_property PACKAGE_PIN V19 [get_ports {LED[3]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[3]}]
set_property PACKAGE_PIN W18 [get_ports {LED[4]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[4]}]
set_property PACKAGE_PIN U15 [get_ports {LED[5]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[5]}]
set_property PACKAGE_PIN U14 [get_ports {LED[6]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[6]}]
set_property PACKAGE_PIN V14 [get_ports {LED[7]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[7]}]
set_property PACKAGE_PIN V13 [get_ports {LED[8]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[8]}]
set_property PACKAGE_PIN V3 [get_ports {LED[9]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[9]}]
set_property PACKAGE_PIN W3 [get_ports {LED[10]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[10]}]
set_property PACKAGE_PIN U3 [get_ports {LED[11]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[11]}]
set_property PACKAGE_PIN P3 [get_ports {LED[12]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[12]}]
set_property PACKAGE_PIN N3 [get_ports {LED[13]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[13]}]
set_property PACKAGE_PIN P1 [get_ports {LED[14]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[14]}]
set_property PACKAGE_PIN L1 [get_ports {LED[15]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[15]}]

## Botones
set_property IOSTANDARD LVCMOS33 [get_ports {BTN[4]}]
set_property IOSTANDARD LVCMOS33 [get_ports {BTN[3]}]
set_property IOSTANDARD LVCMOS33 [get_ports {BTN[2]}]
set_property IOSTANDARD LVCMOS33 [get_ports {BTN[1]}]
set_property IOSTANDARD LVCMOS33 [get_ports {BTN[0]}]
## BTNL -> W19
set_property PACKAGE_PIN W19 [get_ports {BTN[4]}] 
## BTND -> U17
set_property PACKAGE_PIN U17 [get_ports {BTN[3]}] 
## BTNR -> T17
set_property PACKAGE_PIN T17 [get_ports {BTN[2]}] 
## BTNU -> T18
set_property PACKAGE_PIN T18 [get_ports {BTN[1]}] 
## BTNC -> U18
set_property PACKAGE_PIN U18 [get_ports {BTN[0]}] 

## Reloj CLK conectado a W5
set_property IOSTANDARD LVCMOS33 [get_ports {CLK}]
set_property PACKAGE_PIN W5 [get_ports {CLK}]
create_clock -add -name sys_clk_pin -period 10.00 -waveform {0 5} [get_ports CLK]

En la simulación de la FSM vemos como evoluciona el circuito desde S0 a S3 según van cambiando las entradas A y B. También podemos observar cómo se activan las entradas Y0 a Y3.

https://youtu.be/doOEfAsoEmg

Relojes MMCM en FPGAs

noreply@blogger.com (Alberto García Serrano) — Sat, 18 Feb 2017 16:55:17 +0000

Cuando se diseñan sistemas digitales síncronos, una de las reglas de oro es usar un sólo dominio de reloj para todo el circuito. Esto quiere decir que sólo usaremos una única señal de reloj conectada a todos los elementos síncronos del circuito. Por lo tanto debemos evitar una práctica bastante habitual que es usar la salida de un elemento del sistema, como por ejemplo un contador, para obtener la señal de reloj de otro elemento cualquiera.

De hecho, aquí lo hicimos (mal) cuando diseñamos nuestro reloj digital (y usamos la salida del contador de 1Hz como señal de reloj para el contador de segundos). Aunque siempre en aras de la simplicidad y la pedagogía :-)

Así que si queremos que nuestro diseño no tenga problemas de sincronización, el mismo reloj debe ser usado por todos y cada uno de los elementos síncronos. Sin embargo, si nuestro circuito es complejo y sobre todo, si ocupa una gran superficie de la FPGA, la señal de reloj puede sufrir pequeños retardos al distribuirse entre elementos distantes dentro de la FPGA y puede haber un pequeño desfase entre los flancos activos. Para minimizar al máximo este problema, las FPGAs disponen de bloques de reloj dedicados que nos ofrecen módulos PLL y MMCM. Los MMCM (Mixed-Mode Clock Manager) usan conexiones especiales para repartir la señal de reloj por toda la superficie de la FPGA. Vamos a crear un diseño bastante simple en el que veremos cómo usar el reloj MMCM y de paso aplicaremos buenas prácticas para diseñar un sistema digital con un sólo dominio de reloj. En este caso usaremos una FPGA Artix 7 de XILINX (en la placa Basys 3 de Digilent) y el entorno de desarrollo Vivado. El resultado será un contador binario de 8 bits que se mostrará en los LEDs de la placa, con una velocidad de 1Hz (una cuenta por segundo).

Utilizando el módulo MMCM

Para empezar, crearemos un nuevo proyecto en Vivado y también nuestro reloj MMCM. Vamos a utilizar un core IP del catálogo de Vivado, así que pulsamos sobre IP Catalog en el Flow Navigation de la izquierda del GUI. Hacemos doble click en Vivado Repository > FPGA Features and Design > Clocking > Clocking Wizard.

Seleccionamos la primitiva MMCM. en la siguiente imagen vemos qué opciones deben estar marcadas.

En la pestaña Output Clocks, en el campo Output Freq (MHz) Requested ponemos el valor 50.00 (es decir, usaremos un reloj de 50MHz). Marcamos las casillas Reset y Locked.

Pulsamos OK. En la siguiente ventana marcamos Out of context per IP y pulsamos generate.

Cuando aparezca la ventana indicando que se ha generado el core, pulsamos OK.

En el proyecto, añadimos el archivos de restricciones siguiente, en el que se han configurado los pines de la FPGA que vamos a usar (en este caso, el fichero se llama config_GPIO.xdc).

## LEDs
set_property PACKAGE_PIN U16 [get_ports {LED[0]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[0]}]
set_property PACKAGE_PIN E19 [get_ports {LED[1]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[1]}]
set_property PACKAGE_PIN U19 [get_ports {LED[2]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[2]}]
set_property PACKAGE_PIN V19 [get_ports {LED[3]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[3]}]
set_property PACKAGE_PIN W18 [get_ports {LED[4]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[4]}]
set_property PACKAGE_PIN U15 [get_ports {LED[5]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[5]}]
set_property PACKAGE_PIN U14 [get_ports {LED[6]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[6]}]
set_property PACKAGE_PIN V14 [get_ports {LED[7]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[7]}]
set_property PACKAGE_PIN V13 [get_ports {LED[8]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[8]}]
set_property PACKAGE_PIN V3 [get_ports {LED[9]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[9]}]
set_property PACKAGE_PIN W3 [get_ports {LED[10]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[10]}]
set_property PACKAGE_PIN U3 [get_ports {LED[11]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[11]}]
set_property PACKAGE_PIN P3 [get_ports {LED[12]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[12]}]
set_property PACKAGE_PIN N3 [get_ports {LED[13]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[13]}]
set_property PACKAGE_PIN P1 [get_ports {LED[14]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[14]}]
set_property PACKAGE_PIN L1 [get_ports {LED[15]}]
set_property IOSTANDARD LVCMOS33 [get_ports {LED[15]}]
##7 segment display
## segmento A -> W7
set_property PACKAGE_PIN W7 [get_ports {CAT[0]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[0]}]
## segmento B -> W6
set_property PACKAGE_PIN W6 [get_ports {CAT[1]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[1]}]
## segmento C -> U8
set_property PACKAGE_PIN U8 [get_ports {CAT[2]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[2]}]
## segmento D -> V8
set_property PACKAGE_PIN V8 [get_ports {CAT[3]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[3]}]
## segmento E -> U5
set_property PACKAGE_PIN U5 [get_ports {CAT[4]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[4]}]
## segmento F -> V5
set_property PACKAGE_PIN V5 [get_ports {CAT[5]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[5]}]
## segmento G -> U7
set_property PACKAGE_PIN U7 [get_ports {CAT[6]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[6]}]
## segmento DP -> V7
set_property PACKAGE_PIN V7 [get_ports {CAT[7]}]
set_property IOSTANDARD LVCMOS33 [get_ports {CAT[7]}]
## Anodo 0 diplay -> U2
set_property PACKAGE_PIN U2 [get_ports {AN[0]}]
set_property IOSTANDARD LVCMOS33 [get_ports {AN[0]}]
## Anodo 1 diplay -> U4
set_property PACKAGE_PIN U4 [get_ports {AN[1]}]
set_property IOSTANDARD LVCMOS33 [get_ports {AN[1]}]
## Anodo 2 diplay -> V4
set_property PACKAGE_PIN V4 [get_ports {AN[2]}]
set_property IOSTANDARD LVCMOS33 [get_ports {AN[2]}]
## Anodo 3 diplay -> W4
set_property PACKAGE_PIN W4 [get_ports {AN[3]}]
set_property IOSTANDARD LVCMOS33 [get_ports {AN[3]}]
## Reloj CLK conectado a W5
set_property IOSTANDARD LVCMOS33 [get_ports {CLK}]
set_property PACKAGE_PIN W5 [get_ports {CLK}]
create_clock -period 10.000 -name CLK -waveform {0.000 5.000} [get_ports CLK]

En la última línea se definen las restricciones para el reloj externo a la FPGA, que es de 100MHz en el caso de la placa Basys 3 (a partir del cual el MMCM generará los 50MHz que hemos configurado).

El código VHDL del proyecto es el siguiente (nombre del fichero mmcm_top.vhd).

library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
use IEEE.NUMERIC_STD.all;
entity mmcm_top is
Generic ( 
        -- reloj a 50Mhz: 50000000 ciclos en un segundo.
        N : NATURAL := 50000000
    );
    
    Port ( CLK : in STD_LOGIC;
        AN : out STD_LOGIC_VECTOR (3 downto 0);
        LED : out STD_LOGIC_VECTOR (15 downto 0));
end mmcm_top;

architecture Behavioral of mmcm_top is    
    -- contador de 1hz
    signal clk_1hz : NATURAL;
    
    -- señal de habilitación para contador binario
    signal cont_ena : STD_LOGIC;
    
    -- señales para el MMCM
    signal reset_in: STD_LOGIC;
    signal MCLK: STD_LOGIC;
    signal locked: STD_LOGIC;
    
    -- señal de reset global
    signal mrst: STD_LOGIC;
    
    -- registro binario 8 bits para contador
    signal reg_bin: UNSIGNED (7 downto 0);
    
    alias leds_contador is LED (7 downto 0); 
     
begin
-- apagamos displays de 7 seg. (no lo usaremos)
    AN <= "1111";
    
    -- mostrar contador en los leds
    leds_contador <= std_logic_vector(reg_bin);
    
    -- Control reset global del sistema.
    -- La señal locked se pone a 1 cuando el reloj es estable.
    -- Es entonces cuando lo podemos usar.
    -- Por eso aquí lo usamos como reset global del sistema.
    mrst <= not locked;
    
    -- reset del reloj mmcm no activo
    reset_in <= '0'; 
    
    -- Módulo MMCM (Core IP)
    MMCM_clk: entity WORK.clk_wiz_0 port map ( 
        clk_in1     => CLK,
        reset       => reset_in,
        clk_out1    => MCLK,
        locked      => locked );
        
  
    -- contador 1Hz
    process(MCLK)
    begin
        if (MCLK'event and MCLK='1') then
            if clk_1hz < N then 
                clk_1hz <= clk_1hz + 1;
            else
                clk_1hz <= 0;
            end if;
        end if;
    end process;
    
    -- señal de fin de cuenta 1Hz. Pulso para contador binario,
    cont_ena <= '1' when clk_1hz=N else '0';
    
    -- contador binario (sincrono)
    process(MCLK)
    begin
        if (MCLK'event and MCLK='1') then
            if (mrst = '1') then
                -- reset síncrono. 
                -- msrt se pone a 0 cuando el reloj es estable.
                reg_bin <= (others => '0');
             elsif (cont_ena = '1') then
                reg_bin <= reg_bin + 1;
            end if;
        end if;    
    end process;
                     
end Behavioral;

Si todo está correcto, en la jerarquía del proyecto debemos ver que MMCM_clk está subordinado a mmcm_top.

Si observamos que hay un error como el de la figura siguiente no le hacemos caso, parece ser un fallo del entorno, ya que todo está correcto.

Analicemos el código. El contador de 1Hz ya nos debe ser familiar, pero ahora lo hemos modificado un poco. La señal de reloj del contador es MCLK, que es generado por el módulo MMCM.

    -- contador 1Hz
    process(MCLK)
    begin
        if (MCLK'event and MCLK='1') then
            if clk_1hz < N then 
                clk_1hz <= clk_1hz + 1;
            else
                clk_1hz <= 0;
            end if;
        end if;
    end process;
    
    -- señal de fin de cuenta 1Hz. Pulso para contador binario,
    cont_ena <= '1' when clk_1hz=N else '0';

Cuando el contador llega a su fin de cuenta (lo que ocurre con una frecuencia de 1Hz) se pone a 1 la señal cont_ena, que es la señal de habilitación que está conectada al contador binario, cuyo código es el siguiente (en realidad debería ser cuándo llegara a N-1, pero no vamos a perdernos en detalles nimios :-) ).

    -- contador binario (sincrono)
    process(MCLK)
    begin
        if (MCLK'event and MCLK='1') then
            if (mrst = '1') then
                -- reset síncrono. 
                -- msrt se pone a 0 cuando el reloj es estable.
                reg_bin <= (others => '0');
             elsif (cont_ena = '1') then
                reg_bin <= reg_bin + 1;
            end if;
        end if;    
    end process;

Nuestro contador binario es totalmente síncrono y, como vemos, usa la misma señal de reloj MCLK que el contador de 1Hz (ambos contadores usan la misma señal de reloj). La señal msrt está conectada mediante un inversor a la salida locked del MMCM. Esta señal se mantiene a 0 mientras se estabiliza el reloj, y cuando éste está estable se pone a 1, así que lo usamos como señal de reset global para el sistema. Así mantenemos el contador binario parado hasta que el reloj se estabilice, y después ya podemos empezar a contar de forma segura. Cada vez que la señal cont_ena (señal de habilitación) se pone en alto, el contador binario suma 1. Recordemos que en el reloj digital que diseñamos en un artículo anterior usábamos la salida del contador de 1Hz como señal de reloj para el segundero. Es decir teníamos dos dominios de reloj diferentes, que como hemos visto, en un circuito real puede ser una fuente de problemas inesperados. El contador binario, que es un registro de 8 bits, está conectado directamente a los LEDs de la placa.

     -- mostrar contador en los leds
    leds_contador <= std_logic_vector(reg_bin);

Así pues, nuestro diseño sólo tiene una fuente de reloj y, además, al utilizar un módulo de reloj MMCM nos aseguramos de que minimizaremos al máximo los problemas derivados de un posible desfase en los flancos activos en las diferentes partes del circuito.

Timers en MSP430

noreply@blogger.com (Alberto García Serrano) — Sun, 22 Jan 2017 16:11:56 +0000

El MSP430 está provisto de temporizadores (el número de ellos depende de la versión del microcontrolador) que permiten medir intervalos de tiempo y generar eventos basados en tiempo.

No hay que confundir los timers con los relojes de tiempo real que también equipa el MSP430.

Hay dos tipos de timers, Los Timers_A y los Timers_B. Ambos son muy similares, y aunque los de tipo B poseen alguna característica más que los A, aquí no vamos a profundizar en ello.

Alguna de las características de los timers son:

Utilizan un contador de 16 bits con cuatro modos de operación.
Pueden seleccionar y usar diferentes fuentes de reloj.
Utilizan registros de captura/comparación para su configuración.
Permiten salidas PWM (Power Width Modulation).
Pueden usar interrupciones.

Nos apoyaremos en el siguiente esquema de un Timer_A para ver cómo se utiliza.

Hay dos tipos de bloques bien diferenciados. Rodeado de verde está el bloque contador, y de azul el bloque de captura/comparación, de los que puede haber hasta siete en cada timer.

Bloque contador

El bloque contador dispone de un contador de 16 bits llamado TAxR (donde x será 0 para el timer TA0 y así sucesivamente). Para controlar este bloque se utiliza el registro TAxCTL.

El campo TASSEL nos permite seleccionar la fuente de reloj.
El campo ID permite dividir la frecuencia del reloj.
El campo MC permite seleccionar el modo de funcionamiento del contador.
El campo TACLR pone a cero el contador TAxR.
El campo TAIE habilita las interrupciones generadas por el timer.
El campo TAIFG indica si se ha producido una interrupción. Generalmente es un campo para consulta.

Como se puede observar en el campo MC, disponemos de cuatro modos de cuenta:

Stop mode (00b): El timer está parado.
Up mode (01b): El timer cuenta hasta el valor almacenado en el registro TAxCCR0 y vuelve a empezar desde 0.
Continuous mode (10b): El timer cuenta hasta el valor fijo 0xFFFF y vuelve a contar desde 0.
Up/Down Mode (11b): El timer cuenta ascendentemente hasta el valor almacenado en el registro TAxCCR0 y después descendentemente hasta 0.

En el esquema vemos que, además del divisor de frecuencia controlado por el campo ID, hay otro más. Este se controla a través del registro TAxEX0.

Bloque de captura/comparación

Este bloque se controla mediante el registro TAxCCTLn.

Vamos a analizar el uso básico de los contadores, así que nos centraremos en el modo comparación. Así pues, el modo captura no lo analizaremos en este artículo.

El campo CAP nos permitirá seleccionar el modo captura o el modo comparación. Por otro lado, el campo CCIFG del registro TAxCCTLn se activa cuando el valor de la cuenta del registro TAxR alcanza el valor del registro TAxCCRn.

El timer puede generar dos tipos de interrupciones. Por un lado, mediante el bit TAIE del registro TAxCTL podemos generar una interrupción cada vez que el contador alcance el valor 0. El segundo tipo de interrupciones se desencadenan cuando el valor de la cuenta del registro TAxR ha alcanzado el valor almacenado en el registro TAxCCRn (TAxCCRn CCIFG), se habilitan mediante el bit CCIE del registro TAxCCTLn.

Código de ejemplo

Para aclarar conceptos los mejor es analizar un código de ejemplo en el que vamos a poner en práctica lo que se ha visto hasta ahora.

#include <msp430.h>
void main(void)
{
  WDTCTL = WDTPW + WDTHOLD;    // Detiene Watchdog
  P1DIR |= BIT0;               // P1.0 salida
  P1OUT |= BIT0;               // Enciende LED  
  // Configuración TIMER_A:
  // TimerA1, ACLK/1,modo up, reinicia TACLR
  TA1CTL = TASSEL_1 | TACLR | MC_1; 
  // ACLK tiene una frecuencia de 32768 Hz
  // Carga cuenta en TA1CCR0 0.1seg TA1CCR=(0,1*32768)-1
  TA1CCR0 = 3276;  
  TA1CCTL0 = CCIE; // Habilita interrupción (bit CCIE)
  // LPM0, habilitación general de interrupciones
  __low_power_mode_0(); 
  __no_operation();
}

// Rutina de interrupción de TIMER1_A0  
#pragma vector=TIMER1_A0_VECTOR
__interrupt void TIMER1_A0_ISR(void){
  P1OUT ^= 0x01;       // conmuta LED en P1.0
}

Este código hace cambiar de estado, cada 0,1 segundos, el LED conectado al puerto P1.0. Para ello se ha programado la interrupción asociada al timer A1. Mediante la siguiente línea se configura el timer A1 para que utilice el reloj ACLK y configuración en modo UP. También se inicializa el contador poniéndolo a cero.

TA1CTL = TASSEL_1 | TACLR | MC_1;

La interrupción se produce cuando el contador llega al valor almacenado en el registro TACCR0, así que sabiendo que, en este caso, ACLK funciona a una frecuencia de 32768Hz podemos calcular qué valor hay que almacenar en dicho registro. Si quisiéramos que el LED cambiara de estado cada segundo lo tendríamos fácil. Simplemente haríamos TACCR0=32768-1 (siempre restaremos uno al valor deseado). Esto es así porque ACLK cambia de estado 32768 veces en un segundo. Para que el LED cambie un valor arbitrario de veces n usamos la formula TACCR0=(tiempo en segundos * freq. de reloj en Hz)-1.

Así pues, para que cambie de estado cada 0,1 segundos hacemos TA1CCR0=(0,1*32768)-1=3276.

TA1CCR0 = 3276;

Finalmente habilitamos la interrupción del timer mediante el bit CCIE.

TA1CCTL0 = CCIE;

Al entrar en modo LPM0 también estamos activando las interrupciones.

En el siguiente vídeo podemos ver el resultado de la ejecución del anterior código.

Modos de bajo consumo en MSP430

noreply@blogger.com (Alberto García Serrano) — Sun, 15 Jan 2017 10:27:06 +0000

En el artículo hicimos un repaso al módulo de reloj, que está íntimamente ligado a los modos de bajo consumo del microcontrolador MSP430. Una de las ventajas de éste sobre otros microcontroladores es la gran flexibilidad que ofrece a la hora de configurar los relojes y, por ende, el modo de bajo consumo.

El MSP430 permite la desactivación de las fuentes de reloj y de los propios relojes, lo que permite desactivar la CPU y los periféricos mientras no sean necesarios. Además, el módulo de gestión de energía (PMM o Power Management Module) hace posible incluso apagar el core (CPU) del microcontrolador.

En la práctica, la mayoría de las aplicaciones no requieren que la CPU o los periféricos estén siempre funcionando. Imaginemos un dispositivo IoT que envía información sobre contaminación atmosférica a un servidor remoto. En este caso, no es necesario estar constantemente consultando el estado de los sensores y enviando información, sino que puede hacerse cada pocos segundos o incluso minutos. Si la operación de lectura de los sensores más en envío de los datos necesitan, por simplificar, 1 segundo y esta operación se realiza cada 60 segundos, aproximadamente el 98% del tiempo el microcontrolador está ocioso, pero consumiendo energía. Si "dormimos" la CPU y los periféricos durante ese tiempo, el consumo cae drásticamente, y eso en un dispositivo autónomo es lo más importante. Además de poder "dormirlo" necesitaremos una forma de "despertarlo", lo que en general conseguiremos mediante el uso de interrupciones. Evidentemente, esta operación de dormir-despertar conlleva una penalización en tiempo que podría no ser viable en un dispositivo que necesita responder en tiempo real.

Modos de bajo consumo

Los modos de bajo consumo se configuran mediante el registro de estado de la CPU (registro SR).

En concreto se utilizan los bits SCG0, SCG1, OSCOFF y CPUOFF. Mediante la activación de estos bits podemos seleccionar alguno de los modos de bajo consumo definidos en la siguiente tabla.

Los modos LPM3.4 y LPM4.5 no están disponibles en todas las versiones del microcontrolador. En el MSP430F5529, que es el que estoy usando, el modo LPM3.5 no está disponible.

Con todos los bits a 0 el microcontrolador se encuentra en estado activo. En este modo, la CPU y las fuentes de reloj están activas, por lo tanto es el modo de mayor consumo. A cambio la respuesta es más rápida.

En el modo LPM0 y LPM1 se para la CPU como consecuencia de la desactivación del reloj MCLK. La diferencia entre ambos es que en el modo LPM0 el bloque FLL se desactiva si no se está usando el DCO y en el modo LPM1 el FLL se desactiva siempre.

En el modo LMP2 y LMP3 además de MCLK (y por ende la CPU) también se desactiva SMCLK. ACLK permanece activo. Se diferencian en que LPM2 mantiene parte del FLL activo lo que facilita una más rápida estabilización del DCO (y también algo más de consumo).

En el modo LPM4 la CPU y todos los relojes de desactivan.

Los modos LPM3.5 y LPM4.5 paran la CPU y todos los relojes, pero además el PMM desactiva la alimentación del core, lo que implica la perdida del contenido de los registros y la memoria RAM dinámica.

Código de ejemplo

Partimos del mismo código de ejemplo que ya utilizamos en el artículo sobre interrupciones. La única diferencia con aquel código es una sola línea dentro del bucle while:

__low_power_mode_3();

Esta línea hace que el microcontrolador entre en modo LPM3 (es fácil adivinar cómo se entra en, por ejemplo, LPM2 ¿no?), lo que implica que el reloj MCLK (por lo tanto también la CPU) y SMLK están parados. Sin embargo ACLK está activo, lo que permite que el sistema mantenga los sistemas periféricos en marcha y por lo tanto, aunque la CPU esté parada, en el momento en que se produzca una interrupción, la rutina de gestión la va a tratar. Cuando esto sucede, obviamente el reloj MCLK se pone en marcha y permite que se ejecute dicha rutina. Al terminar la rutina restaura desde la pila los valores del registro de estados SR, que recordemos, permiten gestionar los modos de bajo consumo. Así pues, tras la ejecución de la rutina de interrupción el microcontrolador vuelve al estado de bajo consumo en el que se encontraba antes de entrar, ya que se restaura el registro SR tal y como estaba antes de la interrupción.

#include <msp430.h>
int main(void)
{
  WDTCTL = WDTPW + WDTHOLD; // Detiene el watchdog
  P1DIR |= BIT0;            // P1.0 (LED) como salida digital
  P1REN |= BIT1;	    // Resistencia en P1.1 (pulsador derecho)
  P1OUT |= BIT1;	    // Modo pull-up
  P1IE |= BIT1;		    // habilita la interrupción de P1.1
  P1IES |= BIT1;	    // Interrupción por flanco de bajada
  P1IFG &= ~(BIT1);	    // baja flag de interrupción de P1.1
  // Activar interrupciones
  // también puede usarse: __enable_interrupt()
  __bis_SR_register(GIE);   
  // bucle infinito.
  while (1) {
    __low_power_mode_3(); // entramos en modo LPM3
    __no_operation();
  }
}
#pragma vector=PORT1_VECTOR
__interrupt
void P1_ISR(void) {
  if (P1IFG & BIT1) {	// comprobar fuente de la interrupción P1.1
    if (P1IN & BIT1) {	// Estado de P1.1: si está a 1 (no pulsado)
      P1OUT &= ~BIT0;   // Apaga LED
      P1IES |= BIT1;	// activación por flanco de bajada
    } else {            // y si no (pulsado),
      P1OUT |= BIT0;    // Enciende LED
      P1IES &= ~BIT1;	// activación por flanco de subida
    }
    P1IFG &= ~BIT1;	// puesta a cero del flag de interrupción P1.1
  }
}

Vamos a medir, polímetro en mano, qué efecto produce la línea que hemos añadido en el consumo del sistema. Podemos medir la corriente consumida usando el jumper de 3.3V de la placa launchpad (al quitar el jumper tenemos dos pines sobre los que medir).

Como se puede observar en la imagen, sin utilizar el modo de bajo consumo (modo activo) tenemos un consumo en stand by (sin tener en cuenta el consumo al pulsar el botón y encender el LED) de alrededor de 1mA. Ahora veamos el consumo en modo LPM3.

Tenemos 0,66mA de consumo en LPM3. Este consumo aun puede mejorarse, pero nos vale como ejemplo práctico del uso de los modos de bajo consumo del MSP430.