Un programa para calcular el índice de Brillouin
Siguiendo con cuestiones relacionadas con la Ecología teórica, hoy nos centraremos en otro índice de biodiversidad histórico, el índice de Brillouin, uno de los más antiguos aún en activo. Su fórmula es relativamente sencilla, pero el índice tiende a estar en desuso porque puede ser bastante exigente en cuanto a capacidad de cálculo y por ello muchas hojas de cálculo son incapaces de manejarlo. Para evitar este problema en torno a los años 50 se creó una versión simplificada, el índice de Shannon, que hemos visto en el post anterior y es hoy uno de los más extendidos por su sencillez de cálculo (ganada a costa de sacrificar parte de la precisión y ventajas estadísticas que ofrece Brillouin).
Afortunadamente los debianitas somos tipos que siempre guardamos un par de recursos bajo la manga y tenemos a quien poder recurrir cuando, como en éste caso, necesitemos realizar operaciones con números cuyo tamaño puede alcanzar fácilmente los 50.000 dígitos de extensión; así que he escrito el siguiente programa capaz de encargarse de calcular todo por nosotros, en nada de tiempo y sin despeinarse apenas, gracias a una de esos aliados poderosos de los que hablaba, una auténtica bestia llamada Python:
#!/usr/bin/python
# -*- coding: utf-8 -*-
# file brillouin.py
import sys, getopt
from math import log, factorial
def main(argv):
try:
opts, args = getopt.getopt(argv, "hci", ["help","copyright","info-brillouin"])
except getopt.GetoptError:
print "\n Opcion desconocida. Para ayuda use --help"
sys.exit(2)
for opt, arg in opts:
if opt in ("-h", "--help"):
print '\nEste programa calcula el indice de biodiversidad de Brillouin.\n\n uso: $ python brillouin archivo_de_datos\n $ ./brillouin archivo_de_datos\n $ python brillouin -h [--help]\n $ python brillouin -c [--copyright]\n\n Para mas info vea tambien: \n $ python brillouin -i [--info-brillouin]'
sys.exit()
elif opt in ("-c", "--copyright"):
print '\nCopyright (c) 2011 P. Valdés.'
print 'License: This software is provided as-is, without any express or implied warranty. In no event will the author be held liable for any damages arising from the use of this software. \n\n Permission is granted to anyone to use this software for any purpose, including commercial applications, and to alter it and redistribute it freely, subject to the following restrictions: \n\n 1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required. \n 2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.\n 3. This notice may not be removed or altered from any source distribution.'
sys.exit()
elif opt in ("-i", "--info-brillouin"):
print '\nÉste programa permite calcular rápidamente el índice de Brillouin HB = (ln(N!) - ln(Σni!)/N, una medida de biodiversidad utilizada comunmente en Parasitología y Ecología.\nSalvo que se trabaje con comunidades muy sencillas de muy pocos individuos (algo que no suele ocurrir en la vida real) determinar su valor suele\n implicar manejar números cuyo tamaño alcanza fácilmente las decenas de millar de dígitos (ej: un 1 seguido de 50.000 ceros), eso hace que muchos programas y hojas de cálculo tradicionales sean incapaces de computar su valor y arrojen un error en su lugar.\n\n El índice se usa para determinar la diversidad de una poblacion bajo dos posibles premisas:\n 1- Toda la poblacion ha sido muestreada (habitual en comunidades de parásitos), o\n 2- El muestreo no es aleatorio (ej: insectos atraidos hacia trampas de luz)\n\nSi lo que tenemos es una muestra aleatoria de parte de la población debe de utilizarse el índice de Shannon (algo menos preciso) en su lugar.\n\n El archivo_de_datos es un simple archivo de texto con líneas de tipo "etiqueta_de_especie: número_de_ejemplares_encontrados_para_la_misma", ej:\n\n Molicola horridus: 210\n Argulis scutiformis: 12\n Lepeophtheirus nordmanni: 2370\n Floriceps sp. A: 670 ... etc\n\n El separador entre la cadena y el nº debe de ser obligatoriamente un signo ":", el programa es flexible frente la presencia de espacios en blanco extras alrededor, pero acabará con error ante la presencia de cualquier otro tipo de separador o si encuentra cantidades ilógicas (que no puede convertir a un nº entero positivo). Asegúrese especialmente de sustituir todos los tabuladores del archivo por espacios antes de usarlo.'
sys.exit()
if __name__ == "__main__":
main(sys.argv[1:])
(denominador, numerador, d) = (1,1,{})
with open(sys.argv[1]) as data: # Necesita indicar un archivo de datos válido
for line in data:
line = line.strip() # chomp
(k, v) = line.split(':',2)
d[k] = int(v) # carga el diccionario
s = len(d)
n = sum(d.itervalues())
maxv = max(d.itervalues()) #valor maximo
print '\n La muestra contiene', n, 'ejemplares de', s, 'especies distintas.\n La especie mas abundante tiene', maxv, 'ejemplares.'
# ni!*...*nj!
for val in d.itervalues():
if val == 0: continue
elif val == 1: continue
elif val == maxv: continue #simplificamos valmax!, mucho mas rapido
elif val < 0:
print '\nERROR: Revise sus datos.\n No puede tener', val, 'individuos para una especie.'
break
else: denominador *= factorial(val)
# n*(n-1)*(n-2)*... *maxval+1
for i in range(maxv+1, n): numerador *= int(i)
hb = (log(numerador/denominador))/n
print '\n Indice de biodiversidad de Brillouin\n para la muestra examinada = ', hb
data.close()
Tras algún tiempo peleando con la versión de Perl, crear la siguiente versión alternativa con Python me ha llevado un tiempo ridículo en comparación. Para su uso y documentación, podemos copiarlo a un archivo al que llamaremos brillouin. Luego desde terminal, nos situamos en el mismo directorio y lanzamos el comando:
$ python brillouin --helpY eso es todo, Brillouin y Shannon toman valores muy similares y comparables, pero cuando o bien nuestro muestreo comprenda a toda la población o bien sospechemos que no ha sido totalmente aleatorio, obtendremos siempre mejores resultados con el primero.
como sigas haciendo artículos de estos, cualquier día me quedo colgado
me parecen muy interesantes, aunque no acabo de comprenderlos del todo.
Buena entrada y gran idea la de escribir tu propio programa para calcular el Brillouin. Pendiente de posibles futuras entradas
El programa simplemente coge una formula (puedes verla en su documentación), llama a un par de funciones y la calcula usando algunos trucos para ser más eficiente y ahorrar CPU. El quíz de la cuestión es que la fórmula alcanza fácilmente números que pueden ser enormes sobre los que aún tiene que operar para obtener el resultado. Me consta, porque lo he probado, que con suficiente tiempo y una CPU casera echando humo durante un rato el programa es capaz de operar con un número como el siguiente: un 1 seguido de un millón doscientos cincuenta mil ceros. Como comparación un trillón es un 1 seguido de 18 ceros. Es un caso extremo usado simplemente como test pero que ilustra bien el nivel de cálculo en el que nos movemos (de hecho "podría" aparecer en un escenario real). Muchos programas simplemente no pueden con ello así que esa es la razón que me llevó a escribir el programa. También es un ejemplo muy gráfico de lo que son capaces de hacer los grandes como Python o Perl
Salvo que trabajes en temas de conservación probablemente no le encontrarás una utilidad directa, pero sí que puede resultar rentable estudiarlo con calma, especialmente los detalles. Es un ejemplo de como crear un programa para terminal que incluye la posibilidad de pasarle opciones y argumentos.
Si cogemos un programa como por ejemplo cp y vemos su documentación con man cp tenemos algo parecido a esto:
SINOPSIS
cp [OPCION] ORIGEN DESTINO
Como sabemos lo que aparece entre corchetes es optativo, mientras que los argumentos ORIGEN y DESTINO son obligatorios.
En Linux OPCION suele tener una forma corta y otra larga, más fácil de recordar. La forma larga usará dos guiones (por poner un ejemplo concreto, en cp podríamos pasarle la opcion --force) y la forma abreviada usa uno solo. El programa cp hará exactamente lo mismo tanto si le pasamos --force como -f. En el programa que pongo arriba dos equivalentes similares entre los disponibles son --help y -h por ejemplo, y leyendo el código con calma veremos como hacer eso.
Otro detalle interesante, si le pasamos una opción que no esté dentro de la lista de opciones reconocidas (por ejemplo -x) lanza un mensaje de error, a poco que trasteemos con el programa y echemos un vistazo a la documentación de python veremos como hacer eso fácilmente también.
Si hay alguna línea que no comprendes bien simplemente pregunta.
Un saludo
Me alegro de que te haya gustado Debish, un saludo a tí también. Ya tengo otro par de posts en el horno que irán saliendo a su tiempo
Si mis escasos conocimientos y google no me han traicionado
log(10^1.000.000)=log(2^x)1.000.000/log2=x
10^y=2
y=0.301
x=3322259.136
3322259.136/8(bits/byte)=415282.392
415282.392/1024=405.549211kb
REVISAR.
un número igual a 1 seguido de un millón de ceros ocupa menos de medio Mb.
¿Estamos hablando de manejar cuántos números de esos?
Lo bueno que tienen estos lenguajes de scripts (interpretados) es que la ausencia de tipos facilita el manejo de números grandes (ignoraremos que existe typedef) pero su ejecución siempre será más lenta que la de un programa compilado escrito en C. El auténtico dios padre creador de todas las cosas XDDD.
Un saludo, maestro.
mmmh... creo que estamos mezclando dos cosas diferentes aqui, de entrada log(x) es un logaritmo neperiano aquí, pero eso es una cuestion accesoria
Un número de 1 millon de caracteres ocupa simplemente 8 bites por cada caracter, más o menos el doble de lo que te sale (casi un mega). Ciertamente no es mucho realmente, hoy uno puede enchufar un disco externo de muchos gigas por poco dinero, luego puede almacenar miles de números como ese en su ordenador, no es problema en absoluto y si imprimimos el número a un archivo podemos ver que su tamaño es bastante manejable
Pero otra cosa es que el procesador de su CPU pueda manejarlo (multiplicarlo por otro número similar, por ejemplo), de hecho donde realmente se nota el cuello de botella al lanzar el programa es en tiempo de procesador, no en la ram. Tiene algo que ver con no-se-que de los megaflops creo.
No soy experto en hardware y la verdad es que ignoro cual es la capacidad tope de un procesador hoy día (Seguramente mucho mayor de lo que me imagino) así que es fácil que diga alguna tontería, pero lo que sí noto es que salvo que tomemos precauciones bastantes programas se bloquean alrededor de 147! por ejemplo, que es un número de 256 cifras. 147 ejemplares es una cantidad fácil de alcanzar en muchos muestreos. "Algo" ahí tiene un tope de capacidad
Si la ALU de la CPU falla supongo que dejará paso a la FPU y si ésta no es lo bastante potente es irrelevante que tengas un Tera disponible de disco duro porque en algun momento llegará a un número más allá del cual no puede aumentar esa cantidad, y una simple operación de súmale 1 a X, fallará.
> un saludo maestro
Otro saludo también a tí cenizoish, pero de "maestro" en estas cuestiones tengo poco, la verdad... XD
bueh, ya me mareo un poco. Estaba pensando si me iban a caber en la RAM unos cientos de números de esos
es.m.wikipedia.org/wiki/Operaciones_de_punto_flotante_por_segundo
es.m.wikipedia.org/wiki/Operaciones_de_punto_flotante_por_segundo
Vaya, que interesante... así a lo tonto nos hemos enterado de cuales son las computadoras más rápidas del mundo: 8 billones de cálculos por segundo, wow... (yo quiero unoo)
Además mientras estaba tratando de hacerme a la idea del asunto he descubierto la velocidad de mi propio cerebro. Puesto frente a una operación compleja mi cerebro típicamente tiene una capacidad de dos flops.
(uno cuando cae al suelo al ver la operación y el segundo al rebotar... flop, flop)
Bueno, hombre, no te subestimes.
Se estima que el poder de cómputo requerido para procesar olores, sabores, tacto, visión y coordinación motora es del orden de 10 PFLOPS (10 veces el poder de cómputo de «Roadrunner»).
Yo opero a velocidad del orden de mFLOPS, arquitecturas obsoletas, ya sabes. Chau, me piro más rápido que el roadrunner ese XDDD.