Estimando con precisión el número de infectados de coronavirus

En una charla en Facebook salió el problema de a cuánta gente habría que hacer tests para conocer con precisión el número de infectados en España; concretamente, con un error del 1%.

Claro, un 1% es demasiada precisión, porque en el momento de escribir esto, el número oficial de infectados crece cada día un 10%. Olvidémonos de este detalle.

El resultado es chocante, porque hace falta una cantidad enorme de tests. Antes de entrar en faena, veamos dos casos extremos. Primero, tenemos un pueblo con 50 habitantes y queremos conocer el número de infectados con un error del 1%; bueno, es sencillo, hay que hacerle el tests a todos y cada uno de los vecinos. Segundo, tenemos una enfermedad rara que posiblemente padece una persona de cada millón; en realidad no lo sabemos, porque sólo se ha tratado a pacientes que iban a la consulta porque su caso era grave. Hacerle el test a un millón de personas no sería suficiente; posiblemente pillarías a un enfermo, pero a lo mejor no pillarías a ninguno; obviamente, esto no sería suficiente para estimar el número de enfermos con un error del 1%.

Vamos al grano. Imaginemos que en España hay 50M personas, y que el 0,4% de la población está infectada (el dato oficial era más o menos la mitad en el momento de formularse la pregunta). Hacemos un test a N personas. El resultado tendrá una media 0,004*N y una desviación típica raiz(N*0,004*0,996). Supongamos que usamos una confianza del 95%; esto quiere decir que aceptamos que el resultado puede ir desde -2 desviaciones típicas hasta 2, es decir, desde 0,004*N-2*raiz(N*0,004*0,996) hasta 0,004*N+2*raiz(N*0,004*0,996). Como queremos que este resultado tenga un error del 1%, 2*raiz(N*0,004*0,996) = 0,01*0,004*N. Y uf, eso quiere decir que N = 10.000.000.

¿Uh?

Sí, en serio. Si hago 10 millones de tests, encontraré un número de infectados que tendrá de media 40.000 y de desviación típica 199; de forma que si confío en que caerá entre 39.600 y 40.400 ya me aseguro el error de 1%.

Un ejemplo un poco menos chocante. Si sólo testeo a 10.000 personas, espero encontrarme 40 infectados, pero la desviación típica será 6,3, así que el número que me saldrá estará entre 40-2*6,3=27 y 53. Claro, 53 es casi el doble de 27; con 10.000 tests estoy muy, muy lejos del error de 1% en el número de infectados.

Esto es llamativo entre otras cosas porque estamos acostumbrados a que los sondeos electorales con 3.000 personas nos den un error del orden del 1%. Pero es que ahí los porcentajes buscados no son pequeños, los partidos minoritarios ni aparecen. Imaginemos que la mitad de la población estuviese infectada; si tomo 3.000 muestras me saldrá un resultado con media 1.500 y desviación típica raiz(3000*0,5*0,5)=27, así que con confianza 95% el resultado estaría entre 1446 y 1554, y tendría un error menor que 3.6%. Nada que ver.

Una pregunta interesante: imaginemos que hacemos 10.000 tests, pero en vez de escoger a 10.000 españoles al azar, usamos algún tipo de criterio geográfico. ¿Podríamos tener una estimación más precisa? La respuesta es que sí, pero quizás no mucho. Me invento un ejemplo.

Seguimos suponiendo que hay 50M españoles y el % = 0,4% están enfermos. Pero ahora España está dividida en dos "sitios"; el Norte tiene 30M habitantes y %n = 0,6% infectados, y el Sur tiene 20M habitantes y %s = 0,1% infectados (esta diferencia de porcentajes es exagerada). Nosotros no conocemos %, %n y %s, pero obviamente 50 * % = 30 * %n + 20 * %s. En vez de tomar una muestra de m = 10.000 españoles, tomaremos dos muestras de mn norteños y ms sureños, donde nosotros podemos escoger mn y ms con la retricción de que mn+ms = m = 10.000.

Al analizar las muestras nos encontraremos con in y is infectados.

La media de in es %n * mn y su varianza es mn * %n * (1-%n)

La media de is es %s * ms y su varianza es ms * %s * (1-%s)

Así que estimaremos el número total de infectados en el norte como (30M/mn) * in; esto tiene media 30M * %n y varianza mn * %n * (1-%n) * (30M/mn)^2.

Cuando sumamos las estimaciones del norte y del sur, nos sale la media correcta (aquí no hay sorpresa) y varianza

%n * (1-%n) * 30M^2 / mn + %s * (1-%s) * 20M^2 / ms

comparado con la encuesta inicial de una sola muestra para toda España, que tiene una varianza de

% * (1-%) * 50M^2 / m = 996M

Vale, ¿cuál es la mejor forma de escoger mn y ms? Pues resulta que la varianza se minimiza para mn = 7856 y ms = 2144, si lo hacemos así la varianza sale 870M. Claro, esto es un poco hacer trampa, porque para precisar tanto hay que conocer %n y %s. Lo que pasa es que tampoco es algo tan descabellado; si sospechamos que hay más enfermos en el norte que en el sur, parece lógico afinar el muestreo en el norte; si no tuviésemos ninguna preferencia, repartiríamos la muestra 6000/4000 para norte/sur por tener poblaciones 30M/20M; si en vez de 7856/2144 hubiésemos usado 8500/1500, la varianza nos saldría 898M, que tampoco es tan diferente. Si se nos fuese la olla y usásemos 9750/250, la varianza se dispararía y subiría hasta 2149M; pero vaya, cualquier matemático que pasase por ahí nos avisaria de que eso es un disparate.

A lo que vamos: si usamos una sola muestra, nuestra estimación tendría una media 50M * 0,004 = 200.000 infectados con una desviación típica = raiz(varianza) = 31.559 , mientras que si usamos dos muestras, lo menos que podemos conseguir es una desviación típica de raiz(870M) = 29.496. La mejora es como mucho un 7%, y corremos el riesgo de meter la pata y hacerlo mal.

Esto parece muy poco, posiblemente porque, de nuevo, comparamos con las encuestas electorales, que sabemos que consiguen mejoras mucho mayores.

Pero claro, para exprimir la información de una encuesta electoral, se usa un montón de información. Se tienen los resultados de elecciones anteriores, y si decides que vas a usar tal sitio paa hacer la encuesta, te puedes informar de si ahí ha cambiado el censo, si se ha construido una fábrica que atraerá obreros, o urbanizaciones con campos de golf, puedes obtener un montón de datos como estadísticas sobre la declaración de la renta, etc.

Con el coronavirus no tenemos esa información, así que no sabemos cómo sacarle tanto jugo a una muestra.

