Etiqueta: .NET

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
C# Insights

Smart boy’s optimizations

Decía el gran Donald Knuth algo así como que premature optimization is the root of all evil. Santificado sea su nombre…

Calidad de código

Una vez citadas las sagradas escrituras, debo reconocer que mi lado hereje cumple otros mandamientos:

  • Peor que la optimización prematura, es no optimizar nunca. Una vez escuché una mala excusa sobre un programa que tardaba 25 horas en cargar un fichero: «es que nadie me dijo que tenía que ser rápido». Los ordenadores existen, precisamente, para hacer las cosas más rápidamente. No sólo porque es interés directo del usuario del programa que éste termine antes, sino que, además, le interesa el ahorro en electricidad y en el desgaste del propio aparato.
  • Si estás escribiendo una aplicación, te puedes permitir el lujo de esperar a que funcione para buscar los puntos críticos de eficiencia. Pero si estás escribiendo una librería, que se va a utilizar en formas que aún no sospechas… mejor que todo vaya como la seda desde el principio.
  • La mayoría de las optimizaciones (yo diría más bien mejoras) caen en una categoría que yo llamo «mejoras de chico listo», y tienen que ver con la calidad del código que cada programador puede generar sin esfuerzo adicional.

Las optimizaciones de chico listo, por supuesto, dependen de la experiencia del programador, de lo bien que le funcione la memoria y de lo bien que se le dé la detección de patrones, por lo que se trata de una categoría difícil de delimitar. Programar es un arte.

Un ejemplo

En cualquier caso, el propósito de esta entrada es mostrarle algunas de las optimizaciones que he aprendido mirando el código fuente de .NET. Yo las tengo ya en mi memoria de trabajo: las aplico automáticamente cuando detecto que son aplicables.

Trabajando con una implementación de la función erf, tropecé con este código, que evalúa un polinomio en un punto, usando los coeficientes de una tabla:

private static double Evaluate(double z, double[] coefficients)
{
    if (coefficients == null)
        throw new ArgumentNullException(nameof(coefficients));
    int n = coefficients.Length;
    if (n == 0)
        return 0;
    double sum = coefficients[n - 1];
    for (int i = n - 2; i >= 0; --i)
    {
        sum *= z;
        sum += coefficients[i];
    }
    return sum;
}

Esta función se ejecuta varias veces, con distintos coeficientes. Un ejemplo de tabla de coeficientes es ésta:

static readonly double[] ErvInvImpAn =
{
    -0.000508781949658280665617, -0.00836874819741736770379,
    0.0334806625409744615033, -0.0126926147662974029034,
    -0.0365637971411762664006, 0.0219878681111168899165,
    0.00822687874676915743155, -0.00538772965071242932965
};

Este método es un método privado de una clase, y una rápida ojeada me confirmó que las tablas que se le pasan son siempre no nulas, y con longitud mayor que cero. ¿A qué vienen las dos comprobaciones iniciales? Respuesta: es uno de los problemas que causa la «modularidad». Escribes software que no sabes cómo se puede usar, y lo proteges de las cosas más inverosímiles. Pero si es un método privado, tanta precaución sobra. Empezamos por esta simplificación, para ir haciendo boca y verlo todo más claro:

private static double Evaluate(double z, double[] coefficients)
{
    int n = coefficients.Length;
    double sum = coefficients[n - 1];
    for (int i = n - 2; i >= 0; --i)
    {
        sum *= z;
        sum += coefficients[i];
    }
    return sum;
}

El siguiente paso seguramente le sorprenderá: sustituyo la tabla de coeficientes, que ahora es un campo estático de sólo lectura, por esto:

static ReadOnlySpan<double> ErvInvImpAn => new[]
{
    -0.000508781949658280665617, -0.00836874819741736770379,
    0.0334806625409744615033, -0.0126926147662974029034,
    -0.0365637971411762664006, 0.0219878681111168899165,
    0.00822687874676915743155, -0.00538772965071242932965
};

Sorprendente, ¿verdad? Es un truco poco conocido, pero que Microsoft usa a diestra y siniestra en el código de .NET Core. Por razones que en parte se me escapan, el compilador de C# y el JIT transforman esta construcción en una zona de datos dentro de los metadatos del código IL. Y el JIT lo maneja más eficientemente. No hay mucha lógica en que tengamos que usar precisamente un ReadOnlySpan<double>, o que haya que convertir el campo en una propiedad de sólo lectura. Se trata de una marca, o un guiño de complicidad, que utilizan el JIT y el compilador para generar código más eficiente.

