# Tercer Parcial
Escribir un programa usando CUDA que permita aproximar el valor de una integral definida por el
método de Monte Carlo. La integral definida será:

$$
\theta = \int_a^b x^{2}dx
$$

Para cambiar los límites a y b se puede usar un cambio de variables.

## Desarrollo del problema
Para aplicar el método de integración de Montecarlo en este problema, se generan números aleatorios uniformemente distribuidos entre 0 y 1 ($U$), los cuales son evaluados en $f(x) = x^2$. Debido a que esto solo funciona para integrales entre 0 y 1, para cualquier rango arbitrario $[a , b]$, hay que aplicar la siguiente transformación:
$$
\theta = \int_a^b x^2dx = \int_0^1h(y)dy
$$
Donde
$$
h(y) = (b - a) \cdot f(a + (b - a)\cdot x)
$$
La anterior transformación nos permite convertir cualquier integral con limites $[a, b]$ en una con limites $[0, 1]$. Por ende, el método consiste en realizar estimaciones de la integral de la siguiente manera:
$$
\text{est} = (b - a) \cdot f(a + (b - a)\cdot U)
$$
Donde $U$ es un número aleatorio uniformemente distribuido. El resultado estará dado por la media de estos valores:
$$
\text{aprox} = \text{mean}(est)
$$
Hay que mencionar que para obtener buenas aproximaciones tener cantidades considerables de iteraciones (a partir de 1000 o más).
## Implementación en serial

#include <stdio.h>
#include <stdlib.h>
/**
 * Para utilizar la función pow
 */
#include <math.h>
/**
 * Para parsear los argumentos desde la linea de comandos
 */
#include <string.h>
/**
 * Para medir el tiempo de ejecución
 */
#include <time.h>

/**
 * Consigue el valor de una flag proveniente de la linea de ocmandos
 *
 * @param argc cantidad de argumentos
 * @param argv array de argumentos
 * @param flag flag a parsear
 * @return int valor de la flag
 */
double get_flag_value(int argc, char **argv, char flag[])
{
    double n = -1, arg = 0;
    for (int i = 0; i < argc; i++)
    {
        if (strcmp(flag, argv[i]) == 0 && !arg)
        {
            n = atof(argv[i + 1]);
            arg = 1;
        }
        else
        {
            continue;
        }
    }

    if (n == -1)
    {
        printf("Missing flag.\n");
    }

    return n;
}

/**
 * función de la cual se quiere obtener la integral
 */
double f(double x) {
    return pow(x, 2);
}

/**
 * Antiderivada de la función original
 */
double F(double x) {
    return pow(x, 3) / 3;
}

/**
 * Función para generar numeros aleatorios entre 0 y 1
 */
double rand_0_1() {
    return (double)rand() / (double)RAND_MAX;
}

/**
 * Función para calcular el método de integración Montecarlo
 */
double montecarlo_int(double a, double b, double n)
{
    double result = 0;
    for (int i = 0; i <= n; i++)
    {
        result += (b - a) * f(a + (b - a) * rand_0_1());
    }
    result /= n;
    return result;
}

int main(int argc, char **argv)
{
    /**
     * Garantiza una semilla distinta para ejecución del programa
     */
    srand(time(NULL));
    char a_flag[] = "-a";
    char b_flag[] = "-b";
    char n_flag[] = "-n";

    double a, b, n;
    double X;

    a = get_flag_value(argc, argv, a_flag);
    b = get_flag_value(argc, argv, b_flag);
    n = get_flag_value(argc, argv, n_flag);

    printf("a = %f, b = %f, n = %f\n", a, b, n);

    double approx, real, error;

    /**
     * Mide el tiempo de ejecución de la aproximación
     */
    clock_t begin = clock();
    /**
     * Calcula la integral con Montecarlo
     */
    approx = montecarlo_int(a, b, n);
    clock_t end = clock();

    double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;

    /**
     * calcula la integral de la manera estándar
     */
    real = F(b) - F(a);
    /**
     * Obtiene el error de la aproximación
     */
    error = fabs((approx - real) / real);

    printf("Approximated value = %f (Tiempo empleado = %fs, %f iteraciones)\n", approx, time_spent, n);
    printf("Real value = %f\n", real);
    printf("Error = %f\n", error);

    return 0;
}

El resultado de correr el anterior programa para: $a = -2$, $b = 1$ y $n = 25000$ se muestra a cotinuación:

## Implementación en paralelo (CUDA)
Los bloques utilizados para los kernels son de tamaño $\text{BLOCK_SIZE}$ (1024), y el tamaño del grid se cálculo con base en la cantidad de iteraciones requeridas a partir de la siguiente expresión:
$$
\text{grid size} = \frac{n + \text{BLOCK_SIZE} - 1}{\text{BLOCK_SIZE}}
$$
Con CUDA se realizaron las estimaciones de la integral, mientras que la obtención del promedio de estos valores se realizó de manera serial con un `for` loop.