Si lo que nos interesase no fuese estimar el número de infectados sino el número de muertes, podríamos poner en nuestra encuesta más personas mayores y hombres, que mueren más de covid19 que las mujeres. Esto lo podríamos hacer de dos formas: o bien encuestamos más a hombre y a mayores en vez de escoger gente al azar (aplicable al problema de toda España o por zonas); o bien escogemos aquellos sitios donde haya más hombre o más mayores, y después seleccionamos al azar dentro de esos sitios escogidos (sólo para el problema por zonas).

¿Cuánto cuesta encontrar el número más probable en una ruleta?

Alguien me ha preguntado recientemente si teniendo los números que han salido en 6.000 partidas de ruleta se podría averiguar a qué número habría que apostar. La pregunta me ha intrigado...

Es sorprendente, pero 6.000 partidas no son muchas. "Muchas", "pocas" o "suficientes" dependerá de cuánto se aparten las probabilidades de 1/37. He googleado un poco, pero no he encontrado ningún sitio donde hablen de esto; los fabricantes afirman que sus ruletas están maravillosamente equilibradas y que usan la teoría del caos para diseñarlas de forma que sean impredecibles, pero por ninguna parte se menciona nada parecido a una tolerancia.

Admito mi ignorancia pero, como estoy intrigado, no me rindo, así que me voy a inventar una tolerancia.

Lo ideal sería que todos los números tuviesen una probabilidad de 1/37; de esta forma, apuestes al número al que apuestes, esperas perder 1/37 de tu apuesta, porque ganas 35 con probabilidad 1/37, pierdes 1 con probabilidad 36/37, y 35*1/37-1*36/37 = -1/37.

Si las probabilidades reales de un número se desvían menos de un 1%, no pasará gran cosa, el casino seguirá ganando bastante.

El casino empezaría a tener problemas si alguno de los números llegase a tener probabilidad 1/36 (que está un 2.7% por encima de 1/37). Si esto ocurre, es posible "jugar gratis" apostando a ese número, porque 35*1/36-1*(1-1/36)=0; es decir, ni se gana ni se pierde. Quizás el casino podría tolerar que un jugador hiciera esto mientras no se lo contase a los demás, o a lo mejor el casino no se daría cuenta porque le sería difícil detectarlo (lo veremos a continuación).

Lo que sí que no creo que ningún casino pueda tolerar es que un número de la ruleta llegue a tener una probabilidad de 38/(36*37), un 5.5% por encima de 1/37, porque en ese momento el jugador que conozca esto jugará con tanta ventaja sobre el casino como la que el casino tiene normalmente sobre los jugadores; 35*38/(36*37)-1*(1-38/(36*37))=1/37. Así que se ganaría un 1/37 de lo apostado, que es lo que normalmente hace el casino. No sólo sería intolerable, es que además me imagino que se notaría mucho si se hiciese mucho tiempo. Aquí "intolerable" no quiere decir nada mafioso: el casino cambiaría la ruleta y dejaría que el jugador siguiese jugando para que perdiese lo ganado anteriormente.

Así que me he inventado el criterio de que las ruletas en un casino legal (las cosas serán diferentes en timbas en furgonetas) pueden llegar a tener, como mucho, una probabilidad de que salga un número de 75/2664, que está a mitad de camino entre 1/36 y 38/(36*37).

Vale, quizás esto es un disparate; de hecho, yo apostaría a que será muy poco probable que en un casino con ruletas modernas haya algún número con probabilidad de salir por encima de 1/36. Pero bueno, imaginemos que precisamente queremos encontrar una ruleta que está estropeada. El caso es que ya tengo un límite con el que echar cuentas.

El problema entonces es el siguiente. Tenemos una ruleta en la que un número (no sabemos cuál) sale con probabilidad 75/2664, y todos los demás números salen con probabilidad 2589/95904. Tengo los resultados de 6.000 partidas. ¿Puedo decir que el número buscado es simplemente el número que aparece más veces en mi lista de 6.000 resultados?

Sorpresa: no; he simulado esto 100.000 veces, y sólo se acierta en el 8.4% de las ocasiones.

Esto puede parecer sorprendente, pero en realidad está claro por qué ocurre. En promedio, la diferencia entre que un número salga en el 1/37 o en el 75/2664 de 6.000 partidas es 162 y 169; pero la desviación típica del número de veces que sale es aproximadamente 12; así que la diferencia que buscamos (7) es menor que la variación de nuestro número debida al azar. Pero es que en realidad vamos a mirar al máximo del número de veces que han los 37 números; que tiene media 190 y desviación típica 6.1. Es decir, que lo normal es que alguno de los 36 números con probabilidad pequeña tenga suerte y salga 21 veces más que el que en principio debería ser el favorito.

Vale, con 6.000 partidas acertamos el 8.4% de las veces.

Pero con 20.000 partidas acertamos el 16% de las veces.

Y con 200.000 partidas acertamos el 82.6% de las veces.

Se me ha ocurrido otro "juego". Volvamos a una ruleta con un número que sale con probabilidad 75/2664. Primero observo sin jugar 6.000 partidas, y entonces empiezo a jugar, pero sigo contando cuántas veces sale cada número, y siempre apuesto al que ha salido con más frecuencia. En algún momento acertaré con el número correcto y empezaré a ganar sistemáticamente, porque mi ordenador va a jugar millones de partidas. La cuestión es, ¿cuál es la última partida en la que voy perdiendo? Para esto no me basta con simplemente encontrar el número correcto, sino que también tengo que recuperarme de las malas rachas que podría haber tenido al principio.

Los resultados son francamente decepcionantes. He simulado este juego 100.000 veces. En 39 ocasiones ha empezado a ganar desde el principio; lo normal es que se llegue a estar perdiendo unas 200.000 partidas; y en el 2,6% de las simulaciones hubo que jugar más de un millón de partidas antes de quedarse permanentemente en verde. Aquí están los resultados:

Intervalo	n
6000-6000	39
6001-6500	4067
6501-7000	2091
7001-8000	1980
8001-10000	2993
10001-20000	9320
20001-50000	15778
50001-100000	14678
100001-200000	15099
200001-500000	20520
500001-1000000	10781
1000001-inf	2654

Todo esto no encaja muy bien con las historias que todos conocemos de gente que ha hecho saltar la banca. En España son famosos los Pelayo, pero hay muchos héroes parecidos en otros países. Se me ocurren tres explicaciones:

1. Quizás cuando nos cuentan estas historias no enfatizan lo suficiente cuantísimas noches se pasaron apuntando números. O a lo mejor somos nosotros los que no escuchamos eso.
2. Podría ser que las ruletas sean mucho menos imparciales de lo que yo me he inventado. Imaginemos que la probabilidad de que salga un cierto número en una ruleta fuese el 4% en vez del 2.7% normal; entonces, con los datos de 6.000 partidas pasadas acertaríamos el 99.9% de las veces. Quizás el truco no está en trabajar mucho buscando el mejor número en una ruleta normal, sino en encontrar rápido una ruleta desastrosa. Lo que pasa es que me cuesta creer que el casino no se dé cuenta de que tiene una ruleta estropeada; vamos, cuando yo he estado en un casino, los números que salen en las ruletas van apareciendo en una pantalla, así que me imagino que algún ordenador debe de ir comprobando que los números que salen son normales. Quizás parte del truco está en buscar casinos negligentes en vez de limitarse a buscar ruletas defectuosas.
3. He estado pensando en el caso de que sólo un número sea más probable que los demás, pero a lo mejor es más fácil detectar si el eje de la ruleta está desviado; entonces la mitad de los números será un poco más probable que la otra mitad, y esto debería ser bastante más fácil de detectar estadísticamente que un número solo. (Por cierto, no importa mucho que la mesa de la ruleta esté inclinada si la ruleta está bien equilibrada; si, es verdad que las bolas acabarán con más frecuencia en la parte baja de la mesa, pero el número que habrá ahí irá cambiando.)

Aquí dejo el programa usado, en C++.

#include "stdafx.h" #include <string> #include <fstream> #include <iostream> #include "math.h" #include <time.h> #include <set> #include <vector> #include <algorithm> using namespace std; const int NDATOS = 6000; // Generación de número aleatorios unsigned int randomz, randomw; void initializeRandomSeq() { srand( (unsigned)time( NULL ) ); randomz = rand(); randomw = rand(); } double random() { // Get a random number with uniform distribution U(0,1) randomz = 36969*(randomz&65535)+(randomz>>16); randomw = 18000*(randomw&65535)+(randomw>>16); return 0.00000000023283064365386962890625*double((randomz<<16)+(randomw&65535)); } int ruletaJusta() { return (int)(37*random()); } // La ruleta desviada devuelve 36 con probabilidad provDesv, // y el resto de los números con probabilidad (1-probDesv)/36 const double probDesv = 75./2664.; int ruletaDesviada() { double r = random(); if (r < probDesv) return 36; return (int)((r-probDesv)*36/(1-probDesv)); } int _tmain(int argc, _TCHAR* argv[]){ clock_t time0 = clock(); initializeRandomSeq(); /* // Cálculo desviación típica número de veces que sale un número en 6000 partidas int * nEn6000 = new int[100000]; int * count = new int[37]; int suma = 0; for (int repeticion = 0; repeticion < 100000; repeticion++) { for (int i = 0; i < 37; i++) count[i] = 0; for (int partida = 0; partida < 6000; partida++) { count[ruletaJusta()]++; } int max = count[0]; for (int i = 1; i < 37; i++) if (count[i] > max) max = count[i]; nEn6000[repeticion] = max; suma = suma + max; } double promedio = suma / 100000.0; double stdev = 0.0; for (int i = 0; i < 100000; i++) { stdev = stdev + (nEn6000[i]-promedio)*(nEn6000[i]-promedio); } stdev = sqrt(stdev/99999); cout << "promedio = " << promedio << ", desviación típica = " << stdev << "\n"; delete[] nEn6000; delete[] count; */ // Cálculo de la probabilidad de acertar el número más probable de la ruleta int * count = new int[37]; int nAciertos = 0; int nFallos = 0; for (int rep = 0; rep < 100000; rep++) { for (int i = 0; i < NDATOS; i++) count[ruletaDesviada()]++; int imax = 0; int max = count[0]; for (int i = 1; i < 37; i++) { if (count[i] >= max) { max = count[i]; imax = i; } } if (imax == 36) { nAciertos++; }else{ nFallos++; } } cout << "se acierta el " << 100.*nAciertos/(nAciertos+nFallos) << " % de las veces\n"; delete[] count; /* // Apostar al número más frecuente en una ruleta no imparcial, y ver cuándo es la última partida // en la que pierdo int * count = new int[37]; for (int i = 0; i < 37; i++) count[i] = 0; int * countIntervalos = new int[12]; for (int i = 0; i < 12; i++) countIntervalos[i] = 0; int maxIntervalo[12] = {6000, 6500, 7000, 8000, 10000, 20000, 50000, 100000, 200000, 500000, 1000000, 2000000000}; for (int repeticion = 0; repeticion < 100000; repeticion++) { for (int i = 0; i < NDATOS; i++) count[ruletaDesviada()]++; // 6000 partidas iniciales int saldo = 0; // lo que llevo ganado int ultimaPartidaPerdiendo = NDATOS; // la última partida en la que perdí int partida = NDATOS; // la partida actual int n = 0; // el número que ha salido con más frecuencia int max = count[0]; // cuántas veces ha salido el número más frecuente for (int i = 1; i < 37; i++) if (count[i] > max) { max = count[i]; n = i; } // buscar el máximo while (partida < ultimaPartidaPerdiendo + 1000000 && partida < 1900000000) { // cuando llevemos un millón de partidas en verde, creemos que ya no perderemos nunca más partida++; int sale = ruletaDesviada(); if (sale == n) { saldo = saldo + 35; // aposté a n porque era el más frecuente, y he ganado }else{ saldo = saldo - 1; } count[sale]++; if (count[sale] > max) { // mirar si el número que acaba de salir se convierte en el más frecuente n = sale; max = count[sale]; } if (saldo < 0) ultimaPartidaPerdiendo = partida; // mirar si sigo perdiendo } cout << repeticion << " -> " << ultimaPartidaPerdiendo << "\n"; int intervalo = 0; while (ultimaPartidaPerdiendo > maxIntervalo[intervalo]) { intervalo++; } countIntervalos[intervalo]++; } for (int intervalo = 0; intervalo < 12; intervalo++) { if (intervalo == 0) cout << "6000-"; else cout << maxIntervalo[intervalo-1]+1 << "-"; if (intervalo == 11) cout << "inf"; else cout << maxIntervalo[intervalo]; cout << " " << countIntervalos[intervalo] << "\n"; } delete[] count; delete[] countIntervalos; */ cout << "time spent = " << clock() - time0 << " ms\n"; cout << "press any key to end\n"; char cinget; cin.get(cinget); return 0; }

Nota añadida: me confirman que los Pelayo detectaban una mitad de la ruleta en la que la bola tenía más probabilidades de acabar, como si tuviese el eje torcido. Esto tiene mucho más sentido que buscar un número ligeramente mejor que el resto.

Cadena más larga de 100 números

Le he estado dando vueltas a un problema curioso en http://simplementenumeros.blogspot.com.es/2013/05/1130-los-numeros-del-1-al-100.html. Se trata de formar la cadena más larga posible con los números del 1 al 100, de forma que cada número en la cadena sea o bien múltiplo o bien divisor de sus vecinos.