Esto me obliga a crear una nueva versión del amigo Evaluate que acepte un ReadOnlySpan<double> como origen de sus coeficientes. Esta es la nueva versión, con dos optimizaciones adicionales:

private static double Evaluate(
    double z, ReadOnlySpan<double> coeffs)
{
    int n = coeffs.Length;
    ref double rd = ref MemoryMarshal.GetReference(coeffs);
    double sum = Unsafe.Add(ref rd, n - 1);
    for (int i = n - 2; i >= 0; --i)
        sum = Math.FusedMultiplyAdd(
            z, sum, Unsafe.Add(ref rd, i));
    return sum;
}

De las dos nuevas mejoras, la más sencilla es el uso de Math.FusedMultiplyAdd: un método de la clase Math que combina la multiplicación y la suma en una sola instrucción de la CPU, y puede darnos más velocidad y precisión. En este caso, además, he medido que realmente sea ventajosa, porque no siempre lo es.

El segundo cambio tiene dos partes. Como el bucle for utilizado no es un bucle convencional, el JIT actual no puede deducir que no habrán referencias fuera de rango para eliminar las comprobaciones de los índices en tiempo de ejecución. El bucle es descendente, y ni siquiera comienza por el último elemento. No le podemos exigir tanto al JIT.

Lo primero que hago es pedir una managed reference a la primera celda de la tabla de coeficientes:

    ref double rd = ref MemoryMarshal.GetReference(coeffs);
    // Equivalente a:
    // ref double rd = ref coeffs[0];

Esto es más o menos parecido a pedir un puntero al inicio de la tabla. En realidad, C# nos permitiría pedir un puntero al inicio de la tabla, pero el precio sería «fijar» la tabla en memoria para que el recolector de basura no vaya a pensar que no la estamos usando. El puntero que conseguimos con esta técnica es uno que el recolector de basura puede identificar y tener en cuenta. Y la forma normal de pedirlo es la que muestro en los comentarios del fragmento. ¿Por qué no la he usado? Pues porque implicaría una comprobación de rango innecesaria: el JIT generaría una comparación y un salto para verificar si la tabla no está vacía. Para evitarlo, uso MemoryMarshal.GetReference, que es otro truco sucio de Microsoft, para conseguir un puntero al inicio de un array sin costes ocultos.

Lo que sigue es más sencillo: utilizo el método Add de la clase Unsafe para llegar a cada una de las celdas que contienen los coeficientes. Sí, todo es un poco enrevesado, pero una vez que te lo aprendes, no te cuesta nada escribir estas cosas de carrerilla. Me siento en el deber de contárselas. Ya usted decidirá si merece la pena o no usarlas en su propio código cuando lo crea necesario. No son cosas para usar en una aplicación que tienes que escribir en tres meses. Pero creo que tienen un lugar en una librería de código.

Y hay más, claro

Hay montones de trucos similares en el código fuente de .NET. Por ejemplo, imagine que hay que tiene que hacer una comprobación de rango de un índice:

if (0 <= index && index < length) ...

Dos comparaciones, y dos saltos. Las comparaciones son lo de menos. Los dos saltos ralentizan todo. ¿Qué hace Microsoft en estos casos?

if ((uint)index < length) ...

La variable index suele ser un entero con signo. No cuesta nada pedir que el compilador la trate, momentáneamente, como un entero del mismo tamaño, pero sin signo. Si el índice fuese negativo, al tratarlo como un entero sin signo, el valor sería inevitablemente superior al de length. Una sola comparación, y un único salto potencial.

Veamos una variante derivada de este truco. El analizador lexical de Austra tiene que comprobar muchas veces si un carácter es un dígito decimal:

if ('0' <= ch && ch <= '9') ...

La forma más eficiente, sin embargo, es la siguiente:

if ((uint)(ch - '0') < 10u) ...

He introducido una resta, que se ejecuta eficientemente, y he quitado un salto potencial.

De todas maneras, una de mis optimizaciones de chico listo preferidas es muy sencilla. En vez de escribir:

x * x - y * y

un servidor prefiere:

(x + y) * (x - y)

Y es que en la segunda expresión hay una suma de más, pero una multiplicación de menos.

Es agradable tener un cerebro cargado, y estar dispuestos a usarlo.

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#

