Experimentos con Arduino

sábado, 29 de febrero de 2020

Modelo SEAIHR del brote MERS en Korea del año 2015

Continuando con el post anterior, aqui tenemos el modelo SEAIHR del brote MERS en la República de Korea durante el año 2015.
El modelo SEAIHR es una evolución del modelo SEIR tradicional , e incluye dos estados (compartimientos) adicionales:

Susceptible -> Expuesto -> Asintomático -> Infectado > Hospitalizado -> Recuperado

En el enlace PLOS|ONE tenemos el artículo con los detalles de este modelo y su aplicación al brote de Korea.

El script Arduino se muestra a continuación

////
// SEAIHR 2015  MERS-CoV in the Republic of Korea 
// https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0144778#pone-0144778-t002
// JLG biolansys 02/2020
//

int   i,j,t ;
float alfa,beta,beta1,beta2,beta3,gamma,zeta,zeta0,sigma,k1,k2,k20,lambda,theta,mu,nu,tfac,R0;
float pop,sus,expo,inf,rec,asi,hos;
float s1,e1,i1,r1,a1,h1;

void SEAIHR(float sus,float expo,float inf,float rec,float asi,float sigma,float lambda,float k1,float k2,float gamma,float zeta,float beta1,float beta2,float beta3,int tfac,int t) {   
  pop=sus+expo+inf+rec+hos+asi;  // Total Population
  sigma=sigma/tfac;
  lambda=lambda/tfac;
  k1=k1/tfac;
  k2=k2/tfac;
  gamma=gamma/tfac; 
  zeta=zeta/tfac; 
  beta1=beta1/pop/tfac;
  beta2=beta2/pop/tfac;
  beta3=beta3/pop/tfac;  
   
  //Serial.print ("Susceptible" ); Serial.print(" ");
  //Serial.print ("Exposed" ); Serial.print(" ");
  //Serial.print ("Asynthomatic" ); Serial.print(" ");
  Serial.print ("Infected" ); Serial.print(" ");
  Serial.print ("Hospitalized" ); Serial.print(" ");
  //Serial.print ("Recovered" ); Serial.print(" ");
  //Serial.print ("Population" ); Serial.print(" ");
  for (i=1;i<=t;i++) {   
      s1= - beta1*sus*asi - beta2*sus*inf  - beta3*sus*hos ;
      e1=   beta1*sus*asi + beta2*sus*inf  + beta3*sus*hos -  ( sigma) *expo;
      a1= (1- gamma)*sigma*expo -  k1 * asi ;
      i1= gamma*sigma*expo - lambda*inf;
      h1= lambda*inf - k2*hos -  zeta*hos;
      r1= k1*asi + k2*hos +  zeta*hos ;
            
      //Serial.print (sus); Serial.print(" ");
      //Serial.print (expo); Serial.print(" ");
      //Serial.print (asi); Serial.print(" ");
      Serial.print (inf); Serial.print(" ");
      Serial.print (hos); Serial.print(" ");
      //Serial.print (rec); Serial.print(" ");
      //Serial.print ( pop);
      Serial.println(" ");
  
      sus=sus + s1;
      expo=expo + e1;
      asi=asi + a1;
      inf=inf + i1;
      hos=hos + h1;    
      rec=rec + r1;          
  }  
}

void setup() {
  Serial.begin(115200);  
}

void loop() {
  sus=49519960.0; expo=22.0; inf=1.0; rec=0.0 ; hos=0.0; asi=20.0;
  tfac=5.0; t=500;  
  sigma=1.0/5.2;  lambda=1.0/5.0;  k1=1.0/5.0;  k2=1.0/7.0; gamma=0.0348 ;zeta=1.0/15.0; beta1=0.8756; beta2=0.7833;  beta3=0.4568;
  /*
  k2=5.0;
  for (j=0;j<5;j++) {
    k20=1.0/k2;
    SEAIHR(sus,expo,inf,rec,asi,sigma,lambda,k1,k20,gamma,zeta,beta1,beta2,beta3,tfac,t);  
    k2= k2 + 1.0;
    delay (1000);
  }
  delay (10000);
  */
  k2=1.0/7.0;
  zeta=1.0;
  for (j=0;j<42;j++) {
    zeta0=1.0/zeta;
    SEAIHR(sus,expo,inf,rec,asi,sigma,lambda,k1,k2,gamma,zeta0,beta1,beta2,beta3,tfac,t);  
    zeta=zeta + 1.0;
    delay (1000);
  }
  delay (10000);
}

En el script tenemos comentadas algunas variables, de tal forma que sólo se muestran los valores de Infectados y Hospitalizados, Adicionalmente, en la sección loop() hay dos rangos de emulación que permiten variar un parámetro del modelo para ver como las variables cambian.

A continuación un gif animado con el segundo rango del parámetro zeta.

domingo, 23 de febrero de 2020

Modelos deterministicos de enfermedades contagiosas

Lamentablemente estos días las noticias relacionadas con la "epidemia" del coronavirus están muy presentes en todos los medios de comunicación, y también en el día a día.

Un ejemplo de las capacidades de un microcontrolador de 8bits (Arduino) es la simulación de modelos determinísticos de enfermedades contagiosas.

Los modelos SIR y SEIR son dos de los más utilizados, y sus nombres se corresponden con los estados (compartimientos) que simula cada modelo. Estos modelos tratan de representar el número de elementos en un determinado momento (usualmente días) en cada estado, es decir, es una foto instantánea de la población a simular y su distribución en los estados definidos en el modelo.

El modelo SIR (Kermack WO, McKendrick AG .1927) es aplicable a enfermedades como la rubeola, el sarampión y las paperas. Este modelo utiliza tres estados posibles :

Susceptible -> Infectado > Recuperado

Se asume que la suma de todos los elementos en cada compartimiento es constante y representa el total de la población a simular.

A continuación el script de este modelo y una imagen de la simulación:

//
// SIR  with vital dynamics
// beta : The parameter controlling how often a susceptible-infected 
//        contact results in a new infection.
// gamma: The rate an infected recovers and moves into the 
//        resistant phase.
// mu and nu represents the birth and death rates .
// This models a population of a constant size. 
//

int   i,t ;
float beta,gamma,mu,nu,tfac;
float pop,sus,inf,rec;
float s1,i1,r1;

