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.

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#

Reuniendo piezas de un vector

Cuando hay que utilizar instrucciones vectoriales, es siempre aconsejable que se cargue en memoria un trozo consecutivo de memoria para crear vectores. Para un programador ya entrenado en AVX/AVX2/AVX512, por ejemplo, es fácil darse cuenta de que el siguiente bucle (un producto escalar de dos vectores) es un candidato ideal para usar instrucciones SIMD:

double sum = 0.0;
for (int i = 0; i < size; i++)
    sum += p[i] * q[i];

Hay dos zonas de memoria, p y q, y el índice se mueve ascendentemente de uno en uno, en ambos casos. Por lo tanto, para convertir este código en instrucciones vectoriales de AVX, hacemos que cada paso del código cargue cuatro entradas de cada vector, en cada paso del bucle:

double sum = 0;
int i = 0;
if (Fma.IsSupported)
{
    var acc = Vector256<double>.Zero;
    for (; i < size - 4; i += 4)
        acc = Fma.MultiplyAdd(
            Avx.LoadVector256(p + i),
            Avx.LoadVector256(q + i),
            acc);
    acc = Avx.HorizontalAdd(acc, acc);
    sum = acc.ToScalar() + acc.GetElement(2);
}
for (; i < size; i++)
    sum += p[i] * q[i];

Gracias a este sencillo truco, se aprovecha mejor la caché de memoria del procesador. Es cierto que p y q pueden estar en direcciones distantes, pero al menos la carga de cada trozo de cuatro valores no provoca invalidaciones de la caché. En cambio, este otro fragmento de código, extraído del código de Austra que realiza la descomposición de una matriz según la factorización LU, no cumple con esta regla:

double m = c[k];
for (int i = k + 1; i < size; i++)
    c[i] -= m * a[i * size + k];

El almacenamiento en c es consecutivo, pero si reunimos cuatro pasos del bucle, las cargas desde el vector a van a venir de zonas dispersas de la memoria: esto es porque la variable de bucle en el indexador está multiplicada por size. Sería una pena no poder utilizar AVX aquí, porque podríamos reducir por cuatro las restas y multiplicaciones del bucle.

Por fortuna, Intel agregó al conjunto de instrucciones de AVX2 una instrucción que, precisamente, sirve para reunir en una misma carga fragmentos de memoria no consecutivos. Se trata de la instrucción vgatherdpd, que .NET implementa mediante el método Avx2.GatherVector:

public static Vector256<double> GatherVector256(
    double* baseAddress,
    Vector128<int> index,
    byte scale);

A primera vista, parece arameo antiguo. A segunda y tercera vistas, se asemeja más al sumerio de chabolas. Pero no es tan complicado explicar cómo funciona, sobre todo si lo acompañamos con un ejemplo. A este método le pasamos un puntero a un array de números de doble precisión, esa es la parte obvia. Además, le pasamos cuatro índices de tipo entero, en el parámetro de tipo Vector128<int>. Como este tipo contiene 128 bits, es evidente que ahí podemos pasar cuatro enteros de 32 bits; ya veremos cómo. ¿Y el factor de escala? Aquí es un poco redundante, y Microsoft podría habernos ahorrado el trabajo: cuando se trabaja con vectores de double, tenemos que pasar un 8, o lo que es igual, sizeof(double).

Juntando las piezas, he aquí como quedaría el código transformado del fragmento original:

double m = c[k];
int i = k + 1;
if (Avx2.IsSupported)
{
    var vm = Vector256.Create(m);
    var vx = Vector128.Create(0, size, 2 * size, 3 * size);
    for (double* p = a + i * size + k;
        i < size - 4;
        i += 4, p += size * 4)
        Avx.Store(c + i,
            Avx.Subtract(Avx.LoadVector256(c + i),
            Avx.Multiply(Avx2.GatherVector256(p, vx, 8), vm)));
}
for (; i < size; i++)
    c[i] -= m * a[i * size + k];