Una solución es la siguiente cadena de 77 números:

76-38-19-95-5-85-17-34-68-4-92-46-23-69-3-87-29-58-2-62-31-93-1-77-11-99-33-66-22-44-88-8-96-32-64-16-48-24-72-36-18-54-27-81-9-63-21-42-84-28-56-14-98-49-7-35-70-10-80-40-20-100-50-25-75-15-45-90-30-60-12-6-78-39-13-26-52

En realidad hay muchas soluciones con 77 números; por ejemplo, en esta cadena se pueden intercambiar las posiciones de 26 y 52, y de 40 y 80.

La cadena tiene la siguiente estructura:

• 4 múltiplos de 19, y el 5;
• 4 múltiplos de 17, y el 4;
• 4 múltiplos de 23 y el 3;
• 3 múltiplos de 29 y el 2;
• 3 múltiplos de 31 y el 1;
• y finalmente una cadena de 54 números múltiplos de 2, 3, 5, 7, 11 y 13.

Es fácil ver que en la solución no deberían aparecer los primos mayores que 50; 53,59,61,67,71,73,79,83,89 y 97; porque sólo se pueden conectar al 1. Ni tampoco aquellos primos entre 33 y 50 junto con sus dobles, 37 y 74, 41 y 82, 43 y 86, 47 y 94; porque sólo se pueden conectar al 1 y al 2. En cambio, los múltiplos de 17, 19, 23, 29 y 31 sí pueden aparecer, porque son 5 grupos con más de tres números que sólo se pueden conectar al 1,2,3,4,5. Si en la solución apareciese por ejemplo el 37, sería porque alguno de estos primos no aparecería, y por tanto se podría alargar la solución quitando el 37 y poniendo el primo que no estaba.

Aparte de estos 18 números (primos > 50 y primos y dobles entre 33 y 50) que sabemos que no pueden aparecer, hay otros 5 que no aparecen en la solución dada: 51, 55, 57, 65 y 91. El 51 y 47 son múltiplos de 17 y 19 y podrían aparecer, pero entonces tendrían que desaparecer el 85 y el 95.

Así que en realidad sólo se podría mejorar la solución en 3 números. Pero esto no es posible. La cadena larga, a partir del 1, es la más larga posible (demostrado por ordenador). Es verdad que podríamos quitar, por ejemplo, los múltiplos de 31 para dejar libre el 1, con lo cual podríamos enganchar algún otro número. Pero es que en la solución ya hay tres múltiplos de 31; incluso si pudiésemos enganchar los tres números que nos faltan, perderíamos otros 3, con lo cual no mejoraríamos la longitud de la cadena.

Total, que queda por ver cómo encontrar la cadena larga. El problema se puede formalizar como buscar el camino más largo posible en un grafo de 57 vértices. Se podría usar backtracking a lo bruto, pero esto tarda demasiado.

Una primera mejora para podar el árbol es, antes de probar a explorar un nuevo vértice, contar el número de vértices en su componente conexa. Si la componente conexa tiene N vértices, la cadena que ya he explorado tiene longitud L, y de momento el récord de cadena más larga que he encontrado es R, entonces no tiene sentido explorar si N+L <= R, no se podrá batir el récord. Esto mejora una burrada el tiempo de búsqueda, pero sigue siendo lento.

La clave para conseguir explorar todo el árbol de caminos es darse cuenta de que los grafos que quedan por explorar son bastante "fibrosos". Al principio hay unos cuantos números muy conectados, pero tienen a desaparecer pronto y dejar un grafo deshilachado. Por ejemplo, uno que aparece relativamente pronto y que estuve depurando en detalle es éste:

Se entra por el 35. Aunque tiene 21 vértices, está claro que los dos caminos más largos posibles son de sólo 11. Una forma de afinar la cota sería darse cuenta de que en ese grafo hay 6 hojas (vértices con un único vecino) y como sólo podremos usar como mucho uno de ellos (porque ahí se acaba el camino), aunque el grafo tenga 21 vértices podemos acotar los vértices útiles a 21-6+1.

Pero claro, podemos repetir esta observación... en un grafo como éste, si primero quitamos las hojas 21,25,35, nos queda otro grafo con hojas, así que le quitamos 3 y 5 y finalmente nos queda un grafo con 1 vértice sin hojas (15 no es hoja porque está conectado al principio del camino). Así que en ese grafo el camino más largo que podremos recorrer es de 1+2 vértices (1 es el número de vértices en el grafo final y 2 es el número de veces que hemos eliminado hojas).

Acotando de esta forma la longitud del camino máximo dentro de un subgrafo, se puede explorar todo el grafo de 57 vértices en unas 24 horas.

Historias y probabilidades

He estado leyendo un libro, "El andar del borracho: cómo el azar gobierna nuestras vidas". Aparte de la novelesca biografía de Cardano en el capítulo 3, hay una cosa que me ha intrigado.

Imagínate que te pregunto "¿qué es más probable, que a mi empresa le vaya bien el año que viene, o que a mi empresa le vaya bien el año que viene a causa de la recuperación económica mundial?"

Vale, la pregunta es un poco ambigua, y además habría que adornarla un poco y poner otras opciones para que la trampa no sea tan descarada, pero incluso con esta presentación tan chusca, esta pregunta ha desconcertado a varios de mis compañeros de trabajo, varios de los cuales escogieron la segunda alternativa.

Bueno, pues no, la más probable es la primera; siempre que ocurra que a mi empresa le vaya bien _Y_ haya recuperación económica mundial _Y_ la causa de que a mi empresa le vaya bien sea la recuperación económica mundial (tres condiciones), ocurrirá que a mi empresa le irá bien (una condición).

Esta observación debería de ser obvia, pero la gente cae en ella con facilidad en cuanto la trampa se disimula un poco. Cuando quieres convencer a alguien de algo improbable, en vez de analizar las cosas, lo que tienes que hacer es contarle una historia. Esto lo saben perfectamente los abogados que tienen que convencer a un jurado de algo que normalmente resultaría obviamente falso; se inventan una historia en la que cada paso conecte con el siguiente, y así consiguen esquivar la capacidad de la gente de estimar si algo es creíble o no. Eso sí, no intentes esto con el juez: el profesional pasará de tus milongas y mirará el informe de balística. Así es como funcionan también las conspiranoias.

Otro problema chulo de grafos

Hace tiempo leí otro de esos problemas chulos de grafos sin datos, vagamente parecido al que ya conté de los señores Mancha, aunque este es más fácil.

Resulta que en una fiesta los invitados cumplen la siguiente propiedad más bien especialísima: para cada par de invitados, existe exactamente un invitado al que conocen los dos. Pregunta: ¿cómo es la fiesta?

