Autor: Ian Marteens

Categorías
Insights

La mentalidad predominante

Intentaré ser breve.

Hace ya más de diez años, a una persona muy inteligente y a mí nos pusieron a trabajar a destajo en dos proyectos paralelos en Java cuyo núcleo técnico era bastante parecido: había que enviar precios en una red local a toda leche. Por entonces, Java iba sólo por la versión 7, y había una estructura de datos en Java que se había puesto de moda: el «patrón Disruptor». No pongo enlaces porque no se lo merece. La idea del famoso disruptor era sustituir la típica cola bloqueante (o el Channel que ahora tiene .NET Core) por una cola circular, sin bloqueos… y con un hilo dedicado, girando como un demonio constantemente, calentando la CPU y robando recursos. Teóricamente, este desperdicio se traduciría en menos latencia entre procesos… lo cual incluso puede ser cierto.

A esta idea del bucle en perpetuo movimiento se le añadía un poco de polvo de hadas: la cola circular estaba pensada para evitar pedir memoria dentro de lo posible (un problema que provoca el propio Java), y estaba todo calculado para evitar cosas con nombres tan feos como el false sharing (que es un problema real cuando hay concurrencia). Es decir, por estar hecho en Java, el disruptor era ya la leche… aunque terminó teniendo versiones en C# y supongo que hasta en Python. Toda la mierda que los gatos tiran al suelo termina cayendo en una librería de Python y ahí se hace eterna.

La cruda realidad

Como soy incrédulo, sobre todo de las cosas que no me gustan, lo que hicimos (la persona inteligente y yo) fue crear una demo. Probamos el disruptor contra la cola bloqueante de serie de Java, y ganaba el disruptor, aunque no por KO. El ganador absoluto, sin embargo, fue una variante de la cola bloqueante, con el añadido de recuperar e insertar los elementos en lotes. Es una idea estúpidamente simple y efectiva. Si miras la cola y ves que tienes cuatro elementos pendientes de recuperar, aprovecha y tráetelos todos. En Java tienes que ser cuidadoso con estas cosas, y no crear un array cada vez que te traigas n objetos, pero esto tiene una solución tan sencilla como adosar un array a la cola bloqueante en lotes, y reaprovecharlo todo el tiempo.

¿Cuál era el problema del disruptor? Pues que la idea de quemar un hilo empieza pareciendo atractiva, hasta que te das cuenta de que, si esta filosofía la sigues aplicando al resto de la aplicación, te quedas sin CPU en un plisplás. Creo que la única aplicación del disruptor que llegó a usarse en mi empresa de entonces fue en un proceso de logging sin bloqueos. Me dio un ataque de risa cuando me lo contaron, y todavía creo que quien me lo contó bromeaba.

¡Más Java!

Hay un meme rondando Internet y sacado de una peli barata, sobre un productor musical que todo lo resolvía añadiendo more cowbell, es decir, más cencerro, a la pista de percusión de cualquier canción. En el mundo de la programación, el equivalente es añadir más Java a Java. Me explico:

La otra gran decisión a tomar en el par de proyectos paralelos de hace más de diez años era qué íbamos a utilizar para leer y escribir en sockets. Como uno es un idiota iletrado que conoce poco Java, mi razonamiento es que, a no ser que montase una librería de muy bajo nivel en código nativo, cualquier librería basada en sockets que no implementase algún protocolo experimental no iba a ser mejor que la que ya venía en Java. Todas las librerías de sockets de terceros, y había para escoger, lo único que aportaban eran «abstracciones» como superestructura, que presuntamente podían simplificar la programación, pero jamás iban a conseguir que todo fuese más rápido. Por eso, yo me quedé con los sockets de toda la vida, y mi amigo, que es un experto en Java, se decidió por una librería de la famosa fundación Apache.

