1

Programación y Computación paralela (4)

Integración, Monte Carlo,Metropolis

Glen D. Rodríguez R.

Basado en material de J. Demmel

2

Links

• Capítulo 4 del libro de Pacheco

• Ver http://www.old-npac.org/users/gcf/slitex/CPS615NI95/index.html

• Ver http://www.old-npac.org/projects/cps615spring00/presentations-html/cps615nic95/

• Ver http://www.old-npac.org/projects/cpsedu/summer98summary/examples/mpi-c/mpi-c.html

3

Integración

• Se puede integrar en n dimensiones, pero por ahora discutiremos sólo en una dimensión.• Para n mayor de 2, el método Monte Carlo y similares tiene

muchas ventajas

• Se desea calcular la integral de la función f(x) de x=ahasta x=b

• Implícitamente se asume que f(x) es relativamente suave, se determina un conjunto de puntos de “interpolación” o de malla xi en la región a ≤ x ≤ b y se calcula la integral en término de los valores de f(x) y los puntos de la malla.

∫I = f(x) dx

a

b

4

Regla trapezoidal

• Se asume que f)x= es “o constante o lineal” en un intervalo entre dos puntos

• I = 0.5*(b-a)*(f(a)+f(b)), o si no (menos preciso)• I = (b-a)*f((a+b)/2)

x0 xX=a X=b

f(a)

f(b)

5

Regla de Simpson

• Se asume que f(x) es una función cuadrática

• I = (b-a)*(f(a)+ 4f(xc)+f(b))/6

xX=a X=b

f(a) f(b)

xc = (a+ b)/2

f(Xc)

6

Reglas iterativas

• Hay ciertas reglas (Newton-Cotes y Gaussiana) que hacer aproximaciones de mayor orden de f(x)• Son difíciles de programar y no muy “robustas” (funcionan mal

si f(x) varía muy rápido de forma no polinomial)

• En la práctica se usa Simpson o Trapezoidal iterativo –Uar un número impar de puntos para aplicar simpson de forma iterativa

ej.: Aplicar Simpson sobre cada uno de los 18 sub intervalos indiados por

Las flechas abajo. Trapezoidal iterativo involucraría 36 integrales

x0 x4 x8 x12 x16 x20 x24 x28 x32 x36

7

Método Monte Carlo

• El Monte Carlo más sencillo escoje puntos al azar• La Integral es la suma de los valores de la función f(xRandom-i) en

los valores escogidos al azar xRandom-i multiplicados por el intervalo de integración (b-a) y dividido por el número depuntos generados.

• Esta da un método robusto (y sin preferencias) que es fácil de usar en integrales de n-dimensiones o integrales de funciones problemáticas

• Se puede adaptar en límites complicados• Se integra sobre regiones más grandes de lo necesario• Se define f(x) =0 fuera del intervalo [a b]

ba

8

Hallar pi por integración con Monte Carlo

Otra forma de Monte Carlo: g(x,y) es un punto aleatorio en R2.Ambas son equivalentes

9

1

1

Integral ≈ A(n / N)

N= total de puntos(cruces y equis)n= puntos dentrode la integral (cruces)A= área de la simulación(área amarilla) = 1

x

yx2 + y2 =1

x

++

+

+

+

+

+

+

+

+

+

++

x

x

xx

x

x

+

Puntos al azar (x,y), donde x e y tienen distribución uniforme

10

Errores

• Para una integral con N puntos

• Monte Carlo tiene un error del orden de 1/N0.5

• Trapezoidal iterativo : 1/N2

• Simpson iterativo: 1/N4

• Pero en d dimensiones, todos menos el Monte Carlodeben usar una malla de N1/d puntos por lado; no funciona bien para N>3• Monte Carlo aún tiene error 1/N0.5

• Error de Simpson es 1/N4/d

• Cuando usar Montecarlo? Cuando error de Montecarloes menor que error de Simpson para los mismos N puntos? 1/N0.5 < 1/N4/d � 0.5>4/d , o sea d>8

11

Trapezoidal Iterativo y Simpson iterativo

• Trapezoidal: Si tenemos N puntos N xi con