void setup() {
  Serial.begin(19200);
   
  sus=49995.0;        // Initial Susceptible
  inf=5.0;            // Initial Infected
  rec=0.0 ;           // Initial Recupered
  pop=sus+inf+rec;    // Population simulated

  tfac=10.0;          // X axis factor
  beta=0.7/pop/tfac;
  gamma=0.2/tfac;
  mu=0.0/tfac;
  nu=0.0/tfac;
  t=500;              // Dyas simulated

  //Serial.print ("Days" ); Serial.print(" ");
  Serial.print ("Susceptible" ); Serial.print(" ");
  Serial.print ("Infected" ); Serial.print(" ");
  Serial.print ("Recovered" ); Serial.print(" ");
  Serial.print ("Population" ); Serial.print(" ");
  for (i=1;i<=t;i++) {
    
      s1= mu*rec -beta*sus*inf - nu*sus;
      i1= beta*sus*inf - gamma*inf - nu*inf;
      r1= gamma*inf - mu*rec;
       
      //Serial.print (i); Serial.print(" ");
      Serial.print (sus); Serial.print(" ");
      Serial.print (inf); Serial.print(" ");
      Serial.print (rec); Serial.print(" ");
      Serial.println ( pop);
  
      sus=sus + s1;
      inf=inf + i1; 
      rec=rec + r1;   
        
  }
}

void loop() {
  // Nothing to do
}

El modelo SEIR es aplicable a enfermedades como la malaria, ébola y SARS. Este modelo utiliza cuatro estados (compartimientos) posibles :

Susceptible -> Expuesto -> Infectado > Recuperado

Se asume que la suma de todos los elementos en cada compartimiento es constante y representa el total de la población a simular.

//
// SEIR  with vital dynamics
// beta  The parameter controlling how often a susceptible-infected 
//       contact results in a new exposure.
// gamma The rate an infected recovers and moves into the resistant 
//       phase.
// sigma The rate at which an exposed person becomes infective.
// mu  The natural mortality rate (this is unrelated to disease).   
// nu  The natural birth rate  .
// This models a population of a constant size. 
// JLG biolansys 02/2020
//

int   i,t ;
float alfa,beta,gamma,sigma,theta,mu,nu,tfac;
float pop,sus,expo,inf,rec;
float s1,e1,i1,r1;

void setup() {
  Serial.begin(115200);
   
  sus=99999.0;           // Initial Susceptible
  expo=0.0;              // Initial Exposed
  inf=1.0;               // Initial Infected
  rec=0.0 ;              // Initial Recupered
  pop=sus+expo+inf+rec;  // Total Population

  tfac=5.0;             //x axis factor
 
  beta=0.9/pop /tfac;
  gamma=0.5/tfac;
  sigma=0.2/tfac;
  mu=0.01/tfac;
  nu=0.005/tfac;
      
  t=500;  // Days simulated

  //Serial.print ("Days" ); Serial.print(" ");
  Serial.print ("Susceptible" ); Serial.print(" ");
  Serial.print ("Exposed" ); Serial.print(" ");
  Serial.print ("Infected" ); Serial.print(" ");
  Serial.print ("Recovered" ); Serial.print(" ");
  Serial.print ("Population" ); Serial.print(" ");
  for (i=1;i<=t;i++) {
    
      s1= mu*(pop-sus) - beta*sus*inf -nu*sus ;
      e1= beta*sus*inf -  ( mu + gamma) *expo;
      i1= gamma*expo -  (mu + sigma ) *inf  ;
      r1= sigma*inf - mu*rec + nu*sus;
       
      //Serial.print (i); Serial.print(" ");
      Serial.print (sus); Serial.print(" ");
      Serial.print (expo); Serial.print(" ");
      Serial.print (inf); Serial.print(" ");
      Serial.print (rec); Serial.print(" ");
      Serial.println ( pop);
  
      sus=sus + s1;
      expo=expo + e1;
      inf=inf + i1; 
      rec=rec + r1;   
        
  }
}

void loop() {
  // Nothing to do
}

En ambos modelos se consideran los efectos de nacimientos y muertes (no relacionadas con la enfermedad) en la población estudiada.

Los gráficos mostrados se generaron mediante la opción Serial Plotter del IDE de Arduino. Se utiliza la variable tfac para ajustar el gráfico en el eje X.

domingo, 20 de enero de 2019

Red Neuronal en Arduino (machine learning)

El tema de Machine Learning está de moda y como es natural , desde hace algún tiempo ya hay implementaciones de ejemplo para Arduino.

Algunos enlaces de interés relacionados con Arduino:

Una implementación interesante de una Red Neuronal en Arduino esta descrita en el enlace http://robotics.hobbizine.com/arduinoann.html

He modificado ligeramente el script original :

La velocidad de conexión con el terminal es 57600 bps.
Hay una nueva variable (PrintDetail) para indicar si queremos el detalle de los cálculos.
La variable CycleReport indica cada cuantas iteracciones se muestra el resultado detallado.
Al seleccionar detalle de cálculos se incluye la impresión de los valores de los nodos Hidden.

/****************************************************************
 * ArduinoANN - An artificial neural network for the Arduino
 * All basic settings can be controlled via the Network 
 * Configuration section.
 * See robotics.hobbizine.com/arduinoann.html for details.
 ****************************************************************/

#include <math.h>

/****************************************************************
 * Network Configuration - customized per network 
 ****************************************************************/

const int PatternCount = 10;
const int InputNodes = 7;
const int HiddenNodes = 8;
const int OutputNodes = 4;
const int CycleReport = 100;
const float LearningRate = 0.3;
const float Momentum = 0.9;
const float InitialWeightMax = 0.5;
const float Success = 0.0004;

const byte Input[PatternCount][InputNodes] = {
  { 1, 1, 1, 1, 1, 1, 0 },  // 0
  { 0, 1, 1, 0, 0, 0, 0 },  // 1
  { 1, 1, 0, 1, 1, 0, 1 },  // 2
  { 1, 1, 1, 1, 0, 0, 1 },  // 3
  { 0, 1, 1, 0, 0, 1, 1 },  // 4
  { 1, 0, 1, 1, 0, 1, 1 },  // 5
  { 0, 0, 1, 1, 1, 1, 1 },  // 6
  { 1, 1, 1, 0, 0, 0, 0 },  // 7 
  { 1, 1, 1, 1, 1, 1, 1 },  // 8
  { 1, 1, 1, 0, 0, 1, 1 }   // 9
}; 