Formas de montar una escena

Un amigo mío que es director de películas me contó intrigado el siguiente problema. En un libro sobre composición, donde querían enfatizar lo complicado que puede ser montar una escena, contaban que si tienes n tomas, entonces puedes montarlas de N=e·n!-1 formas diferentes para construir la escena, lo cual implica que habrá tantas posibilidades de hacerlo que no será nada fácil decidir entre ellas.

Un ejemplo para aclarar el problema. Quieres mostrar una fiesta en una habitación y tienes cinco tomas; en ellas, la cámara enfoca la puerta, un sofá, el balcón, la mesa con la comida, y la pista de baile. Podrías desechar algunas tomas, o usar sólo un trozo de alguna toma, pero claro, tienes que usar al menos una toma.

¿De cuántas formas se puede hacer esto? Se podrían usar las 5 tomas, mostrándolas ordenadas de 5!=120 formas diferentes. También se podrían usar cuatro tomas; se podrían escoger de comb(5,4)=5 formas, y para cada elección, las tomas seleccionadas se podrían mostrar ordenadas de 4!=24 formas diferentes. O tres, o...

En términos matemáticos, se quiere el número de permutaciones de subconjuntos no vacíos de un conjunto de n elementos.

En total,

Así que la fórmula es rara pero correcta, al redondear el resultado al entero más cercano dará la respuesta correcta.

Ese -1 es peculiar. Sería raro que te preguntasen qué permutaciones hay en el término e·n!. Sin embargo, está muy claro cuál es la permutación que no hay en el -1: la del conjunto vacío.

El problema de las n reinas

La noticia mal contada es que ofrecen un millón de dólares por resolver el problema de las 8 reinas en el tablero de ajedrez. En realidad, el problema es que en un tablero de n x n ya tienes colocadas algunas reinas, y se trata de poner las que faltan. Bueno, no basta con hacerlo, hay que demostrar si existe o no un algoritmo que en tiempo polinómico decida si es posible hacerlo o no. Como han demostrado que este problema es NP completo, en realidad han subido el premio por resolver P = NP de uno a dos millones de dólares.

Pobrecillos, me pregunto cuántas cartas habrán recibido ya de gente reclamando el premio. Aunque lo mismo el propósito era hacerse publicidad aprovechándose de que los periodistas tienen que escribir artículos en cinco minutos.

http://www.hispantv.com/noticias/deporte/352348/premio-million-problema-ajedrez-ocho-reinas

http://jair.org/media/5512/live-5512-10126-jair.pdf

Los colores viven en un cono

El modelo de colores HSV es un plano proyectivo, donde el origen corresponde a los grises y las clases de equivalencia corresponden a "tonos de color" dentro de las que varía la saturación, junto con una dimensión adicional que es el brillo.

Vaya frikez; la próxima vez que me pregunten para qué sirven las geometrías exóticas, ya sé qué responder.

https://en.wikipedia.org/wiki/HSL_and_HSV

https://es.wikipedia.org/wiki/Modelo_de_color_HSV

Algoritmos para el problema de los vales

(Esta entrada es continuación de ésta).

Recordatorio: el problema consiste en que queremos comprar unos productos de precios p_i y nos dan un vale de descuento por cada compra de C euros o mayor; así que tenemos que dividir los productos en varias compras, de forma que consigamos el mayor número posible de vales de descuento. En realidad, el planteamiento original era que el precio tenía que ser mayor que C.

Vale, ya sabemos que es imposible encontrar siempre rápidamente la mejor solución. ¿Pero cómo encontramos rápidamente una solución razonablemente buena?

Mi amigo sugirió este algoritmo:

Yo sugerí este otro. Tiene la pega de que hace falta suponer que todos los precios son múltiplos de un céntimo para poder calcular rápidamente la lista de todos los posibles precios.

• Mientras el precio de los artículos que quedan por pedir exceda C:
  • Averiguar el precio de pedido más pequeño que se puede hacer que exceda C
  • Encontrar un conjunto de artículos que tenga ese precio
  • Hacer un pedido con ese conjunto
• Añadir los artículos que sobran al último pedido, si se ha hecho algún pedido, o hacer con ellos el primer pedido.

Los dos algoritmos se ejecutan en tiempo polinómico, así que o bien P=NP y lo hemos resuelto, o bien estos algoritmos no siempre encuentran la solución óptima. He aquí dos ejemplos en los que fallan:

• C=20 y los precios son {6,6,6,7,8,9}. En este caso, el primer algoritmo falla, porque al meter los dos precios mayores en la misma compra se queda sin encontrar la solución óptima, que es {9,6,6}, {8,7,6}. El segundo algoritmo sí que acierta, porque es posible hacer una compra de 21, luego la otra también tiene que ser de 21.
• C=20 y los precios son {7,7,7,15,15,15}. Ahora el primer algoritmo acierta, porque encuentra la solución {7,15}, {7,15}, {7,15}, pero el segundo algoritmo falla, porque el precio más bajo superior a 20 es 21, así que empieza haciendo la compra {7,7,7}, y luego no puede conseguir dos vales más.

El problema de los vales

He estado pensando con un amigo sobre el siguiente problema.

Quiero comprar en un supermercado n productos de precios estrictamente positivos {p_i}. Por cada compra que haga de más de C euros me dan un vale de descuento. ¿Cómo tengo que particionar los productos que quiero comprar en varias compras, de forma que el número de vales de descuento que consiga sea máximo?

Este problema tiene toda la pinta de ser NP-completo, y tras darle varias vueltas, hemos conseguido una demostración bastante fea, reduciendo a él el problema del bin packing, que sabemos que es NP-completo.

Empezaremos formulando los dos problemas formalmente:

1. Vales. Dados un número de productos n, sus precios {p_i>0}_i=1ⁿ, y la cantidad mínima C con la que obtener un vale, encontrar una partición {S_i} de {p_i} con el número máximo de conjuntos tales que Σ_Sip_i > C. Dada una solución, llamaré p_ki a los precios de los productos de la compra k.
2. Bin packing. Dada una colección ilimitada de contenedores de volumen V, un número de objetos n, y sus volúmenes {a_i<V}_i=1ⁿ, encontrar una partición {S_i} de {a_i} con el número mínimo de conjuntos tales que Σ_Sia_i ≤ V. Dada una solución, llamaré a_ki a los volúmenes de los objetos del contenedor k. (Estos objetos sólo tienen tamaño, no forma; por ejemplo, como al guardar ficheros en un disco duro.) (Podría haber objetos de volumen a_i=V, pero como ellos solos llenarían un contenedor, los quitaría del problema.)

La idea de la demostración es sencilla, consiste simplemente en tomar p_i = V-a_i, pero se complica bastante por el camino.

