Escribir un archivo binario desde Python y leerlo en C

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

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

C:

Máquina de Estados Finitos en C [PIC18F45K50]

Para esta implementación en software de una máquina de estados finitos (FSM: Finite State Machine) haremos uso de punteros a funciones [click en el enlace para ver un video-tutorial]. Usaremos también  el compilador XC8 en MPLAB X y un PIC18F45K50 con las mismas configuraciones que en la entrada anterior.

En una gran cantidad de problemas a resolver mediante un sistema embebido es necesario llevar a cabo un proceso estrictamente secuencial. En estos casos es muy útil y elegante estructurar nuestros programas como una máquina de estados. Quizá muchos de ustedes estén tentados a usar una estructura case para elegir un segmento de código cuando se cumple una condición. Si el numero de estados es pequeño, quizá se una opción fácil. Desafortunadamente el uso de case's es muy ineficiente cuando la máquina de estados y su tabla de transiciones es muy compleja. La implementación que se mostrará en esta entrada parecerá un tanto rebuscada pero tiene la ventaja de tener buen despeño a la hora de escalarla a numero grande de estados.

El diagrama de estados de la FSM de este ejemplo es el siguiente:

Diagrama hecho en LaTeX con TIKz

Las entradas 0/1 representan los estados de un switch conectado a RB0. Los estados se harán visibles enviando el numero correspondiente al Puerto D [esta acción será definida en la función asociada a cada estado]. Obsérvese que está máquina describe a un contador ascendente o descendente según la entrada. Se agregará un retardo de un segundo para poder visualizar las transiciones.

Movámonos al programa. Primero, debemos definir un tipo enumerado para codificar nuestro conjunto de estados:

//-- Definición de tipo de estados
typedef enum {
    Estado_0,      // Estado inicial
    Estado_1,
    Estado_2,
    Estado_3,
    Estado_4,      
    No_de_Edos
} ESTADOS;

Creamos con esto una variable enumerada global que usaremos para almacenar el estado actual de la máquina:

//-- Variable global de estado
ESTADOS Estado_Actual = Estado_0;

El siguiente paso es definir una estructura que nos permitirá construir un objeto para nuestra FSM:

//-- Estructura de Maquina de Estados (FSM))
typedef struct {
    ESTADOS Estado; // Estado miembro del conjunto enumerado ESTADOS
    void (*func)(int); // Acciones asociadas al estado         
} FSM;

Con lo anterior podemos ahora crear una instancia que sera nuestra FSM:

//-- Prototipos de funciones de estado
void Codigo_EDO_0(int);
void Codigo_EDO_1(int);
void Codigo_EDO_2(int);
void Codigo_EDO_3(int);
void Codigo_EDO_4(int);

//-- Creación de la estructura global de la FSM
FSM Maquina_de_Estados[] = {
    {Estado_0,Codigo_EDO_0},
    {Estado_1,Codigo_EDO_1},
    {Estado_2,Codigo_EDO_2},
    {Estado_3,Codigo_EDO_3},
    {Estado_4,Codigo_EDO_4},
};

Es importante que los prototipos de las funciones asociadas a cada estado aparezcan antes de la creación de la estructura (estas funciones serán definidas más adelante en el código). El segmento de código principal queda reducido elegantemente de la siguiente manera:

void main(void) {
    int sw = 0;
    setup();
   
    //-- FSM ---
    while(1){
        sw = PORTBbits.RB0;
        (*Maquina_de_Estados[Estado_Actual].func)(sw);
    }
    return;
}

Código completo:
Circuito utilizado:

Programación multithread en C (pthread.h)

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

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

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

int pthread_create(thread_ID,att,rutina_concurrente,argumento) 

thread_ID: apuntador al identificador de hilo (tipo pthread_t)

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

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

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

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

Ejemplo: Multiplicación matricial 

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

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

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

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

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

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

  pthread_t tid[num_of_threads];
  long rank;

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

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

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

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

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

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

El código completo en C es el siguiente:

Compilamos y ejecutamos desde terminal con:

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

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

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

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

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

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

Instalación de GNU Scientific Library (Debian)

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

Para instalar GSL en Debian y derivados:

sudo apt-get install libgsl0ldbl libgsl0-dev 

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

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

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

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

Compilamos y ejecutamos desde terminal con:

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

Programa de ejemplo en C para MicroBlaze

 

(Tarjeta de desarrollo Symbhia)

#include <xparameters.h>
#include <xgpio.h>

XGpio LEDS;
XGpio SW;

u8 sws; // Variable tipo uint8
int i;

int main(){
//Inicialización
  XGpio_Initialize(&LEDS,XPAR_GPIO_1_DEVICE_ID);
  XGpio_Initialize(&SW,XPAR_GPIO_0_DEVICE_ID);

while(1){
//Lectura de switches
  sws = XGpio_DiscreteRead(&SW,1); //(puntero de instancia,canal)

//Escritura de la variable sws en LEDS
  XGpio_DiscreteWrite(&LEDS,1,sws); // (puntero de instancia,canal,dato)
   }

return 0;

}

Escribir y compilar un programa en C con gcc desde terminal

Abrimos una consola, asegurándonos de que estamos en nuestra carpeta personal (puedes hacerlo con el comando pwd) y escribimos:

$ mkdir Programas  # Creamos una nueva carpeta
$ cd Programas       # Nos movemos a la carpeta Programas
~/Programas$ nano ejemplo.c # Creamos nuevo archivo .c

~/Programas$ gcc -o ejemplo ejemplo.c # Compilamos y creamos código objeto con el nombre ejemplo.o
~/Programas$ ./ejemplo # Ejecutamos código objeto

Función promedio para lectura analógica desde PIC (MikroC)

Cómo muchos hemos notado, la conversión A/D en PIC es un tanto ruidosa. El método de promediación no es la solución definitiva al ruido del ADC, pero mejora considerablemente los resultados:

unsigned int promediolectura(void){

  int i =0;
  unsigned int a = 0;

  for (i=0; i <16;i++){a+ = ADC_Read(0);}
  a = a >>4;

  return a;
}