const byte Target[PatternCount][OutputNodes] = {
  { 0, 0, 0, 0 },  
  { 0, 0, 0, 1 }, 
  { 0, 0, 1, 0 }, 
  { 0, 0, 1, 1 }, 
  { 0, 1, 0, 0 }, 
  { 0, 1, 0, 1 }, 
  { 0, 1, 1, 0 }, 
  { 0, 1, 1, 1 }, 
  { 1, 0, 0, 0 }, 
  { 1, 0, 0, 1 } 
};

/****************************************************************
 * End Network Configuration
 ****************************************************************/


int i, j, p, q, r;
int ReportEvery;
int PrintDetail;  
int RandomizedIndex[PatternCount];
long  TrainingCycle;
float Rando;
float Error;
float Accum;


float Hidden[HiddenNodes];
float Output[OutputNodes];
float HiddenWeights[InputNodes+1][HiddenNodes];
float OutputWeights[HiddenNodes+1][OutputNodes];
float HiddenDelta[HiddenNodes];
float OutputDelta[OutputNodes];
float ChangeHiddenWeights[InputNodes+1][HiddenNodes];
float ChangeOutputWeights[HiddenNodes+1][OutputNodes];

void setup(){
  //Serial.begin(9600);
  Serial.begin(57600);
  randomSeed(analogRead(3));
  ReportEvery  = 1;
  PrintDetail = 1;
  for( p = 0 ; p < PatternCount ; p++ ) {    
    RandomizedIndex[p] = p ;
  }
}  

void loop (){
/****************************************************************
* Initialize HiddenWeights and ChangeHiddenWeights 
****************************************************************/
  for( i = 0 ; i < HiddenNodes ; i++ ) {    
    for( j = 0 ; j <= InputNodes ; j++ ) { 
      ChangeHiddenWeights[j][i] = 0.0 ;
      Rando = float(random(100))/100;
      HiddenWeights[j][i] = 2.0 * ( Rando - 0.5 ) * InitialWeightMax ;
    }
  }
/****************************************************************
* Initialize OutputWeights and ChangeOutputWeights
****************************************************************/

  for( i = 0 ; i < OutputNodes ; i ++ ) {    
    for( j = 0 ; j <= HiddenNodes ; j++ ) {
      ChangeOutputWeights[j][i] = 0.0 ;  
      Rando = float(random(100))/100;        
      OutputWeights[j][i] = 2.0 * ( Rando - 0.5 ) * InitialWeightMax ;
    }
  }
  Serial.print("Initial/Untrained Outputs (Random): ");
  toTerminal();
/****************************************************************
* Begin training 
****************************************************************/

  for( TrainingCycle = 1 ; TrainingCycle < 2147483647 ; TrainingCycle++) {    

/****************************************************************
* Randomize order of training patterns
****************************************************************/
    for( p = 0 ; p < PatternCount ; p++) {
      q = random(PatternCount);
      r = RandomizedIndex[p] ; 
      RandomizedIndex[p] = RandomizedIndex[q] ; 
      RandomizedIndex[q] = r ;
    }
    Error = 0.0 ;
/****************************************************************
* Cycle through each training pattern in the randomized order
****************************************************************/
    for( q = 0 ; q < PatternCount ; q++ ) {    
      p = RandomizedIndex[q];

/****************************************************************
* Compute hidden layer activations
****************************************************************/

      for( i = 0 ; i < HiddenNodes ; i++ ) {    
        Accum = HiddenWeights[InputNodes][i] ;
        for( j = 0 ; j < InputNodes ; j++ ) {
          Accum += Input[p][j] * HiddenWeights[j][i] ;
        }
        Hidden[i] = 1.0/(1.0 + exp(-Accum)) ;
      }

/****************************************************************
* Compute output layer activations and calculate errors
****************************************************************/

      for( i = 0 ; i < OutputNodes ; i++ ) {    
        Accum = OutputWeights[HiddenNodes][i] ;
        for( j = 0 ; j < HiddenNodes ; j++ ) {
          Accum += Hidden[j] * OutputWeights[j][i] ;
        }
        Output[i] = 1.0/(1.0 + exp(-Accum)) ;   
        OutputDelta[i] = (Target[p][i] - Output[i]) * Output[i] * (1.0 - Output[i]) ;   
        Error += 0.5 * (Target[p][i] - Output[i]) * (Target[p][i] - Output[i]) ;
      }

/****************************************************************
* Backpropagate errors to hidden layer
****************************************************************/

      for( i = 0 ; i < HiddenNodes ; i++ ) {    
        Accum = 0.0 ;
        for( j = 0 ; j < OutputNodes ; j++ ) {
          Accum += OutputWeights[i][j] * OutputDelta[j] ;
        }
        HiddenDelta[i] = Accum * Hidden[i] * (1.0 - Hidden[i]) ;
      }

/****************************************************************
* Update Inner-->Hidden Weights
****************************************************************/

      for( i = 0 ; i < HiddenNodes ; i++ ) {     
        ChangeHiddenWeights[InputNodes][i] = LearningRate * HiddenDelta[i] + Momentum * ChangeHiddenWeights[InputNodes][i] ;
        HiddenWeights[InputNodes][i] += ChangeHiddenWeights[InputNodes][i] ;
        for( j = 0 ; j < InputNodes ; j++ ) { 
          ChangeHiddenWeights[j][i] = LearningRate * Input[p][j] * HiddenDelta[i] + Momentum * ChangeHiddenWeights[j][i];
          HiddenWeights[j][i] += ChangeHiddenWeights[j][i] ;
        }
      }

/****************************************************************
* Update Hidden-->Output Weights
****************************************************************/

      for( i = 0 ; i < OutputNodes ; i ++ ) {    
        ChangeOutputWeights[HiddenNodes][i] = LearningRate * OutputDelta[i] + Momentum * ChangeOutputWeights[HiddenNodes][i] ;
        OutputWeights[HiddenNodes][i] += ChangeOutputWeights[HiddenNodes][i] ;
        for( j = 0 ; j < HiddenNodes ; j++ ) {
          ChangeOutputWeights[j][i] = LearningRate * Hidden[j] * OutputDelta[i] + Momentum * ChangeOutputWeights[j][i] ;
          OutputWeights[j][i] += ChangeOutputWeights[j][i] ;
        }
      }
    }

/****************************************************************
* Every 1000 cycles send data to terminal for display
****************************************************************/
    ReportEvery  = ReportEvery  - 1;
    if (ReportEvery == 0)
    {
      Serial.println(); 
      //Serial.println(); 
      Serial.print ("TrainingCycle: ");
      Serial.print (TrainingCycle);
      Serial.print ("  Error = ");
      Serial.print  (Error, 5);
      Serial.print ("  Success = ");
      Serial.print   (Success, 4);

      toTerminal();

      if (TrainingCycle==1)
      {
        ReportEvery = CycleReport -1;
      }
      else
      {
        ReportEvery  = CycleReport;
      }
    }    

/****************************************************************
* If error rate is less than pre-determined threshold then end
****************************************************************/

    if( Error < Success ) break ;  
  }
  //Serial.println ();
  Serial.println(); 
  Serial.print ("TrainingCycle: ");
  Serial.print (TrainingCycle);
  Serial.print ("  Error = ");
  Serial.print  (Error, 4);
  Serial.print ("  Success = ");
  Serial.print   (Success, 4);

  toTerminal();

  Serial.println ();  
  Serial.println ();
  Serial.println ("Training Set Solved! ");
  Serial.println ("--------"); 
  Serial.println ();
  Serial.println ();  
  ReportEvery = 1;
}