La la, how their life goes on

Si insisto tanto en estos temas de instrucciones vectoriales, es parar preparar a mis lectores potenciales para que usen estas técnicas cuando lo crean necesario. Naturalmente, las primeras explicaciones van a ser del tipo heroico: todo hay que hacerlo a mano, razonando desde los primeros principios. Pero, como decía Molly Jones, la mujer de Desmond Jones, la la, how their life goes on: la vida sigue su curso, y lo normal en la vida de un programador es usar las técnicas que todos conocemos y apreciamos para ahorrarnos esfuerzo.

Lo primero para que todos nos ahorremos trabajo, por supuesto, consiste en usar estas cosas a través de una librería bien probada. Ya existen unas cuantas de estas, pero estoy creando Austra porque hay un nicho muy concreto, y porque la librería tiene méritos propios. La idea es hacer de Austra un proyecto open source, por supuesto. En este momento está en GitHub, pero no es por ahora un repositorio público porque la aplicación de «pruebas» está basada en WPF y DevExpress. O me compro personalmente una licencia comercial, o cambio los controles por algo gratuito. No es sencillo. Hay cosas como Avalonia o Uno Platform, que además, son multiplataformas, pero son bastante pobres en componentes. Y no quiero ni oír hablar de .NET MAUI: es un fracaso, de momento. Antes que usar .NET MAUI, prefiero «degradar» la aplicación a Windows Forms (tengo un editor de código bastante bueno) y currarme los gráficos y las cosas que falten. Pero a eso vamos: a que Austra esté disponible gratuitamente como open source, que cualquier que quiera pueda colaborar, que haya un conjunto de tests y benchmarks exhaustivo, etc, etc.

Pero no se escribe un post para explicar lo obvio. Mi verdadero objetivo ahora es explicar un par de técnicas muy tontas que bajan el listón del heroísmo al escribir este tipo de código, y que todos conocemos, pero que nos puede dar miedo usar a primeras, porque no sabes cómo va a interferir en la generación de código de .NET 7 y .NET 8 (que tiene sus cagadas, como todo).

Primer ejemplo: resulta que mucho de estos algoritmos algebraicos calculan el producto escalar de dos zonas de memoria en muchos casos. Hay una clase Vector con su correspondiente producto escalar, e incluso un ComplexVector, pero es mejor tener una rutina que maneje directamente zonas de memoria con un tamaño predeterminado. Austra tiene una clase estática CommonMatrix que precisamente implementa estas cosas. Ésta es la última versión del método en cuestión:

/// <summary>
/// Computes the dot product of two double arrays.
/// </summary>
/// <param name="p">Pointer to the first array.</param>
/// <param name="q">Pointer to the second array.</param>
/// <param name="size">Number of items in each array.</param>
/// <returns>A sum of products.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public unsafe static double DotProduct(
    double* p, double* q, int size)
{
    double sum;
    int i = 0;
    if (Avx.IsSupported)
    {
        Vector256<double> acc = Vector256<double>.Zero;
        for (int top = size & AVX_MASK; i < top; i += 4)
            acc = acc.MultiplyAdd(p + i, q + i);
        sum = acc.Sum();
    }
    else
        sum = 0;
    for (; i < size; i++)
        sum += p[i] * q[i];
    return sum;
}

A primera vista, lo único raro aquí es el atributo que fuerza un inlining agresivo. Pero hay un par de cosillas más que no son métodos habituales de AVX. Los muestro aparte aquí, porque son métodos de extensión de la misma clase CommonMatrix:

/// <summary>Sums all the elements in a vector.</summary>
/// <param name="v">
/// A intrinsics vector with four doubles.
/// </param>
/// <returns>The total value.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static double Sum(this Vector256<double> v)
{ 
    v = Avx.HorizontalAdd(v, v);
    return v.ToScalar() + v.GetElement(2);
}

Resulta que sumar los cuatro elementos de un vector de cuatro dobles no es moco de pavo. Hay teorías que pululan por la Internet sobre cuál es la forma eficiente de lograrlo. Yo no estoy seguro aún de que la mía sea la más eficiente, pero gracias al encapsulamiento en un método separado, si descubro algo mejor, tendré que tocar el código en un único lugar. Obvio, ¿no? Pero hasta comprobar que el JIT de .NET Core no se volvía loco, no me atreví a dar este paso. Ya he comprobado que no pasan cosas raras.

