C++ -> Fundamentos de Programación: Range-v3: Una introducción a la biblioteca (Parte VII)

Artículos de la serie:

Este nuevo post de la serie dedicada a la biblioteca Range-v3 centra su atención en la vista cartesian_product, destinada a iterar productos cartesianos de rangos. Analizaremos el modo de lograr un grado de expresividad en código similar al facilitado por Python, comparando las diferencias de rendimiento con dicho lenguaje.

Ejemplo 11: cartesian_product

Última actualización del artículo: 11 de diciembre de 2021.

En teoría de conjuntos, el producto cartesiano de dos conjuntos A y B, denotado A ⨯ B, se define como el conjunto de parejas ordenadas (a, b) con a ∈ A y b ∈ B. Con carácter general, el producto cartesiano de un número finito de conjuntos A₁, ..., A_n viene dado por el conjunto de n-tuplas para las que cada elemento i-ésimo pertenece a A_i (1 ≤ i ≤ n):

A₁ ⨯ ... ⨯ A_n := {(a₁, ..., a_n) | a₁ ∈ A₁, ..., a_n ∈ A_n}.

La biblioteca Range-v3, objeto de estudio de esta serie de artículos, modela esta importante construcción matemática a través de la vista ranges::views::cartesian_product disponible en la cabecera <range/v3/view/cartesian_product.hpp> [1]. Esta vista enumera el producto cartesiano de n rangos, pudiendo generar bajo demanda todas las n-tuplas (a₁, ..., a_n), con a₁ un elemento del primer rango, a₂ un elemento del segundo rango, etcétera. Existe una propuesta de estandarización, firmada por Sy Brand (Microsoft), encaminada a la futura inclusión de esta vista en el estándar C++23 [2].

A modo de ejemplo, consideremos un breve programa destinado a imprimir en la terminal todos los k-meros de longitud k = 3, es decir, todas las posibles cadenas de bases nitrogenadas pertenecientes al conjunto N = {A, C, G, T} que posean dicha longitud [3]. Esta tarea podría ser completada sin esfuerzo mediante una sucesión de bucles for-each anidados que generen todas las 3-tuplas del producto cartesiano N⨯N⨯N:

   auto const nucleobases = std::array{'A', 'C', 'G', 'T'};

   for (auto&& n1 : nucleobases) {
      for (auto&& n2 : nucleobases) {
         for (auto&& n3 : nucleobases) {
            fmt::print("{} {} {}\n", n1, n2, n3);
         }
      }
   }

La vista cartesian_product resulta especialmente conveniente en un caso como éste, al permitir simplificar el bloque de código anterior en la forma:

   namespace rn = ranges;
   namespace rv = rn::views;
  
   auto const nucleobases = std::array{'A', 'C', 'G', 'T'};      
  
   for (auto&& kmer : rv::cartesian_product(nucleobases, nucleobases, nucleobases)) {
      fmt::print("{}\n", fmt::join(kmer, " "));
   }

Aquí, las 3-tuplas de caracteres char obtenidas al iterar la vista son descompuestas por la función fmt::join de la biblioteca {fmt}lib con el fin de imprimir las bases A, C, G o T separadas por espacios en blanco. De forma alternativa, podríamos emplear una declaración structured binding de la forma:

   for (auto&& [n1, n2, n3] : rv::cartesian_product(nucleobases, nucleobases, nucleobases)) {
      fmt::print("{} {} {}\n", n1, n2, n3);
   }

En cualquiera de los casos anteriores, el listado de los cuatro primeros 3-meros así obtenidos sería:

A A A
A A C
A A G

A A T

A C A

. . .

Al generar la vista cartesian_product, resulta sin duda inconveniente el que debamos repetir el nombre del conjunto nucleobases tantas veces como corresponda a la longitud del k-mero deseado. Es posible simplificar dicha sintaxis sin más que introducir una función repeat_arg que, dado un entero sin signo N, un objeto invocable f y un argumento arg, realice la llamada f(arg,arg,...) con arg remitido N veces:

   template<std::size_t N, typename F, typename Arg>
   constexpr auto repeat_arg(F&& f, Arg&& arg) {
      return /* IILE */ [&]<std::size_t... Is>(std::index_sequence<Is...>){
                return std::invoke(std::forward<F>(f), (static_cast<void>(Is), arg)...);
             }(std::make_index_sequence<N>{});
   }

El uso de static_cast<void> evita warnings ante la existencia de valores no empleados. Observemos, asimismo, que N debe quedar establecida como expresión constante en tiempo de compilación. Nuestro código de impresión de 3-meros puede entonces ser reescrito como:

   auto const nucleobases = std::array{'A', 'C', 'G', 'T'};

   for (auto&& kmer : repeat_arg<3uz>(rv::cartesian_product, nucleobases)) {
      fmt::print("{}\n", fmt::join(kmer, " "));
   }

Soluciones similares pueden proponerse para otros lenguajes de tipado estático como Rust [4]. Comparemos ahora este programa desde la perspectiva de un lenguaje de tipado dinámico, especificamente Python 3, apoyándonos en el uso del iterador product de la biblioteca estándar itertools:

   # .py code
   import itertools

   nucleobases = ['A', 'C', 'G', 'T']

   for kmer in itertools.product(nucleobases, repeat = 3):
      print(*kmer)

