Programación científica (V): Comprobando la ley de Zipf

Artículos de la serie:

1. Representando el atractor de Lorenz con C++20, Boost.Numeric.Odeint y Dlib
2. Análisis dimensional con mp-units, regresiones no-lineales con Dlib.Optimization y visualización de datos con Matplot++
3. Rectas de mejor ajuste por mínimos cuadrados (mp-units, Matplot++)
4. Ceros de una función
5. Comprobando la ley de Zipf
6. Instalación de Matplot++ y mp-units con Mingw-w64

Nube de palabras para los sonetos de Shakespeare

En este artículo nos plantearemos un interesante problema de análisis estadístico de textos. Consideremos una obra literaria de cierta extensión --obtenida, por ejemplo, a través del proyecto Gutenberg [1]. Deseamos:

1. Crear un mapa que asocie a cada frecuencia de ocurrencia en el texto un listado de aquellas palabras que se repitan dicho número de veces. Para ello, deberemos tokenizar el texto adecuadamente.
2. Imprimir en la terminal, en orden decreciente, las cinco mayores frecuencias registradas junto a sus listas de términos asociadas.
3. Registrar el conjunto de valores (ranking, frequency) para las palabras distintas del texto. La primera entrada indicará la posición en el ranking de una palabra, mientras que la segunda entrada contendrá su frecuencia de ocurrencia en el texto. Por ejemplo, si "the" fuese el término más empleado con una frecuencia de uso de 1345, deberemos registrar el punto (1, 1345) ; si "of" fuese el segundo término en el ranking con una frecuencia de 703, registraríamos (2, 703), etcétera. Si dos o más términos compartieran una misma frecuencia de uso, les asignaríamos posiciones de ranking consecutivas y la misma frecuencia.

Utilizaremos los datos obtenidos en el punto 3 para representar gráficamente la relación log10(frequency) vs log10(ranking) y verificar el cumplimiento de la ley de Zipf [2]. Es decir, comprobaremos que la secuencia de puntos obtenida sigue una tendencia de tipo potencia de la forma

frequency = a·(ranking)^b   [ecuación 1]

donde a > 0 y b < 0 son constantes reales, siendo |b| ligeramente superior a uno.

En primer lugar, incluiremos las cabeceras relevantes para nuestro proyecto, tanto del estándar C++20 como de las bibliotecas Boost [3], {fmt}lib [4], Matplot++ [5] y Range-v3 [6]:

   #include <cctype>    #include <cmath>    #include <cstdlib>    #include <map>    #include <string>    #include <vector>    #include <boost/math/statistics/linear_regression.hpp>    #include <boost/tokenizer.hpp>    #include <fmt/format.h>    #include <matplot/matplot.h>    #include <range/v3/core.hpp>    #include <range/v3/action/remove_if.hpp>    #include <range/v3/action/sort.hpp>    #include <range/v3/action/transform.hpp>    #include <range/v3/view/group_by.hpp>    #include <range/v3/view/indices.hpp>    #include <range/v3/view/intersperse.hpp>    #include <range/v3/view/map.hpp>    #include <range/v3/view/take.hpp>