Expliquémonos: antes de entrar en el bucle, creamos un vector de números reales a partir del valor de la variable m. Lo que sigue es lo que nos interesa: creamos por adelantado, un Vector128<int> con cuatro offsets o desplazamientos diferentes, múltiplos del tamaño de la fila de la matriz, que viene en size. Y ahora viene el detalle artesanal delicado: al iniciar el bucle, creo una variable adicional de puntero, p, que apunte i * size + k entradas a partir de la dirección inicial, que está en a. Recuerde que leíamos el valor de a[i * size + k] en cada paso del bucle original, y que i es la variable de bucle. Si hubiese querido ser más claro, debería haber escrito:

double* p = a + (k + 1) * size + k

Pero me puede la pereza, y al fin y al cabo, en ese preciso instante, i es equivalente a k + 1, ¿no?

El detalle es que, en cada paso del bucle AVX, la variable p se incrementa en 4 * size entradas (es decir, en cuatro filas enteras de size entradas). Digo que es un detallazo elegante porque la alternativa obvia sería incrementar directamente los cuatro índices enteros que guardamos en el vector vx. ¿Por qué no lo he hecho así? Pues porque, al trabajar con vectores de 128 bits, la instrucción de suma pertenecería al conjunto de instrucciones SSE2… y resulta que el cambio entre AVX y SSE tiene un coste pequeño, pero no desdeñable.

… y este es el momento más apropiado para exclamar ¡pero coño! (con perdón), ¿esto no tendría que hacerlo el compilador, o el JIT, o quien narices sea?. Pues sí. De hecho, teóricamente, Hotspot, el JIT de Java, tendría que ocuparse de estas cosas. En la práctica, la vectorización automática en Hotspot es un truño como un puño. Y el famoso escape analysis tampoco hace lo que se supone que tendría que hacer. O si lo hace, no se nota. Al final, Java ha añadido una librería de vectores basados en instrucciones SIMD, al estilo de .NET. Los compiladores de C/C++ hacen mejor papel en estos casos, pero no son omnipotentes, y a veces tienes que echarles una mano.

NOTA: AVX2, creo recordar, se introdujo en la cuarta generación de Intel, Haswell, en 2013. La primera implementación en hardware de vgatherdpd no era para tirar cohetes, y todavía hay páginas y comentarios en Internet que lo recuerdan. La implementación mejoró notablemente en Skylake, la sexta generación (2015). Yo tuve un Haswell durante muchos años, y empecé a probar AVX2 con esa máquina. Mi penúltimo portátil fue un Skylake. Ahora mismo tengo un PC con CPU de oncena generación, Tiger Lake, que es un pepino. Lo menciono porque, por casualidad, Tiger Lake es la última generación de procesadores de «consumo» que incluye las instrucciones AVX512: con la llegada de las CPU con núcleos heterogéneos, estas instrucciones se deshabilitan porque los núcleos «eficientes» no las implementan. Ahora mismo, si quieres AVX512, tienes que irte a Xeon, o a AMD.

Categorías
FinTech

Austra

Otro pequeño interludio: si os interesan estas cosas, echad un vistazo a la ayuda online de Austra. Es un lenguaje de expresiones orientado al manejo de series financieras y modelos econométricos. Está todavía en desarrollo, pero ya va tomando forma. Entre otras cosas interesantes, permite usar funciones lambda, y la librería, Austra.Library, utiliza AVX2 para optimizar las operaciones sobre matrices, vectores y series temporales. Hasta 200/300 dimensiones, es generalmente más eficiente que las alternativas, incluso las que usan Intel MKL. A partir de ese punto, naturalmente, Intel MKL es más eficiente. En cualquier caso, tengo la posibilidad de mejorar la implementación para ese volumen de información. Y como última alternativa, puedo incorporar Intel MKL como una opción más; opcional, por supuesto, porque esta librería tiene licencias de pago.

Mi idea es, en el momento en que esté suficientemente maduro, convertir parte del código en open source. Aparte de añadir más algoritmos a la librería básica, quiero dar alternativas al analizador del lenguaje. Ahora mismo, el lenguaje genera árboles de expresiones, que a su vez generan código nativo sobre la marcha. Para una línea de comandos, sin embargo, puede que esto sea excesivo. No voy a destruir el compilador actual, pero quiero subir, antes de abrir el código, un intérprete más modesto que, al mismo tiempo, seguramente será más eficiente para el propósito de la herramienta.