void toTerminal()
{
  for( p = 0 ; p < PatternCount ; p++ ) { 
    if ( PrintDetail==1 )
    { 
      Serial.println(); 
      Serial.print ("  Pattern ");
      Serial.print  (p);      
      Serial.print ("  --> Input: ");
      for( i = 0 ; i < InputNodes ; i++ ) {
        Serial.print (Input[p][i], DEC);
        Serial.print (" ");
      }
      Serial.print ("  Target: ");
      for( i = 0 ; i < OutputNodes ; i++ ) {
        Serial.print (Target[p][i], DEC);
        Serial.print (" ");
      }
      Serial.print ("  Hidden: ");
    } 
  
/****************************************************************
* Compute hidden layer activations
****************************************************************/

    for( i = 0 ; i < HiddenNodes ; i++ ) {    
      Accum = HiddenWeights[InputNodes][i] ;
      for( j = 0 ; j < InputNodes ; j++ ) {
        Accum += Input[p][j] * HiddenWeights[j][i] ;
      }
      Hidden[i] = 1.0/(1.0 + exp(-Accum)) ;
    }

/****************************************************************
* Compute output layer activations and calculate errors
****************************************************************/

    for( i = 0 ; i < OutputNodes ; i++ ) {    
      Accum = OutputWeights[HiddenNodes][i] ;
      for( j = 0 ; j < HiddenNodes ; j++ ) {
        Accum += Hidden[j] * OutputWeights[j][i] ;
      }
      Output[i] = 1.0/(1.0 + exp(-Accum)) ; 
    }
    if ( PrintDetail==1 )
    {
      for( i = 0 ; i <  HiddenNodes ; i++ ) {       
        Serial.print ( Hidden[i], 4);
        Serial.print (" ");
     }
      Serial.print ("  Output: ");
      for( i = 0 ; i < OutputNodes ; i++ ) {       
        Serial.print (Output[i], 4);
        Serial.print (" ");
     }
   }
  }

}

sábado, 12 de septiembre de 2015

Simulador e IDE Arduino en la web

Las aplicaciones que normalmente utilizamos en nuestros proyectos de Arduino son programas que instalamos localmente en nuestro ordenador.

Sin embargo, el auge de aplicaciones cooperativas en la web ha permitido que se hayan hecho desarrollos muy interesantes que nos permiten realizar proyectos sin necesidad de instalar programas localmente.

Dos ejemplos de ello son 123d circuits y codebender.

123d circuits es una página web de la empresa Autodesk que permite hacer simulaciones de circutos electrónicos, diseñar placas, etc. Una característica muy interesante de esta página es que permite realizar circuitos utilizando placas Arduino como uno de sus componentes. Además del diseño del circuito podemos cargar código en la placa y simular su ejecución. El registro es gratis aunque hay modalidades de pago que ofrecen mayores funcionalidades.

codebender es otro ejemplo de una aplicación cooperativa en la web para trabajar con Arduino, en este caso es un IDE que nos permite desarrollar códígo y cargarlo a traves del puerto USB conectado a nuestra placa. La página esta muy bien elaborada y cuenta con muchas funcionalidades.

Es necesario instalar plugins para los navegadores (firefox y chrome). Adicionalmente en la página se recomienda instalar los drivers de comunicación para los puertos USB. En mi caso no los instalé y las pruebas básicas de comunicación con una placa Leonardo fueron satisfactorias.

Al igual que en el caso anterior el registro es gratuito pero tienen opciones de pago con funcionalidades adicionales.

No todo el mundo es dado a utilizar aplicaciones en la web, a pesar de la tendencia cada vez mayor a utilizar la web como nuestro "gran repositorio". Sin embargo hay ocasiones en que esta tecnología ofrece ventajas que no ese encuentran en aplicaciones instaladas localmente. Lo bueno es que están ahi, disponibles para cuando se necesiten (y gratis !!!).

domingo, 1 de marzo de 2015

Números primos grandes en Arduino