Gracias a la función fileread de la biblioteca Matplot++ podremos cargar el texto completo a analizar en un objeto std::string. El archivo de origen será indicado por el usuario a través del segundo argumento de la línea de comandos. El texto será convertido a minúscula y tokenizado convenientemente con el fin de construir un mapa que posea como clave una frecuencia de ocurrencia determinada y, como valor asociado, un vector con las palabras repetidas dicho número de veces (punto 1 del enunciado del problema). La biblioteca Range-v3 nos permitirá realizar estas operaciones adoptando un estilo funcional de programación (véase esta serie de artículos anterior para más detalles):

   auto main(int argc, char* argv[]) -> int    try {       namespace mp = matplot;       namespace rn = ranges;       namespace ra = rn::actions;       namespace rv = rn::views;       if (argc != 2) {          fmt::print(stderr, "invalid number of command-line arguments\n");          return EXIT_FAILURE;       }       auto const freq_tokens_map = /* IILE: mapa frecuencia -> vector de palabras */ [&]{          auto const tokens = /* IILE: vector de tokens ordenados lexicográficamente */ [&]{             // texto completo en minúscula:             auto const text = mp::fileread(argv[1])                             | ra::transform([](unsigned char c){ return std::tolower(c); });             return boost::tokenizer{text}                  | rn::to<std::vector<std::string>>                  | ra::sort                  | ra::remove_if(rn::empty);          }();          auto res = std::map<int, std::vector<std::string>, rn::greater>{};          for (auto tkn_grp : tokens | rv::group_by(rn::equal_to{})) {             res[rn::distance(tkn_grp)].push_back(*tkn_grp.begin());          }          return res;       }();

En el código anterior, tokens es un vector de strings que recoge, una a una, las distintas palabras empleadas en el texto, convertidas a minúscula y ordenadas lexicográficamente. Los separadores empleados por defecto por boost::tokenizer son los espacios en blanco (std::isspace(ch)!=0) [7] y las puntuaciones (std::ispunct(ch)!=0) [8]. La reducida longitud de la mayoría de las palabras así obtenidas activará la optimización de strings cortos propia del estándar C++11 y posteriores. El bucle for toma el vector de tokens ordenado lexicográficamente y lo somete a una vista group_by que identifica el inicio y el final de cada subgrupo tkn_grp de palabras idénticas. Dichos subgrupos son entonces iterados por el bucle. Observemos que la longitud de cada subgrupo se corresponde con la frecuencia de ocurrencia de la palabra correspondiente en el texto. Dicha longitud es, pues, empleada como clave del mapa, remitiendo al fondo de su vector asociado una copia del primer representante del subgrupo, *tkn_grp.begin().

Habiendo registrado las claves del mapa en orden decreciente mediante la política de ordenación ranges::greater, resulta inmediato imprimir en la terminal las cinco mayores frecuencias de ocurrencia junto a sus correspondientes listados de palabras, separadas por comas (punto 2 del enunciado del problema):

      // imprimimos en la terminal las cinco primeras parejas clave-valor del mapa:       for (auto const& [freq, tokens] : freq_tokens_map | rv::take(5)) {          fmt::print("Frequency {:<5}: ", freq);          for (auto const& tk : tokens | rv::intersperse(", ")) {             fmt::print("{}", tk);          }          fmt::print("\n{:_^30}\n", "");       }

Aquí, la vista intersperse introducirá una coma y un espacio entre una palabra y otra de existir más de un token para una frecuencia determinada, pero no tras la última palabra del vector.

Procederemos a continuación a generar dos secuencias numéricas cuyas entradas estarán en correspondencia uno-a-uno:

• Una vista de nombre ranks, que dará cuenta de las posiciones de ranking de las palabras registradas en el mapa (las primeras posiciones corresponderán a los términos más utilizados).
• Un vector de nombre freqs, que registrará la frecuencia de ocurrencia asociada a cada una de dichas palabras.

Recordemos que, tal y como explicamos al inicio del post, dos palabras con idéntica frecuencia de uso recibirán posiciones consecutivas en el ranking:

      auto const freqs = /* IILE: vector de frecuencias */ [&]{          auto res = std::vector<double>{};          for (auto const& [freq, tokens] : freq_tokens_map) {             res.insert(res.end(), tokens.size(), freq);          }          return res;       }();             // vista de posiciones de ranking asociadas a las frecuencias anteriores:       auto ranks = rv::closed_indices(1ull, freqs.size());

En relación al vector de frecuencias freqs, observemos que éste podría ser inicializado de forma alternativa empleando una transformación de tipo flat map:

      auto const freqs = freq_tokens_map                        | rv::for_each([](auto const& p){                             auto const& [freq, tokens] = p;                             return rv::repeat_n(freq, tokens.size());                          })                        | rn::to<std::vector<double>>;

Las herramientas estadísticas de la biblioteca Boost.Math nos permitirán realizar un ajuste por mínimos cuadrados al conjunto de puntos (ranking, frequency) según la relación lineal ln(frequency) = ln(a) + b·ln(ranking), denotando ln el logaritmo natural. Esta fórmula se obtiene de la ecuación 1 proporcionada más arriba sin más que aplicar el logaritmo a sus miembros izquierdo y derecho. La ordenada en el origen y_interc = ln(a), la pendiente b y el coeficiente de correlación al cuadrado r_sqr de dicha recta de mejor ajuste se calculan como:

      // ajuste por mínimos cuadrados ln(freq) = ln(a) + b·ln(rank):       namespace bstat = boost::math::statistics;       auto ln = [](double d){ return std::log(d); };       auto const [y_interc, b, r_sqr] = bstat::simple_ordinary_least_squares_with_R_squared(                                            ranks | rv::transform(ln) | rn::to<std::vector>,                                            freqs | rv::transform(ln) | rn::to<std::vector>);

Tan sólo resta representar en escala log10-log10 tanto la curva de mejor ajuste de tipo potencia frequency = a·(ranking)^b, con a = exp(y_interc), como la nube de puntos (ranking, frequency). Sirviéndonos nuevamente de la biblioteca Matplot++, codificaríamos (punto 3 del enunciado del problema):

      // representación gráfica:       mp::gcf()->quiet_mode(mp::on);       auto const a = std::exp(y_interc);       auto const x = mp::logspace(0, 5);       mp::loglog(x, mp::transform(x, [&](double r){ return a*std::pow(r, b); }));       mp::hold(mp::on);       mp::loglog(ranks | rn::to<std::vector>, freqs, "o");       mp::title(fmt::format("freq = {:.1f}·pow(rank,{:.3f}) | R^2 = {:.3f}", a, b, r_sqr));       mp::xlabel("ranking position");       mp::ylabel("frequency");       mp::show();    }    catch (std::exception const& e) {       fmt::print(stderr, "exception -> what(): {}\n", e.what());    }

A modo de ejemplo, de emplear la obra Ulysses de James Joyce disponible en [1], obtendríamos las siguientes cinco mayores frecuencias de ocurrencia:

Frequency 15047: the

______________________________

Frequency 8259 : of

______________________________

Frequency 7215 : and

______________________________

Frequency 6530 : a

______________________________

Frequency 5033 : to

Deberemos descender hasta la frecuencia 829 para encontrar un vector de palabras de longitud superior a la unidad:

Frequency 829  : my, up

La gráfica obtenida sería, en tal caso (a = 30862.9 y b = -1.032):

Así, como cabía esperar, la palabra a la cabeza del ranking es aproximadamente dos veces más frecuente que la segunda palabra del ranking, tres veces más frecuente que la tercera palabra del ranking, etcétera.

Se proporciona a continuación un fichero CMakeLists.txt para la construcción del proyecto anterior, asumiendo que su código esté contenido en un fichero de nombre main.cpp:

   cmake_minimum_required(VERSION 3.20)    project(zipf       VERSION 0.1.0       DESCRIPTION "code provided by dgvergel.blogspot.com"       LANGUAGES CXX    )    add_executable(${PROJECT_NAME} main.cpp)    target_compile_features(${PROJECT_NAME} PRIVATE cxx_std_20)    set_target_properties(${PROJECT_NAME} PROPERTIES       CXX_EXTENSIONS OFF       RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_BINARY_DIR}/debug       RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_BINARY_DIR}/release    )    #-------------------------------------------------------------------    # integración de bibliotecas de terceros:    find_package(Boost 1.77.0 REQUIRED)    find_package(fmt 8.0 REQUIRED)    find_package(Matplot++ 1.1.0 REQUIRED)    find_package(range-v3 0.11.0 REQUIRED)    target_link_libraries(${PROJECT_NAME} PRIVATE       Boost::headers       fmt::fmt       Matplot++::matplot       range-v3    )    #-------------------------------------------------------------------    # política de avisos:    if (CMAKE_CXX_COMPILER_ID MATCHES "MSVC")       target_compile_options(${PROJECT_NAME} PRIVATE /W3 /WX)    elseif (CMAKE_CXX_COMPILER_ID MATCHES "Clang|GNU")       target_compile_options(${PROJECT_NAME} PRIVATE -Wall -Wextra -pedantic -Werror)    endif()    #-------------------------------------------------------------------    # mensajes a emitir en construcción:    message("-- CMake Generator: ${CMAKE_GENERATOR}")      get_target_property(CFEATURES ${PROJECT_NAME} COMPILE_FEATURES)    message("-- Target compile features: ${CFEATURES}")      get_target_property(COPTIONS ${PROJECT_NAME} COMPILE_OPTIONS)    message("-- Target compile options: ${COPTIONS}")    

---

Referencias bibliográficas

1. Proyect Gutenberg - https://www.gutenberg.org/
2. Zipf's law - https://en.wikipedia.org/wiki/Zipf%27s_law
3. Boost - https://www.boost.org/
4. {fmt}lib - https://github.com/fmtlib/fmt
5. Matplot++ - https://github.com/alandefreitas/matplotplusplus
6. Range-v3 - https://github.com/ericniebler/range-v3
7. std::isspace - https://en.cppreference.com/w/cpp/string/byte/isspace
8. std::ispunct - https://en.cppreference.com/w/cpp/string/byte/ispunct