Es curioso, sin embargo, que el mínimo y el máximo de un vector se calcule más eficientemente con un método como el siguiente:

[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static double Max(this Vector256<double> v)
{
    var x = Sse2.Max(v.GetLower(), v.GetUpper());
    return Math.Max(x.ToScalar(), x.GetElement(1));
}

En este caso, no parece que la transición temporal a una instrucción SSE haga mucho daño al código que la rodea. De todas maneras, observe que para la reducción final al escalar hay que ir a por el segundo elemento directamente. Cosas de Intel.

Éste es otro método que usa el código original, omitiendo esta vez los comentarios XML:

internal unsafe static Vector256<double> MultiplyAdd(
    this Vector256<double> summand,
    double* multiplicand,
    double* multiplier) =>
    Fma.IsSupported
        ? Fma.MultiplyAdd(
            Avx.LoadVector256(multiplicand),
            Avx.LoadVector256(multiplier),
            summand)
        : Avx.Add(summand, Avx.Multiply(
            Avx.LoadVector256(multiplicand),
            Avx.LoadVector256(multiplier)));

Hay un par de variantes más que usan vectores directamente en vez de direcciones, y variantes adicionales para MultiplySubtract y MultiplyAddNegated, que sons métodos FMA. Estas aparentes tonterías ahorran líneas y líneas de código. En mi caso, como tengo memoria fotográfica, prefiero que el código fuente sea lo más pequeño posible para poder recordarlo más adelante.

Categorías
C#

Blazor

A flaming Blazor

Empecé mis experimentos con Blazor en los lejanos tiempos de .NET 5. Hace un par de años, quiero decir.

Con un poco de contexto me explicaré mejor: prefiero una buena «aplicación de escritorio» antes que cualquier página web. Puede ser porque el tipo de proyectos en los que he invertido más tiempo han sido siempre aplicaciones con alta densidad de información, como las llaman ahora. Pero también porque he podido comparar la productividad de equipos trabajando para la web y para aplicaciones de escritorio de toda la vida. La pérdida de productividad, según mi experiencia personal, se agrava cuando es Angular la herramienta de front-end. Conozco TypeScript en profundidad, y he hecho proyectos con Angular y con Svelte. Svelte es muchísimo más productivo, e infinitamente menos pesado luego en producción. Pero la industria tiene cierta fijación con el puñetero Angular.

Las primeras pruebas con Blazor las hice con la idea de usarlo para crear prototipos rápidos. Funcionó de maravillas. Pero no estaba muy puesto en Bootstrap por entonces, y el diseño visual de los ejemplos de Microsoft era, y sigue siendo, abominable. Me preocupaba, además, la falta de controles de terceros. El cliente, por ejemplo, estaba empeñado en usar un gráfico de radar en la página principal del proyecto… y los componentes que manejábamos no traían el dichoso gráfico. Al final, el front-end terminó haciéndose en React, pero el prototipo de Blazor sigue existiendo y usándose internamente, porque en varios aspectos, es más rápido y potente.

Una de las lecciones más importantes de mi carrera como desarrollador, es que si te toca crear una librería, un servidor o algo que normalmente no tenga una interfaz de usuario… es mejor que te crees una interfaz propia, aunque no te la paguen ni te la agradezcan. Es tu seguro de vida para que no te culpen si el front-end o la aplicación hecha por otros es lenta. La vida es dura.

Ahora mismo, acabo de terminar un proyecto grande con Blazor y ASP.NET. Bueno, realmente el prototipo, porque los «sabios» de la «industria» siguen obcecados con Angular y Spring Boot, y ahora habrá que rehacerlo todo. Blazor es ya una herramienta de desarrollo madura. Yo he aprendido un poco más de CSS y HTML (desastres ambos, en mi opinión) y la parte de la fealdad ya está más o menos superada. Los componentes de terceros han evolucionado también. Y hay mejores libros y blogs sobre cómo funciona Blazor, que no siempre es lo que aparenta. La velocidad de desarrollo sigue siendo incomparable. He terminado la funcionalidad necesaria del proyecto en un mes. Calculamos que rehacerlo todo en Angular nos costará tres veces más tiempo. La estimación de tiempo no es mía: mi papel ha sido rebajarla, porque al fin y al cabo, ya sabemos cómo hacer todo.

Me siento un poco como Sísifo. Pero uno termina por acostumbrarse a estas cosas.