Recordando algo de matemáticas podemos decir que un número primo es un número entero que sólo es divisible entre 1 y si mismo. Los números primos son importantes para el común de los mortales porque una gran parte de la seguridad en internet se basa en utilizar un algoritmo muy dificil de descifrar ya que depende de dos números primos muy grandes.
Para mortales no tan comunes , los matemáticos, los números primos son especialmente interesantes ya que no se ha encontrado un algoritmo eficaz para determinar si un número es o no primo. El único algoritmo totalmente determinístico es la criba de eratóstenes , pero este procedimiento requiere recursos inmensos de tiempo y poder de cálculo a medida que los números son mas grandes.
Para otros mortales, como los que tenemos algo que ver con informática, desde muy temprano hemos tenido que desarrollar un programa que emule la criba de eratóstenes. Inclusive se utilizan ciertos algoritmos para determinar el poder de cálculo de un procesador (benchmark).
Es por ello que siempre encontraremos a alguien gastando el tiempo en desarrollar algún programa de búsqueda de números primos para cualquier procesador que tenga entre manos. Y claro, el Arduino no va a quedar huerfano en esta tan encomiable labor !!!
Buscando en la web me encontré con la siguiente entrada en el foro de Arduino http://forum.arduino.cc/index.php?topic=192209.0 que trata sobre la búsqueda de números primos grandes en la plataforma Arduino.
Para quien tenga unos quince minutos le recomiendo que la revisen aunque sea por encima, ya que entre otras cosas podemos ver como evoluciona el programa desarrollado por el Sr. Nick Gammon (por cierto, muy prolífico en temas de electronica y específicamente Arduino).
Como a mi esto de la teoría de números también me atrae, quise hacer un pequeño aporte.
Haciéndole unas pequeñas modificaciones al programa del Sr. Gammon logré que el programa se ejecutara en un Atmega644 a 20Mhz con una LCD donde se muestran los números primos que van siendo identificados. El programa guarda en el eeprom interno la información cada 1000 números primos, de tal forma que continúa despues de un reset o apagado. Este número puede ser cambiado en el programa.
La primera que vez que se graba el programa en la memoria flash del micro debemos eliminar los caracteres de comentarios en las líneas 43,44 y 45 , una vez que el programa ha encontrado los primeros 1000 números primos comentamos de nuevo estas líneas y grabamos de nuevo el programa en el micro.

// 64 bits primes
PROGMEM prog_uchar primes[] = {
    3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 
   37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 
   83, 89, 97,101,103,107,109,113,127,131,137,
  139,149,151,157,163,167,173,179,181,191,193,
  197,199,211,223,227,229,233,239,241,251
};
const byte nPrimes = sizeof(primes)/sizeof(primes[0]);
#include <EEPROM.h>
#include <LiquidCrystal.h>
LiquidCrystal lcd(7, 6, 5, 4, 3, 2); //Configuración de pines del LCD

#include "Pseudoprimes.h";
const uint16_t nPSP = (sizeof(psps)/sizeof(psps[0])) >> 3;
  
uint64_t candidate;
uint64_t leido;
uint64_t Barr_m;
uint8_t Barr_k;

long start;
long timerStart, timerRuntime;
long found = 0; // Number of primes found
int led = 0;
int Direccion = 0;
int Informacion;
char prime_eee[16];
int  grabar_cada=1000; //Cada 1000 primos se graba la información en el eeeprom
int  grabar_van =0;
float  tiempo_min;
float  tiempo_min0;

void setup() 
{ 
  lcd.begin(16, 2);
   
  
  leido=EEPROMReadu64(0);
  found=EEPROMReadFloat(30);
  tiempo_min0=EEPROMReadInt(40);
  
  //leido =  1999999999999999ULL;
  //found=0;
  //tiempo_min0=0;
  
  lcd.setCursor(0, 0);
  printll(leido ); 
  candidate=leido;
  lcd.setCursor(0, 1);
  lcd.print((int)tiempo_min0 ); 
  lcd.setCursor(9, 1);
  lcd.print(found); 
  delay(5000);
   
  pinMode(led, OUTPUT);   
}  // end setup()

void loop() 
{
  uint64_t nextPSP;
  int nextPSPindex;
  for (int i=0;i<nPSP;i++) {
    nextPSP = fetchPSP(i);
    if ((nextPSP == candidate) || (nextPSP > candidate)) {
      nextPSPindex = i;
      break;
    }
    nextPSP = 0;
    nextPSPindex = 0;
  }
  start = millis();
  
  if ((candidate & 1) == 0) {
    candidate++;
  }
  
  for ( ; candidate < 9999999999999999ULL; candidate += 2) {
    if (nextPSP != candidate) {
      if (isPrime(candidate)) {
        found++;
        showPrime (candidate);
        timerRuntime= 0;
      }
    }
    else {
      nextPSPindex++;
      if (nextPSPindex < nPSP) {
        nextPSP = fetchPSP(nextPSPindex);
      }
    }
  }
  found = 1;
  showPrime(1);
  found++;
  showPrime(2);
  candidate = 3;  // back to start
}  // end loop()

uint64_t mulmodBarrett(uint64_t a, uint64_t b, uint64_t c) {
  // Calculate a*b %c

  uint64_t ah; 
  uint64_t al;
  uint64_t ph; 
  uint64_t pm;
  uint64_t pl;

  multll(a, b, &ph, &pl);
  ah = ph;
  al = pl;

  if (ah == 0) {
    if (al < c) {
      return al;
    }
    else {
      if (al == c) {
        return 0;
      }
      else {
        al %= c;
        return al;
      }
    }
  }
  
  mult4lx2l(ah, al, Barr_m, &ph, &pm, &pl);
  uint8_t Bk = Barr_k << 1;
  if (Bk > 64) {
    if (Bk > 128) {
      Bk -= 128;
      pl = ph >> Bk;
    }
    else {
      Bk -= 64;
      pl = (pm >> Bk) | (ph << (64 - Bk));
    }
  }
  else {
    pl = (pl >> Bk) | (pm << (64 - Bk));
  }      
      
  multll(pl, c, &ph, &pm);
  pm = al - pm;
  if (pm > ah) {
    ph--;
  }
  ph = ah - ph;
  if (pm > c) {
    pm -= c;
  }
  return pm;
} // end mulmodBarrett()

void mult4lx2l(uint64_t ah, uint64_t al, uint64_t b, uint64_t *prodh, uint64_t *prodm, uint64_t *prodl) {
  // Multiply a and b
  multll(ah, b, prodh, prodm);
  uint64_t prodt;
  multll(al, b, &prodt, prodl);
  *prodm += prodt;
  if (*prodm < prodt) {
    (*prodh)++;
  }
}  