Y llegó el día del estreno. Mi aplicación llegaba a leer hasta 250.000 mensajes por cada hilo o canal que habilitase en el proceso. La otra aplicación reventó por un límite absurdo que imponía la librería de terceros. El problema se solucionó, por supuesto, tras algo de trabajo adicional. Pero yo, que soy un idiota, me quedé con una moraleja: añadir más Java a Java no suele aportar nada de valor.

Tiempos modernos

Como sospechará el lector, nadie saca una historia de hace diez años si no ha pasado algo parecido recientemente. En efecto. Una librería escrita en C# resuelve un problema de negocios muy interesante. Pero es una librería a secas, y el equipo que la mantiene no da abasto para conectarla a todas las diferentes fuentes de datos posibles. El tiempo de respuesta, además, es sumamente importante, en este caso.

Entonces apareció un Fervoroso Creyente de la Religión Verdadera, que es Java (y Martin Fowler es uno de sus profetas). La solución propuesta fue añadir Java, y polvos de unicornio. Y, como Java iba a estar obligatoriamente en otro proceso, traer un equivalente moderno del disruptor: una librería de comunicaciones que «resuelve» el problema de la «latencia» poniendo un proceso intermedio que gira como un púlsar al que se le ha ido la olla, y un par de hilos en cada extremo haciendo lo mismo. Es decir, la idea del hilo eterno del disruptor multiplicada por tres. Lo probaron en una máquina razonable, con sólo esos tres componentes: un publicador, el intermediario, y un subscriptor. El invento añadió sólo una decena de microsegundos, o eso me contaron, a la transmisión. Eso ya me habría hecho dudar, si lo hubiese sabido a tiempo. El problema es que la aplicación necesita transmitir datos de unos ocho canales. Tres por ocho, aquí y en Javalandia, siempre ha sido igual a veinticuatro. Estoy escribiendo en la máquina a la que llamo El Pepino, y que es un Core i9 con 8 núcleos físicos y 16 núcleos virtuales por hyperthreading. Es decir, menor estrictamente que veinticuatro.

Creo que todavía están apretado y soltando tuercas para intentar que la cosa funcione. Best of the luck, my friends.

Prometí ser breve, pero no lo conseguí.

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
Austra

List comprehensions in Austra

Hay quien traduce el término inglés list comprehension literalmente como «comprensión de listas» o, aún peor, «listas de comprensión». Lo interesante es que David Turner, el autor del lenguaje funcional Miranda, llamó inicialmente a estas expresiones Zermelo-Frankel expressions, pero alguien lo convenció para llamarlas list comprehensions, que en inglés no suena tan mal. Yo, porque soy un tipo caprichoso, cuando traduzca el término en el contexto de Austra, las llamaré constructores de secuencias, a secas. Menos pretencioso, y más comprensible, creo.

Un truco para escribir menos

¿Que es un constructor de secuencias? Pues es un truco sencillo para eso mismo: construir secuencias, pero usando menos código, con el añadido de que el resultado es normalmente más sencillo de leer (no siempre). Imagine que tenemos un vector de números reales, y queremos quedarnos con los que son enteros divisibles por dos, para elevarlos al cuadrado. En Austra, hasta ahora, haríamos esto, suponiendo que tenemos el vector ya almacenado en una variable global v:

v.filter(x => x % 2 = 0).map(x => x^2)

Con el mecanismo nuevo de construcción de secuencias, la expresión anterior se reduciría a esto:

[x in v : x % 2 = 0 => x^2]

No hace falta contar caracteres para ver que realmente hemos escrito menos. Nos hemos ahorrado los paréntesis de los métodos, y los parámetros de las funciones lambdas se han reducido a uno solo, reflejando el hecho de que los valores que fluyen por el constructor son casi siempre (casi) del mismo tipo. El resultado del constructor, en este caso, es un vector, y puedo seguir aplicando métodos y operadores tras el corchete de cierre. Por ejemplo, puedo hacer algo algo estúpido como añadir uno a cada valor (lo podía haber hecho al elevar al cuadrado):

