Etiqueta: Austra

Categorías
Austra

Destrucción creativa

El compilador de Austra sigue mejorando. Y la librería también, que todo hay que decirlo.

Finalmente, he podido añadir generación de valores aleatorios con AVX2, además de la que ya había implementado con AVX512. La dificultad residía en la necesidad de rotar circularmente un entero largo. AVX512 te da la instrucción. Pero es fácil ver que en AVX2 es posible hacer lo mismo con dos desplazamientos lógicos y una conjunción lógica.

El otro avance importante es la vectorización del algoritmo de Box-Muller para generar distribuciones normales. Ya teníamos un logaritmo natural para vectores de 256 y de 512 bits. Ahora añadimos una función que calcula simultáneamente el seno y el coseno para estos dos tipos de vectores. De hecho, hay una variante especial que es la que usa Box-Muller, en la que todos los valores están entre cero y dos veces pi, y es un poco más rápida. Añadí esta función, la utilicé en el código del generador normal y, finalmente, se ha usado en las clases de vectores, matrices y secuencias.

Operaciones in situ

El compilador intenta ahora optimizar expresiones como la siguiente:

vector1 + vector2 + vector3

¿Qué hay aquí para optimizar? A primera vista, parecería que esa expresión es candidata para reescribirse mediante un combinador lineal, que de hecho existe y se usa para algunas optimizaciones ya:

vec([1, 1, 1, vector1, vector2, vector3)

Se podría utilizar, pero este constructor añadiría tres multiplicaciones escalares innecesarias. Lo que sí nos ahorraríamos sería la creación de dos búferes temporales como resultado de las dos sumas de la expresión original. Lo que ha aprendido a hacer el compilador es a modificar la expresión original a este otro patrón, que es más general:

(vector1 + vector2).InplaceAdd(vector3)

InplaceAdd no se puede llamar directamente desde el lenguaje: es una instrucción peligrosa, porque sobrescribe el búfer del primer operando. El compilador tiene que detectar que el primer operando es de tipo temporal, y no va a utilizarse en el resto de la fórmula, o en el peor de los casos, en el resto de la sesión. Pero en el ejemplo mostrado, la primera suma es evidentemente una construcción temporal, independientemente del sitio de donde estemos sacando las tres variables usadas.

Naturalmente, podemos hacerlo mejor si ya sabemos qué son vector1 y sus dos amigos. Si son variables de sesión, no hay nada más que podamos hacer. Lo mismo ocurre si son parámetros de una función definida por el usuario o variables locales, introducidas por una cláusula let. Pero supongamos que la expresión original fuese la siguiente:

vec::random(10) + vector2 + vector3

En ese caso, el primer operando se crea en un constructor y sólo se usa en esa expresión. El compilador puede entonces crear el siguiente código equivalente:

vec::random(10).InplaceAdd(vector2).InplaceAdd(vector3)

No importa de dónde salgan vector2 y vector3. El búfer del primer operando puede ser reutilizado, y eso es lo que hacemos. El resultado de InplaceAdd, dicho sea de paso, también se puede sobrescribir con total seguridad.

Este tipo de optimizaciones es mejor que las haga un compilador, y no la librería. Es el compilador quien tiene toda la información sobre el uso de un operando. De momento, nos movemos sobre terreno seguro. Hay algunas operaciones que podrían optimizarse con este truco, pero de momento no lo hacemos. Por ejemplo, si una variable local sólo se utiliza una vez, no hay memoria compartida que tengamos que respetar. De momento, no contamos el número de usos de cada variable local.

Finalmente, hemos ampliado las optimizaciones que hacíamos sobre vectores reales a vectores complejos. Es un poco más de código, pero no se hace más lenta la compilación.

Se supone que la señorita de la imagen es Freya bailando sobre calaveras. Eran Kali y Shiva quienes practicaban este tipo de danza, pero estaba seguro de que la tía que la IA iba a generar para Freya iba a ser más guapa. ¿Y las calaveras? Pues están debajo. Pero los generadores de Inteligencia Artificial siguen teniendo problemas para generar manos humanas. De hecho, la imagen que utilizo para Austra la llamo en mi cabeza «la chica del eccema». Llevo meses haciendo retoques a sus manos, y todavía tienen problemas. Pero me gusta el tono del pelo, el óvalo facial y la postura. No se puede tener todo.

Categorías
C#

Xoshiro256**

God and Einstein, playing dice at Heaven

Aunque parezca que el título de la entrada lo tecleó mi gato, es en realidad el nombre de un algoritmo de generación de números aleatorio, inventado por Sebastiano Vigna y David Blackman (enlaces al final).

El nombre es una combinación de las operaciones principales del algoritmo: xor, shift, rotate. En realidad, hay toda una familia de algoritmos similares, con nombres que se parecen mucho. El xoshiro256** es, simplemente, el que ha adoptado .NET Core desde la versión 6 (implementación aquí).

Xoshiro256** es muy rápido, es robusto (aunque no te recomiendan que lo uses en aplicaciones de criptografía), y utiliza muy poca memoria. Como se puede ver en la implementación de .NET 8, sólo necesita cuatro variables de tipo ulong. Es decir, 32 bytes por cada instancia del generador.

Xoshiro vectorial

Austra utiliza generadores de números aleatorios a diestra y siniestra. ¿Es fácil escribir una versión de este algoritmo usando instrucciones AVX? La respuesta es: no, a no ser que uses directamente AVX512F. El problema es que este algoritmo realiza una rotación. Es curioso que los lenguajes modernos de programación no te den directamente esta operación como parte de los operadores de bit. C# tiene un >> y un <<, por ejemplo, pero las rotaciones tienes que buscarlas en la clase BitOperations. De todos modos, el problema es que AVX/AVX2 no tiene una operación SIMD de rotación. Podríamos simularla con shifts y máscaras, pero perderíamos parte de las ventajas de la vectorización. Hay implementaciones en GitHub de generadores aleatorios con AVX2, pero no dan una ganancia de velocidad destacable.

Hay otro problema con el que hay que tener sumo cuidado. Pongamos como ejemplo el constructor de la clase DVector que genera un vector con valores aleatorios:

/// 
/// Creates a vector filled with a uniform distribution generator.
/// 
/// Size of the vector.
/// A random number generator.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public DVector(int size, Random rnd);

Omito la implementación para simplificar. Lo que quiero subrayar es que este constructor recibe una instancia explícita de la clase Random. ¿Por qué? Pues porque al cliente de la librería puede necesitar resultados repetibles. Es decir, puede que para que un vector «aleatorio» tenga siempre los mismos resultados, podría ser que la instancia de Random se inicializase siempre con la misma semilla. Si no nos interesa usar una semilla, casi siempre utilizaremos Random.Shared, que es una propiedad estática de Random, que es única para cada hilo y que se crea por demanda.

Podríamos complicarnos la vida al escribir la versión vectorizada de este constructor, y averiguar cómo recuperar la semilla del generador que nos pasa el cliente para crear un generador vectorial con esa misma semilla. Pero:

  1. No tengo claro que recuperar esa semilla sea tarea fácil
  2. El generador aleatorio vectorial que he programado no tiene, de momento, la posibilidad de inicializarse con una semilla (y no me parece tampoco sencillo).

Por lo tanto, las rutinas que se han acelerado con el generador vectorial comprueban si el generador que han recibido del cliente es Random.Shared. Eso lo interpretamos como que al cliente no le interesa la repetibilidad de los resultados. Naturalmente, hay que verificar antes si AVX512F está disponible en el ordenador.

Cuando usamos el generador vectorial, los tiempos de ejecución se reducen a una cuarta o quinta parte. No se reducen a la octava parte mecánicamente porque en la mayoría de estas rutinas pesa más la asignación de memoria y el posterior llenado de la misma. Pero bajar el tiempo de ejecución a la cuarta parte me parece meritorio.

Enlaces

Estos son los enlaces a la clase de Austra.Library que genera ocho números de golpe por llamada, y a la página de Sebastiano Vigna, el autor del algoritmo que hemos vectorizado. Mi implementación, con toda seguridad, está abierta a mejoras. Por ejemplo, no me convence la forma en que tengo que convertir un valor ulong en otro de tipo double, porque no es vectorial. Debe existir algo más eficiente, pero de momento no lo he encontrado. De todas maneras, ahí tiene mi implementación, por si le es útil en alguno de sus proyectos.

Queda también pendiente la vectorización de números aleatorios provenientes de una distribución normal. Austra utiliza el método de Box-Muller, pero este método exige el cálculo de un logaritmo y de senos y cosenos. No es tarea imposible, pero tengo primero que proporcionar una rutina que genere vectorialmente senos y cosenos y, quizás antes, decidir si me interesa realmente el método de Box-Muller, o si merece la pena ir directamente a una implementación vectorial del algoritmo del zigurat. Todo, en su debido momento.

Categorías
Austra

Funciones definidas por el usuario

Vamos al grano, o, como diría Haskell B. Curry, «let’s cut to the chase». Esto ya se puede hacer en AUSTRA (en cuanto libere la próxima versión):

let mcd(a, b: int): int =
    let m = a % b in iff(m = 0, b, mcd(b, m)) in
        mcd(80, 140)

Esta es una versión recursiva del máximo común divisor, calculado no con restas, sino con el módulo de la división. Observaciones importantes:

  • Estoy declarando la función como una función local del script. Me falta permitir ahora la declaración de funciones como si fuesen definiciones paramétricas. No lo he hecho todavía porque la implementación de este tipo de funciones se realiza mediante lambdas, que se asignan a una variable local. La recursividad es posible porque la variable local está disponible, con toda la gloria de su prototipo, cuando se compila el cuerpo de la función. Cuando se trate de una definición, probablemente use un truco parecido, pero no es tan inmediato (internamente) como cuando defino una función como parte de una cláusula let.
  • Como se trata de una función recursiva, observe que he definido su tipo de retorno explícitamente. Si no, cuando el compilador encuentre la llamada recursiva a mcd va a tener que volverse loco infiriendo cuál es el tipo de retorno. Ese tipo de inferencias es posible, pero ahora mismo el compilador no está preparado para ello. En mi defensa, recuerde que incluso el gran F# necesita el modificador rec para declarar funciones recursivas. Y F# sí es un lenguaje funcional con todas las de la ley.
  • En este caso, hay una cláusula let anidada dentro de la definición de función. Austra tenía una regla para «aplanar» siempre estas cláusulas en el nivel superior, pero aquí me interesa violar la regla. Eso, o tengo que calcular dos veces el módulo de los dos parámetros.

Otro ejemplo de función recursiva, que utiliza también una cláusula let anidada, aunque esta vez, de tipo diferente:

let fact(n: int) =
    let f(n, acc: int): int = iff(n <= 1, acc, f(n - 1, n * acc)) in
        f(n, 1);
fact(10);
  • Esta es la archifamosa función factorial, pero en vez de definirla en la forma más simple, la defino con una función auxiliar que permite recursividad "por la cola".
  • Aunque el factorial es normalmente recursivo, esta vez no se llama a sí mismo, y no hay que declarar explícitamente el tipo de retorno.
  • Por el contrario, la función interna f sí se llama a sí misma, y hay que tener cuidado con su prototipo.

Claro está que podemos declarar funciones sin tanta fanfarria. Por ejemplo, el factorial es mejor programarlo así:

let fact(n: int) = [x in 2..n].prod;
fact(10);

Recuerde que empezamos con una secuencia de enteros, disfrazada de constructor de secuencias, por lo que no necesitamos poner los valores de la secuencia en un vector, ocupando memoria.

Por cierto, mire lo que podemos hacer ahora con las secuencias:

let collatz(n: int) =
    iseq::unfold(1000000, n, x => iff(x % 2 = 0, x / 2, 3x + 1))
    .until(x => x = 1);
collatz(137)

Esta función genera la secuencia de la conjetura de Collatz para el número entero suministrado como parámetro. Es un problema interesante, y aparentemente fácil, pero que aún no está resuelto. La novedad es que ahora Austra tiene un método Until y un método While que sirven para este tipo de acrobacias. Nuestra función mcd podría haberse también programado sin recursión, mezclando unfold con uno de estos métodos.

Categorías
C# FinTech

Secuencias

Esto es una perogrullada: en un lenguaje de programación funcional, las cosas se suelen hacer de forma diferente. Y probablemente escribir «perogrullada» sea una pedantería. Dejo el párrafo vivo por pereza, pero no me lo tenga en cuenta.

Factorial

Austra no pretende ser un lenguaje de programación funcional. Ni siquiera pretende ser «Turing-completo». En algún punto del camino, he pensado en diseñar un lenguaje «de verdad», e incluso en saltarme la maquinaria de .NET para compilarlo a código nativo directamente. Pero falta un poco para eso.

De todas maneras, hay que intentar que se puedan hacer todas las cosas posibles en Austra. Por ejemplo, ¿cómo calculas un factorial, si no tienes todavía funciones recursivas, y no quieres introducir bucles? Hasta ahora, la solución era parecida a ésta:

vec(10, i => i + 1).prod

No estamos definiendo una función, sino que estamos usando un parámetro a dedo, pero funciona. Creamos un vector de 10 entradas, lo llenamos con los números del 1 al 10, y multiplicamos todos los elementos del vector. Incluso podemos «optimizar» un poco la expresión:

vec(8, i => i + 2).prod * 2

Es absurdo dejar el 1 en el vector, y ya que tenemos diez elementos, en vez de dejarlos nueve, quito también el 2. Eso lo hago porque sé que internamente el producto mete los ocho elementos en dos registros AVX y no tengo que manejar el elemento nuevo aparte. El problema es que, de todas maneras, estamos pidiendo memoria para el vector.

En Austra 2.0 (que ya tiene soporte para .NET 8 y AVX512, dicho sea de paso), ya están implementadas las secuencias de valores reales, y puedes hacer esto otro, que es más natural:

seq(2, 10).prod

La nueva clase se llama seq, y es una copia descarada del diseño de enumerables y LINQ (o de los streams de Java, o de las «list comprehensions» de tantos otros lenguajes funcionales). Ahora tenemos un seq, pero luego vendrán iseq y cseq, para secuencias de enteros y complejos. Observe también que he acortado vector a vec, para usar cvec en vez de complexvector, y para que un futuro vector de enteros se pueda llamar ivec a secas.

Las secuencias pueden hacer casi todas las cosas que hacen los vectores, y en la mayoría de los casos, generarlas es más sencillo: con el vector, usamos una función lambda. También permiten ahorrar código porque aplican los mismos trucos que LINQ para objetos. Por ejemplo:

-- Esta es una secuencia que divide el intervalo
-- de 0 a 2*pi en 360 trozos.
-- Multiplicamos cada número por dos...
seq(0, τ, 360) * 2
-- .. pero internamente, lo transformamos en esto:
seq(0, 2τ, 360)

Recursividad eficiente

Y esto último todavía no está implementado, pero es bueno pensar en estas cosas antes de lanzarse a programarlas. ¿Cómo definiríamos una función para el factorial que fuese medianamente eficiente, pero usando recursividad? Lo que se suele hacer en un lenguaje funcional, pero con una notación al estilo Austra:

def fact(n: int) =
  let f = (n, acc: int) =>
    if n = 1 then acc else f(n - 1, acc * n) in
      f(n, 1)

El truco está en usar una función auxiliar, que al estilo Austra sería una lambda en una cláusula let. La función auxiliar es recursiva por la cola, por lo que el propio JIT de .NET la puede convertir en un bucle. La recursividad en lambdas va a ser problemática, porque estoy usando directamente Linq.Expressions para generar código, pero hay un truco sencillo, que es el que voy a usar en las primeras pruebas: declarar una variable de tipo lambda, asignarle un puntero vacío, y a continuación, asignarle una expresión lambda que ya podrá usar la variable recursivamente. La alternativa es usar un combinador Y. Es un truco que viene, precisamente, del cálculo lambda, sobre el que escribiré en algún momento. No es que sea cómodo o eficiente, pero es interesante.

De todas maneras, me he dado cuenta, revisando la documentación, que en F# hay que anunciar cuando una función va a ser recursiva. Una alternativa que tengo que pensar es generar un «método dinámico» en estos casos, con el esfuerzo adicional de generar directamente el código IL. Pero tengo ya experiencia de generar IL con Freya, y puedo reutilizar código.

Categorías
C# FinTech

Ostara

Todavía es work-in-progress, pero ya tenemos una versión open-source de una aplicación de escritorio, en WPF, para Austra. El código está ya en GitHub.

El editor de código es AvalonEdit. Reconoce la sintaxis del lenguaje de fórmulas e implementa una versión bastante decente de completamiento de código. Esto es lo que aporta la versión WPF respecto a la aplicación de consola: es mejor para aprender a usar el lenguaje. Los gráficos están hechos con OxyPlot, de momento.

Faltan cosas, tanto en la librería como en la aplicación, para darme por satisfecho, pero todo es cuestión de tiempo. Me gustaría, sobre todo, terminar de definir un mecanismo genérico de carga de series desde fuentes de datos externas.

Quiero también actualizar el código de la función de autocorrelación, e implementarla usando la transformada de Fourier que ya viene incluida. E incluir, de una puñetera vez, el modelo ARIMA completo y algo de GARCH y familia. De momento, sólo está incluido el modelo de series autoregresivas. También está pendiente una implementación con AVX de generación de números aleatorios.

Nota: autocorrelación ya actualizada. Como hay que añadir ceros al final de la serie para evitar que se cuelen las correlaciones cíclicas, la FFT se aplica a un número de muestras que es potencia de dos, y se utiliza el sub-algoritmo más eficiente.

Categorías
C#

Austra en GitHub

Tras muchos nervios y algún titubeo, ya he subido AUSTRA a GitHub:

Está el código completo de la librería (proyecto Austra.Library), el compilador del lenguaje de fórmulas (Austra.Parser), una aplicación de consola para ejecutar fórmulas en el lenguaje (Austra.REPL), y un par de proyectos más, para tests y para benchmarks.

Todo esto es work in progress, por supuesto. Me queda pendiente subir los dos packages (librería y compilador) a NuGet. Lo que me urgía, no obstante, era dar acceso al código fuente. Está subido con licencia MIT, que es la que entiendo que es más permisiva. Recuerde que el objetivo final de todo esto es docente: tener una referencia online de cómo usar código SIMD, punteros y esas cosas, y servir como base para mi futuro libro (Unsafe C#).

Categorías
C#

Austra PolySolve

C’est la vie: elle est dure et souvent courte.

Es imposible escribir sobre ecuaciones algebraicas sin mencionar a Évariste Galois, quien no sólo cerró una larga historia de intentos de solución de este tipo de ecuaciones, sino que además tuvo una vida corta y trágica.

Dicen que los egipcios y los babilonios eran capaces de resolver las ecuaciones de segundo grado y, por supuesto, las lineales. Las de tercer y cuarto grado tuvieron que esperar al Renacimiento Italiano. Y luego, la teoría se estancó: nadie era capaz de resolver una ecuación de quinto grado general; sólo algunas versiones restringidas.

Lagrange estuvo a punto de demostrar que las ecuaciones de quinto grado y superiores no tenían una solución general. Fue Ruffini quien lo consiguió, aunque con algunos pequeños errores, que más tarde corrigió Abel. De todos modos, la teoría que Galois puso por escrito en 1830, cuando sólo tenía 18 años, tenía un alcance mucho mayor, y ofrecía una estructura más completa y versátil para estudiar las ecuaciones algebraicas. Mientras que el teorema de Ruffini-Abel se centraba en la solubilidad de ecuaciones algebraicas por medio de funciones elementales (como exponentes, logaritmos, etc.), la criatura que inventó nuestro héroe introdujo los llamados grupos de Galois para analizar la solubilidad y las simetrías en un marco más general. No solo abordó la solubilidad de ecuaciones, sino que la teoría es aplicable también en otras áreas de las matemáticas, como la teoría de números y la geometría algebraica.

Por desgracia, el artículo presentado en 1830 por nuestro héroe no tuvo el éxito que merecía. Cauchy le pasó la patata caliente a Poisson, y Poisson no entendió ni papa del tema. El rechazo cabreó a Galois, pero aprovechó para rescribir la demostración, y fue esta modificación la que finalmente fue reconocida, entre otros, por el propio Cauchy. Eso sí: tras la muerte de Galois…

Galois tenía una cabecita muy loca, le pirraba la política, y para colmo, estaba algo deprimido por la muerte de su padre. 1832 fue un año difícil para el chico. Estuvo en prisión un par de veces, por ciscarse en Louis Philippe, el penúltimo rey de Francia. Al salir de la cárcel por segunda vez, se enredó en un duelo absurdo, teóricamente por una coquette, aunque no es descartable que todo fuese una trampa de sus enemigos políticos. Se pasó la noche anterior al duelo escribiendo una carta sobre sus últimos avances matemáticos. Al día siguiente, interpuso su abdomen en la trayectoria de una bala, y sus contrincantes lo dejaron tirado sobre la hierba como a un chien. Un transeúnte lo vio y lo llevó al hospital, pero al día siguiente se reunió con su Creador, probablemente por culpa de una peritonitis.

And all the king’s horses and all the king’s men, couldn’t put Évariste together again.

Raíces reales

Mi curiosidad por estos temas viene de cuando tenía unos diez u once años: encontré la solución razonada de las ecuaciones de segundo grado, en un libro de electrónica, y me dio por intentar resolver por mi cuenta el problema de las cúbicas. No lo conseguí. Tropecé por casualidad con la sustitución de Vieta, pero no conseguí algo mucho más sencillo: cómo eliminar el término cuadrático, que suele ser el primer paso de la solución. Pero compré un libro que explicaba la fórmula cúbica y la cuártica, y me convertí en un friqui de las mates.

Volví a enredar con ecuaciones algebraicas en 2005, cuando me dio por probar si se podía escribir un ray tracer decente en C#. Es bastante frecuente tener que resolver ecuaciones de tercer y cuarto grado para calcular intersecciones entre rayos de luz y determinados tipos de objetos. La particularidad es que, en este contexto, sólo se necesitan las soluciones reales. Cuando las cosas se ponen feas, existe una técnica para encontrar las raíces reales de cualquier polinomio utilizando las secuencias de Sturm. Naturalmente, este algoritmo es una sólo una aproximación iterativa.

Raíces complejas, todas

Cuando estás escribiendo una librería como AUSTRA, te interesa resolver el problema más general, que es encontrar todas las raíces, ya sean complejas o reales, de un polinomio arbitrario. ¿Se acuerda de los valores propios? El método que utilizo en AUSTRA está basado en ellos.

Supongamos que queremos resolver la ecuación:

$$c_0 + c_1x + c_2x^2 + \cdots + c_{n-1}x^{n-1} + x^n$$El término de mayor grado está normalizado para que su coeficiente sea la unidad. Ahora formamos la siguiente matriz, conocida como «matriz de Frobenius»:

$$F=\pmatrix{0&0&0&\cdots&0&-c_0\cr
1&0&0&\cdots&0&-c_1\cr
0&1&0&\cdots&0&-c_2\cr
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\cr
0&0&0&\cdots&1&-c_{n-1}}$$Nos planteamos entonces encontrar los valores propios de $F$, que deben cumplir esta igualdad:

$$F\vec{v} = \lambda\vec{v}$$donde $\vec{v}$ es uno de los vectores propios. Si reordenamos los términos, nos encontramos con esto:

$$(F-\lambda I)\vec{v}=0$$donde $I$ es la matriz identidad. Para que esta igualdad se cumpla, el determinante de $(F-\lambda I)$ debe ser igual a cero. Y resulta que el determinante de $(F-\lambda I)$ es, precisamente, la ecuación original. Qué listo era Frobenius.

AUSTRA tiene un método muy eficiente para calcular valores propios, incluso en casos como estos, en los que la matriz no es simétrica. Por lo tanto, para resolver un polinomio primero lo normalizamos, luego creamos su matriz de Frobenius, y finalmente calculamos sus valores propios. La función global polySolve es la que se encarga de la implementación, en el lenguaje funcional de AUSTRA. En la aplicación de consola, podemos teclear lo siguiente:

> set v = [5, 4, 3, 2, 1]
ans ∊ ℝ(5)
5  4  3  2  1
> polysolve(v)
ans ∊ ℂ(4)
<0,137832; 0,678154>   <-0,537832; 0,358285>
<0,137832; -0,678154>  <-0,537832; -0,358285>

polySolve puede recibir tanto un vector con los coeficientes, como los coeficientes sueltos. En este caso, estamos resolviendo la ecuación de cuarto grado $5x^4+4x^3+3x^2+2x+1=0$, y el resultado son cuatro números complejos, conjugados a pares.

¿Quiere comprobar que las raíces son realmente soluciones de la ecuación? Hagamos esto entonces:

> polysolve(v).map(c => polyeval(c, v))
ans ∊ ℂ(4)
<-1,33227E-15; -7,77156E-16>   <-1,33227E-15; 4,44089E-16>
 <-1,33227E-15; 7,77156E-16>  <-1,33227E-15; -4,44089E-16>

polyEval sirve para evaluar un polinomio para un argumento complejo o real, y el método map crea un nuevo vector complejo calculando sus entradas con una función lambda, al estilo del método Select de LINQ. Incluso tenemos una función poliDerivative que, con los mismos argumentos que polyEval, evalúa la derivada del polinomio que le pasamos en la coordenada que le digamos. Esto, a su vez, es muy conveniente para buscar raíces reales con el método de Newton-Raphson… que también ofrece AUSTRA (función solve, a secas).

¿Librería o lenguaje?

Por supuesto, todo esto sería igual de fácil, eficiente y elegante, o quizás un poco más, si simplemente enchufásemos el package Austra.Library a un proyecto en .NET Core y utilizásemos directamente las clases. Pero he querido mostrar este ejemplo en el lenguaje de fórmulas de AUSTRA como demostración de un caso de uso importante para el lenguaje: es una forma rápida y sencilla de poner a prueba la funcionalidad de la librería.

Y hay más casos de uso, que explicaré más adelante.

Categorías
C#

El Gran Secreto de los Complejos

Al grano: el Gran Secreto de los Números Complejos es que, si quieres utilizar instrucciones AVX para acelerar los cálculos, la mejor forma de representarlos no es la que todos imaginamos: la parte real y, a continuación, la parte imaginaria.

Partamos de una regla básica de las instrucciones vectoriales:

  • Es mejor manejar estructuras de arrays que arrays de estructuras.

Observe que ésta es una píldora difícil de tragar en la Programación Orientada a Objetos. Complex es una clase que ya está (bien) definida en System.Numerics, pero para simplificar la explicación, voy a fingir que la definimos nosotros. Con la POO en mente, comenzaríamos definiendo la estructura, junto con sus métodos, y la haríamos probablemente implementar algunas interfaces, por completitud, en este plan:

public readonly struct Complex
{
    public double Real { get; }
    public double Imaginary { get; }

    public Complex(double re, double im) =>
        (Real, Imaginary) = (re, im);

    // Y así, sucesivamente…
}

Si quisiéramos entonces un vector de números complejos, haríamos algo parecido a esto:

public readonly struct ComplexVector
{
    private readonly Complex[] values;

    public unsafe ComplexVector(Complex[] values) =>
        this.values = values;

    // Y así, sucesivamente…
}

Pues bien, ahora llego yo (o la sacrosanta Realidad, si lo prefiere hacer menos personal) y le digo que la mejor forma de programar un vector de complejos, al menos si queremos acelerarlo con AVX, es la siguiente:

public readonly struct ComplexVector
{
    private readonly double[] re;
    private readonly double[] im;

    // Omito verificaciones de igual longitud
    // para simplificar el ejemplo.
    public ComplexVector(double[] re, double[] im) =>
        (this.re, this.im) = (re, im);

    public unsafe ComplexVector(Complex[] values)
    {
        this.re = new double[values.Length];
        this.im = new double[values.Length];
        fixed (double* p = re, q = im)
            for (int i = 0; i < values.Length; i++)
                (p[i], q[i]) = values[i];
    }

    // Y así, sucesivamente…
}

Podemos dejar la estructura Complex original: nos es útil. Pero al representar la lista de complejos, es mejor que cada campo vaya en su propia lista. Es cierto también que deberíamos utilizar AVX para convertir un array tradicional de complejos en un vector: existe la posibilidad, pero no lo voy a mostrar aquí, para simplificar. He omitido también un método de extensión que he añadido en una clase estática para poder "deconstruir" fácilmente un complejo en sus componentes. No tiene mucha trascendencia, pero ahí va, para que no haya tantos espacios en blanco:

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Deconstruct(
    this Complex cmp, out double re, out double im) =>
    (re, im) = (cmp.Real, cmp.Imaginary);

¿Cómo se da cuenta uno?

¿Cómo se da cuenta uno de estas cosas? StackOverflow está lleno de consejos de este tipo, escritos por personas que ya lo han sufrido en sus carnes. Pero a uno no se le enciende la bombilla hasta que se pega con su propia pared. En este caso, fue intentando acelerar una Transformada Discreta de Fourier para AUSTRA. El código sin acelerar usaba los complejos como se han usado casi toda la vida: cada real con su parte imaginaria. A priori, cuando uno no conoce bien AVX, se imagina que será relativamente sencillo manejar dos números complejos dentro de un vector de cuatro valores reales, y que el conjunto de instrucciones va a estar de tu parte. Craso error.

Al final, tuve que limitarme a acelerar algunas partes que, oportunistamente, "se dejaban", casi siempre con código SSE y vectores de sólo 128 bits. Lo normal, cuando he acelerado otros algoritmos, ha sido reducir los tiempos de ejecuciones de cuatro hasta incluso ocho o diez veces. En este caso, la mejora sólo ha sido la mitad, en la mayoría de los casos, y en algunos, la tercera parte. Como resultado, tengo pendiente replantearme todo el asunto de la Transformada Discreta de Fourier, pero utilizando listas paralelas para los reales y los imaginarios.

Quiero que vea, de todas maneras, lo sencillo que es usar la técnica de estructura de listas en vez de listas de estructuras. El siguiente método es el producto escalar de dos vectores complejos, representados como Dios manda:

public static unsafe Complex operator *(
    ComplexVector v1, ComplexVector v2)
{
    if (v1.Length != v2.Length)
        throw new VectorLengthException();
    fixed (double* pr = v1.re, pi = v1.im,
                   qr = v2.re, qi = v2.im)
    {
        double sumRe = 0, sumIm = 0;
        int i = 0, size = v1.Length;
        if (Avx.IsSupported)
        {
            Vector256<double> accRe = Vector256<double>.Zero;
            Vector256<double> accIm = Vector256<double>.Zero;
            for (int top = size & ~3; i < top; i += 4)
            {
                var vpr = Avx.LoadVector256(pr + i);
                var vpi = Avx.LoadVector256(pi + i);
                var vqr = Avx.LoadVector256(qr + i);
                var vqi = Avx.LoadVector256(qi + i);
                accRe = Avx.Add(accRe,
                    Avx.Multiply(vpr, vqr)
                       .MultiplyAdd(vpi, vqi));
                accIm = Avx.Add(accIm,
                    Avx.Multiply(vpi, vqr)
                       .MultiplySub(vpr, vqi));
            }
            sumRe = accRe.Sum();
            sumIm = accIm.Sum();
        }
        for (; i < size; i++)
        {
            sumRe += pr[i] * qr[i] + pi[i] * qi[i];
            sumIm += pi[i] * qr[i] - pr[i] * qi[i];
        }
        return new(sumRe, sumIm);
    }
}

Si no lo recuerda del álgebra, cuando se trata de vectores complejos, el producto escalar usa la conjugada del segundo operando. Esto es: la parte imaginaria del segundo operando invierte su signo. Esto es lo que permite que el producto escalar de un vector consigo mismo sea un valor real.

El código vectorial tiene una correspondencia uno a uno con el código escalar que maneja la parte de los arrays que puede sobrar al final. Si repasa la entrada anterior, la del algoritmo de Welford, notará el mismo patrón: a pesar de que el algoritmo escalar es bastante oscuro, es relativamente sencillo convertir esa parte en código vectorial. La parte más complicada de la entrada pasada era cuando había que mezclar los cuatro acumuladores individuales. Es el mismo problema que hemos visto ya unas cuantas veces cuando calculamos un producto escalar: acumular es sencillo. Lo complicado es sumar después los cuatro acumuladores.

El Misterioso Constructor Postergado

Para no dejar demasiados cabos sueltos, aquí tiene una posible implementación del código que reparte partes reales e imaginarias a sus respectivos arrays. Mi primer impulso fue utilizar Avx2.GatherVector: leer cuatro partes reales saltándome las imaginarias, y luego leer cuatro partes imaginarias. Pero, por desgracia, el tiempo de ejecución del constructor se disparó al doble. No hay forma humana de predecir estas cosas, que no sea prueba, benchmark y error.

La versión que funciona, y que reduce casi a la mitad el tiempo de la versión de más arriba, lee cuatro complejos en dos vectores de 256 bits. Lo primero que hace es usar Avx.Shuffle para "barajar las cartas" y juntar todas las partes reales en un mismo vector de 256 bits, y las imaginarias en otro. No me pida que calcule estas cosas de memorias. Cuando me tocan estos marrones, tengo todavía que pillar un cuaderno y lápiz, e irme a páginas como ésta, para repasar los diagramas. He visto también que el Shuffle se puede sustituir también por llamadas a UnpackHigh/UnpackLoad. Es probable que estas llamadas den mejores tiempos, pero no me ha dado tiempo a hacer la prueba.

El problema de Shuffle (y de las alternativas mencionadas) es que te deja los números en el orden [1, 3, 2, 4]. Si no es importante respetar el orden, se pueden quedar así. Pero si hay que reordenar los elementos, hay que usar Avx2.Permute4x64 para ello. En general, AVX intenta, dentro de lo posible, no pasar valores de una mitad a la otra mitad del vector. Hay que usar cosas introducidas en AVX2 para conseguirlo. Por ese motivo, el constructor verifica si Avx2.IsSupported antes de lanzarse al río:

public unsafe ComplexVector(Complex[] values)
    : this(values.Length)
{
    fixed (double* p = re, q = im)
    fixed (Complex* r = values)
    {
        int i = 0;
        if (Avx2.IsSupported)
        {
            for (int top = values.Length & ~7; i < top; i += 4)
            {
                var v1 = Avx.LoadVector256((double*)(r+i));
                var v2 = Avx.LoadVector256((double*)(r+i+2));
                Avx.Store(p + i, Avx2.Permute4x64(
                    Avx.Shuffle(v1, v2, 0b0000), 0b11011000));
                Avx.Store(q + i, Avx2.Permute4x64(
                    Avx.Shuffle(v1, v2, 0b1111), 0b11011000));
            }
        }
        for (; i < values.Length; i++)
            (p[i], q[i]) = r[i];
    }
}

Números: en mi máquina, un i9-11900K, crear un ComplexVector directamente a partir de un array de 1024 complejos, tardaba más o menos un milisegundo. Con las mejoras AVX2, tarda 650 microsegundos. Casi la mitad. Y lo mejor, para mi gusto, es que no he tenido que usar paralelismo con tareas. El usuario de la librería ya usará ese paralelismo cuando lo considere necesario, y tendrá las manos más libres.

Como regalo, le dejo la conversión inversa: de vector complejo a array de complejos:

public unsafe static explicit
    operator Complex[](ComplexVector v)
{
    Complex[] result = new Complex[v.Length];
    fixed (double* p = v.re, q = v.im)
    fixed (Complex* r = result)
    {
        int i = 0;
        if (Avx2.IsSupported)
        {
            for (int top = v.Length & ~3; i < top; i += 4)
            {
                var vr = Avx.LoadVector256(p + i);
                var vi = Avx.LoadVector256(q + i);
                Avx.Store((double*)(r + i),
                    Avx2.Permute4x64(Avx.Permute2x128(
                    vr, vi, 0b0010_0000), 0b11_01_10_00));
                Avx.Store((double*)(r + i + 2),
                    Avx2.Permute4x64(Avx.Permute2x128(
                    vr, vi, 0b0011_0001), 0b11_01_10_00));
                }
            }
        for (; i < result.Length; i++)
            r[i] = new(p[i], q[i]);
    }
    return result;
}

El tiempo de ejecución se reduce "sólo" a las tres cuartas partes, pero yo creo que merece la pena.

Categorías
C#

El algoritmo de Welford

En la entrada sobre la varianza, vimos que podíamos tener problemas de estabilidad numérica si intentábamos calcular la varianza en una sola pasada sobre los datos, usando inocentemente la definición matemática. La solución de entonces fue usar un algoritmo de dos pasos: calcular la media en el primer paso, y en el segundo, calcular la varianza de la muestra menos la media. Había otra posibilidad: usar el primer valor de la secuencia como estimado malo de la media, y restar ese valor a las sucesivas muestras.

¿Podríamos hacer algo mejor si corrigiésemos el estimado de la media sobre la marcha? Resulta que se puede, y el primero en darse cuenta fue Welford, allá por el 1962. Donald Knuth incluyó el algoritmo en el segundo tomo de The Art of Computer Programming. El algoritmo original sólo calculaba la media y la varianza sobre la marcha, pero Timothy Terriberry, en 2007, lo amplió para que calculase momentos superiores. Este algoritmo está implementado, por ejemplo, en Math.Net Numerics, aunque la implementación es mejorable.

¿Qué nos da?

En AUSTRA, la clase que implementa este algoritmo se llama Accumulator. Hay también una clase simplificada, SimpleAccumulator, que sólo calcula los dos primeros momentos, con el beneficio evidente de tener que ejecutar menos trabajo.

La definición de Accumulator, junto con su constructor principal y los campos y propiedades de almacenamiento, es la siguiente:

/// Calculates statistics by adding samples.
public sealed class Accumulator
{
    ///Minimum value.
    private double min = double.PositiveInfinity;
    ///Maximum value.
    private double max = double.NegativeInfinity;
    ///Estimated mean.
    private double m1;
    ///Accumulated second moment.
    private double m2;
    ///Accumulated third moment.
    private double m3;
    ///Accumulated fourth moment.
    private double m4;

    ///Gets the total number of samples.
    public long Count { get; private set; }

    ///Creates an empty accumulator.
    public Accumulator() { }

    /* … */
}

La información que nos interesa se obtiene a través de estos campos, por medio de propiedades calculadas:

///Returns the minimum value.
public double Minimum => Count > 0 ? min : double.NaN;
///Returns the maximum value.
public double Maximum => Count > 0 ? max : double.NaN;
///Gets the sample mean.
public double Mean => Count > 0 ? m1 : double.NaN;
///Gets the unbiased variance.
public double Variance =>
    Count < 2 ? double.NaN : m2 / (Count - 1);
///Gets the unbiased standard deviation.
public double StandardDeviation =>
    Count < 2 ? double.NaN : Sqrt(m2 / (Count - 1));
///Gets the unbiased population skewness.
public double Skewness =>
    Count < 3
    ? double.NaN
    : Count * m3 * Sqrt(m2 / (Count - 1))
        / (m2 * m2 * (Count - 2)) * (Count - 1);
///Gets the unbiased population kurtosis.
public double Kurtosis =>
    Count < 4
    ? double.NaN
    : ((double)Count * Count - 1)
        / ((Count - 2) * (Count - 3))
        * (Count * m4 / (m2 * m2) - 3 + 6.0 / (Count + 1));

He omitido, por brevedad, el cálculo de propiedades como PopulationVariance, PopulationSkewness y demás. De todas maneras, están disponibles en el código de Austra, y en el propio código de Math.NET Numerics.

¿Qué le damos?

Para alimentar al acumulador, hay que pasarle las muestras, al menos en principio, de una en una, por medio de un método que hemos nombrado Add:

/// Adds a sample to this accumulator.
/// The new sample.
public void Add(double sample)
{
    ++Count;
    double d = sample - m1, s = d / Count;
    double t = d * s * (Count - 1);
    m1 += s;
    m4 += (t * s * (Count * (Count - 3 + 3)
        + 6 * s * m2 - 4 * m3) * s;
    m3 += (t * (Count - 2) - 3 * m2) * s;
    m2 += t;
    if (sample < min) min = sample;
    if (sample > max) max = sample;
}

Hay algunas pequeñas mejoras en el código anterior, respecto al original. Hay algunas multiplicaciones menos, y está todo preparado por si quisiéramos usar alguna instrucción de fusión de multiplicación y suma. No las he usado porque tengo algunas dudas sobre la eficiencia en .NET Core. Es cierto que siempre tienes la ventaja de la mayor exactitud, pero ya veremos dónde sí se usan (en pocas palabras: donde realmente importan).

Combinando acumuladores

Lo mejor de todo es que podemos combinar los valores en dos acumuladores independientes y generar un acumulador de los datos conjuntos. Esto nos permitiría, por ejemplo, dividir una muestra grande en cuatro partes, calcular cuatro acumuladores en paralelo, y luego mezclarlos en el resultado final.

public static Accumulator operator +(
    Accumulator a1, Accumulator a2)
{
    if (a1.Count == 0) return a2;
    if (a2.Count == 0) return a1;

    long n = a1.Count + a2.Count, n2 = n * n;
    double d = a2.m1 - a1.m1, d2 = d * d;
    double d3 = d2 * d, d4 = d2 * d2;
    double m1 = (a1.Count * a1.m1 + a2.Count * a2.m1) / n;
    double m2 = a1.m2 + a2.m2 + d2 * a1.Count * a2.Count / n;
    double m3 = a1.m3 + a2.m3
        + d3 * a1.Count * a2.Count * (a1.Count - a2.Count) / n2
        + 3 * d * (a1.Count * a2.m2 - a2.Count * a1.m2) / n;
    double m4 = a1.m4 + a2.m4 + d4 * a1.Count * a2.Count
            * (a1.Count * (a1.Count - a2.Count)
                + a2.Count * a2.Count) / (n2 * n)
        + 6 * d2 * (a1.Count * a1.Count * a2.m2
            + a2.Count * a2.Count * a1.m2) / n2
        + 4 * d * (a1.Count * a2.m3 - a2.Count * a1.m3) / n;
    return new() {
        Count = n,
        m1 = m1, m2 = m2, m3 = m3, m4 = m4,
        min = Min(a1.min, a2.min),
        max = Max(a1.max, a2.max),
    };
}

El código es complicado, y es fácil equivocarse copiando y pegando (ya me ha pasado). De todas maneras, es una clase que es fácil de testear.

El primer impulso es dejarlo aquí, y confiar en el paralelismo con tareas para cuando queremos acelerar el código. Mi problema con esto es que, cuando escribo código para una librería, prefiero realizar la aceleración básica con código vectorial (AVX o lo que esté disponible). ¿Por qué? Pues porque por experiencia, el programador que usa luego la biblioteca prefiere tener la opción del paralelismo por tareas para su propio código. Es cierto que las tareas se combinan más o menos bien en .NET, gracias al thread-pool que obtienes del entorno de ejecución, sin esforzarte demasiado; en Java, todo es más complicado con sus malditos executors.

Drum roll, please...

Prefiero, por lo tanto, ganar todo lo que pueda en paralelismo en una librería a golpe de instrucciones SIMD. Y esto es precisamente lo que hacemos en Accumulator con el siguiente método y algunos más que lo llaman:

public unsafe void Add(double* samples, int size)
{
    int i = 0;
    if (Avx.IsSupported && size >= 16)
    {
        var vMin = Vector256.Create(double.PositiveInfinity);
        var vMax = Vector256.Create(double.NegativeInfinity);
        var vM1 = Vector256<double>.Zero;
        var vM2 = Vector256<double>.Zero;
        var vM3 = Vector256<double>.Zero;
        var vM4 = Vector256<double>.Zero;
        var v3 = Vector256.Create(3.0);
        var v4 = Vector256.Create(4.0);
        var v6 = Vector256.Create(6.0);
        long c = 0;
        for (int top = size & CommonMatrix.AVX_MASK;
             i < top; i += 4)
        {
            c++;
            var vSample = Avx.LoadVector256(samples + i);
            vMin = Avx.Min(vMin, vSample);
            vMax = Avx.Max(vMax, vSample);
            var vd = Avx.Subtract(vSample, vM1);
            var vs = Avx.Divide(vd,
                Vector256.Create((double)c));
            var vt = Avx.Multiply(Avx.Multiply(vd, vs),
                Vector256.Create((double)(c - 1)));
            vM1 = Avx.Add(vM1, vs);
            var t1 = Avx.Multiply(Avx.Multiply(vt, vs),
                Vector256.Create((double)(c * (c - 3) + 3)));
            var t2 = Avx.Multiply(Avx.Multiply(vs, vM2), v6);
            var t3 = Avx.Multiply(v4, vM3);
            vM4 = vM4.MultiplyAdd(Avx.Subtract(
                Avx.Add(t1, t2), t3), vs);
            t1 = Avx.Multiply(vt,
                Vector256.Create((double)(c - 2)));
            t2 = Avx.Multiply(vM2, v3);
            vM3 = vM3.MultiplyAdd(Avx.Subtract(t1, t2), vs);
            vM2 = Avx.Add(vM2, vt);
        }
        var acc01 = Mix(c,
            vM1.ToScalar(), vM2.ToScalar(),
            vM3.ToScalar(), vM4.ToScalar(),
            vM1.GetElement(1), vM2.GetElement(1),
            vM3.GetElement(1), vM4.GetElement(1));
        var acc23 = Mix(c,
            vM1.GetElement(2), vM2.GetElement(2),
            vM3.GetElement(2), vM4.GetElement(2),
            vM1.GetElement(3), vM2.GetElement(3),
            vM3.GetElement(3), vM4.GetElement(3));
        var a = Mix(c + c,
            acc01.m1, acc01.m2, acc01.m3, acc01.m4,
            acc23.m1, acc23.m2, acc23.m3, acc23.m4);
        if (Count == 0)
            (Count, m1, m2, m3, m4)
                = (4 * c, a.m1, a.m2, a.m3, a.m4);
        else
        {
            long acCnt = 4 * c, n = Count + acCnt, n2 = n * n;
            double d = a.m1 - m1, d2 = d * d;
            double d3 = d2 * d, d4 = d2 * d2;

            double nm1 = (Count * m1 + acCnt * a.m1) / n;
            double nm2 = m2 + a.m2 + d2 * Count * acCnt / n;
            double nm3 = m3 + a.m3
                + d3 * Count * acCnt * (Count - acCnt) / n2
                + 3 * d * (Count * a.m2 - acCnt * m2) / n;
            m4 += a.m4 + d4 * Count * acCnt
                    * (Count * (Count - acCnt) 
                        + acCnt * acCnt) / (n2 * n)
                + 6 * d2 * (Count * Count * a.m2
                    + acCnt * acCnt * m2) / n2
                + 4 * d * (Count * a.m3 - acCnt * m3) / n;
            (m1, m2, m3, Count) = (nm1, nm2, nm3, n);
        }
        min = Min(min, vMin.Min());
        max = Max(max, vMax.Max());
    }
    for (; i < size; ++i)
        Add(samples[i]);

    static (double m1, double m2, double m3, double m4) Mix(
        long c,
        double a1, double a2, double a3, double a4,
        double b1, double b2, double b3, double b4)
    {
        long n = c + c, n2 = n * n;
        double d = b1 - a1, d2 = d * d, d4 = d2 * d2;
        return (
            (a1 + b1) / 2,
            a2 + b2 + d2 * c / 2,
            a3 + b3 + 3 * d * (b2 - a2) / 2,
            a4 + b4 + d4 * c / 8 + 3 * d2 * (b2 + a2) / 2
               + 2 * d * (b3 - a3));
    }
}

Observe que la precondición para aprovechar las instrucciones vectoriales es tener todo un array de muestras a nuestra disposición. Si nos diesen un IEnumerable<double>, tendríamos que hacer maniobras como materializar las muestras en grupos de cuatro, en un array, y alimentar así al animalito vectorial.

El código es relativamente sencillo, si miramos con atención. La parte AVX prácticamente repite el código del Add escalar. Por cada campo de Accumulator hay un vector de doble precisión. La excepción es la propiedad Count, y la tratamos diferente porque para los cuatro acumuladores virtuales que maneja el método, la cantidad de muestras es siempre la misma.

Esto es una ventaja cuando tenemos que mezclar los resultados de los cuatro acumuladores. La función interna estática Mix aprovecha la igualdad de los contadores para simplificar algebraicamente algunas fórmula. Observe, por ejemplo, que la fórmula para el m3 combinado es más sencilla, al anularse uno de los términos.

Una vez que hemos mezclado los cuatro acumuladores parciales, mezclamos el resultado, a su vez, con los valores que pueda haber ya en el propio acumulador (si los hubiera). Aquí no podemos simplificar tanto, porque los contadores nuevos y antiguos pueden ser muy diferentes, aunque en el caso en el que el acumulador inicial no tuviese muestras, es todo más simple.

Si quiere hacerse una idea de cuánto mejora este tipo de procesamiento vectorial, los benchmarks que he ejecutado me dan casi cinco veces más velocidad. Es extraño, porque yo esperaría una mejora de 4x, pero puede deberse a que aquí sí hacemos uso de las instrucciones FMA vectoriales, cuando están disponibles. Las instrucciones FMA están escondidas en los métodos de extensión MultiplyAdd que presenté en esta entrada.

Por cierto, la niña de la imagen de la entrada tiene poco que ver con el algoritmo, pero estoy usando imágenes generadas por AI, entre otros motivos, para evitar problemas de derechos de autor. En este caso, le pedí a la AI que generase una niña perdida e indefensa en un universo digital simulado. En parte, la AI me hizo caso; en parte, ignoró la petición. Pero el resultado me gusta, y ahí lo tiene.

Categorías
C#

Safe indexers

El lenguaje de AUSTRA es un sencillo lenguaje de fórmulas, inspirado mayormente en la Programación Funcional. Esto lo hace fácil de usar, y sobre todo, lo hace bastante seguro: no nos deja estropear los datos de una serie que hemos obtenido, por ejemplo, desde una fuente de pago. Al mismo tiempo, nos obliga a ser «creativos» para resolver problemas que serían más sencillos en un lenguaje tradicional.

Por ejemplo, digamos que quiero crear un vector de 1024 elementos, con la serie de números cuadrados. Eso es sencillo, si usamos el constructor (o «método de clase») apropiado. En el lenguaje de AUSTRA, se hace así:

vector::new(1024, i => (i + 1)^2)

AUSTRA soporta constructores con nombres: a diferencia de C#, en los que todos los constructores se definen con el nombre de la clase, AUSTRA nos permite distinguir entre constructores por medio de sus nombres. Por ejemplo, estas son maneras alternativas de construir una matriz:

-- Una matriz de 10x10, con una función lambda para las celdas.
matrix::new(10, 10, (fila, columna) => fila * 10 + columna)
-- Una matriz de 2x4, a la que le pasamos dos vectores.
matrix::rows([1, 2, 3, 4], [5, 6, 7, 8])
-- Una matriz de covarianza de cuatro series temporales.
matrix::cov(aapl, msft, dax, esx)

En un lenguaje como el de MATLAB, casi seguramente, tendríamos tres funciones globales para conseguir lo mismo. Con el truco de AUSTRA, evitamos identificadores globales, que terminan colisionando entre ellos y teniendo que adoptar nombres crípticos. Nos evitamos también las cábalas que hay que hacer en lenguajes como C# o Java para averiguar cuál es el constructor adecuado de acuerdo a los parámetros que estamos pasando. En el fondo, un constructor de AUSTRA termina llamando indistintamente a un constructor de C# o a un método estático. Lo importante es que el lenguaje nos abstrae de estos detalles de implementación.

Volviendo al ejemplo del vector, hemos utilizado un constructor que recibe el tamaño deseado, y una función lambda que devuelve los valores de cada celda. Digamos ahora, para complicarlo, que lo que queremos es la secuencia de Fibonacci. La solución más sencilla que se me ha ocurrido es permitir que la función lambda pueda tener un parámetro adicional que apunte al propio vector que estamos construyendo. Este ejemplo sigue siendo idéntico al anterior, pero ya estamos pasando un parámetro adicional en la función lambda, aunque no lo utilicemos de momento:

vector::new(1024, (i, v) => (i + 1)^2)

¿Es esto limpio y seguro? Por supuesto. El parámetro v nunca va a ser nulo, y de hecho, ya nos llega medio cocido: la primera vez que se llama la función lambda, todos sus elementos están a cero. La siguiente vez, estará asignado el primero elemento. Y así sucesivamente. Tenemos un contrato que nos garantiza el orden en que se van a inicializar los elementos del vector. Tenga presente que una posible implementación, para un vector suficientemente grande, podría usar paralelismo para rellenar segmentos concurrentemente. Pero cuando usamos esta variante del constructor, tenemos la garantía de que la inicialización va a ser secuencial.

Ahora veamos una primera versión de Fibonacci:

vector::new(1024, (i, v) =>
  if i <= 1 then 1 else v[i-1] + v[i-2])

Esto funciona perfectamente. El if de marras no es una instrucción: es una expresión condicional ternaria, como la interrogación y los dos puntos en C/C# y familia. En el caso más habitual, sumamos los dos elementos anteriores, que ya estarán inicializados. Si nos diese por hacer referencia a un elemento posterior, no sería un problema, excepto que su valor sería cero. Y la expresión condicional nos ahorra una excepción de acceso fuera de rango.

AUSTRA tiene un mecanismo más potente para estos casos, sin embargo:

vector::new(1024, (i, v) =>
  if i = 0 then 1 else v{i-1} + v{i-2})

Esta vez, tenemos un solo caso especial: el del primer elemento del vector. El cambio está en la forma en que accedemos a los elementos de v: con llaves, en vez de corchetes. Esto es lo que he llamado un safe indexer, o indexador seguro. Si el valor del índice está en rango, todo procede como de costumbre. En caso contrario, la expresión devuelve un cero. Es como si tuviésemos una memoria infinita, en la que una región de la misma puede tener valores distintos de cero, pero fuera de esa región, todo es cero. Como se trata de un lenguaje funcional, además, no existe la posibilidad de escribir en la región "externa" de la memoria del vector: no podemos hacer asignaciones ni a v[0] ni a v{v.length}, por ejemplo, aunque la primera expresión se refiera a un elemento que realmente existe.

Otras aplicaciones

En Econometría, se conoce como serie temporal autoregresiva de orden p a la que se construye de acuerdo a esta fórmula:
$$X_t = \sum_{i=1}^{p}{\phi_i X_{t - i}} + \epsilon_t$$
El símbolo $\epsilon_t$ se refiere a variables aleatorias con una distribución normal estándar. Los coeficientes $\phi_i$ son números reales que definen el comportamiento de la serie. Se trata, por supuesto, de una serie que combina ruido blanco con una retroalimentación limitado de valores anteriores. Por ejemplo, en high-frequency trading, los precios de ejecución de un activo suelen poderse representar con una serie autoregresiva, en muchos casos.

Este es el caso de uso que me ha obligado a implementar estos indexadores seguros. La forma de construir una serie autoregresiva en AUSTRA es la siguiente:

let r=vector::nrandom(1024) in
    vector::new(r.length, (i, v) =>
        r[i] + 0.7*v{i-1} + 0.1*v{i-2})

Primero construimos un vector aleatorio usando el constructor vector::nrandom. La cláusula let nos permite crear una variable local a la fórmula, a la que podemos hacer referencia en lo que queda de fórmula. El resto es fácil de adivinar: no hace falta un indexador seguro para acceder a r, pero sí a los elementos anteriores del vector que estamos construyendo.

Implementación en C#

La implementación del indexador seguro en C# puede resultar interesante, por las técnicas de "bajo nivel" que utiliza. Esta es la función de la clase Vector que implementa dicha funcionalidad:

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public double SafeThis(int index) =>
    (uint)index >= values.Length
    ? 0.0
    : Unsafe.Add(
        ref MemoryMarshal.GetArrayDataReference(values),
        index);

Vector es, en AUSTRA, una estructura muy ligera que se limita a encapsular un campo values de tipo double[]. La implementación podría ser una consulta de rango, seguida de un acceso "normal" a los datos del vector:

public double SafeThis(int index) =>
    index >= 0 && index < values.Length ? values[index] : 0.0;

Pero AUSTRA es una librería que se va a utilizar en cosas que a priori no te puedes imaginar, y en estos casos, optimizar bien el código no puede calificarse de "optimización prematura".

En este caso, he usado varios trucos, que son bastante usados por el propio código de .NET Core:

  • Primero, está el truco de convertir el índice en un entero sin signo, para ahorrarnos una comparación y un salto. Si nos pasan un índice negativo, el valor reconsiderado como entero sin signo va a ser inevitablemente mayor que la longitud del array.
  • Luego está el truco indescriptiblemente sucio de usar MemoryMarshal.GetArrayDataReference para obtener la dirección de memoria donde comienza el array. Este truco sólo puede fallar si value fuese un puntero nulo, pero no es el caso.
  • Finalmente, el método Unsafe.Add suma el índice a esa posición inicial y devuelve el valor del elemento. Este es un truco más tolerable.

Al final, el compilador y el JIT generan más o menos el mismo código que para un acceso normal "verificado", pero con la particularidad de devolver cero si el índice está fuera de rango, que es lo que queríamos. No se añade código innecesario.