bool Miller(uint64_t source)
{
  if(source <= 3) {
    return true;
  }
  if((source & 1) == 0) {
    return false;
  }

  uint64_t d = source - 1;
  int s = 0;
  while((d & 1) == 0)
  {
    d >>= 1;
    s ++;
  }

  uint8_t a;
  // si se incluyen los demas primos se incrementa en 500ms 
  const byte bases[] = {2, 3, 5, 7 };// 11, 13, 17, 19, 23 };
  const byte nbases = sizeof(bases)/sizeof(bases[0]);
  for(int i = 0; i < nbases; i++)
  {
    a = bases[i];
    uint64_t x = powmodll(a, d, source);
   if((x == 1) || (x == (source - 1))) {
      continue;
    }
    bool stillOK = false;
    for(int r = 1; r < s; r++)
    {
      x = mulmodBarrett(x, x, source);
      if(x == 1)
        return false;
      if(x == (source - 1)) {
        stillOK = true;
        break;
      }
    }
    if(!stillOK)
      return false;
  }
  return true;
}  // end Miller()

bool isPrime(uint64_t n) {
  if (!check(n)) {
    return false;
  }
  Barr_k = Barrett_k(candidate);
  Barr_m = Barrett_m(candidate,Barr_k);
  if (n > 0x10000) {
    if (!Miller (n)) {
      return false;
    }
  }
  return true;
}

void showPrime (uint64_t n) 
{
  digitalWrite(led, HIGH);
  printll(n);
  grabar_van++;
  tiempo_min=tiempo_min0+(millis() - start)/60000.0;
  if (grabar_van==grabar_cada) {
    EEPROMWriteu64(0, n);  
    EEPROMWriteFloat(30, found);
    EEPROMWriteInt(40, (int)tiempo_min );
    lcd.setCursor(15, 1);
    lcd.print("*");
    grabar_van=0;
  }  else {
    lcd.setCursor(15, 1);
    lcd.print(" ");
  }
  lcd.setCursor(0, 1);
  lcd.print((int)tiempo_min );
  lcd.setCursor(4, 1);  lcd.print("m Van ");  
  lcd.setCursor(10, 1);
  lcd.print(found ); 
  digitalWrite(led, LOW);
  //start = millis();
}  // end showPrime()

bool check(uint64_t n) {
  if (n <= 3) {
    return true;
  }
  if (!(n & 1)) {
    return false;
  }
  uint8_t i = 1;
  uint8_t f = pgm_read_byte_near(primes);
  while ((i <= nPrimes) && (f < n)) {
    if ((n % f) == 0) {
      return false;
    }
    f = pgm_read_byte_near(primes+i);
    i++;
  }
  return true;
} // end check()

uint64_t fetchPSP(int i) {
  i <<= 3;
  uint64_t n = 0;
  for (int8_t k=0;k<sizeof(n);k++) {
    n = n <<= 8;
    n |= pgm_read_byte_near(psps+i+k);
  }
  return n;
} // end fetchPSP()

void printll(uint64_t n) {
  // print an unsigned long long
  if (n == 0) {  
    return;
  }
  const byte DIGITS = 21;
  char s[DIGITS+1];
  for (byte i=0;i<DIGITS;i++) {
    s[i] = 'X';
  }
  s[DIGITS] = 0;
  byte i = DIGITS-1;
  while (n > 0) {
    s[i] = (n % 10) | '0';
    n = n / 10;
    i--;
  }
  i++;
  char *z = &s[i];
  
  lcd.clear();
  lcd.setCursor(0, 0);
  lcd.print(z );
   
  // Palindrome
  
} // end printll()

uint8_t Barrett_k(uint64_t n) {
  uint8_t k = 0;
  uint8_t t;
  while (n > 0xFF) {
    n >>= 8;
    k += 8;
  }
  t = n;
  while (t > 0) {
    t >>= 1;
    k++;
  }
  return k;
}

uint64_t Barrett_m(uint64_t n, uint8_t k) {
  uint8_t t;
  uint64_t iml, imh;
  
  if (k<32) {
    t = 2 * k;
  }
  else {
    t = 2 *k - 64;
  }
  iml = 1;
  for (uint8_t i=8;i<=t;i+=8) {
    iml <<= 8;
  }
  t &= 7;
  for (uint8_t i=0;i<t;i++) {
    iml <<= 1;
  }
  imh = 0;
  if (k>=32) {
    imh = iml;
    iml = 0;
  }
  return quotll(n, imh, iml);
}
  
uint64_t powmodll(uint8_t a, uint64_t b, uint64_t c) {
  // Calculate a**b % c,
  uint64_t x = 1;
  uint64_t y = a;

  while(b > 0) {
    if(b & 1) {
        x = mulmodBarrett(x, y, c);
      }
      y = mulmodBarrett(y, y, c);
      b >>= 1;
    }
    return x;
} // end powmodll()

void multll(uint64_t a, uint64_t b, uint64_t *imh, uint64_t *iml) {
  // Multiply a and b

  *imh = (a >> 32) * (b >> 32);
  uint64_t imm = (a & 0xffffffffUL) * (b >> 32);
  *iml = (a >> 32) * (b & 0xffffffffUL);
  imm += *iml;
  if (imm < *iml) {     // Check for carry
    *imh += 1ULL << 32;
  }
  *iml = (a & 0xffffffffUL) * (b & 0xffffffffUL);
  imm += *iml >> 32;
  if (imm < (*iml >> 32)) {     // Check for carry
    *imh += 1ULL << 32;
  }
  *iml = (*iml & 0xffffffffUL) + (imm << 32);
  *imh += imm >> 32;
}

uint64_t quotll(uint64_t c, uint64_t imh, uint64_t iml) {
  uint64_t thisQuotient;
  
  if (imh == 0) {
    return (iml / c);
  }
  uint64_t quotient = 0;
  for (uint8_t i=0;i<8;i++) {
    if ((imh >> 56) != 0) {
      if (imh > c) {
        thisQuotient = (imh / c);
        quotient |= thisQuotient;
        imh = imh - (thisQuotient * c);
      }
      else {
        if (imh == c) {
          quotient |= 1;
          imh = 0;
        }
      }
    }
    imh = imh << 8;
    imh |= iml >> 56;
    iml = iml << 8;
    quotient = quotient << 8;
  }
  if (imh > c) {
    quotient |= (imh / c);
  }
  return quotient;
}

void printllHex(uint64_t n) {
  char t[19] = "0x0000000000000000";
  char s[] = "0123456789ABCDEF";

  for (int8_t i=17;i>1;i--) {
    t[i] = s[n & 0xf];
    n = n >> 4;
  }
} 

void startTimer() {
  timerStart = micros();
}

void accTimer() {
  timerRuntime += micros() - timerStart;
}