Vamos al grano. Supongamos que puedo resolver los problemas de vales en tiempo polinómico, y que tengo un problema de bin packing.

Para este problema de bin packing tendré un ε>0 definido de la siguiente manera:

ε < min( min{a_i}, (min{Σs_i : {s_i}⊂{a_i}, Σs_i>V} - V) / n)

Es decir, ε es menor que cualquier a_i, y de todas las sumas de volúmenes mayores que V, me quedo con la menor, resto V, y divido entre n. Este ε tiene la propiedad de que si Σ(a_ki-ε)<V, entonces Σa_ki≤V; porque, de lo contrario, Σa_ki>V y estaría dentro del conjunto de la definición de ε, luego ε≤(Σa_ki-V)/n, y por tanto V≤Σa_ki-n·ε; pero, dado que como mucho hay n términos en la suma, V≤Σ(a_ki-n), y de aquí la contradicción.

Claro, esta definición plantea un problema filosófico; dado que en principio hay un número exponencial de posibles sumas, ¿cómo calcular ε en tiempo polinómico? Es una buena pregunta, que afea esta demostración. En realidad, no hay que calcular ése epsilon, vale cualquier número menor. Por ejemplo, si los volúmenes de todos los objetos fuesen múltiplos de 1 centímetro cúbico, entonces podría tomar ε = 1cc / n, sin preocuparme de si el mínimo es 1cc o 2cc o 3cc o lo que sea.

Sigamos. Para un m arbitrario (que no tiene por qué ser el número mínimo de contenedores, es simplemente un número) puedo considerar el siguiente problema de vales (nótese que ε no depende de m):

1. Hay n·m productos;
2. el precio de los n primeros productos es V-a_i+ε;
3. el precio de los siguientes n·(m-1) productos, que llamaré comodines, es V;
4. para conseguir un vale tengo que hacer una compra de más de C = (n-1)·V.

A continuación demostraré cuatro proposiciones sobre este problema que depende de m.

Primera, todos los productos en este problema tienen un precio p_i≤V, y los únicos que tienen precio V son los comodines. Porque si V-a_i+ε>V, entonces ε>a_i.

Segunda, si el problema de vales con m tiene soluciones con m compras, entonces todas sus compras tienen exactamente n productos. Porque si una compra tuviese n-1 productos, la suma de sus precios sería como mucho (n-1)·V = C, de forma que no conseguiría un vale, así que cada compra contiene al menos n productos. Pero en total hay m compras y n·m productos, así que cada compra tiene que contener exactamente n productos.

Tercera, si el problema de bin packing tiene una solución con m contenedores en la que ningún contenedor está vacío, entonces el problema de vales tiene al menos una solución con m vales, en la que ninguna compra consiste sólo en comodines. Esta solución se obtiene a partir de la solución del problema de bin packing, convirtiendo cada contenedor en una compra, mandando los j_k objetos del contenedor k con volumen a_ki a productos de la compra k con precio V-a_ki+ε, y luego añadiendo n-j_k comodines de precio V hasta tener n productos; se consigue un vale por cada compra porque, como Σa_ki≤V, entonces

Σp_ki = Σ(V-a_ki+ε) + (n-j_k)·V
= j_k·V - Σa_ki + j_k·ε + n·V - j_k·V
= n·V + j_k·ε - Σa_ki
≥ n·V + ε - V = C + ε > C

donde sabemos que j_k≥1 por la condición de que hay al menos un objeto en cada contenedor, que se traduce en que en la compra k hay al menos un producto que no es un comodín.

Cuarta, si el problema de vales tiene una solución con m vales, entonces el problema de bin packing tiene una solución con m contenedores (a diferencia de la tercera proposición, aquí se permiten contenedores vacíos y compras de sólo comodines). Como antes, generamos la solución del problema de bin packing a partir de la solución del problema de vales, mandando la compra k al contenedor k; como Σp_ki>C=(n-1)·V, y en esa suma hay n-j_k comodines,

Σ(V-a_ki+ε) + (n-j_k)·V > (n-1)·V

j_k·V - Σ(a_ki-ε) + (n-j_k)·V > (n-1)·V

V > Σ(a_ki-ε)

Al definir ε se demostró que esa última desigualdad implica Σa_ki ≤ V, y por tanto se satisfacen las condiciones para ser una solución del problema de bin packing; además, como había al menos un producto que no era un comodín, tiene que haber al menos un objeto de volumen >0 en el contenedor k, que por tanto no está vacío.

Juntando todo esto, el siguiente algoritmo encuentra la solución del problema de bin packing en tiempo polinómico:

• Si no hay objetos, devolver m = 0.
• Para m = 1, 2, ...
•     Si el problema de vales con m tiene una solución, devolver m.

Para empezar, mientras m sea tan pequeño que el problema de bin packing no tenga soluciones porque m contenedores no bastan, el problema de vales tampoco la tendrá, por la cuarta proposición. Pero justo cuando alcancemos el m mínimo para que bin packing tenga una solución sin contenedores vacíos, entonces el problema de vales tendrá una solución, con m compras donde ninguna compra consistirá sólo en comodines, por la tercera proposición. Aquí hay que hacer la observación de que con el m mínimo no puede haber soluciones con contenedores vacíos... porque entonces eliminaríamos un contenedor vacío y tendríamos una solución con m-1 contenedores, de forma que m no sería mínimo.

Finalmente, ese algoritmo se ejecuta en tiempo polinómico porque el problema de bin packing contiene al menos n datos (los volúmenes de los n objetos), de forma que repetir algo n veces es polinómico en el tamaño del problema. Pero metiendo cada objeto en un contenedor diferente obtengo una solución trivial con n contenedores, así que el m mínimo tiene que ser menor o igual que n. Como cada ejecución del bucle se hace en tiempo polinómico, el tiempo total también lo es.

Quizás todo quede más claro con un ejemplo. Quisiera resolver el problema de bin packing consistente en guardar n=3 objetos de volúmenes 4, 6 y 7 en contenedores de volumen V=10. Como todos los volúmenes son múltiplos de 1, puedo tomar ε = 1/n = 1/3. Hay más de 0 objetos, luego no me basta con 0 contenedores. Pruebo m=1; es decir, resuelvo el problema de los vales con precios 10-a_i+ε = 6.33, 4.33 y 3.33, y me dan un vale por compras de más de (n-1)·V = 20; obviamente no hay solución, porque el precio de todos los productos juntos, 6.33+4.33+3.33=18, no llega a 20. Pruebo ahora con m=2; a los productos de antes añado tres comodines de precio V=10, así que tengo 6 productos de precios 6.33, 4.33, 3.33, 10, 10 y 10. Ahora sí que hay soluciones, en concreto una: las dos compras son {6.33, 4.33, 10}>20 y {10, 10, 3.33}>20; por tanto, hay soluciones del problema de bin packing, que obtengo eliminando los comodines y restando del volumen y ε, para quedarme con {4,6}≤10 y {7}≤10. También puedo mirar qué habría ocurrido para m=3; en este caso podría encontrar una solución como {6.33, 4.33, 10}>20, {10, 10, 3.33}>20 y {10, 10, 10}>20, que me indicaría que m es demasiado grande porque contiene una solución con n=3 comodines; pero también podría encontrar una solución como {10, 10, 3.33}>20, {10, 10, 4.33}>20 y {10, 10, 7.33}>20, que me indicaría que 3 contenedores son suficientes, pero me dejaría con la duda de si hay soluciones con menos contenedores (como el algoritmo prueba todos los m en orden creciente, nunca llego a tener la duda, siempre habré probado m=2 antes de probar m=3).

