- Representando el atractor de Lorenz con C++20, Boost.Numeric.Odeint y Dlib
- Análisis dimensional con mp-units, regresiones no-lineales con Dlib.Optimization y visualización de datos con Matplot++
- Rectas de mejor ajuste por mínimos cuadrados (mp-units, Matplot++)
- Ceros de una función
- Comprobando la ley de Zipf
- Instalación de Matplot++ y mp-units con Mingw-w64
Introducción
Última actualización del post: 24 de octubre de 2021.
En este post, continuación natural de un artículo previo, estudiaremos el empleo de la biblioteca de unidades mp-units [2] como medio para lograr un análisis dimensional en tiempo de compilación, sin que ello represente sobrecoste alguno en tiempo de ejecución. A modo de aplicación, el artículo analizará un experimento real de laboratorio conducente a la determinación de la aceleración de la gravedad terrestre mediante tiros parabólicos. Su tratamiento de datos requerirá el empleo del algoritmo iterativo de Levenberg-Marquardt facilitado por la biblioteca Dlib [3] para la resolución de problemas de mínimos cuadrados no-lineales. Se hará uso, asimismo, de las funcionalidades de representación de datos ofrecidas por la biblioteca Matplot++ [4].
El problema de las unidades
Consideremos una partícula que, partiendo del reposo, inicie un movimiento 1-dimensional con aceleración a constante. El desplazamiento neto Δx experimentado desde su partida tras un intervalo de tiempo Δt vendrá dado por a·Δt2/2. De acuerdo con ello, podríamos codificar una función que calculase el desplazamiento a partir de la aceleración y el intervalo temporal en la forma siguiente:
auto uarm_displacement(double a, double t) -> double { return a*t*t/2.0; }
Incluso un ejemplo tan sencillo como éste exhibe varios de los peligros que deseamos mitigar en nuestro estudio y que derivan de emplear un mismo tipo primitivo básico (en este caso, double) para representar magnitudes diferentes (distancia, tiempo, velocidad, etcétera) o medidas de una misma magnitud expresadas en unidades distintas (metro vs milla, por ejemplo). En efecto, imaginemos una llamada a la función que, por error, invirtiese el orden de sus argumentos:
auto const a = 0.5; // implícitamente en metros por segundo cuadrado
auto t = 12.0; // implícitamente en segundos
auto x = uarm_displacement(t, a); // oops!
Esta incidencia pasaría fácilmente desapercibida, corrompiendo el cálculo numérico de distancias. De forma más señalada aún, la función anterior no protege al programador frente al uso indebido de unidades de medida incoherentes:
auto const a = 0.5; // en metros por segundo cuadrado
auto t = 1.0; // en minutos
auto x = uarm_displacement(a, t); // oops!
Podría darse el caso, incluso, de que la fórmula incluida en la función de desplazamiento fuese dimensionalmente incorrecta, de forma que el valor retornado no correspondiese con la magnitud esperada (en nuestro caso, una longitud):
auto uarm_displacement(double a, double t) { return a*t/2.0; } // oops, retorna velocidades
A continuación, analizaremos cómo una biblioteca de unidades puede ayudarnos a evitar estos y otros problemas similares en el ámbito de la programación científica.
Biblioteca mp-units
Entre las principales características de la biblioteca de análisis dimensional mp-units (versión 0.8.0) caben resaltar su diseño moderno adaptado a los últimos estándares del lenguaje --en particular, la introducción de módulos en C++20--, así como los mensajes de error de compilación comprensibles alcanzados gracias al uso de conceptos. Esta herramienta ha sido propuesta para su estandarización de cara a C++23/26 [5].
La biblioteca interpreta una cantidad como un objeto poseedor de:
- Dimensión. Así, por ejemplo, el Sistema Internacional (SI) reconoce siete dimensiones básicas (longitud, tiempo, masa, corriente eléctrica, temperatura, cantidad de sustancia e intensidad lumínica) a partir de las cuales poder construir, mediante combinaciones de potencias de las primeras, las dimensiones derivadas (como la aceleración, la fuerza o la energía).
- Unidad de medida, de nombre y notación aceptadas por convenio (como el amperio A para la intensidad de corriente o el tesla T para la inducción magnética dentro del SI).
- Magnitud o valor numérico, representada mediante un escalar (como int y double, o incluso tipos no-primitivos definidos por el programador).
La plantilla de clase units::quantity presenta tres parámetros asociados a los conceptos anteriores, siendo double el tipo de representación subyacente adoptado por defecto:
namespace units {
template<Dimension D, UnitOf<D> U, Representation R = double>
class quantity {
public:
using dimension = D;
using unit = U;
using rep = R;
// ...
};
}
La signatura de la función uarm_distance de la sección anterior podría ahora codificarse de modo que los argumentos pertenecieran a tipos diferentes --expresados en las unidades que les correspondan en el Sistema Internacional--, si bien permitiendo el mismo álgebra de operaciones que con el tipo primitivo double:
namespace uisq = units::isq;
namespace usi = uisq::si;
using len_si = usi::length<usi::metre>;
using tm_si = usi::time<usi::second>;
using spd_si = usi::speed<usi::metre_per_second>;
using acc_si = usi::acceleration<usi::metre_per_second_sq>;
auto uarm_displacement = [](acc_si a, tm_si t) -> len_si { return a*t*t/2.0; };
auto const a = acc_si{0.5};
auto t = tm_si{12.0};
auto x = uarm_displacement(a, t); // usi::length<usi::metre> de valor 36.0
auto y = uarm_displacement(t, a); // error: dimensión incorrecta de argumentos
La biblioteca permite la conversión implícita de cantidades correspondientes a una misma dimensión, siempre que ello no implique una pérdida de precisión. Se trata de una decisión de diseño que otras bibliotecas de unidades optan por no adoptar [5]. A modo de ejemplo, y definiendo previamente la macro UNITS_REFERENCES para la biblioteca:
using namespace usi::time_references;
auto const a = acc_si{0.5};
auto t = 1.0*min; // usi::time<usi::minute>
// convertimos implícitamente el tiempo t de minutos a segundos:
auto x = uarm_displacement(a, t); // usi::length<usi::metre> de valor 900.0
Un fallo de codificación en la fórmula del desplazamiento, análogo al señalado en la sección anterior, produciría ahora un error de compilación:
auto uarm_displacement(acc_si a, tm_si t) -> len_si
{
return a*t/2.0; // error: el valor retornado es una velocidad
}
De forma notable, mp-units integra conceptos --en el sentido definido por C++20-- asociados a una amplia variedad de magnitudes físicas, lo que permite generalizar de forma directa funciones como la anterior. En efecto, consideremos el siguiente cálculo de velocidades medias:
auto average_speed_1 = [](len_si s, tm_si t) -> spd_si { return s/t; };
auto average_speed_2 = [](uisq::Length auto s, uisq::Time auto t) -> uisq::Speed auto
{ return s/t; };
using namespace usi::international::references;
auto x = 100.0*mi; // usi::length<usi::international::mile>
auto t = 2.0*h; // usi::time<usi::hour>
auto v_1 = average_speed_1(x, t);
fmt::print("{}\n", v_1); // output: 22.352 m/s
auto v_2 = average_speed_2(x, t);
fmt::print("{}\n", v_2); // output: 50 mi/h
La función average_speed_2 se encuentra templatizada, estando los parámetros de la plantilla restringidos por los conceptos Length, Time y Speed. Dicha función admite --al contrario que average_speed_1-- longitudes y tiempos en unidades diversas, generando resultados expresados en unidades coherentes con aquéllas. Ello elimina las conversiones intermedias de los argumentos a metros y segundos que tienen lugar al invocar average_speed_1. De necesitar expresar la velocidad v_2 en metros por segundo, la conversión podría posponerse hasta el momento oportuno mediante una operación quantity_cast:
fmt::print("{}", units::quantity_cast<spd_si>(v_2)); // output: 22.352 m/s
Como podremos comprobar, mp-units sigue, en este y otros aspectos, la filosofía de la biblioteca <chrono>.
Caso de ejemplo: Tiro parabólico en el laboratorio
Supongamos que se ha procedido al lanzamiento del proyectil para una serie de ángulos distintos φ (medidos en grados sexagesimales), anotándose para cada tiro la velocidad vexp (en millas por hora) y el alcance r (en centímetros). Las ternas (φ, vexp, r) se encuentran registradas como líneas de un fichero CSV de nombre parabolic_throwing_data.csv. Éstas, separadas en dos columnas con el fin de reducir espacio de presentación, contienen los valores siguientes:
0, 5.68, 61.5 0, 7.99, 87.55, 5.66, 64.5 5, 7.99, 99.510, 5.64, 72.0 10, 7.99, 113.015, 5.61, 77.3 15, 7.94, 124.520, 5.59, 82.0 20, 7.85, 135.525, 5.55, 85.5 25, 7.83, 145.030, 5.50, 87.5 30, 7.90, 152.535, 5.50, 89.3 35, 7.81, 157.040, 5.53, 89.0 40, 7.85, 158.545, 5.48, 86.0 45, 7.70, 156.550, 5.53, 82.5 50, 7.83, 152.055, 5.48, 77.0 55, 7.83, 141.560, 5.50, 70.0 60, 7.83, 128.565, 5.50, 60.5 65, 7.74, 111.570, 5.53, 51.5 70, 7.85, 95.575, 5.50, 41.0 75, 7.78, 72.0
Dichos datos corresponden a d = 3.3 cm y h = 29.8 cm.
Codificaremos un programa que, de forma dimensionalmente correcta, calcule la aceleración de la gravedad g en el laboratorio a través de la fórmula anterior empleando el algoritmo iterativo para regresiones no-lineales de Levenberg-Marquardt. Para ello, usaremos la implementación de dicho método proporcionada por la biblioteca Dlib [3].
Procederemos, en primer lugar, a incluir las cabeceras estándar a emplear en el programa (para obtener el código completo de esta sección, bastará concatenar los bloques de código siguientes de fondo negro):
#include <algorithm>
#include <cmath>
#include <concepts>
#include <cstdlib>
#include <fstream>
#include <functional>
#include <numbers>
#include <sstream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <vector>
Nos serviremos, asimismo, de las siguientes cabeceras de las bibliotecas Dlib, Matplot++, mp-units y Range-v3:
#include <dlib/optimization.h>
#include <matplot/matplot.h>
#include <range/v3/core.hpp>
#include <range/v3/view/transform.hpp>
#include <units/format.h>
#include <units/math.h>
#include <units/generic/angle.h>
#include <units/isq/si/acceleration.h>
#include <units/isq/si/length.h>
#include <units/isq/si/speed.h>
#include <units/isq/si/international/speed.h>
Incluiremos, finalmente, la biblioteca terminal.hpp facilitada en un artículo anterior, con el fin de facilitar la lectura segura de datos por la terminal:
#include "terminal.hpp"
Consideremos una línea CSV almacenada en un objeto std::string. Definiremos una función variádica parse capaz de parsear dicha línea y retornar, en base a los tipos units::quantity<> que indiquemos como parámetros de la plantilla, una tupla con las cantidades físicas correspondientes:
namespace un = units;
namespace uisq = un::isq;
namespace usi = uisq::si;
template<typename T>
concept Readable = requires (std::istream& is, T& t){ is >> t; };
template<un::Quantity... Qs>
requires (Readable<typename Qs::rep> and ...)
[[nodiscard]] auto parse(std::string const& csv_line, char delim = ',') -> std::tuple<Qs...>
{
auto isstr = std::istringstream{csv_line};
auto set_value = [&](un::Quantity auto& q) -> void {
using Q = std::remove_reference_t<decltype(q)>;
if (auto str_data = std::string{}; std::getline(isstr, str_data, delim)) {
if (auto val = typename Q::rep{}; std::istringstream{str_data} >> val) {
q.number() = val;
}
else throw std::runtime_error{"format error in CSV entry"};
}
else throw std::runtime_error{"missing entries in CSV line"};
};
auto res = std::tuple<Qs...>{};
std::apply([&](auto& ...q){ (..., set_value(q)); }, res);
return res;
}
El concepto units::Quantity cuidará de que los parámetros de la plantilla sean, en efecto, cantidades en el sentido definido por la biblioteca mp-units. Por su parte, el concepto Readable garantizará que los tipos de representación subyacentes en dichas cantidades (simplemente de tipo aritmético en este artículo) puedan ser extraídos a partir de un objeto std::istream. Observemos, finalmente, que la función std::apply [6] de la cabecera <functional>, empleada en el código anterior en combinación con un fold expression [7], nos permite iterar los tipos almacenados en la tupla y almacenar en ella, según el orden en que aparezcan en la línea CSV, las distintas cantidades físicas. A modo de ejemplo, una llamada a la función parse de la forma
auto const csv_line = std::string{"15.7, 9.0"};
auto const [x, t] = parse<usi::length<usi::metre>, usi::time<usi::second>>(csv_line);
assert(x == usi::length<usi::metre>{15.7});
assert(t == usi::time<usi::second>{9.0});
parseará dos cantidades en el Sistema Internacional: una distancia x en metros y un tiempo t en segundos, de ser éstos los valores almacenados en el string csv_line.
La siguiente función permitirá la lectura segura por la terminal de una cantidad física cualquiera, sometida a un predicado unario a nuestra elección. En particular, observemos que la función añadirá al mensaje facilitado como primer argumento las unidades de medida que correspondan:
template<un::Quantity Q, std::predicate<typename Q::rep> Cond = std::identity>
[[nodiscard]] auto read(std::string_view message, Cond cond = Cond{}) -> Q
{
auto res = Q{};
auto const mssg_with_units = fmt::format("{} in {:%q}: ", message, res);
terminal::prompt(mssg_with_units, res.number(), cond);
return res;
}
A continuación, proporcionaremos una serie de alias convenientes tanto para las cantidades almacenadas en el CSV como para la aceleración de la gravedad a calcular. El resto del código se adaptaría de forma transparente ante cualquier futura modificación en las unidades de longitud y/o velocidad:
auto main() -> int
try {
// ángulo en grados sexagesimales:
struct degree : un::named_unit<degree, "deg", usi::prefix> { };
using ang_type = un::quantity<un::dim_angle<degree>, degree, double>;
auto deg_to_rad = [](ang_type const& a){
return un::angle<un::radian>{2*std::numbers::pi*a.number()/360};
};
using len_type = usi::length<usi::centimetre>;
using spd_type = usi::speed<usi::international::mile_per_hour>;
using acc_type = usi::acceleration<usi::metre_per_second_sq>;
using csv_info_type = std::tuple<ang_type, spd_type, len_type>;
Seguidamente, solicitaremos al usuario que indique la pequeña distancia de vuelo d desde el lanzamiento del proyectil hasta su paso por la barrera de luz, así como la altura inicial h del tiro parabólico:
auto gtr0 = []<std::floating_point T>(T v){ return v > T{}; };
auto const d = read<len_type>("distance d", gtr0);
auto const h = read<len_type>(" height h", gtr0);
Los mensajes anteriores aparecerán en la terminal acompañados de la unidad de medida asociada a len_type (en este caso, el centímetro), como así muestra la siguiente captura de pantalla:
Llegados a este punto, almacenaremos las ternas (φ, vexp, r) del fichero CSV dentro de un vector estándar de objetos csv_info_type de nombre data. Tendremos presente que el primer valor de la terna contiene el ángulo de lanzamiento en grados sexagesimales, el segundo su velocidad experimental de salida en millas por hora y el tercero el alcance horizontal del tiro medido en centímetros:
auto const data = /* IILE: vector de ternas (φ [deg], v_exp [mph], r [cm]) */ []{
auto csv = std::ifstream{"../../parabolic_throwing_data.csv"};
if (!csv) throw std::ios::failure{"unable to open the CSV file"};
return ranges::getlines(csv)
| ranges::views::transform([](std::string const& ln) -> csv_info_type {
return parse<ang_type, spd_type, len_type>(ln);
})
| ranges::to<std::vector>;
}();
El bloque de código anterior adopta un estilo esencialmente funcional gracias a la biblioteca Range-v3 (consúltese esta serie de artículos para más detalles). Las siglas IILE indican la invocación automática de la expresión lambda (véase este artículo acerca de esta técnica).
La siguiente expresión lambda auxiliar será empleada por el algoritmo de Levenberg-Marquardt. Dado un punto experimental (φ, vexp) y un valor de aceleración de la gravedad g (el parámetro a optimizar), la función retornará el alcance horizontal teórico r calculado de acuerdo con la fórmula indicada al inicio de esta sección:
auto hrange = [&](ang_type phi, spd_type v_exp, acc_type g) -> len_type {
auto const t = deg_to_rad(phi).number(),
c = std::cos(t),
s = std::sin(t);
auto const v0_sq = un::pow<2>(v_exp) + 2*g*d*s;
un::Dimensionless auto const chi = g*h/v0_sq;
return v0_sq*(s + un::sqrt(un::dimensionless<un::one>{s*s} + 2*chi))*c/g;
};
La ejecución de la siguiente expresión lambda proporcionará el valor óptimo de la aceleración de la gravedad g que minimiza el error estándar de los residuos (RSE):
auto const [g, rse] = /* IILE: algoritmo Levenberg-Marquardt */ [&]{
using parameter_type = dlib::matrix<acc_type::rep, 1, 1>;
auto residual = [&](csv_info_type const& info, parameter_type const& param){
auto const& [phi, v_exp, r] = info;
auto const g = acc_type{param(0)};
return (hrange(phi, v_exp, g) - r).number();
};
auto param = parameter_type{10.0};
auto const rss = 2.0*dlib::solve_least_squares_lm(
dlib::objective_delta_stop_strategy{1.e-7},
residual, dlib::derivative(residual), data, param);
return std::tuple{acc_type{param(0)}, len_type{std::sqrt(rss/(data.size() - 1))}};
}();
El modelo distingue a g como único parámetro, de tipo parameter_type. La lambda residual es la encargada de calcular los residuos para el problema de mínimos cuadrados. Ésta toma una terna experimental (φ, vexp, r) contenida en un objeto de tipo csv_info_type y compara su alcance horizontal r con aquél generado por nuestro modelo hrange evaluado en (φ, vexp, g). La biblioteca Dlib requiere que residual retorne una valor de tipo coma flotante [8], por lo que las operaciones que conlleven la pérdida temporal de unidades --algo considerado inseguro en este contexto-- han sido cuidadosamente tomadas en consideración. Finalmente, y ante la complicada expresión adoptada por la derivada parcial del residuo respecto a g, aproximamos numéricamente su valor a (residual(info,g+eps)-residual(info,g-eps))/(2*eps) --siendo eps=1.e-7-- a través del método dlib::derivative.
Así, en el caso de los valores numéricos proporcionados para el archivo CSV al inicio de esta sección --con d = 3.3 cm y h = 29.8 cm--, una operación final de impresión
fmt::print("Inferred g = {:%.2Q %q}\n", g);
fmt::print("Residual standard error = {:%.2Q %q}\n", rse);
nos proporcionará una aceleración de la gravedad g de valor 9.75 m/s2 con un RSE final de 2.17 cm.
Estos resultados se encuentran en excelente acuerdo con el valor esperado en el laboratorio, teniendo en cuenta que las medidas de los alcances r se encontraban afectadas de un error de ± 3 cm.
Finalmente, sirviéndonos de la biblioteca Matplot++, podremos visualizar la superficie de mejor ajuste r = hrange(φ, vexp, g) sin más que añadir unas pocas líneas de código adicionales:
// ángulos (p), velocidades (v) y alcances (r) numéricos sin unidades:
auto to_number = [&]<std::size_t N>(){
return data
| ranges::views::transform([](auto const& t){ return std::get<N>(t).number(); })
| ranges::to<std::vector>;
};
auto const p = to_number.operator()<0>();
auto const v = to_number.operator()<1>();
auto const r = to_number.operator()<2>();
auto const [pmin, pmax] = std::ranges::minmax_element(p); // iteradores a mín y máx
auto const [vmin, vmax] = std::ranges::minmax_element(v);
namespace mp = matplot;
mp::gcf()->quiet_mode(mp::on);
mp::xlabel(fmt::format("φ [{:%q}]", ang_type{}));
mp::ylabel(fmt::format("vexp [{:%q}]", spd_type{}));
mp::zlabel(fmt::format("r [{:%q}]", len_type{}));
auto f = [&](double a, double s){ return hrange(ang_type{a}, spd_type{s}, g).number(); };
mp::fsurf(f, {*pmin, *pmax}, {*vmin, *vmax})->face_alpha(.7);
mp::colorbar();
mp::hold(mp::on);
mp::stem3(p, v, r, "filled")->color("r").line_width(1.5).marker_face_color("w");
mp::hold(mp::off);
mp::show();
mp::view(90, 0);
mp::show();
mp::view(180, 0);
mp::show();
}
catch (std::exception const& e) {
fmt::print("exception caught: {}\n", e.what());
return EXIT_FAILURE;
}
Referencias bibliográficas
- Wikipedia - Mars Climate Orbiter - https://en.wikipedia.org/wiki/Mars_Climate_Orbiter
- mp-units - https://mpusz.github.io/units/
- Dlib - Optimization - http://dlib.net/optimization.html
- Matplot++ - https://alandefreitas.github.io/matplotplusplus/
- P1935R2 - A C++ Approach to Physical Units - http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2020/p1935r2.html
- Cppreference - std::apply - https://en.cppreference.com/w/cpp/utility/apply
- Cppreference - Fold expressions - https://en.cppreference.com/w/cpp/language/fold
- Dlib - Least Squares Abstract - http://dlib.net/dlib/optimization/optimization_least_squares_abstract.h.html#solve_least_squares
No hay comentarios:
Publicar un comentario