I= (b-a) (f(x0)+2f(x1)+2f(x2)+ ……+2f(x(N-2))+f(x(N-1))/(2(N-1))

• O con Simpson:

I= (b-a) (f(x0)+4f(x1)+2f(x2)+ ……+4f(x(N-2))+f(x(N-1))/(3(N-1))

• Se suman los f(xi) con diferentes pesos• Note que ambos aproximan su cálculo en una región

tamaño (b-a)/(N-1) (el doble de eso para Simpson) que tiende a cero si N es grande, así que la aprox. Es mejor si N se agranda.

12

Las funciones de peso

• Se suman varios wi f(xi) y las diferentes reglas sólo con

diferentes wi

• Trapezoidal: para índices = 0, 1, 2, 3, ..., N-1Se usa el patrón

wi ∝ 1, 2, 2, … 2, 2, 1

• Para los mismos índices, Simpson con N impar tieneelpatrón

wi ∝ 1, 2, 4, 2, 4, … 2, 4, 2, 4, 1

13

Computación paralela

• Diferentes procesadores calculan independientemente diferentes puntos en la integral y luego se suman los resultados de cada procesador

• Lo más natural es el escenario de la figura de arriba, con cada procesador calculando una integral sobre una subregión

• x0 al x8 en proc. 0 … x24 to x32 en proc. 3 etc.

• La integral total es la suma de las sub-integrales calculadas en cada procesador

• Pero podría asignarse PUNTOS y no REGIONES a los procesadores• x0 x4 x8 .. x32 al proc. 0

• x3 x7 x11 .. x31 al proc. 3

x0 x4 x8 x12 x16 x20 x24 x28 x32

10 2 3

14

Integración con reglas y Paralelismo

• Es un ejemplo clásico de “master worker” donde los procesadores individuales calculan arte de la integral.

• Se necesita usar “MPI rank” para determinar quien se encarga de cada rango de la integración

• Luego se suma la contribución de cada procesador• Suma global acumulativa (uno tras otro), por árbol o por

reducción al “master”.

15

Ejemplo 2: Monte Carlo con dist.uniforme#include <stdio.h>#include <stdlib.h>#include <time.h>#include <math.h>#include "mpi.h"

int main(int argc, char *argv[]){

int nprocs, myproc;int i, n, acum, sum_acum;double pi_25= 3.141592653589793238462643;double pi, w, x, y, error;double t0, t1;

/* MPI initialization */MPI_Init( &argc, &argv );MPI_Comm_size( MPI_COMM_WORLD, &nprocs );MPI_Comm_rank( MPI_COMM_WORLD, &myproc );

/* Determine the number of divisions to use */if( argc == 2 ) { /* Take n from the command line argument */

sscanf( argv[1], "%d", &n );} else {

if (myproc == 0){printf("Enter the number of random points: (0 quits) ");fflush(stdout);scanf("%d",&n);}

}

16

Sigue….if( myproc == 0 ) printf("Calculating pi using %d points\n", n);

MPI_Bcast( &n, 1, MPI_INT, 0, MPI_COMM_WORLD); /* ??? Is this needed ??? *//* Start the timer after a barrier command */

MPI_Barrier( MPI_COMM_WORLD );t0 = MPI_Wtime();

pi = 0.0;sum_acum=0;acum=0;srand( (unsigned)time( NULL )+myproc );w = 1.0 / n; for( i=myproc; i<n; i+=nprocs ) {

x = (double) rand()/RAND_MAX;y = (double) rand()/RAND_MAX;if (((x*x)+(y*y)) <= 1.0) {acum=acum+1;}

}

MPI_Allreduce( &acum, &sum_acum, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );pi= ((double) sum_acum)/(double) n;pi = pi*4;error = fabs( pi - pi_25 );

t1 = MPI_Wtime();if( myproc == 0 ) {

printf("The calculated pi = %f (error = %f)\n", pi, error);printf("The calculation took %f seconds on %d nodes\n", t1-t0, nprocs);

}MPI_Finalize();

}

17

Cómo mejorar la precisión de Monte Carlo

• Error absoluto de regla rectangular (peor que

trapezoidal), usando 10000 intervalos = 0.000199

• Error absoluto de regla trapezoidal usando 10000 intervalos = 0.000001

• Error absoluto promedio de Monte Carlo usando 10000 puntos y 10 simulaciones = 0.013521

• Como mejorar: no usar distribución uniforme. Se

necesita una "guía" al muestreo aleatorio que enfatiza el muestreo donde la función es grande o varia

rápidamente.

• Ver ejemplo

• En ese caso, el error absoluto promedio de Monte Carlo

usando 10000 puntos y 10 simulaciones = 0.006479

18

1

1

Integral ≈ ????

N= total de puntos(cruces y equis)n= puntos dentrode la integral (cruces)A= área de la simulación(área amarilla)= 0.857

x

yx2 + y2 =1

++

+

+

+

+

+

+

+

+

+

++

x

x

x

+

Puntos al azar (x,y), donde x, y no tienen dist. uniforme

19

Muestreo por importancia

• Quiero integrar en 1-D

• Se mete la función peso w

• Que debe satisfacer

• Cambiando la variable, con dy = w(x) dx, y(0) = 0, and y(1) = 1.

20

Muestreo por importancia

• Eso nos da

• La función w se escoge de tal forma que f/w sea cercana a una constante. El muestreo uniforme en y se modifica (vuelve no uniforme) al ser afectado por w. O dicho de otra forma, la integral anterior es equivalente a la inicial, pero transformada por el cambio de variable en una función con menos variaciones.

21

Programa Ejemplo 3

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

int main(int argc, char *argv[]){

int nprocs, myproc;int i, n;double pi_25= 3.141592653589793238462643;double pi, w, x, y, error, acum, sum_acum;double t0, t1;

/* MPI initialization */MPI_Init( &argc, &argv );MPI_Comm_size( MPI_COMM_WORLD, &nprocs );MPI_Comm_rank( MPI_COMM_WORLD, &myproc );

/* Determine the number of divisions to use */if( argc == 2 ) { /* Take n from the command line argument */

sscanf( argv[1], "%d", &n );} else {

if (myproc == 0){printf("Enter the number of random points: (0 quits) ");fflush(stdout);scanf("%d",&n);}

}

22

Sigue…

if( myproc == 0 ) printf("Calculating pi using %d points\n", n);

/* Broadcast the number of divisions to all nodes */

MPI_Bcast( &n, 1, MPI_INT, 0, MPI_COMM_WORLD); /* ??? Is this

needed ??? */

/* Start the timer after a barrier command */

MPI_Barrier( MPI_COMM_WORLD );

t0 = MPI_Wtime();

pi = 0.0;

sum_acum=0.0;

acum=0.0;

srand( (unsigned)time( NULL )+myproc );

w = 1.0 / n;

/* Each processor starts at a different value and computes its

* own contributions to pi. */

23

Sigue…

for( i=myproc; i<n; i+=nprocs ) {x = (double) rand()/RAND_MAX;y = (double) rand()/RAND_MAX;if ((x>0.5) && (x<=0.8) && (y<=1.5-x) && (((x*x)+(y*y)) <= 1.0)) acum=acum+0.857;if ((x>0.5) && (x<=0.8) && (y>1.5-x)) i=i-nprocs;if ((x>0.8) && (y<=1.77-1.4*x) && (((x*x)+(y*y)) <= 1.0)) acum=acum+0.857;if ((x>0.8) && (y>1.77-1.4*x)) i=i-nprocs;if ( (x<=0.5) && (((x*x)+(y*y)) <= 1.0) ) acum=acum+0.857;

}

MPI_Allreduce( &acum, &sum_acum, 1, MPI_DOUBLE, MPI_SUM,

MPI_COMM_WORLD );

pi= ((double) sum_acum)/(double) n;pi = pi*4;error = fabs( pi - pi_25 );

t1 = MPI_Wtime();if( myproc == 0 ) {

printf("The calculated pi = %f (error = %f)\n", pi, error);printf("The calculation took %f seconds on %d nodes\n", t1-t0,

nprocs);}MPI_Finalize();

}

24

Paradigmas de Monte Carlo paralelo

• Master-worker: se genera los números

aleatorios en el nodo 0 y se pasan a los demás.

Todos los nodos computan.

25

Paradigmas de Monte Carlo paralelo

• Cliente-servidor: es un servidor de números

aleatorios, no procesa, sólo crea y envía los

randoms

26

Paradigmas de Monte Carlo paralelo

• Full workers: cada nodo genera sus propios

números aleatorios. El nodo 0 sólo cuida que

se siembren diferentes semillas o que se use la

misma semilla pero diferentes rangos.

27

Resolviendo Ec. de Poisson con Montecarlo

• Sea la ecuación de Poisson:

• Con fuente sinusoidal : -2π2 sen(π x) sen(π y)

• Y cero en todas las fronteras

• En diferencias finitas:

28

Solución con “random walks”

• De un punto (i,j) comenzar el “random walk”.La

probabilidad de moverse a cualquiera de los 4

puntos vecinos es ¼.

• Generar un número aleatorio para escoger el

vecino.

• Añadir g(x,y) en la nueva posición

• Repetir hasta llegar a una frontera

• Esto fue solo UNA “random walk”

• Después de N “random walks” el estimado de

u(i,j) es:

29

Paralelización de esta solución

• Hacer el siguiente proceso para todos los puntos

• Comenzar en una frontera• De afuera hacia adentro

• Fila por fila

• Actualizar (intercambiar) data en el límite entre procesadores

• Actualizar el walk dentro de los puntos de un procesador

• Si el walk sale de un procesador, informarle al otro procesador

30

Metrópolis

• Se usa para integrales con funciones de peso.

• Genera un random walk pero guiado por una función

de peso w(x). Ejemplo en 2-D comenzando de un punto (xi,yi):

1. Escoger δ, el step size

2. Generar dos aleatorios R1, R2 uniformes en el rango [-1,1]

3. El nuevo punto será1. xT

i+1 = xi + δR1

2. yTi+1 = yi + δR2

4. Evaluar el ratio de w en el punto actual vs. el punto anterior1. r= w(xT

i+1 ,yTi+1) / w(xi,yi)

31

Metrópolis

5. Si r>1 aceptar el nuevo punto y regresar a (2)1. xi+1 = xT

i+1 ; yi+1 = yTi+1

6. Si r<1 aceptar el nuevo punto con probabilidad r. O sea, generar un aleatorio uniforme η en [0, 1] y

aceptar el punto solo si r>η (hacer lo mismo que en el paso 5)

7. Caso contrario, rechazar el nuevo punto y volver a (2)

El step size no debe ser ni muy chico ni muy grande.

Los puntos obtenidos se usan luego para integrar, como

en MonteCarlo (Algo*n/N).

La paralelización es similar a la discutida para las otras

integrales con Montecarlo