Ahora un ejemplo en el que el truco con ε resulta crucial. Tendremos n=2, volúmenes 4 y 6, y C=10. La solución es, obviamente, {4, 6}≤10. Pero si no fuese por el ε, el problema de vales no tendría solución, ya que los precios serían 10-4=6 y 10-6=4, pero claro, 6+4≯C=(n-1)·V=10. En este caso ε tiene que ser menor que 4 (no hay sumas mayores que 10), así que si tomo ε=3, los precios son 10-4+ε=9 y 10-6+ε=7. Ahora sí que existe una solución para m=1, {9, 7}>10, y por tanto {10-9+ε, 10-7+ε} = {4, 6} ≤ 10 es una solución del problema de bin packing.

Es una reducción muy fea, lo admito...

Dos demostraciones de que existen infinitos números primos

Mientras pensaba otro problema se me han ocurrido estas dos demostraciones de que existen infinitos números primos. No es que éste sea un resultado novedoso, ¿no?, pero como son medio chulas, pues las pongo aquí.

Primera demostración. Sean r₁, r₂, ... , r_n n números impares relativamente primos entre sí. Entonces el número

N=3·(2^{r₁·r₂·...·r_n}-1)

tiene al menos n+1 factores impares relativamente primos entre sí. Lo que implica que existen al menos n+1 números primos, y por inducción, existen tantos primos como se quiera.

Los n+1 factores son 3 y s_i = 2^n_i-1. Es fácil ver que 3 divide a N pero no a s_i; también se ve que s_i divide a N porque para todos a y b impares

2^ab-1 = (2^a-1)·(2^a(b-1)-2^a(b-2)+2^a(b-3)-...+2^2a-2^a+1)

Falta demostrar que los s_i son relativamente primos entre sí; esto es un caso particular de la siguiente proposición: dos números a y b son relativamente primos entre sí si y sólo si lo son 2^a-1 y 2^b-1. Como

2^a+b-1 = 2^a(2^b-1)+2^a-1

tenemos que

resto ( 2^a+b-1 , 2^b-1 ) = resto ( 2^a-1 , 2^b-1 )

y

mcd(2^a+b-1, 2^b-1) = mcd(2^a-1, 2^b-1)

y por tanto se puede aplicar el algoritmo de Euclides para calcular el máximo común divisor entre 2^a-1 y 2^b-1 simplemente al número de unos. Finalmente, a y b son primos relativos si y solo si el algoritmo de Euclides aplicado a a y b da 1, si y solo si el algoritmo de Euclides aplicado a 2^a-1 y 2^b-1 da 1, si y solo si 2^a-1 y 2^b-1 son primso relativos. Quizás un ejemplo lo aclare; para a=7 y b=17:

mcd(7, 17)	mcd(27-1, 217-1) en binario
mcd(7 ,17) =	mcd(11111112, 111111111111111112 =
mcd(7 ,3+7*2) =	mcd(11111112, 1112 + 11111112 · 100000010002) =
mcd(7 ,3) =	mcd(11111112, 1112) =
mcd(3 ,7) =	mcd(1112, 11111112) =
mcd(3 ,1+3*2) =	mcd(1112, 12 + 1112 · 100102 =
mcd(3 ,1) =	mcd(1112, 12) =
1	12

Segunda demostración. Sea p un número primo, y q un divisor primo de 2^p-1. Entonces 2^p=1 (mod q), luego p es un múltiplo del orden de 2 en N/qN. Pero p es primo y por tanto sólo tiene dos divisores; como el orden de 2 no puede ser 1, entonces el orden de 2 tiene que ser p. Pero en N/qN todos los órdenes son menores que q, por tanto p < q. De forma que si p es primo, todos los divisores primos de 2^p-1 son mayores que p, y por tanto, para cada primo p existe otro primo mayor, que se puede encontrar "simplemente" factorizando 2^p-1, de aquí que existan infinitos primos.

Lo injusto de los sorteos por letras

Todos hemos visto alguna vez algún sorteo en el que se escoge una letra aleatoriamente y los ganadores se seleccionan por orden alfabético, empezando por la primera persona cuyo apellido tenga esa letra como inicial. Todos nos hemos dado cuenta de que es un sistema injusto; por ejemplo, es un chollo llamarse Abad, porque probablemente conseguirás aquello que se sortee tanto si sale la W como si sale la X, Y, Z o la A.

La mayoría de nosotros, incluido yo, hemos pensado que la diferencia de probabilidades sería pequeña... bueno, pensar, lo que se dice pensar, no lo hemos pensado, simplemente lo hemos dado por supuesto.

Bueno, pues no, la diferencia no es pequeña. Yo me he dado cuenta después de que mi amigo Raúl Corvillo me llamase la atención sobre dos blogs, "Un dato vale más que mil palabras" y "La ciencia para todos". También está este artículo de una mujer cabreantemente apellidada Grima.

El problema es que las iniciales de los apellidos en España están distribuidas de una forma más caprichosa de lo que se podría esperar; estos datos están sacados del primer blog:

	frecuencia	% sobre el total
A	2.884.390	6,7%
B	2.263.664	5,2%
C	3.969.992	9,2%
D	1.747.696	4,0%
E	781.910	1,8%
F	1.877.528	4,3%
G	4.857.351	11,2%
H	992.297	2,3%
I	424.730	1,0%
J	722.854	1,7%
K	55.885	0,1%
L	2.250.441	5,2%
M	5.291.515	12,2%
N	699.534	1,6%
O	803.973	1,9%
P	3.042.595	7,0%
Q	185.195	0,4%
R	3.565.620	8,2%
S	3.201.882	7,4%
T	1.425.424	3,3%
U	171.705	0,4%
V	1.631.083	3,8%
W	48.578	0,1%
X	14.690	0,0%
Y	92.553	0,2%
Z	269.539	0,6%
TOTAL	43.272.624	99,8%

Por si algún día le pasase algo al blog de Eduardo, también hago una copia de su histograma, que ilustra maravillosamente lo irregular que es la distribucion de las iniciales:

Para echar las cuentas, he escrito este programita:

#include "stdafx.h" #include <iostream> using namespace std; int _tmain(int argc, _TCHAR* argv[]){ int i, winner, letter, imax, imin; double freq[26] = {6.7, 5.2, 9.2, 4.0, 1.8, 4.3, 11.2, 2.3, 1.0, 1.7, 0.1, 5.2, 12.2, 1.6, 1.9, 7.0, 0.4, 8.2, 7.4, 3.3, 0.4, 3.8, 0.1, 0.0, 0.2, 0.6 }; double pfirst[26]; double plast [26]; double sum = 0.0, p, noreach=0.0, percent, max, min, nochance = 0.0; for (i=0; i<26; i++) { sum += freq[i]; pfirst[i] = 0; plast [i] = 0; } cout << "sum = " << sum << "\n"; cout << "escriba la probabilidad de ganar como porcentaje:\n"; cin >> p; cout << "\n"; for (winner=0; winner<26; winner++) { percent = p; letter = winner; if (freq[letter]>p) nochance += freq[letter]-p; while (percent > freq[letter]) { pfirst[letter] += 100.0/26.0; plast [letter] += 100.0/26.0; percent -= freq[letter]; letter = (letter+1)%26; } pfirst[letter] += 100.0/26.0; } max = 0; min = 10000; for (i=0; i<26; i++) { if (pfirst[i]>max) { max = pfirst[i]; imax = i; } if (plast[i]<min) { min = plast[i]; imin = i; } } cout << "% de poblacion que no puede ganar: " << nochance << "\n"; cout << "la probabilida de que gane la primera " << (char)(65+imax) << " es " << pfirst[imax] << "\n"; cout << "la probabilida de que gane la ultima " << (char)(65+imin) << " es " << plast [imin] << "\n"; }

Empecemos imaginando que se sortea algo que va a ganar el 1% de los participantes; por ejemplo, se sortean 100 plazas de campamentos para niños y se presentan 10.000. Resulta que el 79% de los participantes no tiene ninguna posibilidad de conseguir una plaza, mientras que los 100 primeros solicitantes cuyo apellido empiece por la A tienen nada más y nada menos que el 19,23% de probabilidades de conseguirla.

Hasta aquí no hay ninguna sorpresa, era obvio que el sistema era especialmente injusto si la probabilidad de ganar era pequeña.

La auténtica sorpresa empieza a aparecer cuando descubrimos que, al aumentar el número de premios, el sistema no se hace justo rápidamente. Por ejemplo, si la probabilidad "global" de ganar es un razonable 15%, hay gente que tiene un 31% de posibilidades de ganar (las primeras A) mientras que otros tienen sólo un 4% (las últimas G).

¿Y si el porcentaje de premiados fuese el 50%? Ahora las primeras C ganan con probabilidad 62%, mientras que las últimas M sólo tienen un 38% de probabilidad de ganar. Es decir, Cabaretero tiene un 60% más de probabilidades de ganar que Mutante.

Si vamos al extremo de dar premios al 90% de los solicitantes, resulta que las primeras A, G y M tienen garantizado el premio, mientras que las últimas M sólo tienen un 69,23% de probabilidades de conseguirlo.

Supongamos que estuviésemos dispuestos a aceptar que un sorteo es "tolerablemente injusto" (he aquí un oxímoron) si al más beneficiado le da sólo un 25% más de probabilidades de ganar que al más perjudicado. ¿Cuándo es tolerable un sorteo basado en las iniciales de los apellidos? Sorpresa, sólo si al menos el 97% de los participantes van a conseguir el premio. No era esto lo que yo me esperaba.

Es previsible que este sistema funcione incluso peor en pueblos pequeños donde algunos apellidos se repiten mucho.

Francamente, me parece un sistema inaceptable, especialmente cuando hay una solución obvia y al alcance de cualquier notario del siglo XXI: se numeran las solicitudes (ya sea alfabéticamente, usando otro criterio, o sin usar ningún criterio), se escoge un número al azar, y se otorgan los premios a partir de ahí.

Simulación de olas

Le estamos poniendo olas al simulador de aerogeneradores de Gamesa para ver qué pasa cuando el oleaje aporrea la base de la torre de un aerogenerador en el mar. En esta simulación, las olas llegan a tener hasta 7 metros de altura y el agua tiene 20 metros de profundidad. Las flechas indican la velocidad del agua. La mayor parte del trabajo es de David Campillo.

He intentado publicar aquí un video en formato .avi, pero no he conseguido que funcione, así pongo un enlace en Facebook.

Sintaxis bizarra

En el trabajo, al hacer un interfaz entre MATLAB y una librería nuestra en C++, hemos seguido unas instrucciones quizás un poco anticuadas, y nos hemos encontrado con unas líneas bizarras de C++. Tras una investigación sintáctica, las rarezas se pueden simplificar a esto:

double a[2][2][2][2];
double (*data)[2][2][2][2];
data = new double[110][2][2][2][2];
data = (double (*)[2][2][2][2]) new double[32];

La única línea que resulta familiar es la primera; declara que a es una variable de un tipo peculiar, está formada por 16 doubles en un array unidimensional (sin los cuatro niveles usuales de indirección) al que se accede a través de no un índice, sino 4; es decir, el decimocuarto double no es a[13] sino a[1][1][0][1] (o quizás a[1][0][1][1], no estoy seguro ahora).

La segunda línea, que ya empieza a ser rara, simplemente declara que data es un puntero a un objeto del tipo descrito antes.

La tercera línea construye un objeto del tipo anterior pero de dimensiones [110][2][2][2][2]. Obsérvese que esto es posible porque todas las dimensiones son estáticas, y por tanto el compilador puede saber qué tamaño tendrá el objeto construido. La gracia está en que, aunque este objeto no tiene las dimensiones de los objetos a los que apunta data, en realidad se puede ver como un array de 110 objetos del tipo de data, luego el puntero del new se puede convertir al tipo de puntero de data. O algo parecido.

La cuarta línea construye un array de 32 doubles de los de toda la vida, y entonces convierte el puntero al tipo de puntero adecuado para poder asignarlo a data, de forma que los 32 doubles se ven como dos objetos del tipo double[2][2][2][2]. Una de las cosas llamativas de esta sintaxis es que, si le quitas los paréntesis al (*), el compilador no te entiende.

Para quien se haya quedado intrigado; Visual Studio es capaz de compilar esta línea, pero ¿cómo la interpreta?

double** b = (double**)(double* (***)[2][2][10]) new (double* (**)[2][50][2]) ();