Algunas notas sobre “Scaling Genetic Algorithms using MapReduce”

En esta publicación, los autores (Abhishek Verma, Xavier Llor`, David E. Goldberg y Roy H. Campbell) porponen la utilización del framework MapReduce (en particular con la implementación de Hadoop) para implementar GA paralelos. Fijan posición en que la implementación tradicional basada en MPI requiere mucho conocimiento sobre la máquina donde se monta, y que, por el contrario, MapReduce permite un nivel de escalabilidad que no introduce ningún cuello de botella.

En particular, se enfocan en las variantes tradicionales y sencilla de los GA para brindar una implementación de referencia.

De los pasos del algoritmo clásico genético, la evaluación de los fitness tiene una analogía directa con la función Map, la cuál evaluará el fitness del individuo y guardará track (en un counter) del mejor individuo, finalmente persistiéndolo en el HDFS, el proceso cliente (el que inició la tarea MapReduce) chequea recurrentemente en este archivo si se ha alcanzado el criterio de convergencia.

Para poder implementar la selección distribuida, hay que tener control o mantener estado global, lo cual es un tanto complejo en MapReduce, por lo que el único componente que puede realizar este tipo de tarea es el Partitioner. El problema puntual por el cual no conviene utilizar la implementación por defecto es que durante la etapa de Shuffling (el momento en el que el framework reparte, mezcla y distribuye las emisiones de los mappers a los reducers) es porque utiliza una función de hash sobre la key (módulo la cantidad de reducers) para distribuir la carga, generando inconvenientes de convergencia del GA. Por otro lado, mientras el algoritmo genético pogrese, muchos individuos serán iguales, por lo que serán tratados por el mismo reducer, con la reducción de performance que esto implicaría. Es por eso que los autores utilizaron un partitioner especial que en vez de utilizar un hash, utiliza un número pseudoaleatorio cuya semilla es (mapperId * time).

El mecanismo de selección que han utilizado es el de torneo, ya que es el único que puede implementarse con facilidad en MapReduce, ya que la ruleta implicaría el mantenimiento de un estado global y un par de iteraciones más. Para esto han decidido tomar individuos al azar y ver quien es el ganador n veces (n es el tamaño de la población) del torneo para someterlo a la cruza. La operación de crossover se realiza en conjunto con la selección de dos individuos.

El tipo de crossover que utilizaron es el uniforme.

Cuando las poblaciones son muy grandes, la inicialización serial de los individuos tomó mucho tiempo, por lo que decidieron agregar una tanda de procesamiento map reduce en el que el map genera un individuo al azar y el reduce es trivial.

Las medidas de performance que realizaron fueron hechas con la implementación de un problema trivial utilizado para benchmark llamado OneMax Problem que consiste en encontrar la máxima cantidad de 1’s en los cromosomas de los individuos, es decir, un individuo tal que todos sus alelos sean 1’s. El análisis realizado arrojó que:

  • Análisis de convergencia: Para 104 variables, converge en 220 iteraciones con un promedio de 149 segundos por iteración.
  • Escalabilidad con carga constante por nodo: Incrementando la cantidad de recursos, no cambia el tiempo de procesamiento por iteración.
  • Escalabilidad con carga constante general: Fijaron la cantidad de variables en 50.000 e incrementaron el número de mappers, logrando una reducción del tiempo de iteración. Aunque el speed-up general queda limitado por la ley de Amdahl (“la mejora obtenida en el rendimiento de un sistema debido a la alteración de uno de sus componentes está limitada por la fracción de tiempo que se utiliza dicho componente”) dado el overhead de Hadoop (unos 10 segundos para iniciar y terminar un job).
  • Escalabilidad con incrementos en el tamaño del problema: Utilizaron el máximo de recursos y fueron incrementando el número de variables, logrando un incremento razonable con incrementos de población superlineales (n log n).

El gran inconveniente de trabajar con MPI en clústers standard es la poca tolerancia a fallos del framework. MapReduce no tiene este inconveniente. Introducen ciertas críticas a la implementación llamada MRPGA respecto a la escalabilidad.

Otros enfoques han tratado de modificar MapReduce para que se parezca a algo más consumible desde los GA, lo interesante de este trabajo es que no han intentado martillar MapReduce, sino todo lo contrario.

Como línea de investigación proponen comparar rendimientos con MPI y demostrar la importancia de GA escalables en las aplicaciones prácticas.

Por otro lado cabe destacar que no se proveen detalles de la implementación a bajo nivel, de modo que mucho se puede conjeturar sobre la implementación específica que han realizado, y que no se han aventurado a explorar alternativas un poco más complejas que el GA básico.

Algunas notas sobre “A Survey of Parallel Genetic Algorithms”

En esta publicación de 1998, Erick Cantú-Paz hace hincapié en la facilidad de implementar GA paralelos y en sus ganancias en rendimiento. Para eso intenta colectar, organizar y presentar las publicaciones más representativas en este aspecto. Sin embargo se va a centrar en las investigaciones en donde los GA tengan múltiples poblaciones independientes.

Haciendo un paneo por la historia de los enfoques paralelos en GA, destaca el trabajo de Betke de 1976, en donde describió una implementación paralela de GA convencionales. Otro estudio interesante fue el de Grefenstette de 1981 que propuso mapear las arquitecturas paralelas existentes con GA, trabajando en 4 prototipos distintos, siendo uno de ellos el primer prototipo del enfoque de poblaciones múltiples, siendo el prototipo que más tardíamente se implementó (1987).

Luego de una breve introducción describiendo someramente el mecanismo de los GA, pasa a clasificar los GA paralelos:

    1. Global single-population master-slave GA: Donde existe una única población, pero el procesamiento de la evaluación del fitness es altamente distribuida. Cada individuo puede competir y cruzarse con cualquier otro. Se dice que el algoritmos es sincrónico si debe esperar la evaluación de todos los fitness (procesados en paralelo). Hay posibilidades de implementación asincrónica, pero en esos casos el mecanismo de GA no funciona de manera tradicional. La mayoría de las implementaciones son sincrónicas y reducen la velocidad de evaluación al tiempo de evaluación más lento de cada generación. Pueden ser implementados en modelos de memoria compartida, como de memoria distribuida. Cabe destacar que para paralelizar el proceso de selección habrá obstáculos por el overhead de comunicación, aunque Branke en 1997 paralelizó distintos tipos de selección global en una grilla 2D de procesadores mostrando que sus algoritmos eran óptimos para la topología utilizada (O(\sqrt{n}) en una grilla de \sqrt{n} x \sqrt{n}). La cruza y la mutación pueden ser fácilmente paralelizables utilizando la misma idea que para la evaluación del fitness. En conclusión, son fácilmente implementables y pueden ser un método eficiente y se puede aplicar la teoría disponible para los GA tradicionales.
    2. Single-population fine-grained GA: Consiste en una población espacialmente estructurada, en donde la selección y cuza se restringen a una pequeña vecindad del individuo. Manderick y Spiessens (1989) observaron que el rendimiento se degradaba a medida que el tamaño de la vecindad aumentaba. También pudieron demostrar que para un tamaño de subpoblaciones s y un largo de cromosoma l, la complejidad de tiempo del algoritmo es O(s+l) u O(s(log s)+l) pasos dependiendo del esquema de selección utilizado. Una conclusión interesante es que no se puede determinar una regla general para comparar el método fine-grained con el coarse-grained, aunque en 1992 Gordon, Whitley y Böhm demostraron que el camino crítico de un GA fine-grained es más corto que el de GA de población múltiple, aunque no incluyeron consideraciones de comunicaciones.
    3. Multiple-population (multiple-deme) coarse-grained GA: Son los más sofisticados, dado que consisten en varias subpoblaciones que intercambian individuos entre si eventualmente (migración). Cabe destacar que los efectos de la migración todavía no han sido comprendidos enteramente. Introduce cambios en el esquema de funcionamiento de los GA tradicionales. Tienen una analogía con el enfoque llamado “modelo de isla” de la genética de poblaciones. Una observación interesante sobre los estudios realizados en este ambiente relevan ciertas preguntas de interés: ¿Cuál es el ancho de banda de comunicaciones necesario para hacer que se comporten como GA global?¿Cuál es el costo de estas comunicaciones?¿Es el costo de las comunicaciones lo suficientemente pequeño como para hacer de esta una alternativa viable para el diseño de GA paralelos? Otro de los parámetros estudiados (particulamente en la investigación de Tanese en donde comparaba el rendimiento de un GA serial con respecto a la implementación paralela con y sin costo de comunicación) es el de la frecuencia de migración de individuos. En la misma investigación encontró que una implementación de múltiple población (r poblaciones de n individuos, sin costo de comunicación) encontraba soluciones de la misma calidad que la implementación serial con una población de r.n individuos. Los parámetros que controlan la migración de individuos son:
  • La topología que define las conexiones entre las subpoblaciones.
  • Una tasa de migración que controla cuántos individuos deben migrar
  • Un intervalo de migración que afecta a la frecuencia de las migraciones

Una de las conclusiones tiene una implicancia biológica notable: hay cambios favorables en los cromosomas que se extienden más rápido en poblaciones más pequeñas que en las más grandes.

  1. GA Paralelos Jerárquicos: Es una combinación de múltiples poblaciones y procesamiento maestro-esclavo. Dentro de esta rama hay investigadores que han combinado métodos. Este enfoque puede reducir el tiempo de ejecución más que cualquiera de las otras técnicas por si sola.

Probablemente el primer intento de fundamentar teóricamente el rendimiento de un GA paralelo fue una publicación de Pettey y Leuze de 1989, en donde se describe una derivación del teorema de los esquemas para GA paralelos. Una cuestión teórica muy importante es poder determinar si (y bajo qué condiciones) un GA paralelo puede encontrar una mejor solución que un GA serial. Hay dos observaciones en el trabajo de Starkweather (et al.):

  • Las subpoblaciones relativamente aisladas pueden converger a diferentes soluciones, y la migración y recombinación pueden combinar soluciones parciales. Se especula que por esto los GA paralelos con baja tasas de migración pueden encontrar mejores soluciones para funciones separables.
  • Si la recombinación de soluciones parciales resulta en individuos con fitness más pobre (quizás porque la función no sea separable), entonces los GA seriales pueden tener una ventaja.

La migración se considera sincrónica cuando en un determinado momento migran todos los individuos que se esperaba que migren. Hay un enfoque asincrónico. Hay enfoques que utilizan la migración asincrónica tanto cuando las subpoblaciones tienden a converger como cuando convergen totalmente. De esta manera surge una pregunta muy interesante: ¿Cuál es el mejor momento para migrar?

Topología de las comunicaciones: Es un factor fundamental en los GA paralelos. Si la topología tiene una conectividad densa, se supone que las buenas soluciones se expandirán más rápido a lo largo de las poblaciones. Por otro lado, si la topología es más rala, lo harán más lento, permitiendo la aparición de diferentes soluciones. La línea general de implementación dice que la topología de comunicación debe ser estática a lo largo de toda la ejecución. El enfoque de topología dinámica es utilizado cuando puede determinarse a qué subpoblación conviene más que migre un individuo.

Entre los avances más recientes para la fecha de publicación, el autor señala:

  • Se concluyo que la cantidad óptima de workers para el enfoque GA paralelo Global es S*=\sqrt{\dfrac{nT_{f}}{T_{c}}} donde n es el tamaño de la población, Tf es el tiempo que toma hacer una evaluación de fitness y Tc es el tiempo de comunicación. El speed-up óptimo es 0.5 S*.
  • Cuando las subpoblaciones están aisladas, el speed-up no es tan significante.
  • El speed-up es mucho mejor cuando las subpoblaciones están comunicadas.
  • Se puede determinar un número óptimo de subpoblaciones (asociado al tamaño de las mismas) que maximiza el speed-up.
  • Hay varios problemas que todavía no han sido resueltos, entre ello:
      • Determinar la tasa de migración que hace que las subpoblaciones distribuidas se comporten como una única población global.
      • Determinar la topología de comunicaciones adecuada que permita la combinación de buenas soluciones, pero que no resulte en un costo excesivo de comunicación.
      • Encontrar si hay un número óptimo de subpoblaciones que maximice la confiabilidad del GA.

Algunas notas sobre “Genetic Algorithms: Concepts and Applications”

En este artículo publicado en 1996 en una publicación de la IEEE, los autores, K.F. Man, K.S. Tang y S.Kwong se propusieron dar una base conceptual de algoritmos genéticos para presentárselo a ingenieros. El enfoque ingenieril sobre los GA no sólo se basa en las aplicaciones industriales posibles, sino que también, en su fundamento, intenta conformar un cuerpo de conocimiento técnico que sirva como marco de trabajo para herramientas de diseño en la ingeniería industrial.

Se hace hincapié en el trabajo fundacional de J.H.Holland de la década de 1970, enfocándose en los problemas intratables bajo otros mecanismos.

Para dar un breve pantallazo sobre los conceptos básicos, se toma el algoritmo genético clásico, en el que los cromosomas están formados por arreglos de bits, divididos en genes que expresan las características de la solución que está siendo codificada por el individuo/cromosoma.

Un valor de fitness es calculado para cada individuo utilizando la función de fitness asignada al problema, que indica cuán buena es la solución evaluada para resolver el problema en cuestión. El modelo que utiliza para representar la selección es el de la ruleta de probabilidades.

A través de la evolución genética, los individuos tienen la tendencia a mejorar su desempeño a lo largo de las generaciones.

El ciclo básico del algoritmo genético se expresa en el siguiente pseudocódigo:

Standard Genetic Algorithm ()
{
  // start with an initial time
  t := 0;
  // initialize a usually random population of individuals
  initpopulation P (t);
  //evaluate fitness of all initial individuals of population
  evaluate P (t);
  // test for termination criterion (time, fitness, etc.)
  while not done do
    // increase the time counter
    t := t + 1;
    // select a sub-population for offspring production
    P' := selectparents P (t);
    //recombine the ”genes”of selected parents
    recombine P (t);
    // perturb the mated population stochastically
    mutate P (t);
    //evaluate it’s new fitness
    evaluate P (t);
    // select the survivors from actual fitness
    P := survive P,P‘ (t);
  od
}

La selección de los parámetros de probabilidad de mutación y crossover es en si un problema complejo, un problema de optimización no lineal. De hecho, su configuración es dependiente de la función objetivo. En este punto se cita a K.A. DeJong y W.M. Spears, en “An analysis of the interacting roles of population size and crossover in genetic algorithms” y a J.J. Grefenstette en “Optimization of control parameters for genetic algo-
rithms”, donde se sugiere que para una población grande, la probabilidad de cruza debiera ser del orden de 0,6 y la de mutación, de 0,001; siendo, para una población más pequeña, la probabilidad de cruza, de 0,9 y la de mutación, 0,01.

En el ámbito teórico, los autores citan al gran trabajo de Michalewicz (“Genetic Algorithms + Data Structures = Evolution Programs”), de allí toman la teoría de los esquemas. Se llama esquema a una estructura de máscara sobre el arreglo de bits del que están formados los cromosomas. Estas máscaras funcionan como los comodines de las expresiones regulares, o de los comandos para explorar los sistemas de archivos de los sistemas operativos.

El orden de estos esquemas va a estar dado por la cantidad de posiciones fijas (con 0 o 1) en el cromosoma (las que no tienen asterisco). Es decir, que el orden del esquema determina el grado de especificidad del mismo.

Se define la longitud del esquema como la diferencia entre la posición fija más alejada del mismo y la posición fija más cercana a la primer posición (es como restar el índice de un array).

Se denota como \zeta(S,t) a la cantidad de cromosomas matcheados por el esquema S en la generación t de la corrida del algoritmo genético en cuestión. Como f(S,t) se va a denotar al fitness promedio de los individuos que matchean con el esquema S en la generación t. Se llama L a la longitud de un cromosoma.

Finalmente, se va a denotar F(t) a la sumatoria de los fitness de todos los individuos de la generación t. Y es aquí donde creemos que el artículo incluye un pequeño error, dado que dice que F(t) denota el número promedio de ocurrencias de S, lo que no tendría mucho sentido, dado que ya fue estudiado.

Siguiendo adelante en el razonamiento, y tomando las operaciones de crossover en un punto y mutación, se puede obtener una ecuación de crecimiento de un esquema de la siguiente manera:

\zeta(S,t+1) >= \zeta(S,t) \dfrac {f(S,t)}{F(t)} . [1-p_{c} \dfrac{\delta (S)}{L-1} - o(S) p_{m}]

Y es aquí donde encontramos otra diferencia con el trabajo de Michalewicz, ya que en el mismo se utiliza la notación de valor medio o promedio para F(t), es decir la sumatoria de todos los fitness de los individuos de la generación t dividido su cantidad, llamémosla N (tamaño de la población). Por tanto la ecuación de crecimiento, debiera ser escrita como:

\zeta(S,t+1) >= \zeta(S,t) \dfrac {f(S,t).N}{F(t)} . [1-p_{c} \dfrac{\delta (S)}{L-1} - o(S) p_{m}]

Cabe aclarar que F(t)>0 y que L>1

Allí se cita el teorema de los esquemas que dice que los esquemas cortos, de bajo orden, que están por encima del promedio, reciben un crecimiento exponencial en la cantidad de matcheos en las subsiguientes generaciones de un algoritmo genético.

Luego se plantea la hipótesis de los bloques de construcción, en donde se afirma que un algoritmo genético busca rendimientos cercanos al óptimo de un problema a través de la yuxtaposición de esquemas cortos, de bajo orden y alto rendimiento, llamados bloques de construcción (building blocks).

Los building blocks me remiten mentalmente al concepto de meme, acuñado por Dawkins en “The Selfish gene”, a pesar de que no creo que la naturaleza sea egoísta, siguiendo a Chaitin, en mi humilde opinión la naturaleza es más indiferente que egoísta.

Una conclusión interesante es la siguiente: dado que el mecanismo de los GA no es gobernado ni por ecuaciones diferenciales, así como tampoco se comporta como una función continua, puede concebirse que esta estructura sencilla pueda ser utilizada como un optimizador en muchas aplicaciones prácticas. Tampoco cabe ninguna duda acerca de que los GA puedan ser utilizados en problemas en los que las técnicas de hill-climbing o basadas en gradientes presentan inconvenientes.

Aunque también se hace énfasis en las limitaciones de los GA y se sugieren varias líneas de investigación, como otras estructuras de representación de cromosomas por fuera de los arreglos de bits.

Se aclara que las predicciones de los GA pueden carecer de sentido o llevar a resultados erróneos en determinados problemas, refiriéndose al problema que se da cuando se abordan temáticas con funciones deceptivas.

Como modificaciones estructurales a los algoritmos genéticos se proponen diversidad de cambio en la representación de los cromosomas. También se proponen cambios sobre la escala de valores de la función de fitness (y función objetivo), nombrando como alternativas:

  • Escala lineal
  • Escala de potencia
  • Sigma Truncation

Respecto a los mecanismos de selección, además de la ruleta se porpone Stochastic Universal Sampling (SUS) donde la complejidad algorítmica es similar, esquemas de rankings y otros.

Para la operación de crossover se propone la cruza multipunto, la cruza uniforme (por máscara) y la PMX (Partially matched crossover).

Una forma diferente de incrementar el rendimiento del algoritmo genético es intentar reordenar los genes dentro del cromosoma para encontrar patrones que tengan un mejor potencial evolutivo.

Las estrategias de reinserción en la población también son tema de discusión, generalmente se utiliza el reemplazo generacional en donde las descendencias reemplazan totalmente a los individuos que los generaron. Existe una manera elitista de reemplazar individuos, en donde permanecen los dos mejores entre “padres” e “hijos”.

Otro punto calificado de controversial es el ajuste de las probabilidades de mutación y cruza.

En ciertos procesos dinámicos es difícil mantener estático el “entorno”. En el enfoque tradicional de algoritmos genéticos existen complicaciones para variar el entorno, es por eso que se han desarrollado extensiones como la posibilidad de mantener individuos diploides con una representación tri-alélica para decidir dominancia de genes, de modo de poder ser más flexible a la hora de cambiar el enterno en el que son evaluados los individuos. Otras técnicas para el cambio de entorno incluyen la migración de individuos para generar mayor diversidad, o la introducción en determinado momento de un mecanismo de hipermutación (un cambio radical en la tasa de mutación).

Respecto de la paralelización, se dice que la performance serial de los GA no es buena respecto de otros mecanismos de optimización o búsqueda, dado que no se basa en heurística alguna. Para eso se han desarrollado modelos de procesamiento en paralelo, siendo los principales:

  • Global: Hay una sola población y se procesa en paralelo, siendo más apto para los enfoques de paralelismo en memoria compartida.
  • Migración de individuos en donde hay islas de evolución que corren en paralelo y cada tanto migran individuos.
  • Difusión en donde cada individuo se le asigna una ubicación en una plantilla de 2 dimensiones y sólo puede cruzarse con individuos “cercanos”, este enfoque de difusión ha evolucionado al enfoque celular.

Uno de los nichos de utilización de GA es el de optimización multiobjetivo, dado que matemáticamente es muy complejo resolver este tipo de problemas. Se ha demostrado que los GA son particularmente aptos para encontrar conjuntos de soluciones Pareto-óptimas.

Los problemas más evidentes que se detectan en los GA son:

  • Deception: Existen funciones que son GA-deceptivas, en donde la combinación de building blocks positivos evolutivamente hablando no generan una descendencia mejor que sus antecesores.
  • Estancamiento en subóptimos: Dado que los GA no buscan óptimos globales, existe la posibilidad de que la evolución se estanque en un óptimo local y no pueda salir de allí. Existen diversas técnicas para evitar este tipo de comportamientos.
  • Tiempo real: Los GA no son una técnica aplicable fácilmente a problemas de tiempo real.

Finalmente se describen situaciones prácticas en las que los GA han sido exitosamente implementados, describiendo particularidades de cada implementación:

  • Parámetros e Identificación de Sistema (IIR)
  • Procesos de Control (por ejemplo proceso de neutralización de pH)
  • Robótica (navegación, cálculo de trayectorias)
  • Reconocimiento de patrones (imágenes, tanto en el proceso de extracción, como en el de reconocimiento en si)
  • Reconocimiento del habla
  • Diseño e ingeniería (VLSI)
  • Planificación y Scheduling (travelling salesman, pump scheduling, job-shop scheduling, etc)
  • Sistemas de clasificación (detección de pérdidas de gas y control de cañerías).

Tecnología emergente a partir de los algoritmos genéticos:

  • Integración con las redes neuronales: en donde los GA pueden ser buenos seleccionando las reglas de aprendizaje o los parámetros de las redes neuronales, pudiendo funcionar como buenos supervisores de aprendizaje.
  • Integración con lógica difusa: este punto es más complejo porque no hay una relación de feedback tan evidente como con las reglas neuronales.

Integrales por suma de Riemann: MPI & MapReduce

La Suma de Riemann es una técnica de integración numérica para calcular valores de integrales definidas (bajo ciertas condiciones, como la continuidad). Visualmente lo que estamos haciendo es calcular el área por “debajo” de la curva definida por la función en cuestión.

Para calcular esa superficie, este método divide el continuo del eje x en “barras” o rectángulos, como se ve en la imagen:

Riemann_Integration_3Estos rectángulos van desde el eje x hasta la curva, dependiendo el enfoque que se tome (sumas izquierdas, derechas, máximas, mínimas, etc.)

El error en este tipo de sumas es alto respecto a otros métodos, pero mientras las bases de los rectángulos sean más chicas (tiendan a cero), más converge la suma hacia el resultado de la integral.

Un programa en lenguaje C que hace ese cálculo para la función x2 en el intervalo [0;10] con 100 rectángulos (pasos):

#include <stdio.h>
#include <stdlib.h>

#define XINICIAL 0
#define XFINAL 10
#define PASOS 100

float calcularFuncion(float x);

int main (int argc, char *argv[])
{
      float resultado = 0.0;        
      float incremento = ((float)XFINAL-(float)XINICIAL)/(float)PASOS;
      float semiIncremento = incremento/2;
      float x = XINICIAL;
      int i;
      for (i=0;i<PASOS;i++)
      {
         resultado += incremento * (calcularFuncion(x)+calcularFuncion
                                             (x+incremento))/(float)2;
         x+=incremento;
      }

      printf("El resultado de es %f \n",resultado);
}

float calcularFuncion(float x)
{
  return x*x;
}

Un código realmente simple. Si quisiéramos procesarlo en un cluster utilizando MPI (OpenMP) en C, tenemos dos opciones:

  1. Que el nodo master le pase a los workers cada posición x del rectángulo a calcular, el worker calcule todos sus rectángulos asignados, sume sus superficies y devuelva un valor al nodo master.
  2. Que el nodo master divida la curva a asignar en tantas partes como workers disponibles tenga el cluster y les asigne una posición de inicio y final para que cada worker calcule con este método las superficies parciales y devuelva un valor al nodo master.

Veamos una implementación posible en C (con MPI) de la primer opción:

#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"

#define XINICIAL 0
#define XFINAL 10
#define PASOS 100
#define MASTER 0
#define TAG_X 1
#define TAG_RESULT_PARCIAL 2

float calcularFuncion(float x);

int main (int argc, char *argv[])
{
      int rank, value, size, dest,i,workers;  
      float resultado,incremento,x,resultParcial;
      MPI_Status status; 

      incremento = ((float)XFINAL-(float)XINICIAL)/(float)PASOS;

      MPI_Init( &argc, &argv );

      MPI_Comm_rank( MPI_COMM_WORLD, &rank );
      MPI_Comm_size( MPI_COMM_WORLD, &size ); 

      workers=size-1;

      if (size<2)
      {
    printf("No se puede procesar con %i nodos\n",size);
    return 0;
      }

      if (rank==MASTER)
      {
    dest=1;
    x=XINICIAL;
    // Hacemos un Round Robin
    for (i=0;i<PASOS;i++)
    {
      MPI_Send(&x, 1, MPI_FLOAT, dest, TAG_X, MPI_COMM_WORLD); 
      dest++;
      if (dest==size) dest=1;
      x+=incremento;
    }
    for (i=MASTER+1;i<size;i++)
    {
      MPI_Recv(&resultParcial, 1, MPI_FLOAT, i, TAG_RESULT_PARCIAL, 
                                          MPI_COMM_WORLD, &status); 
      resultado +=resultParcial;
    }

    printf("El resultado de es %f \n",resultado);

 }
 else
 {
    for (i=0; i<PASOS-(((int)(PASOS/workers))*workers)>=rank?
                                          (int)(PASOS/workers)+1:
                                          (int)(PASOS/workers); i++) 
    {
      MPI_Recv(&x, 1, MPI_FLOAT, MASTER, TAG_X, MPI_COMM_WORLD, 
                                                          &status); 
      resultParcial += incremento * (calcularFuncion(x)+
                            calcularFuncion(x+incremento))/(float)2;

    }
    MPI_Send(&resultParcial,1, MPI_FLOAT, MASTER, TAG_RESULT_PARCIAL,
                                                    MPI_COMM_WORLD);
      }      
}

float calcularFuncion(float x)
{
  return x*x;
}

Veamos una implementación posible en C (con MPI) de la segunda opción:

#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"

#define XINICIAL 0
#define XFINAL 10
#define PASOS 1000
#define MASTER 0
#define TAG_INICIO 1
#define TAG_FIN 2
#define TAG_RESULT_PARCIAL 3

float calcularFuncion(float x);

int main (int argc, char *argv[])
{
      int rank, value, size, i,workers;  
      float resultado,incremento,x,resultParcial,xfinal;
      MPI_Status status; 

      incremento = ((float)XFINAL-(float)XINICIAL)/(float)PASOS;

      MPI_Init( &argc, &argv );

      MPI_Comm_rank( MPI_COMM_WORLD, &rank );
      MPI_Comm_size( MPI_COMM_WORLD, &size ); 
      resultParcial=0;
      workers=size-1;

      if (size<2)
      {
    printf("No se puede procesar con %i nodos\n",size);
    return 0;
      }

      if (rank==MASTER)
      {
    for (i=MASTER+1;i<size;i++)
    {
      x=XINICIAL+(i-1)*XFINAL/(size-1);
      MPI_Send(&x, 1, MPI_FLOAT, i, TAG_INICIO, MPI_COMM_WORLD); 
      xfinal=XINICIAL+i*XFINAL/(size-1);
      MPI_Send(&xfinal, 1, MPI_FLOAT, i, TAG_FIN, MPI_COMM_WORLD); 
    }
    for (i=MASTER+1;i<size;i++)
    {
      MPI_Recv(&resultParcial, 1, MPI_FLOAT, i, TAG_RESULT_PARCIAL, 
                                            MPI_COMM_WORLD, &status); 
      resultado +=resultParcial;
    }
    printf("El resultado de es %f \n",resultado);    
}
else
{
     MPI_Recv(&x, 1, MPI_FLOAT, MASTER, TAG_INICIO, MPI_COMM_WORLD, 
                                                          &status); 
     MPI_Recv(&xfinal, 1, MPI_FLOAT, MASTER, TAG_FIN, MPI_COMM_WORLD, 
                                                          &status); 
     while (x<xfinal)
     {
        resultParcial += incremento * (calcularFuncion(x)+
                           calcularFuncion(x+incremento))/(float)2;
        x+=incremento;
     }
     MPI_Send(&resultParcial,1,MPI_FLOAT, MASTER, TAG_RESULT_PARCIAL
                                                 , MPI_COMM_WORLD);
 }             
}

float calcularFuncion(float x)
{
  return x*x;
}

En términos de consumo de red, parece ser mejor la segunda opción.

A continuación vamos a ver un planteo de este mismo problema implementado con MapReduce de Hadoop en su versión 1. Por motivos de simplicidad no se utilizó el File System distribuido, que es una de las ventajas que tiene este enfoque:

La interfaz IFuncion es la responsable de darles una “cara” consumible a las funciones:

package ar.edu.rdipasquale;

public interface IFuncion {
    public abstract Double calcularY(Double x);
}

La Función x2 se implementa de la siguiente manera:

package ar.edu.rdipasquale;

public class FuncionXCuadrado implements IFuncion{

    public FuncionXCuadrado() {
    }

    @Override
    public Double calcularY(Double x)
    {
        return x*x;
    }
}

La clase que se encarga de configurar el job y ejecutarlo se llama Riemann. Por simplicidad, no utilizará el file system distribuido y generará un archivo con una clave y un valor. La clave será la función (el id de la función) a integrar y el valor será el comienzo del rectángulo en el eje x:

package ar.edu.rdipasquale;

import java.io.PrintWriter;
import org.apache.hadoop.conf.Configuration;
import org.apache.hadoop.conf.Configured;
import org.apache.hadoop.fs.Path;
import org.apache.hadoop.io.DoubleWritable;
import org.apache.hadoop.io.Text;
import org.apache.hadoop.mapred.FileInputFormat;
import org.apache.hadoop.mapred.FileOutputFormat;
import org.apache.hadoop.mapred.JobClient;
import org.apache.hadoop.mapred.JobConf;
import org.apache.hadoop.mapred.TextInputFormat;
import org.apache.hadoop.mapred.TextOutputFormat;
import org.apache.hadoop.util.Tool;
import org.apache.hadoop.util.ToolRunner;

public class Riemann extends Configured implements Tool {

    private static final Double XINICIAL=(double)0;
    private static final Double XFINAL=(double)10;
    private static final int PASOS=100;
    private static final int FUNCIONES_A_INTEGRAR=1;

    public int run(String[] args) throws Exception {

        JobConf conf = new JobConf(Riemann.class);
          conf.setJobName("Riemann");      
          conf.setOutputKeyClass(Text.class);
          conf.setOutputValueClass(DoubleWritable.class);

          Double x=XINICIAL;

          // Habría que aprovechar el HDFS en vez del local... 
          PrintWriter out = new PrintWriter("funcionesaintegrar.txt");
          for (int i=0;i<FUNCIONES_A_INTEGRAR;i++)
              for (int j=0;j<PASOS;j++)
              {
                  out.println(String.valueOf(i+1)+";"+
                                                   String.valueOf(x));
                  x=(j+1)*XFINAL/PASOS;
              }
          out.flush();
          out.close();

          conf.setMapperClass(MapRiemann.class);
          conf.setCombinerClass(ReduceRiemann.class);
          conf.setReducerClass(ReduceRiemann.class);

          conf.setInputFormat(TextInputFormat.class);
          conf.setOutputFormat(TextOutputFormat.class);

          FileInputFormat.setInputPaths(conf, 
                                  new Path("funcionesaintegrar.txt"));
          FileOutputFormat.setOutputPath(conf, new Path("salida"));

          conf.set("riemann.funcion", 
                               "ar.edu.rdipasquale.FuncionXCuadrado");
          conf.set("riemann.incremento", 
                             String.valueOf((XFINAL-XINICIAL)/PASOS));
          JobClient.runJob(conf);        
          return 0;
    }

    public static void main(String[] args) throws Exception {
        int res = ToolRunner.run(new Configuration(), new Riemann(), 
                                                               args);
        System.exit(res);
    }
}

El Mapper será la Clase MapRiemann que tomará cada línea del archivo en cuestión y emitirá la superficie de cada rectángulo:

package ar.edu.rdipasquale;

import java.io.IOException;
import java.util.StringTokenizer;
import org.apache.hadoop.io.DoubleWritable;
import org.apache.hadoop.io.LongWritable;
import org.apache.hadoop.io.Text;
import org.apache.hadoop.mapred.JobConf;
import org.apache.hadoop.mapred.MapReduceBase;
import org.apache.hadoop.mapred.Mapper;
import org.apache.hadoop.mapred.OutputCollector;
import org.apache.hadoop.mapred.Reporter;

public class MapRiemann extends MapReduceBase implements 
                   Mapper<LongWritable, Text, Text, DoubleWritable> 
{
    private IFuncion funcion;
    private Double incremento;

    @Override
    public void configure(JobConf job) {
        try {
            funcion = (IFuncion)Class.forName(
                          job.get("riemann.funcion")).newInstance();
            incremento=Double.parseDouble(
                                     job.get("riemann.incremento"));
        } catch (Exception e) {
            e.printStackTrace();
        }
    }

    @Override
    public void map(LongWritable arg0, Text value,
            OutputCollector<Text, DoubleWritable> output, 
            Reporter arg3)
            throws IOException {

        String line = value.toString();
        StringTokenizer tokenizer = new StringTokenizer(line,";");
        String funcionAIntegrar=tokenizer.nextToken();
        Double valor=Double.parseDouble(tokenizer.nextToken());
        output.collect(new Text(funcionAIntegrar), 
               new DoubleWritable(incremento*
                           (funcion.calcularY(valor)+
                           funcion.calcularY(valor+incremento))/2));
    }
}

El Reducer está implementado en la clase ReduceRiemann y se encargará de sumar todas las superficies calculadas por los mappers (para cada función, obviamente, no mezclará peras con manzanas):

package ar.edu.rdipasquale;

import java.io.IOException;
import java.util.Iterator;
import org.apache.hadoop.io.DoubleWritable;
import org.apache.hadoop.io.Text;
import org.apache.hadoop.mapred.MapReduceBase;
import org.apache.hadoop.mapred.OutputCollector;
import org.apache.hadoop.mapred.Reducer;
import org.apache.hadoop.mapred.Reporter;

public class ReduceRiemann extends MapReduceBase implements 
                  Reducer<Text, DoubleWritable, Text, DoubleWritable> 
{
    @Override
    public void reduce(Text key, Iterator<DoubleWritable> values,
            OutputCollector<Text, DoubleWritable> output, 
            Reporter arg3)
            throws IOException {

      double sum = 0;
      while (values.hasNext()) 
          sum += values.next().get();
      output.collect(key, new DoubleWritable(sum));                
    }   
}

Como combiner se utiliza la misma clase utilizada como reducer para poder optimizar la performance realizando algunas agregaciones locales, lo que puede observarse en la configuración del job en la clase Riemann.

Sería tonto hacer una comparativa de performance entre los dos enfoques, ya que son herramientas que se utilizan para fines distintos.

El nivel de performance en el cálculo en la plataforma C+MPI es probablemente inalcanzable.

La plataforma MapReduce (en Hadoop) se va a destacar en trabajos orientados a grandes volúmenes de datos, donde minimizará el costo de red frente a la plataforma MPI.

Por otro lado la plataforma MPI se aplica a lenguajes de programación general, y si bien MapReduce puede utilizarse en C++, Java, Python y otros, los algoritmos deben expresarse en términos de concatenación de Mappers y Reducers, lo cual agrega una complejidad, entre las que está la no existencia de un estado, es decir, que son implementaciones stateless (esto en algoritmos de grafos puede llevar a dolores de cabeza, o a soluciones por debajo del óptimo en términos de performance).

A favor del enfoque MapReduce vamos a decir que es tolerante a fallos, en una buena instalación, en un cluster lo suficientemente grande y bien administrado, se nos puede romper un nodo en la mitad del proceso y el mismo finalizará sin ningún problema. Eso es muy difícil de alcanzar con C-MPI. Asimismo el nivel de seguimiento de las tareas está administrado por la herramienta provista por Hadoop, de modo que el monitoreo ya viene de fábrica.

En un clúster heterogéneo, las soluciones que brindamos en C serían tan performantes como el nodo menos performante del cluster, ya que asignamos la carga de manera lineal. En cambio en MapReduce la carga es administrada por el framework, de modo que asignará a cada uno según sus necesidades y les pedirá según sus posibilidades (Ahí está el verdadero Socialismo del siglo XXI en los clusters de computadoras).