([x in v : x % 2 = 0 => x^2] + 1).plot

Pero también podía haber transformado el vector resultante con una matriz, o cualquier otra cosa posible con un vector.

También podía haber creado una secuencia en el constructor, aunque los datos viniesen de un vector, o de una matriz:

[x in seq(v) : x % 2 = 0 => x^2];
[x in seq(v1^v2) : x % 2 = 0 => x^2];
[x in iseq(1, 100) : x % 2 = 0 => x^2];
[x in 1..100 : x % 2 = 0 => x^2];

El primer ejemplo usa una secuencia basada en un vector. El segundo, una secuencia basada en una matriz que se genera a partir de dos vectores. En el tercero, simplemente uso una secuencia de enteros construida a partir de un range. Y en el último ejemplo uso más «syntatic sugar» para generar la secuencia directamente a partir de un rango.

Ojo con las series

Con las series, hay que tener un poco de cuidado, porque el método que filtra una serie recibe como parámetro una lambda del tipo Func<Point<Date> bool>, mientras que el método Map usa una de tipo Func<double, double>. El constructor de secuencias lo tiene en cuenta, pero tenemos que recordarlo:

let mean = MSFT.mean in
    [x in MSFT : x.date >= jan2020 and x.value >= mean];

Observe que el filtro utiliza tanto la fecha como el valor de los puntos de la serie. Además, he omitido la transformación. Podemos omitir el filtro, la transformación (o proyección, en terminología SQL y C#) o incluso ambos.

Cuantificadores lógicos

De todas maneras podemos ir un poco más lejos que Python y Haskell. Vamos a comenzar por algo sencillo. ¿Es el 97 un número primo? Vamos a preguntarlo usando constructores de listas:

[all x in 2..96: 97 % x != 0]

La expresión anterior no devuelve una lista, sino un valor de tipo lógico. Será verdadero si alguno de los números entre 2 y 96 divide al número 97. Podríamos haber usado una cota superior más baja, por supuesto, pero no quiero complicar la explicación con detalles innecesarios.

all y any, en AUSTRA, no son palabras reservadas, pero en este caso se consideran palabras reservadas contextuales. La expresión anterior es equivalente a esta otra:

iseq(2..96).all(x => 97 % x != 0)

Esto, naturalmente, aunque sea ligeramente interesante, es sólo un rodeo hacia nuestro objetivo. ¿Qué tal si quiero todos los números primos del 2 al 1000?

[y in 2..1000 : all x in 2..96: y % x != 0]

La presencia de dos caracteres : nos está indicando que hay dos expresiones entre los corchetes. De hecho, estamos usando una función lambda anidada dentro de otra, y la más interna está «capturando» el parámetro de la más externa:

iseq(2, 1000).filter(y => iseq(2, y - 1).all(x => y % x != 0))

Y si quisiéramos elevar cada primo al cuadrado, añadiríamos una función de proyección al engendro que hemos creado:

[y in 2..1000 : all x in 2..96: y % x != 0 => y^2]

¡Chúpate esa, Python…!

Categorías
Música

Babelle

¡Tengo una sobrina artista! Naturalmente, el tío opina que la chica lo hace muy bien. Pero aquí le dejo el enlace a su última canción, para contrastar opiniones:

Es más guapa, por cierto, de lo que aparenta la foto, pero esa también es una opinión personal.

Categorías
Austra

Unfold

¿Cómo puedo calcular la serie de Fibonacci en Austra, utilizando una secuencia? Con vectores (voy a usar vectores y secuencias reales para evitar problemas de desbordamiento), es relativamente sencillo, si aprovechamos los safe indexers en una expresión lambda:

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

El problema de esta solución es que consume memoria. Tenemos los dos últimos números generados a la izquierda (por así decirlo) del número que estamos generando, pero en realidad, es porque tenemos todo la memoria del vector ya reservada. ¿Y si quisiéramos hacerlo con las nuevas secuencias?

La solución en lenguajes funcionales

Los lenguajes funcionales suelen usar una función unfold para estos casos. Por ejemplo, en F# se define unfold de esta manera:

val unfold : ('State -> ('T * 'State) option) -> 'State -> seq<'T>

La función, en sí, escrita como una secuencia infinita, sería más o menos esto.

let fib = Seq.unfold (fun (lastValue, currentValue)
    -> Some (lastValue,
    (currentValue, lastValue + currentValue))) (1, 1)

Ahora, las explicaciones: necesitaríamos que Austra fuese un lenguaje con tipos genéricos e inferencia automática de tipos… y no es un objetivo mío a corto plazo. Nyx, que es una idea de lenguaje con capacidades numéricas y aceleración que me guardo bajo la manga, tendrá genericidad, pero no estoy del todo convencido que tenga que ser funcional «puro». La función en F#, que sería casi idéntica en Haskell, necesita un «estado» para ir arrastrando en las sucesivas llamadas. En el caso de Fibonacci, el estado son los últimos dos números generados.

De regreso a Austra

No obstante, podemos hacer algo útil para el 99% de los casos. Estos son los tres métodos que se han implementado en C# para soportar unfold en el lenguaje de fórmulas:

public static DSequence Unfold(int size,
  double seed,
  Func<double, double> unfold);
public static DSequence Unfold(int size,
  double seed,
  Func<int, double, double> unfold);
public static DSequence Unfold(int size,
  double first, double second,
  Func<double, double, double> unfold);

La idea es proporcionar sobrecargas para los casos más frecuentes. El primer caso es para cuando el estado es un solo número real. El segundo caso es más sutil: el estado es un número real… más la posición del elemento que se va a generar. Si estuviésemos en un lenguaje funcional, esa posición tendría que estar explícitamente representada en el estado. Y en el tercer caso, el estado son dos números reales, como necesitamos para la secuencia de Fibonacci.

Veamos ahora cómo se usa esto en el lenguaje de fórmulas:

-- Potencias de 2, desde 2 a 1024.
seq::unfold(10, 2, x => 2x);
-- Serie de Maclaurin para exp(1).
1 + seq::unfold(100000, 1, (n, x) => x / (n + 1)).sum;

He tenido que sumar explícitamente un uno a la serie de Maclaurin, porque era lo más sencillo. Si quiere practicar, piense en cómo calcular ln(2) con una secuencia de estas.

La secuencia de Fibonacci se calcula así:

seq::unfold(50, 1, 1, (x, y) => x + y);

Con números enteros, haríamos esto:

iseq::unfold(30, 1, 1, (x, y) => x + y);

No es la solución perfecta, pero es fácil de entender, y encima se ejecuta con rapidez.

Algo curioso relacionado con el diseño de lenguajes: imagine un lenguaje que ofrece un vec<double> y un vec<complex>, en vez del vec y el cvec del lenguaje de fórmulas de Austra. ¿Se ha dado cuenta de que el vector de complejos no podría soportar todas las operaciones del vector de números reales? El caso más evidente: no se puede ordenar el vector de complejos, ni calcular la entrada con el menor o mayor valor. Los números complejos no soportan un orden total que tenga sentido. ¿Cómo haría usted para que una clase genérica tuviese en cuenta estos detalles? Ahora mismo, si quieres hacer esto en C#, tienes que definir una clase genérica «recortada», en plan Vec<T>, y las clases finales serán clases que ya no serán genéricas, pero que añadirían la funcionalidad no común. Algo así es que lo hace Math.NET.

Ese es el problema que intento resolver: un lenguaje de programación con auto-vectorización, en el que programar una librería como la de Austra no sea tan complicado, y que me deje definir tipos matemáticos como el vector de complejos y el vector de reales utilizando una base común genérica. Hay más cosas en las matemáticas que no se ajustan a la idea de la programación de que «lo derivado tiene más operaciones».

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# 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# 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.