void EEPROMWriteu64(int address, uint64_t value)
       {
       byte eight = (value & 0xFF);
       byte seven = ((value >> 8)  & 0xFF);
       byte six   = ((value >> 16) & 0xFF);
       byte five  = ((value >> 24) & 0xFF);
       byte four  = ((value >> 32) & 0xFF);
       byte three = ((value >> 40) & 0xFF);
       byte two   = ((value >> 48) & 0xFF);
       byte one   = ((value >> 56) & 0xFF);

       EEPROM.write(address, eight);
       EEPROM.write(address + 1, seven);
       EEPROM.write(address + 2, six);
       EEPROM.write(address + 3, five);
       EEPROM.write(address + 4, four);
       EEPROM.write(address + 5, three);
       EEPROM.write(address + 6, two);
       EEPROM.write(address + 7, one);
       }

uint64_t EEPROMReadu64(long address)
       {
       uint64_t eight = EEPROM.read(address);
       uint64_t seven = EEPROM.read(address + 1);
       uint64_t six   = EEPROM.read(address + 2);
       uint64_t five  = EEPROM.read(address + 3);
       uint64_t four  = EEPROM.read(address + 4);
       uint64_t three = EEPROM.read(address + 5);
       uint64_t two   = EEPROM.read(address + 6);
       uint64_t one   = EEPROM.read(address + 7);

       return ((eight << 0)  & 0xFF) + 
              ((seven << 8)  & 0xFFFF) + 
              ((six << 16)   & 0xFFFFFF) + 
              ((five << 24)  & 0xFFFFFFFF) +
              ((four << 32)  & 0xFFFFFFFFFFULL) + 
              ((three << 40) & 0xFFFFFFFFFFFFULL) +
              ((two << 48)   & 0xFFFFFFFFFFFFFFULL) + 
              ((one << 56)   & 0xFFFFFFFFFFFFFFFFULL)
       ;
       }
       
void EEPROMWriteFloat(int p_address, long p_value)
{
  byte Byte1 = ((p_value >> 0) & 0xFF);
  byte Byte2 = ((p_value >> 8) & 0xFF);
  byte Byte3 = ((p_value >> 16) & 0xFF);
  byte Byte4 = ((p_value >> 24) & 0xFF);

  EEPROM.write(p_address, Byte1);
  EEPROM.write(p_address + 1, Byte2);
  EEPROM.write(p_address + 2, Byte3);
  EEPROM.write(p_address + 3, Byte4);
}       

float EEPROMReadFloat(int p_address)
{
  byte Byte1 = EEPROM.read(p_address);
  byte Byte2 = EEPROM.read(p_address + 1);
  byte Byte3 = EEPROM.read(p_address + 2);
  byte Byte4 = EEPROM.read(p_address + 3);

  long firstTwoBytes = ((Byte1 << 0) & 0xFF) + ((Byte2 << 8) & 0xFF00);
  long secondTwoBytes = (((Byte3 << 0) & 0xFF) + ((Byte4 << 8) & 0xFF00));
  secondTwoBytes *= 65536; // multiply by 2 to power 16 - bit shift 24 to the left

  return (firstTwoBytes + secondTwoBytes);
}

void EEPROMWriteInt(int p_address, int p_value)
      {
      byte lowByte = ((p_value >> 0) & 0xFF);
      byte highByte = ((p_value >> 8) & 0xFF);

      EEPROM.write(p_address, lowByte);
      EEPROM.write(p_address + 1, highByte);
      }


unsigned int EEPROMReadInt(int p_address)
      {
      byte lowByte = EEPROM.read(p_address);
      byte highByte = EEPROM.read(p_address + 1);
      return ((lowByte << 0) & 0xFF) + ((highByte << 8) & 0xFF00);
      }

En la imágen podemos ver la información mostrada en la pantalla del lcd : en la primera línea tenemos el último número primo encontrado, en la segunda línea tenemos durante cuantos minutos ha estado ejecutándose el programa y el número de primos encontrados.

Una mejora al programa sería lograr que se grabaran los primos encontrados en una memoria eeprom externa o en una tarjeta de memoria. Ya veremos ...

Actualización del 09/03/2015
En la imágen podemos ver que después de 9768 minutos ( 163 horas) se han encontrado 731.696 números primos.

viernes, 27 de febrero de 2015

Leer y escribir variables del tipo uint64_t en el eeeprom ( unsigned long long )

La libreria eeeprom original de Arduino permite leer y escribir solamente variables del tipo byte . En caso de querer utilizar variables de otro tipo tenemos que segmentar estas variables en bytes.
En la página arduino.cc hay código adicional que permite utilizar variables de diferentes tipos.Uno de estos programas , muy sencillo por cierto, permite trabajar con variables del tipo long.
A partir de este código se pueden generar subrutinas para trabajar con variables de 64 bits (unsigned long long , o variables de 64 bits ).

En el recuadro tenemos un programa sencillo con un lazo que escribe y lee en la eeeprom utilizando dichas subrutinas:

#include <EEPROM.h>

uint64_t candidate;
uint64_t leido;

void setup() {
  Serial.begin(115200); 
  while (!Serial) { }
  candidate= 1999999999999999ULL; 
}


void loop() {
   candidate=candidate + 1;   
   EEPROMWriteu64(0, candidate); 
   leido=EEPROMReadu64(0);
   printll(leido );  
   delay(5000);
}

void printll(uint64_t n) {
  // print an unsigned long long
  if (n == 0) {   
    return;
  }
  const byte DIGITS = 21;
  char s[DIGITS+1];
  for (byte i=0;i<DIGITS;i++) {
    s[i] = 'X';
  }
  s[DIGITS] = 0;
  byte i = DIGITS-1;
  while (n > 0) {
    s[i] = (n % 10) | '0';
    n = n / 10;
    i--;
  }
  i++;
  char *z = &s[i];
   Serial.println (z);
} // end printll()