A diferencia del código en C++, el número de repeticiones del rango de nucleobases se resuelve en Python de forma dinámica (argument unpacking). Se trata de una distinción crucial entre ambos lenguajes. El hecho de que el número de argumentos de una función, así como sus tipos, deba ser determinado estáticamente en C++ [5] es fuente, en el ejemplo considerado, de claras limitaciones: el parámetro N de la plantilla repeat_arg no puede ser establecido en tiempo de ejecución. Sin embargo, ello es también fuente de beneficios. Así, la siguiente línea de código

   fmt::print("{}\n", rn::distance(repeat_arg<14uz>(rv::cartesian_product, nucleobases)));

será capaz de imprimir el número de 14-meros (i.e., 4¹⁴ = 268'435'456) invirtiendo un tiempo nulo de ejecución en calcularlo: la longitud total de la vista podrá ser deducida en tiempo de compilación, sin necesidad de iterarla. Así lo atestigua la siguiente sección de instrucciones generadas por el compilador GCC 11.2 activando el nivel de optimización -O3:

   mov   QWORD PTR [rsp], 268435456
   call  fmt::v7::vprint(fmt::v7::basic_string_view<char>, fmt::v7::format_args)

Una de las principales características de la biblioteca Range-v3 consiste en la posibilidad de componer operaciones tales como vistas y acciones. Así, por ejemplo, supongamos que deseamos determinar el número de 14-meros que contienen al menos una vez la secuencia {AGCGGT} (en ese orden). Este número puede calcularse fácilmente de forma exacta gracias a que la secuencia escogida no puede auto-superponerse [6]:

9⨯4⁸ - 6⨯4² = 589'728,

donde el término sustraído evita que contabilicemos dos veces los 14-meros que presenten la secuencia repetida. Interesados por motivos ilustrativos en calcular este número forzando la iteración de una vista, codificaríamos (supuesta la definición de las nucleobases en los códigos anteriores):

   auto const seq = std::array{'A', 'G', 'C', 'G', 'G', 'T'};

   auto const k = 14uz;
   auto k_mers_with_seq = repeat_arg<k>(rv::cartesian_product, nucleobases) 
                        | rv::filter([&](auto&& tup){
                             auto const kmer = tuple_to_array(tup);
                             return not rn::search(kmer, seq).empty();
                          });
         
   fmt::print("number of {}-mers containing ({}) = {}\n",
              k, fmt::join(seq, ", "), rn::distance(k_mers_with_seq));

donde, con el fin de facilitar la búsqueda de la secuencia {AGCGGT} en un 14-mero dado, hemos definido la función tuple_to_array para transformar una n-tupla de valores de un tipo común en un array de longitud n:

   template<typename T = void, typename Tuple>
   constexpr auto tuple_to_array(Tuple&& tuple)
   {
      return std::apply(
                []<typename... Args>(Args&&... vs){
                   using U = std::conditional_t<std::is_same_v<T, void>,
                                                std::common_type_t<Args...>, T>;
                   return std::array<U, sizeof...(vs)>{std::forward<Args>(vs)...};
                },
                std::forward<Tuple>(tuple)
             );
   }

El programador usuario puede especificar explícitamente el tipo a contener en el array (parámetro T) o permitir que éste sea deducido como el tipo común existente en la tupla (std::common_type_t<Args...>), si existe.

Nuevamente, comprobamos que el código anterior en C++ moderno adopta una estructura similar a la que emplearíamos de forma idiomática en Python 3:

   # .py code
   seq = 'AGCGGT'

   k = 14
   num_k_mers_with_seq = sum(1 for kmer in itertools.product(nucleobases, repeat = k)
                                        if seq in ''.join(kmer))

   print(f'number of {k}-mers containing {seq} = {num_k_mers_with_seq}')

Sin embargo, hagamos notar que, dada por ejemplo una máquina Intel(R) Core(TM) i3-7020U CPU a 2.30GHz y 8GB de memoria RAM, la versión del programa escrita en C++ sería capaz de iterar la vista filtrada de 14-meros en aproximadamente 3 segundos, mientras que su equivalente en Python podría llegar a demorarse hasta dos minutos. Ello invitaría a investigar optimizaciones adicionales en este último caso, posiblemente mediante el uso de bibliotecas alternativas como NumPy.

Nota: Mi agradecimiento a José Manuel López (Universidad Europea) por sus valiosos comentarios acerca de soluciones idiomáticas y optimizaciones de código en Python.

Referencias bibliográficas:

Range-v3 - User manual - https://ericniebler.github.io/range-v3/
Sy Brand - P2374R1 - views::cartesian_product - http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p2374r1.html
Wikipedia - k-mer - https://en.wikipedia.org/wiki/K-mer
stack overflow - https://stackoverflow.com/questions/44139493/in-rust-what-is-the-proper-way-to-replicate-pythons-repeat-parameter-in-iter
stack overflow - https://stackoverflow.com/questions/56026842/python-like-dynamic-argument-unpacking-in-c
StackExchange - https://math.stackexchange.com/questions/3454913/probability-of-having-a-specific-k-mer-in-a-randomised-sequence-of-length-13

C++ -> Fundamentos de Programación

Range-v3: Una introducción a la biblioteca (Parte VII)

Ejemplo 11: cartesian_product

No hay comentarios:

Publicar un comentario