#include <stdio.h>
#include <stdlib.h>
/**
 * Para utilizar la función pow
 */
#include <math.h>
/**
 * Para parsear los argumentos desde la linea de comandos
 */
#include <string.h>
#include <cuda.h>
#include <curand.h>
#include <curand_kernel.h>
#include <chrono>


/**
 * Función para obtener los parametros pasados desde consola
 */
double get_flag_value(int argc, char **argv, char flag[])
{
    double n = -1, arg = 0;
    for (int i = 0; i < argc; i++)
    {
        if (strcmp(flag, argv[i]) == 0 && !arg)
        {
            n = atof(argv[i + 1]);
            arg = 1;
        }
        else
        {
            continue;
        }
    }

    if (n == -1)
    {
        printf("Missing flag.\n");
    }

    return n;
}

/**
 * función de la cual se quiere obtener la integral
 */
__device__ double f(double x)
{
    return pow(x, 2);
}

/**
 * Antiderivada de la función original
 */
double F(double x)
{
    return pow(x, 3) / 3;
}

/**
 * Kernel para calcular el método de integración montecarlo
 */
__global__ void montecarlo_gpu(double *estimations, curandState *state, double a, double b, int N)
{
    unsigned i = threadIdx.x + (blockDim.x * blockIdx.x);

    curandState localState = state[i];

    if (i < N)
    {
        estimations[i] = (b - a) * f(a + (b - a) * curand_uniform(&localState));
        state[i] = localState;
    }
}

/**
 * Kernel para inicializar los estados de los números aleatorios a generar
 */
__global__ void Init(curandState *state)
{
    unsigned idy = threadIdx.x + (blockDim.x * blockIdx.x);
    curand_init(1234, idy, 0, &state[idy]);
}

int main(int argc, char **argv)
{
    char a_flag[] = "-a";
    char b_flag[] = "-b";
    char n_flag[] = "-n";

    double a, b;
    int n, BLOCK_SIZE = 1024;

    a = get_flag_value(argc, argv, a_flag);
    b = get_flag_value(argc, argv, b_flag);
    n = floor(get_flag_value(argc, argv, n_flag));

    double *cuda_random_estimations, *random_estimations;
    /**
     * Array con los estados para gen los números aleatorios
     */
    curandState *array_states;

    /**
     * inicializa los eventos
     */ 
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    /**
     * Reserva memoria para las estimaciones que se calculan en la GPU y que luego se copian en la CPU.
     */
    cudaMalloc((void **)&cuda_random_estimations, n * sizeof(double));
    cudaMalloc((void **)&array_states, n * sizeof(curandState));
    random_estimations = (double *)malloc(n * sizeof(double));

    printf("a = %f, b = %f, n = %d\n", a, b, n);

    double approx, real;

    int grid_size = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;

    /**
     * Empieza a medir el tiempo que le toma al programa
     */
    cudaEventRecord(start);

    /**
     * Lanza el kernel para incializar cuRAND
     */
    Init<<<grid_size, BLOCK_SIZE>>>(array_states);

    /**
     * Lanza el kernel para calcular la integración Montecarlo
     */
    montecarlo_gpu<<<grid_size, BLOCK_SIZE>>>(cuda_random_estimations, array_states, a, b, n);

    /**
     * Copia los resultados calculados en la GPU (cuda_random_estimations) en la CPU (random_estimations)
     */
    cudaMemcpy(random_estimations, cuda_random_estimations, n * sizeof(double), cudaMemcpyDeviceToHost);
    /**
     * Sincroniza los kernels para que todos hayan acabado antes de continuar
     */
    cudaDeviceSynchronize();


    /**
     * Calcula el promedio de las estimaciones, el cual corresponde a la aproximación de la integral
     */
    for (int i = 0; i < n; i++)
    {
        approx += random_estimations[i];
    }

    approx /= n;

    /**
     * calcula la integral de la manera estándar
     */
    real = F(b) - F(a);

    /**
     * Termina de medir el tiempo
     */
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);

    printf("Approximated value = %f (Tiempo empleado = %fms, %d iteraciones)\n", approx, milliseconds, n);
    printf("Real value = %f\n", real);

    /**
     * Libera la memoria
     */
    cudaFree(array_states);
    cudaFree(cuda_random_estimations);
    free(random_estimations);

    return 0;
}

El resultado de correr el anterior programa para: $a = -2$, $b = 1$ y $n = 25000$ se muestra a cotinuación:

## Conclusiones
Como se puede apreciar de los resultados de ambos programas, el realizado en CUDA obtiene una aproximación más precisa, pero el tiempo empleado en realizar los calculos es mayor. Lo anterior puede deberse a que con la cantidad de iteraciones realizadas, el poder de la parelelización con CUDA no sale a flote.

Sin embargo, la razón de peso por la que el programa en paralelo corre en un mayor tiempo, es debido a que este se ejecutó utilizando Google Colab, el cual agrega cierto overhead por parte de la red y la máquina utilizada.