void EEPROMWriteu64(int address, uint64_t value)
       {
       byte eight = (value & 0xFF);
       byte seven = ((value >> 8)  & 0xFF);
       byte six   = ((value >> 16) & 0xFF);
       byte five  = ((value >> 24) & 0xFF);
       byte four  = ((value >> 32) & 0xFF);
       byte three = ((value >> 40) & 0xFF);
       byte two   = ((value >> 48) & 0xFF);
       byte one   = ((value >> 56) & 0xFF);

       EEPROM.write(address, eight);
       EEPROM.write(address + 1, seven);
       EEPROM.write(address + 2, six);
       EEPROM.write(address + 3, five);
       EEPROM.write(address + 4, four);
       EEPROM.write(address + 5, three);
       EEPROM.write(address + 6, two);
       EEPROM.write(address + 7, one);
       }

uint64_t EEPROMReadu64(long address)
       {
       uint64_t eight = EEPROM.read(address);
       uint64_t seven = EEPROM.read(address + 1);
       uint64_t six   = EEPROM.read(address + 2);
       uint64_t five  = EEPROM.read(address + 3);
       uint64_t four  = EEPROM.read(address + 4);
       uint64_t three = EEPROM.read(address + 5);
       uint64_t two   = EEPROM.read(address + 6);
       uint64_t one   = EEPROM.read(address + 7);

       return ((eight << 0)  & 0xFF) + 
              ((seven << 8)  & 0xFFFF) + 
              ((six << 16)   & 0xFFFFFF) + 
              ((five << 24)  & 0xFFFFFFFF) +
              ((four << 32)  & 0xFFFFFFFFFFULL) + 
              ((three << 40) & 0xFFFFFFFFFFFFULL) +
              ((two << 48)   & 0xFFFFFFFFFFFFFFULL) + 
              ((one << 56)   & 0xFFFFFFFFFFFFFFFFULL)
       ;
       }  <br>

Se incluye otra subrutina que permite imprimir en el puerto serial las variables leidas (original de Nick Gammon).

domingo, 5 de octubre de 2014

Simulación de controlador PID en micro 8051

Los controladores PID (Proporcional/Integral/Derivativo) (wikipedia) son muy utilizados para controlar procesos de todo tipo.

Existe mucha documentación en la web donde se explica la teoría asociada y desarrollos tanto analógicos como discretos.

Simplificando podemos decir que la idea es medir un valor y compararlo con el valor deseado. En base a la diferencia de ambos valores (error) se genera una señal de ajuste para lograr que el valor medido se aproxime al valor deseado. El valor de ajuste se calcula mediante tres componentes :

Proporcional : Basado en el error actual (presente)
Integral : Basado en el historial de errores (pasado)
Derivativo : Basado en la proyección del error (futuro)

Vajuste = Kp * Proporcional + Ki * Integral + Kd * Derivativo

Las constantes de cada factor de la fórmula se seleccionan para determinar la curva de respuesta mas adecuada al proceso a controlar.

Cuando queremos utilizar un micro como controlador PID se utiliza un algoritmo discreto (se mide a ciertos intervalos de tiempo el valor de salida, se calcula el ajuste y se aplica la corrección).

A continuación tenemos un programa en pascal que permite simular el algoritmo PID en un micro de 8bits:

//   PID controller simulation

//  Serial configuration
//  Baud  Bits Parity Stop FlowControl
//  xxxx  8    N      1    None 
//  Led conected to port P0.0
//  20MHz -> 7200*2 = 14400 bauds

Program PID;
{$O+} 

Const
  Osc = 20000000; 
  BaudRateTimerValue = Byte (- Osc div 12 div 32 div 7200);
  UMAX= 0.2;   //Superior range
  UMIN=-0.2;   //Inferior range 
 
Var 
    SerialPort: Text;
   
    DelayTimer: Word; Volatile;
    control,set_point,measured_value,actual_error,error_previous:real;
    Pc,Ic,Dc:real;
    Kp,Ki,Kd:real;
    r:real;
    
Procedure Delay (DelayTime: Word);
var i,j:word;
begin
  for i:= 0 to DelayTime do   begin
    for j:= 0 to 32000 do begin end;
  end;
end;    

Procedure WriteToSerialPort; Assembler;
Asm
  CLR   TI
  MOV   SBUF, A
  @WaitLoop:
  JNB   TI, @WaitLoop
end;
 
Function ReadFromSerialPort: Char;
Var ByteResult: Byte absolute Result;
begin
  While not RI do;
  RI := False;
  ByteResult := SBUF;
  { Echo character }
  Asm
    MOV    A, Result
    LCALL  WriteToSerialPort
  end;
end;
 
Procedure Init;
begin
  TL1  := BaudRateTimerValue;
  TH1  := BaudRateTimerValue;
  TMOD := 100001;     
  SCON := %01010000;     
  TI   := True;         
  TR1  := True;  
  PCON := PCON or $80; //SMOD=1  por lo que dobla el baud de la fórmula inicial
  Assign (CurrentIO, ReadFromSerialPort, WriteToSerialPort);         
end;
 
{$DefaultFile On }

begin
  Init;
  Randomize;
  
  Writeln ('Discrete PID simulation');
  Kp:=10;
  Ki:=1;
  Kd:=1;
  set_point:=10.5;      //Desired Value
  measured_value:=5.5;  //Initial value
  error_previous:=set_point - measured_value;  //Initial error
  Ic:=0;
  r:=0;
  Writeln ('Desired Value:',set_point:10:5);
  Repeat
    actual_error   := set_point - measured_value; 
   
    Pc := actual_error;   //Current error
    Ic := Ic + error_previous;  //Sum of previous errors
    Dc :=  error_previous - actual_error ;  //Difference with previous error
    control:= Kp*Pc + Ki*Ic + Kd*Dc; 
    if (control > UMAX) then control := UMAX;    
    if (control < UMIN) then control := UMIN;
    error_previous := actual_error;  //error_previous holds the previous error
    writeln ('Measured Value=',measured_value:10:5,' Control=',control:7:3, 
             '  Pc:',Pc:10:5,' Ic:',Ic:10:5,' Dc:',Dc:10:5,' R:',r:10:5);
    delay(2);
    r:=random(100)/500;   // random interference
    measured_value:= measured_value + control + r ;
  until False;
end.

En la imagen se muestra la salida de la simulación donde se puede observar como el valor medido se va incrementando:

En la siguiente imagen podemos ver que el valor medido ya se aproximó al valor deseado, desplazándose hacia arriba y hacia abajo dependiendo del valor de interferencia r .

Es interesante ver como los factores Proporcional/Integral/Derivativo varían a lo largo de la simulación.

Algunos enlaces interesantes:

Páginas