lunes, 10 de diciembre de 2018
No enseñes C... si lo que quieres es enseñar C++
Aquí dejo una muy interesante charla (en inglés) pertenenciente a la conferencia de C++ 2015 sobre cómo enfocar la docencia de un curso introductorio de C++. Me toca de cerca ahora, y me ha resultado muy útil y gratificante ver que, si bien a un nivel muy básico, coincide plenamente con el enfoque que creo debería plantearse en la enseñanza de C++ para Biología Computacional. Espero que la disfrutéis:
sábado, 3 de febrero de 2018
Información y apareamiento no aleatorio
Os dejo un enlace a mi último trabajo en el que formalizo la desviación del apareamiento al azar en términos de la teoría de la información. La desviación del apareamiento al azar tiene importancia evolutiva en la medida en que como consecuencia de esa desviación se alteran, de manera sistemática, las frecuencias fenotípicas y genéticas de la población.
El apareamiento aleatorio implica información cero, las desviaciones del apareamiento al azar se computan en cantidad de información. El índice de información general se pude expresar como la suma de componente dentro de sexo y entre sexos.
Aviso a navegantes de que el trabajo es bastante matemático (al menos desde el punto de vista del biólogo promedio, esperemos que esto irá cambiando). Dejo un enlace al artículo completo aquí (acceso libre provisto por la editorial hasta 25 Marzo, también hay un preprint en biorxiv o al autor directamente por email):
https://authors.elsevier.com/a/1WVHB9EXxAG-p
El apareamiento aleatorio implica información cero, las desviaciones del apareamiento al azar se computan en cantidad de información. El índice de información general se pude expresar como la suma de componente dentro de sexo y entre sexos.
Aviso a navegantes de que el trabajo es bastante matemático (al menos desde el punto de vista del biólogo promedio, esperemos que esto irá cambiando). Dejo un enlace al artículo completo aquí (acceso libre provisto por la editorial hasta 25 Marzo, también hay un preprint en biorxiv o al autor directamente por email):
https://authors.elsevier.com/a/1WVHB9EXxAG-p
sábado, 25 de noviembre de 2017
Miríadas...
... de contrastes estadísticos. El problema de la corrección de múltiples test, algo que jamás pudieron imaginar los estadísticos de comienzos del siglo XX, pero que tampoco lo pudieron hacer los de mediados de siglo que a lo sumo, pensaron en la repetición de unas decenas de tests, pero ¿miles? ¿decenas de miles? ¿centenares de miles?? ¡qué locura!.. pero ahí estamos. Aquí mi reciente aportación al asunto: https://doi.org/10.1093/bioinformatics/btx746
que espero poder convertir en una entrada divulgativa sobre el problema de los contrastes múltiples, también llamado coloquialmente 'corrección multitest' en algún futuro no muy lejano (aunque quién sabe, dada la extraña costumbre del tiempo de fluir y desaparecer sin despedirse).
que espero poder convertir en una entrada divulgativa sobre el problema de los contrastes múltiples, también llamado coloquialmente 'corrección multitest' en algún futuro no muy lejano (aunque quién sabe, dada la extraña costumbre del tiempo de fluir y desaparecer sin despedirse).
martes, 12 de julio de 2016
Botas, batas y bits
La verdad es que hace no tanto tiempo que si alguien decía -soy bióloga- lo primero que se pensaba era, ¡de bota! y se podía imaginar a la intrépida investigadora, amiga de los animales y las plantas, recorriendo el monte, identificando, clasificando y recogiendo muestras de todo tipo. Hacia la última década del siglo XX algunos biólogos comenzaron a matizar -eh que yo soy de bata, no de bota- y se les podía ver cuán jóvenes CSI captando la pequeña muestra en el escenario adecuado que en el laboratorio llevaría hasta el ADN y de ahí a la prueba o identificación certera... Lo cierto es que, ya en pleno siglo XXI, los biólogos seguimos, por supuesto, siendo amigos de plantas y animales, pero ya no sólo somos de bota o de bata. Este es el siglo de la tercera B biológica: de bota, de bata y... ¡de bit!
Como en tantos otros campos los ordenadores han supuesto una revolución también en la biología. Pero esta es más profunda de lo que muchos piensan. Al hablar de biología e informática, quizás el término bioinformática nos viene a la mente. Este es un concepto actualmente muy ligado a la genética y el análisis molecular del ADN. Los ordenadores permitieron en primera instancia mejorar muchísimo la eficiencia del proceso de secuenciación. En consecuencia la cantidad de datos disponibles de ADN de muchas especies, poblaciones e individuos pertenecientes a ellas es actualmente enorme, y sigue creciendo. Cuando hablamos de ADN sabemos que estamos tratando con moléculas de ácido desoxirribonucleico en las que subyace la información genética de los organismos. Una de las informaciones más interesantes que podemos extraer de estas moléculas es la secuencia de bases nitrogenadas (Adenina, Citosina, Guanina y Timina) que caracterizan los diferentes genes y que pueden variar entre especies, poblaciones e incluso, a veces, entre individuos de la misma población. Para representar estas moléculas una vez secuenciadas sólo necesitaríamos papel y lápiz, o un editor de texto en un ordenador, y anotar la secuencia de bases nitrogenadas. Por ejemplo, si hemos secuenciado un fragmento de 5 bases, esta secuencia podría ser AATCA. Pero si en vez de una secuencia de 5 bases hablamos de varias secuencias con muchos miles, cientos de miles o incluso millones de bases entonces los investigadores necesitan utilizar pequeños programas (scripts) para manejarlas de modo que las pueden copiar de unos ficheros a otros, contar el número de bases de un determinado tipo, buscar una secuencia o patrón concreto de bases, comparar unas secuencias con otras, etcétera.
Sin embargo, como ya dijimos, la revolución del bit en biología, es más profunda y por tanto no sólo se trata de los programas y bases de datos para manejar las secuencias. A comienzos del siglo XX, cuando la integración de la genética mendeliana y la teoría evolutiva dio lugar a la disciplina conocida como genética de poblaciones, ello incentivó importantes desarrollos estadísticos para la resolución de los problemas que la nueva disciplina planteaba. Actualmente, en vez de tratar con unos pocos genes tratamos con los genomas enteros. La integración de estos datos genómicos con los modelos poblacionales y ecológicos, plantea nuevos problemas y requiere el desarrollo de nuevos métodos y modelos.
Si estudiamos un proceso biológico que cambia a lo largo del tiempo (dinámico) podemos considerar al menos 3 tipos de modelos. En primer lugar, modelos de tipo biológico (por ejemplo, datos de apareamientos entre dos poblaciones vecinas de una especie), que consisten en información (datos, mediciones,…) que está sometida a ruido. Este ruido proviene del muestreo, errores de medida o simplemente de eventos aleatorios de la naturaleza. El segundo tipo, el modelo matemático, sería una representación formal del modelo biológico en la cual se establecen relaciones entre las características de interés. En estos modelos los números no están sometidos a ruido. Si iteramos el modelo matemático podríamos resolver las variables fundamentales y obtener soluciones predictivas para el mismo. Si el ajuste entre el modelo matemático y el biológico es aceptablemente bueno entonces estas predicciones podrán ser muy útiles. En nuestro ejemplo, podríamos relacionar frecuencias de apareamiento con determinadas características de los individuos y de estos con su composición genética y así estaríamos en posición de conocer la evolución de los fenotipos y genotipos de esas poblaciones tanto hacia adelante como hacia atrás en el tiempo. El problema de estos modelos matemáticos es que, o bien son muy simples (no se ajustan al modelo biológico) o bien no somos capaces de resolverlos cuando el modelo biológico subyacente es muy complejo. Entra aquí en juego el tercer tipo, los modelos de simulación, los cuales resuelven los cálculos de modelos dinámicos matemáticos complejos. Es decir, los modelos de simulación nos permiten plantear y resolver modelos matemáticos complejos de modo que mejore el ajuste respecto a los modelos biológicos.
Y hemos llegado al mensaje clave de esta entrada. La necesidad cada vez más acuciante de desarrollar simulaciones de ecosistemas evolutivos a la luz de los datos genómicos y poblacionales. Si bien es cierto que este tipo de modelos llevan ya un tiempo aplicándose en biología, por ejemplo los llamados modelos basados en individuos (IBM por sus siglas en inglés) que se utilizan en ecología y genética de poblaciones, es ahora cuando emerge una extensión de los mismos combinando información genómica, ecológica y demográfica. Hablamos de modelos eco-genéticos o mejor aún eco-genómicos. Estos nos permitirán avanzar en el conocimiento teórico de aspectos muy diversos e importantes de la ecología y la evolución de las especies a la vez que realizar y contrastar predicciones. Estos modelos se aplican tanto en conservación, como en biomedicina, también para poner a punto estimadores de estadísticos de interés, en estudios de impacto de actividades humanas sobre la diversidad genética, pesquerías, etc., etc.
A pesar de que existen algunas herramientas de simulación eco-genética, pocas veces se ajustan exactamente a nuestras necesidades específicas. Hace tiempo escuché una gran verdad que personalmente he experimentado muchas veces -entendemos realmente algo cuando sabemos programarlo, en caso contrario sólo creemos entenderlo-. Surge entonces la cuestión clave. ¿Qué podemos hacer si somos biólogos y queremos programar nuestros propios bio-simuladores? Como todas las disciplinas que se mueven en la espumeante cresta de la ola del conocimiento, el campo de la simulación biológica, que enmarcaré dentro de la biología computacional, es una ciencia multidisciplinar donde se echa mano de conocimientos de genética, ecología y evolución pero también de programación, estadística y matemáticas. Respecto a los lenguajes de programación, aunque mi favorito es C++ quizás Python (https://www.python.org/) es el que mejor combina simplicidad con buenas prestaciones en términos de estructura y posibilidades del lenguaje. Dispone de funciones predefinidas de utilidad en biología siendo relativamente sencillo realizar los primeros programas útiles. Existen además infinidad de tutoriales, comunidades y foros donde la persona interesada encuentra ayudas de todo tipo para dar los primeros pasos y avanzar en su conocimiento “pythonico”. Así que finalizo este post con una exhortación para los nuevos biólogos de bits, os digo: ¡simulad, simulad malditos!
Para saber más
Haddock, S. H. D., and C. W. Dunn. 2011. Practical computing for biologists. Sinauer, Sunderland, Mass.
Hartmann, A. K. 2015. Big Practical Guide to Computer Simulations. World Scientific.
Hoban, S. 2014. An overview of the utility of population simulation software in molecular ecology. Molecular Ecology 23:2383-2401.
Perkel, J. M. 2015. Programming: Pick up Python. Nature 518:125-126
Como en tantos otros campos los ordenadores han supuesto una revolución también en la biología. Pero esta es más profunda de lo que muchos piensan. Al hablar de biología e informática, quizás el término bioinformática nos viene a la mente. Este es un concepto actualmente muy ligado a la genética y el análisis molecular del ADN. Los ordenadores permitieron en primera instancia mejorar muchísimo la eficiencia del proceso de secuenciación. En consecuencia la cantidad de datos disponibles de ADN de muchas especies, poblaciones e individuos pertenecientes a ellas es actualmente enorme, y sigue creciendo. Cuando hablamos de ADN sabemos que estamos tratando con moléculas de ácido desoxirribonucleico en las que subyace la información genética de los organismos. Una de las informaciones más interesantes que podemos extraer de estas moléculas es la secuencia de bases nitrogenadas (Adenina, Citosina, Guanina y Timina) que caracterizan los diferentes genes y que pueden variar entre especies, poblaciones e incluso, a veces, entre individuos de la misma población. Para representar estas moléculas una vez secuenciadas sólo necesitaríamos papel y lápiz, o un editor de texto en un ordenador, y anotar la secuencia de bases nitrogenadas. Por ejemplo, si hemos secuenciado un fragmento de 5 bases, esta secuencia podría ser AATCA. Pero si en vez de una secuencia de 5 bases hablamos de varias secuencias con muchos miles, cientos de miles o incluso millones de bases entonces los investigadores necesitan utilizar pequeños programas (scripts) para manejarlas de modo que las pueden copiar de unos ficheros a otros, contar el número de bases de un determinado tipo, buscar una secuencia o patrón concreto de bases, comparar unas secuencias con otras, etcétera.
 |
| Créditos de la imagen: Flickr/CC BY-SA 2.0, www.yourgenome.org,
www.jax.org, en.wikipedia.org, www.foodsafetynews.com
|
Si estudiamos un proceso biológico que cambia a lo largo del tiempo (dinámico) podemos considerar al menos 3 tipos de modelos. En primer lugar, modelos de tipo biológico (por ejemplo, datos de apareamientos entre dos poblaciones vecinas de una especie), que consisten en información (datos, mediciones,…) que está sometida a ruido. Este ruido proviene del muestreo, errores de medida o simplemente de eventos aleatorios de la naturaleza. El segundo tipo, el modelo matemático, sería una representación formal del modelo biológico en la cual se establecen relaciones entre las características de interés. En estos modelos los números no están sometidos a ruido. Si iteramos el modelo matemático podríamos resolver las variables fundamentales y obtener soluciones predictivas para el mismo. Si el ajuste entre el modelo matemático y el biológico es aceptablemente bueno entonces estas predicciones podrán ser muy útiles. En nuestro ejemplo, podríamos relacionar frecuencias de apareamiento con determinadas características de los individuos y de estos con su composición genética y así estaríamos en posición de conocer la evolución de los fenotipos y genotipos de esas poblaciones tanto hacia adelante como hacia atrás en el tiempo. El problema de estos modelos matemáticos es que, o bien son muy simples (no se ajustan al modelo biológico) o bien no somos capaces de resolverlos cuando el modelo biológico subyacente es muy complejo. Entra aquí en juego el tercer tipo, los modelos de simulación, los cuales resuelven los cálculos de modelos dinámicos matemáticos complejos. Es decir, los modelos de simulación nos permiten plantear y resolver modelos matemáticos complejos de modo que mejore el ajuste respecto a los modelos biológicos.
Y hemos llegado al mensaje clave de esta entrada. La necesidad cada vez más acuciante de desarrollar simulaciones de ecosistemas evolutivos a la luz de los datos genómicos y poblacionales. Si bien es cierto que este tipo de modelos llevan ya un tiempo aplicándose en biología, por ejemplo los llamados modelos basados en individuos (IBM por sus siglas en inglés) que se utilizan en ecología y genética de poblaciones, es ahora cuando emerge una extensión de los mismos combinando información genómica, ecológica y demográfica. Hablamos de modelos eco-genéticos o mejor aún eco-genómicos. Estos nos permitirán avanzar en el conocimiento teórico de aspectos muy diversos e importantes de la ecología y la evolución de las especies a la vez que realizar y contrastar predicciones. Estos modelos se aplican tanto en conservación, como en biomedicina, también para poner a punto estimadores de estadísticos de interés, en estudios de impacto de actividades humanas sobre la diversidad genética, pesquerías, etc., etc.
A pesar de que existen algunas herramientas de simulación eco-genética, pocas veces se ajustan exactamente a nuestras necesidades específicas. Hace tiempo escuché una gran verdad que personalmente he experimentado muchas veces -entendemos realmente algo cuando sabemos programarlo, en caso contrario sólo creemos entenderlo-. Surge entonces la cuestión clave. ¿Qué podemos hacer si somos biólogos y queremos programar nuestros propios bio-simuladores? Como todas las disciplinas que se mueven en la espumeante cresta de la ola del conocimiento, el campo de la simulación biológica, que enmarcaré dentro de la biología computacional, es una ciencia multidisciplinar donde se echa mano de conocimientos de genética, ecología y evolución pero también de programación, estadística y matemáticas. Respecto a los lenguajes de programación, aunque mi favorito es C++ quizás Python (https://www.python.org/) es el que mejor combina simplicidad con buenas prestaciones en términos de estructura y posibilidades del lenguaje. Dispone de funciones predefinidas de utilidad en biología siendo relativamente sencillo realizar los primeros programas útiles. Existen además infinidad de tutoriales, comunidades y foros donde la persona interesada encuentra ayudas de todo tipo para dar los primeros pasos y avanzar en su conocimiento “pythonico”. Así que finalizo este post con una exhortación para los nuevos biólogos de bits, os digo: ¡simulad, simulad malditos!
Para saber más
Haddock, S. H. D., and C. W. Dunn. 2011. Practical computing for biologists. Sinauer, Sunderland, Mass.
Hartmann, A. K. 2015. Big Practical Guide to Computer Simulations. World Scientific.
Hoban, S. 2014. An overview of the utility of population simulation software in molecular ecology. Molecular Ecology 23:2383-2401.
Perkel, J. M. 2015. Programming: Pick up Python. Nature 518:125-126
martes, 3 de noviembre de 2015
Generación de variables aleatorias con la librería Random de C++11
A partir de la nueva versión del estándar de C++ conocida como C++11 (o en su última actualización C++14 http://isocpp.org/blog/2014/08/we-have-cpp14) la librería estándar aporta un buen número de clases y tipos relacionados con la generación de números aleatorios y sus distribuciones. Esta es una novedad de interés para la biología computacional puesto que gran cantidad de algoritmos y simulaciones se apoyan en el modelo de Monte Carlo cuyo pilar fundamental es la generación de números aleatorios. Antes de esta inclusión era el propio desarrollador quién debía encargarse de construir su propia librería de funciones aleatorias o importar alguna de distribución libre. Algunas librerías de este tipo eran las antiguas genlib.c y randlib.c, o, más recientemente, ALGLIB http://www.alglib.net/ y la librería científica de la GNU (GSL) http://www.gnu.org/software/gsl/ que incluían, además de otras muchas funciones matemáticas, también las relacionadas con la generación de números aleatorios y distribuciones de probabilidad.
Sin embargo, la nueva librería de números aleatorios incluida en C++11 además de algunas de las más importantes distribuciones, aporta también algunos de los mejores y más modernos generadores de aleatoriedad. Sin entrar en detalles que nos llevarían más allá de la intención de este post, podemos aclarar que en la generación de variables aleatorias distinguimos dos aspectos fundamentales. Primero, necesitamos un mecanismo generador de aleatoriedad (random engine) es decir una fuente que genere patrones "impredecibles" de bits. En segundo lugar, necesitamos funciones que aprovechando la susodicha fuente de aleatoriedad generen números según alguna distribución de probabilidad deseada. De manera un poco confusa, a la combinación de un generador de aleatoriedad y una distribución de probabilidad se le llama generador de números aleatorios (random number generator, RNG).
Debe notarse que el generador de aleatoriedad es en realidad un generador de pseudoaleatoriedad (PRNG) es decir se genera aleatoriedad desde el punto de vista estadístico pero basándose en un modelo causal y determinista. Una simple función matemática que para un valor de entrada (semilla) genera una secuencia (muy larga) de valores de salida. Estos valores deben ser estadísticamente uniformes e independientes. La impredecibilidad radica en que estas funciones generan resultados muy distintos para pequeñas variaciones en los estados iniciales. Por tanto, si se cambia la semilla, el estado inicial del mecanismo generador cambia y la secuencia subsiguiente de números aleatorios será muy distinta a la anterior. La ventaja del uso de procesos pseudoaleatorios es que el desarrollador puede, para una misma semilla, generar exactamente la misma secuencia de números una y otra vez, lo cual es extremadamente útil a la hora de comprobar y depurar código, corregir errores, etc.
La clave para conseguir emular la deseada impredecibilidad de muchos fenómenos estocásticos será conseguir semillas del modo menos predecible posible. Una práctica muy común es usar los milisegundos (o nanosegundos) del reloj del procesador. Esta práctica sigue siendo perfectamente válida para la mayoría de aplicaciones en biología computacional si bien ha quedado obsoleta para algunos usos como la seguridad informática (por ejemplo el hackeo del sitio web de tecnología Hacker News se pudo hacer gracias a que usaban números pseudolaeatorios incializados con una semilla basada en el reloj del sistema en milisegundos: https://blog.cloudflare.com/why-randomness-matters/).
Por otro lado, para conseguir generadores realmente aleatorios se pueden usar procesos basados en hardware como alteraciones debidas al ruido electrónico o efectos cuánticos, sin embargo, desde el punto de vista de la simulación, se pierde eficiencia, reproducibilidad y portabilidad. Por tanto es el uso de generadores pseudoaleatorios el que está más extendido, siendo el elemento fundamental de la simulación Monte Carlo y por tanto una parte importantísima de la biología computacional.
El lector interesado puede encontrar un buen tutorial sobre la generación de números aleatorios en C++11 en el siguiente enlace: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3551.pdf.
En lo que sigue haré un breve resumen sobre los generadores aleatorios en C++11 y algunos ejemplos de su uso en simulación de biología de poblaciones.
Mecanismos generadores de aleatoriedad en la librería estándar de C++11
Según el tipo de algoritmo que se utilice para generar los patrones aleatorios de bits distinguimos entre diferentes clases de generadores. Por ejemplo, los generadores congruentes lineales (http://en.wikipedia.org/wiki/Linear_congruential_generator) y uno de los más usados actualmente, el algoritmo "Mersenne Twister" http://en.wikipedia.org/wiki/Mersenne_Twister.
Recordad que necesitábamos dos componentes para generar variables aleatorias: un mecanismo generador y una distribución. En C++11 vienen predefinidos varios mecanismos generadores de pseudoaleatoriedad, lo mismo que diferentes distribuciones. Un resumen de todo ello se puede consultar en el siguiente enlace http://en.cppreference.com/w/cpp/numeric/random.
Personalmente suelo utilizar uno de los generadores más comunes como es el std::mt19937 o su versión de 64 bits std::mt19937_64 aunque para casos puntuales donde se van a generar pocos valores, como por ejemplo dentro de una función de muestreo sin reemplazamiento, uso el generador por defecto std::default_random_engine.
El siguiente paso es generar la semilla (seed) que inicializa la serie. Una opción es
auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();
que guarda en seed el número de tics de reloj desde el tiempo Unix o epoch (https://es.wikipedia.org/wiki/Tiempo_Unix). Como este número depende de la precisión del reloj del sistema (microsegundos en mi máquina) podemos querer definirlo con una precisión distinta. Si por ejemplo queremos milisegundos haríamos:
auto seed =
chrono::duration_cast < chrono::milliseconds
O nanosegundos:
chrono::duration_cast
Por supuesto, hay mejores (y más complicadas) formas de generar la semilla. Algunas sugerencias pueden consultarse aquí: http://stackoverflow.com/questions/24334012/best-way-to-seed-mt19937-64-for-monte-carlo-simulations.
Distribuciones de probabilidad en C++11
Existen varios tipos genéricos de distribuciones tanto discretas como continuas. Algunos ejemplos son: uniforme (incluye uniforme entera y real), Bernoulli (incluye Bernoulli, binomial, binomial negativa y geométrica), Poisson (incluye Poisson, exponencial y gamma entre otras) y normal (incluye normal, log-normal y chi-cuadrado entre otras).
Hay detalles a los que vale la pena estar atento, como por ejemplo el rango de los valores generado. Así, la uniforme entera genera valores en el intervalo cerrado [a,b] mientras que la uniforme real lo hace en el semicerrado [a,b).
Algunos ejemplos de generación y uso de números aleatorios
Lo primero que hacemos al comienzo de nuestro programa es definir el mecanismo generador deseado, la semilla y las distribuciones que necesitamos:
mt19937_64 gener; // este es el mecanismo generador que usarán las distribuciones
// la semilla:
auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();
gener.seed(seed); // inicializa el generador con la semilla
// uniforme entre 0 (incluido) y 1 (no incluido):
uniform_real_distribution
// la normal de media m y desviación típica s:
double m(0.0), s(1.0);
normal_distribution
// Poisson de media avemut=5.2
double avemut(5.2);
poisson_distribution< int>
Después se utilizarán unas u otras según sea necesario.
Ensayos Bernoulli: supervivencia y fertilidad en una población biológica
Un ensayo Bernoulli con probabilidad p es un experimento con sólo dos resultados posibles donde la probabilidad de uno de ellos es p y la del otro 1 – p. Podemos imaginar muchos procesos biológicos como ensayos Bernoulli. Por ejemplo, la probabilidad de sobrevivir dado un valor de eficacia w o la probabilidad de ser progenitor dado un valor de fertilidad f etc etc. Consideremos la eficacia de un individuo en su aspecto de viabilidad es decir, la probabilidad de sobrevivir hasta alcanzar la madurez sexual. Entonces, la supervivencia de un individuo X cualquiera con valor de eficacia w se puede modelar como una Bernoulli. La probabilidad de que ese individuo sobreviva vendrá dada por
if(X.w > uniform(gener)) ... ; // sobrevive
Nótese que si X.w = 1 el individuo siempre sobrevive (recuerda que el intervalo de la uniforme real era semicerrado) pero si X.w es 0 nunca lo hará.
Inicializar los individuos de una población con un fenotipo normalmente distribuido
Esta es una asunción muy común en biología. Una población P en la que estamos considerando un carácter cualquiera: color, tamaño, altura, etc... se distribuye de manera normal para ese carácter. Si queremos generar una población inicial de media m y desviación s hacemos:
for(auto indiv=0; indiv< P.size(); ++indiv)
P[indiv].caracter= normal(gener);
Generar los individuos que van a mutar en una población
Dado que la ocurrencia de una mutación es un suceso de muy baja probabilidad se suele modelar mediante la distribución de Poisson. Si la probabilidad de mutación en un genoma haploide es 0.01 y tenemos N individuos diploides entonces esperamos 0.02N nuevas mutaciones en cada generación. Podemos forzar la ocurrencia de ese número exacto de nuevos mutantes en cada generación o, si queremos ser más realistas, dejar que varíe según una Poisson definida con el parámetro 0.02N. De este último modo el número de mutantes oscilará en cada generación aunque el valor medio a lo largo de las generaciones tenderá a 0.02N.
mutantes=poisson(gener);
Estos son sólo algunos ejemplos bastante sencillos de entre la multitud de posibilidades que los generadores aleatorios nos ofrecen para simular procesos estocásticos en biología.
Un último apunte respecto a la semilla es que, si estamos realizando varias réplicas de un proceso biológico relativamente largo, puede ser conveniente iniciar cada réplica con una nueva semilla de modo que actualizamos el generador:
// actualiza el generador si se quieren realizar nuevas simulaciones independientes
gen.seed(chrono::high_resolution_clock::now().time_since_epoch().count());
Contrastar la corrección de un generador de números pseudoaleatorios
Es siempre conveniente comprobar que hemos implementado correctamente nuestros generadores. Para ello podemos comparar valores esperados de cada distribución, como su media y varianza, con los obtenidos en una secuencia dada de números pseudoaleatorios.
La media de una uniforme U(0,1) es 1/2 y la varianza 1/12. Si generamos una cantidad n de valores bajo esta distribución, la media y varianza de esos valores deberían ir aproximándose a los esperados según aumenta n.
Lo mismo podemos hacer con otras distribuciones. Por ejemplo si un individuo tiene viabilidad 0.7 esa es su probabilidad de supervivencia. Si repetimos el proceso n veces esperamos que el número de veces que sobreviva se aproxime a 0.7n y la varianza de los valores de sobrevivir (1) o morir (0) se aproximará a 0.21n. Igualmente, si inicializamos una población de fenotipos utilizando una normal de media m y desviación s, el fenotipo medio y la desviación en la población obtenida deberían aproximarse a m y s según aumenta el tamaño de la población.
Para una información más detallada sobre estos y otros tipos de test, puede consultarse http://www.euclides.dia.uned.es/aurquia/Files/Simulacion%20Teoria%202008_09.pdf
Bibliografia
- The C++ Standard Library. A Tutorial and Reference. Autor: Nicolai Josuttis. Ed. Addison-Wesley 2012.
- Random Number Generation in C++11.Document #: WG21 N3551. url: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3551.pdf. Date: 2013-03-12.Revises: None.Project: JTC1.22.32 Programming Language C++.Walter E. Brown.
- Matsumoto, M.; Nishimura, T. (1998). "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator". ACM Transactions on Modeling and Computer Simulation 8 (1): 3–30.
- Urquia: "Simulacion, Texto base de teoria", 200.url: http://www.euclides.dia.uned.es/aurquia/Files/Simulacion%20Teoria%202008_09.pdf
sábado, 14 de marzo de 2015
Una bibliografía C++ comentada para biólogos (computacionales)
Libros de C++ y recursos online hay infinitos. Lo que en este post traslado es sólo mi experiencia sobre una serie muy limitada de recursos bibliográficos y/o electrónicos que a mí personalmente me han sido, y siguen siendo, de gran utilidad en mi desempeño profesional. Algunas de mis recomendaciones pueden ser discutibles dado que a día de hoy cualquier persona interesada dispone de ingente cantidad de información vía blogs, manuales online, foros, etc., etc. Pero, ya digo, es mi experiencia y, en mi opinión, algunos buenos libros, cuando se hace el esfuerzo de consultarlos con detenimiento, siguen conteniendo información que a menudo es difícil, sino prácticamente imposible, encontrar online.
Antes de mencionar cualquier material o recurso vale la pena dedicar unas líneas al entorno de desarrollo a utilizar (IDE, por sus siglas en inglés). En mi caso, por diversas circunstancias, he tenido que programar bajo sistemas Windows, Mac y Linux. Por ello se que, un paso importante a la hora de disponer de un método de trabajo eficiente es encontrar un entorno que nos permita generar no sólo código portable sino proyectos de desarrollo intercambiables entre los diferentes sistemas operativos. Hoy día hay herramientas gratuitas y muy potentes, como pueden ser NetBeans y Eclipse, pero personalmente mis preferencias se inclinan, en el caso de C++, por Code::Blocks el cual puede integrarse con diferentes compiladores y que se distribuye gratuitamente bajo licencia GPL v3.0 l. Lo podemos bajar con el compilador GCC ya incluido por defecto. Esta última opción en el caso de Windows implica la instalación (la hace Code::Blocks de manera automática) de MinGW (Minimalist GNU for Windows) es decir, el sistema GNU (sistema operativo tipo-Unix pero libre) mínimo para poder trabajar bajo Windows. MinGW incorpora la colección de compiladores GNU que son más conocidos bajo el acrónimo GCC (GNU Compiler Collection) incluyendo los compiladores de C, C++, ADA y Fortran. Utilizando Code::Blocks, el que esto escribe, ha compilado y ejecutado sin problemas y sin cambiar ni un punto y coma, el mismo proyecto bajo Windows, Linux y Mac. Es cierto que para proyectos de gran complejidad y posiblemente enfocados a sistemas operativos concretos, puede haber otros IDEs más completos y potentes. Sin embargo, para el desarrollo de proyectos de mediana complejidad en biología computacional me parece que Code::Blocks supone un compromiso más que correcto entre facilidad de configuración y optimización, manejo y eficiencia.
Recursos online
Hay infinidad de manuales disponibles gratuitamente online. Para comprobarlo basta teclear en vuestro buscador favorito (no os voy a engañar, el mío es google): "c++ pdf" y obtendréis inmediatamente un listado de manuales en el formato indicado que podréis bajar sin ningún problema a vuestro ordenador. Pero aparte de todos esos pdfs, si disponemos de una conexión a internet, podemos acceder una serie de páginas que dan información y ejemplos rápidamente.
1.- Aquí se puede ver por encima la estructura de la Biblioteca Estándar de C++
http://es.wikipedia.org/wiki/Biblioteca_est%C3%A1ndar_de_C%2B%2B
2.- La siguiente NO es la mejor referencia de C++
http://www.cplusplus.com/reference/
Pues ha sido criticada por contener algunas inexactitudes conceptuales o código no eficiente en algunos de sus ejemplos pero, sin embargo, es una página útil para consulta rápida de métodos asociados a contenedores como vector o métodos de string o algoritmos. Es ligera y carga relativamente rápido. Útil si deseamos obtener información concreta sobre algunas clases de la biblioteca. En mi caso suelo hacer consultas rápidas sobre, por ejemplo, vector:
http://www.cplusplus.com/reference/vector/vector/
3.- El mítico stackoverflow
http://stackoverflow.com/questions/tagged/c%2B%2B
En realidad cuando busco una forma eficiente, o simplemente una forma, de hacer algo, lo que hago es "googlear" el término de búsqueda junto al término c++ y de la lista de resultados consulto en primer lugar los de stackoverflow. Por ejemplo si googleamos: median c++
http://stackoverflow.com/questions/2114797/compute-median-of-values-stored-in-vector-c
4.- El curso de C++ de Zator Systems
http://www.zator.com/Cpp/index.htm
Es un curso con algunos muy buenos ejemplos y explicaciones. Su resumen de operadores de C++ está en mi pestaña de favoritos.
Libros
1.- C++ básico. Autora: Carmen Fernández. Ed. Starbook (http://www.starbook.es/libros/C-BASICO/2245/978-84-936896-8-1)
En mi caso este es un descubrimiento reciente y es el libro que recomiendo actualmente para comenzar en C++ (no en C) desde 0. Tiene, a mi modo de ver, dos ventajas importantes: la primera, es que es un manual breve, sencillo y directo (de muy buen precio) y la segunda es que, a diferencia de muchos otros libros y manuales, menciona dos elementos fundamentales de la Librería Estándar de C++, como son las clases string y vector. Ambas clases facilitan inmensamente la vida del programador con respecto a las herramientas de funcionalidad similar disponibles en C. Otro aspecto interesante de este libro es que los ejemplos del mismo han sido desarrollados utilizando Code::Blocks de modo que dedica un par de páginas a presentar el IDE, explicar cómo se crea un proyecto etc.
2.- Programación en C/C++. Guias prácticas de Anaya multimedia
¿Quién no ha consultado alguna vez en sus inicios alguna de estas guías? Suelen ser manuales de bajo coste, cómodos y útiles para consultas rápidas.
3.- Enciclopedia del lenguaje C++ (C++ estándar). Autor: Fco. Javier Ceballos. Ed. Ra-Ma
Este autor es muy prolífico y tiene varios libros sobre C++ y en otros lenguajes. A mí me gusta este en particular. Es un libro muy completo (y "gordo" 1000 páginas), muy bien explicado y con buenos ejemplos.
4-. Numerical recipes in C++. The Art of Scientific Computing. Autores: Press et al. Ed. Cambridge University Press
Que yo sepa, no hay traducción al español. La edición en inglés que tengo es la del 2002 en su reimpresión de 2005. En cualquier caso, lo digo así de claro, este libro es una de mis biblias. Ello, a pesar de ser consciente de que ha sido bastante criticado, tanto por su pobre diseño desde el punto de vista de la programación orientada a objetos (es poco más que una traducción de C a C++ aprovechando las ventajas de la librería estándar), como desde el campo de la computación numérica, por presentar varios algoritmos que distan de ser los más eficientes en el estado actual del arte de algunos de los muchos campos que toca. Sin embargo, este es un riesgo a asumir, al menos en mi opinión, ante una obra enciclopédica que, en apenas 1000 páginas, introduce conceptos, explica brevemente los métodos y muestra el código completo en C++ de algoritmos que abarcan cuestiones tan dispares, a la par que útiles, como pueden ser la descomposición de Cholesky de una matriz cuadrada (no os riais que sí que sí que es muy útil, cualquiera que haya tenido que generar series de variables correlacionadas lo sabe), la generación de números aleatorios, ordenación, autocorrelación, momentos de una distribución, cálculo de probabilidades de una Chi cuadrado, la tansformada de Fourier, etc, etc. Es un libro que siempre tengo a mano, a partir del cual he programado infinidad de funciones para calcular probabilidades, medias, varianzas, covarianzas... En otras ocasiones me ha servido como una guía de información rápida sobre determinados conceptos y métodos. Una obra de consulta fundamental.
5.- The C++ Standard Library. A Tutorial and Reference. Autor: Nicolai Josuttis. Ed. Addison-Wesley 2012
Este no es un libro para aprender C++ ni algoritmos numéricos. Como su nombre indica es un tutorial y una referencia sobre la Librería Estándar de C++. Personalmente, la mayoría de mis consultas sobre la Librería Estándar las hago utilizando los recursos online, pero este libro incluye algunas explicaciones sobre determinadas características y ventajas de la librería que es difícil (aunque supongo que no imposible) encontrar online. Es un libro claro y profundo, con explicaciones muy buenas, por ejemplo, sobre las ventajas e inconvenientes de los distintos tipos de contenedores (arrays, deques, lists, vectors...), las nuevas características del nuevo estándar C++11 (incluyendo la nueva y magnífica librería de números aleatorios), la biblioteca de algoritmos, etc. Una obra magistral indicada para aquellos que ya saben programar en C++ y desean sacarle todo el partido a la librería estándar.
Y por ahora esto es todo. Hay, por supuesto, muchos más libros y recursos. Pero, como ya dije, estos me parecen útiles para biólogos computacionales o al menos lo han sido en mi caso, a la hora de programar modelos de simulación, algoritmos genéticos y métodos de contraste estadístico en C++. Otros libros de referencia, más avanzados, que quedan lejos (o muy lejos) de la intención de este post y que por tanto sólo me limito a mencionar serían:
Programming. Principles and Practice Using C++ (Stroustrup)
C++ Templates. The Complete Guide (Vandermonde & Josuttis)
Exceptional C++ (Herb Sutter)
Modern C++ Design (Alexandrescu) (muy lejos)
Antes de mencionar cualquier material o recurso vale la pena dedicar unas líneas al entorno de desarrollo a utilizar (IDE, por sus siglas en inglés). En mi caso, por diversas circunstancias, he tenido que programar bajo sistemas Windows, Mac y Linux. Por ello se que, un paso importante a la hora de disponer de un método de trabajo eficiente es encontrar un entorno que nos permita generar no sólo código portable sino proyectos de desarrollo intercambiables entre los diferentes sistemas operativos. Hoy día hay herramientas gratuitas y muy potentes, como pueden ser NetBeans y Eclipse, pero personalmente mis preferencias se inclinan, en el caso de C++, por Code::Blocks el cual puede integrarse con diferentes compiladores y que se distribuye gratuitamente bajo licencia GPL v3.0 l. Lo podemos bajar con el compilador GCC ya incluido por defecto. Esta última opción en el caso de Windows implica la instalación (la hace Code::Blocks de manera automática) de MinGW (Minimalist GNU for Windows) es decir, el sistema GNU (sistema operativo tipo-Unix pero libre) mínimo para poder trabajar bajo Windows. MinGW incorpora la colección de compiladores GNU que son más conocidos bajo el acrónimo GCC (GNU Compiler Collection) incluyendo los compiladores de C, C++, ADA y Fortran. Utilizando Code::Blocks, el que esto escribe, ha compilado y ejecutado sin problemas y sin cambiar ni un punto y coma, el mismo proyecto bajo Windows, Linux y Mac. Es cierto que para proyectos de gran complejidad y posiblemente enfocados a sistemas operativos concretos, puede haber otros IDEs más completos y potentes. Sin embargo, para el desarrollo de proyectos de mediana complejidad en biología computacional me parece que Code::Blocks supone un compromiso más que correcto entre facilidad de configuración y optimización, manejo y eficiencia.
Recursos online
Hay infinidad de manuales disponibles gratuitamente online. Para comprobarlo basta teclear en vuestro buscador favorito (no os voy a engañar, el mío es google): "c++ pdf" y obtendréis inmediatamente un listado de manuales en el formato indicado que podréis bajar sin ningún problema a vuestro ordenador. Pero aparte de todos esos pdfs, si disponemos de una conexión a internet, podemos acceder una serie de páginas que dan información y ejemplos rápidamente.
1.- Aquí se puede ver por encima la estructura de la Biblioteca Estándar de C++
http://es.wikipedia.org/wiki/Biblioteca_est%C3%A1ndar_de_C%2B%2B
2.- La siguiente NO es la mejor referencia de C++
http://www.cplusplus.com/reference/
Pues ha sido criticada por contener algunas inexactitudes conceptuales o código no eficiente en algunos de sus ejemplos pero, sin embargo, es una página útil para consulta rápida de métodos asociados a contenedores como vector o métodos de string o algoritmos. Es ligera y carga relativamente rápido. Útil si deseamos obtener información concreta sobre algunas clases de la biblioteca. En mi caso suelo hacer consultas rápidas sobre, por ejemplo, vector:
http://www.cplusplus.com/reference/vector/vector/
3.- El mítico stackoverflow
http://stackoverflow.com/questions/tagged/c%2B%2B
En realidad cuando busco una forma eficiente, o simplemente una forma, de hacer algo, lo que hago es "googlear" el término de búsqueda junto al término c++ y de la lista de resultados consulto en primer lugar los de stackoverflow. Por ejemplo si googleamos: median c++
http://stackoverflow.com/questions/2114797/compute-median-of-values-stored-in-vector-c
4.- El curso de C++ de Zator Systems
http://www.zator.com/Cpp/index.htm
Es un curso con algunos muy buenos ejemplos y explicaciones. Su resumen de operadores de C++ está en mi pestaña de favoritos.
Libros
1.- C++ básico. Autora: Carmen Fernández. Ed. Starbook (http://www.starbook.es/libros/C-BASICO/2245/978-84-936896-8-1)
En mi caso este es un descubrimiento reciente y es el libro que recomiendo actualmente para comenzar en C++ (no en C) desde 0. Tiene, a mi modo de ver, dos ventajas importantes: la primera, es que es un manual breve, sencillo y directo (de muy buen precio) y la segunda es que, a diferencia de muchos otros libros y manuales, menciona dos elementos fundamentales de la Librería Estándar de C++, como son las clases string y vector. Ambas clases facilitan inmensamente la vida del programador con respecto a las herramientas de funcionalidad similar disponibles en C. Otro aspecto interesante de este libro es que los ejemplos del mismo han sido desarrollados utilizando Code::Blocks de modo que dedica un par de páginas a presentar el IDE, explicar cómo se crea un proyecto etc.
2.- Programación en C/C++. Guias prácticas de Anaya multimedia
¿Quién no ha consultado alguna vez en sus inicios alguna de estas guías? Suelen ser manuales de bajo coste, cómodos y útiles para consultas rápidas.
3.- Enciclopedia del lenguaje C++ (C++ estándar). Autor: Fco. Javier Ceballos. Ed. Ra-Ma
Este autor es muy prolífico y tiene varios libros sobre C++ y en otros lenguajes. A mí me gusta este en particular. Es un libro muy completo (y "gordo" 1000 páginas), muy bien explicado y con buenos ejemplos.
4-. Numerical recipes in C++. The Art of Scientific Computing. Autores: Press et al. Ed. Cambridge University Press
Que yo sepa, no hay traducción al español. La edición en inglés que tengo es la del 2002 en su reimpresión de 2005. En cualquier caso, lo digo así de claro, este libro es una de mis biblias. Ello, a pesar de ser consciente de que ha sido bastante criticado, tanto por su pobre diseño desde el punto de vista de la programación orientada a objetos (es poco más que una traducción de C a C++ aprovechando las ventajas de la librería estándar), como desde el campo de la computación numérica, por presentar varios algoritmos que distan de ser los más eficientes en el estado actual del arte de algunos de los muchos campos que toca. Sin embargo, este es un riesgo a asumir, al menos en mi opinión, ante una obra enciclopédica que, en apenas 1000 páginas, introduce conceptos, explica brevemente los métodos y muestra el código completo en C++ de algoritmos que abarcan cuestiones tan dispares, a la par que útiles, como pueden ser la descomposición de Cholesky de una matriz cuadrada (no os riais que sí que sí que es muy útil, cualquiera que haya tenido que generar series de variables correlacionadas lo sabe), la generación de números aleatorios, ordenación, autocorrelación, momentos de una distribución, cálculo de probabilidades de una Chi cuadrado, la tansformada de Fourier, etc, etc. Es un libro que siempre tengo a mano, a partir del cual he programado infinidad de funciones para calcular probabilidades, medias, varianzas, covarianzas... En otras ocasiones me ha servido como una guía de información rápida sobre determinados conceptos y métodos. Una obra de consulta fundamental.
5.- The C++ Standard Library. A Tutorial and Reference. Autor: Nicolai Josuttis. Ed. Addison-Wesley 2012
Este no es un libro para aprender C++ ni algoritmos numéricos. Como su nombre indica es un tutorial y una referencia sobre la Librería Estándar de C++. Personalmente, la mayoría de mis consultas sobre la Librería Estándar las hago utilizando los recursos online, pero este libro incluye algunas explicaciones sobre determinadas características y ventajas de la librería que es difícil (aunque supongo que no imposible) encontrar online. Es un libro claro y profundo, con explicaciones muy buenas, por ejemplo, sobre las ventajas e inconvenientes de los distintos tipos de contenedores (arrays, deques, lists, vectors...), las nuevas características del nuevo estándar C++11 (incluyendo la nueva y magnífica librería de números aleatorios), la biblioteca de algoritmos, etc. Una obra magistral indicada para aquellos que ya saben programar en C++ y desean sacarle todo el partido a la librería estándar.
Y por ahora esto es todo. Hay, por supuesto, muchos más libros y recursos. Pero, como ya dije, estos me parecen útiles para biólogos computacionales o al menos lo han sido en mi caso, a la hora de programar modelos de simulación, algoritmos genéticos y métodos de contraste estadístico en C++. Otros libros de referencia, más avanzados, que quedan lejos (o muy lejos) de la intención de este post y que por tanto sólo me limito a mencionar serían:
Programming. Principles and Practice Using C++ (Stroustrup)
C++ Templates. The Complete Guide (Vandermonde & Josuttis)
Exceptional C++ (Herb Sutter)
Modern C++ Design (Alexandrescu) (muy lejos)
domingo, 22 de febrero de 2015
Paso de argumentos por valor y por referencia. Aspectos de interés en simulación biológica.
En C++ cuando pasamos los argumentos en un procedimiento o función decimos que el paso puede ser por valor o por referencia. Antes de continuar vamos a clarificar nuestra terminología recordando que en la definición de una función hablamos de los parámetros que recibe la función. Es decir si defino la función suma en C++
void suma(int x, int y, int z){
z=x+y;
}
entonces x, y, z son los parámetros o argumentos formales de la función en los que además se ha indicado el tipo de dato al que corresponden, que en este caso es un entero (int). Cuando dentro de cualquier programa utilizo la función suma pasándole como argumentos distintos variables o valores como por ejemplo
suma(a,b,c); suma(2,3,c);
tanto a,b,c como 2 y 3 son argumentos de la función que se corresponderán dentro del cuerpo de la función a sus parámetros según el orden indicado. Por tanto a se corresponderá con x, b con y, c con z.
Paso por valor
El paso por valor es el mecanismo por omisión en C++. En el paso por valor, las variables que se pasan en una llamada a la función son copiadas a nuevas variables que se corresponden a los parámetros formales y que existirán durante el tiempo que se ejecuta la función. Por ejemplo en el siguiente código
http://pastebin.com/1z8AJvz9
en la línea 6 la variable c tiene valor 7. Cuando se llama a la función suma pasándole como tercer argumento la variable c, lo que ocurre es que la función suma copia el valor de c en una nueva variable que se corresponde con el parámetro z. Todo lo que se le haga a z le ocurre a la copia de c no a c. Lo mismo podemos decir de las variables a y b. Por tanto al ejecutarse suma(a,b,c) la copia de c que llamamos z pasa a ser la suma de las copias de a y b es decir z = x+y = 0+0 = 0. Pero cuando salimos de la función la copia se destruye y la variable c sigue teniendo su valor 7. Debido a que en la llamada suma(a,b,c) no se pasan las variables a, b y c sino una copia de su valor, hablamos de paso por valor.
Paso por referencia
En este caso lo que se le pasa a la función no es una copia del valor de la variable sino una refeencia a esa misma variable. Decimos que el parámetro se convierte en un alias para el argumento, es decir, es la misma variable pero con otro nombre. Como ahora el valor que toman los parámetros de la función no es una copia sino que corresponden a las variables externas a la función los cambios que ocurran dentro de la función tendrán alcance fuera de ella (algo que no ocurría en el paso por valor). Siguiendo con el ejemplo anterior si cambiamos la manera de expresar los parámetros en la definición de suma de modo que en vez de
void suma(int x, int y, int z)
escribimos
void suma(int x, int y, int& z)
lo que le estamos indicando al compilador de C++ es que el parámetro z se pase por referencia. Si volvemos al ejemplo anterior, al llamar a suma con el argumento c ya no se haría una copia de c sino que sería la misma variable c que dentro dela función tendría un segundo nombre (como un apodo o alias) llamado z y por tanto al tomar dentro de la función el valor z = 0+0=0 al salir de la función c valdría 0.
Biología
Veamos ahora un detalle importante relacionado con la eficiencia computacional en relación al paso por valor o por referencia. Imaginemos que tenemos que modelar una población biológica y deseamos programar una simulación donde los individuos se reproducen de modo que una pareja de padres tienen un hijo. Si quisiéramos programar una función reproducción podríamos a priori pensar en algo como
void reproduce(Individuo padre, Individuo madre, Individuo& hijo){..código de la función..}
Donde el tipo de los parámetros es una clase llamada Individuo. Además hemos sido hábiles y nos dimos cuenta de que podemos pasar el padre y la madre por valor puesto que no necesitan cambiar al pasar por la función. Sin embargo el hijo lo vamos a calcular dentro y presumiblemente querremos conservar el valor calculado al salir de la función, por tanto hemos decidido, sabiamente, pasar el parámetro hijo por referencia. Esto hecho así, es correcto y posiblemente obtengamos el comportamiento deseado.
Hay sin embargo un posible problema relacionado con la eficiencia. Imaginemos que nuestra clase Individuo es una abstracción de un organismo diploide, es decir con dos series de cromosomas y que además de los cromosomas le queremos asignar un valor constante para el carácter que sea (por ejemplo color) y otros valores variables como altura, edad etc. Respecto a los cromosomas, vamos a considerar 10 pares de cromosomas cada uno con un millón de posiciones que pueden tener valor 0 o 1. Es decir, el Individuo padre consta de 20 vectores de longitud un millón de posiciones cada uno además de una constante y algunas variables. Si recordamos lo anteriormente dicho para el paso por valor, cuando llamemos a la función reproduce el primer argumento de la llamada se copiará en el primer parámetro (padre). Esta copia exigirá copiar cada una de las millones de posiciones de los vectores así como las constantes, variables y cualquier procedimiento que contengan los objetos de la clase Individuo. Total para usar, sin cambiarlos, unos valores que ya existían fuera de la función. Incurrimos por tanto en un gasto innecesario de memoria y tiempo de computación. En modelos de simulación con miles de individuos y procesos que se repiten durante muchos miles de generaciones el coste puede ser prohibitivo. ¿Cuál es la solución? Pues obviamente es el paso por referencia. La definición adecuada de la función sería
void reproduce(const Individuo& padre, const Individuo& madre, Individuo& hijo){..código de la función..}
Donde todos los argumentos se pasan por referencia de modo que no se realiza la copia de los millones de posiciones cromosómicas sino que simplemente se pasa una referencia al comienzo de la dirección del objeto Individuo correspondiente. El modificador const delante de padre y madre indica que estos parámetros no se modifican dentro de la función e impiden que se haga por error puesto que el compilador se quejará si intentamos cambiarlo en el código.
Otra alternativa que quizás genera un código más claro sería no pasar ningún argumento que se modifique y hacer que la función devuelva un nuevo individuo que será el hijo:
Individuo reproduce(const Individuo& padre, const Individuo& madre){ Indiv hijo; /* codigo....*/ return hijo; }
Sin embargo esta última opción no será posible si es necesario cambiar más de un objeto (si queremos por ejemplo generar dos hijos de cada pareja de padres) y tampoco es lo más razonable cuando queremos modificar vectores. En estos casos el paso por referencia no constante es más conveniente.
Punteros (paso por dirección)
Una última opción para el paso de argumentos en C++ es que el parámetro de la función sea un puntero. En este caso estamos pasando la dirección de la variable y cualquier cambio en el contenido de esa dirección afectará al valor de la variable. En el paso por dirección obtenemos por tanto un comportamiento igual al del paso por referencia. Pero ahora el nombre del parámetro ya no es un alias de la variable sino que corresponde a la dirección de la variable y por tanto hay que desreferenciarla para acceder al contenido de esa dirección (usando el operador *). Dejamos aquí está discusión pues queda más allá de la intención del post. Para saber más sobre referencias y punteros y cuando sería más conveniente el uso de unos u otros se puede consultar el libro de Stroustrup citado en la bibliografía. Personalmente suelo usar siempre referencias puesto que simplifica el código de la función ya que no es necesario desreferenciar la variable (una referencia equivale a un puntero que se desreferencia automáticamente). El paso de punteros puede generar problemas adicionales si ocurren llamadas a la función usando punteros nulos.
Bibliografía
Lenguajes de programación. Principios y práctica. K.C. Louden.
Programming: Principles and Practice Using C++. B. Stroustrup
Object-Oriented Programming in C++. R. Lafore.
Effective C++. S. Meyers.
void suma(int x, int y, int z){
z=x+y;
}
entonces x, y, z son los parámetros o argumentos formales de la función en los que además se ha indicado el tipo de dato al que corresponden, que en este caso es un entero (int). Cuando dentro de cualquier programa utilizo la función suma pasándole como argumentos distintos variables o valores como por ejemplo
suma(a,b,c); suma(2,3,c);
tanto a,b,c como 2 y 3 son argumentos de la función que se corresponderán dentro del cuerpo de la función a sus parámetros según el orden indicado. Por tanto a se corresponderá con x, b con y, c con z.
Paso por valor
El paso por valor es el mecanismo por omisión en C++. En el paso por valor, las variables que se pasan en una llamada a la función son copiadas a nuevas variables que se corresponden a los parámetros formales y que existirán durante el tiempo que se ejecuta la función. Por ejemplo en el siguiente código
http://pastebin.com/1z8AJvz9
en la línea 6 la variable c tiene valor 7. Cuando se llama a la función suma pasándole como tercer argumento la variable c, lo que ocurre es que la función suma copia el valor de c en una nueva variable que se corresponde con el parámetro z. Todo lo que se le haga a z le ocurre a la copia de c no a c. Lo mismo podemos decir de las variables a y b. Por tanto al ejecutarse suma(a,b,c) la copia de c que llamamos z pasa a ser la suma de las copias de a y b es decir z = x+y = 0+0 = 0. Pero cuando salimos de la función la copia se destruye y la variable c sigue teniendo su valor 7. Debido a que en la llamada suma(a,b,c) no se pasan las variables a, b y c sino una copia de su valor, hablamos de paso por valor.
Paso por referencia
En este caso lo que se le pasa a la función no es una copia del valor de la variable sino una refeencia a esa misma variable. Decimos que el parámetro se convierte en un alias para el argumento, es decir, es la misma variable pero con otro nombre. Como ahora el valor que toman los parámetros de la función no es una copia sino que corresponden a las variables externas a la función los cambios que ocurran dentro de la función tendrán alcance fuera de ella (algo que no ocurría en el paso por valor). Siguiendo con el ejemplo anterior si cambiamos la manera de expresar los parámetros en la definición de suma de modo que en vez de
void suma(int x, int y, int z)
escribimos
void suma(int x, int y, int& z)
lo que le estamos indicando al compilador de C++ es que el parámetro z se pase por referencia. Si volvemos al ejemplo anterior, al llamar a suma con el argumento c ya no se haría una copia de c sino que sería la misma variable c que dentro dela función tendría un segundo nombre (como un apodo o alias) llamado z y por tanto al tomar dentro de la función el valor z = 0+0=0 al salir de la función c valdría 0.
Biología
Veamos ahora un detalle importante relacionado con la eficiencia computacional en relación al paso por valor o por referencia. Imaginemos que tenemos que modelar una población biológica y deseamos programar una simulación donde los individuos se reproducen de modo que una pareja de padres tienen un hijo. Si quisiéramos programar una función reproducción podríamos a priori pensar en algo como
void reproduce(Individuo padre, Individuo madre, Individuo& hijo){..código de la función..}
Donde el tipo de los parámetros es una clase llamada Individuo. Además hemos sido hábiles y nos dimos cuenta de que podemos pasar el padre y la madre por valor puesto que no necesitan cambiar al pasar por la función. Sin embargo el hijo lo vamos a calcular dentro y presumiblemente querremos conservar el valor calculado al salir de la función, por tanto hemos decidido, sabiamente, pasar el parámetro hijo por referencia. Esto hecho así, es correcto y posiblemente obtengamos el comportamiento deseado.
Hay sin embargo un posible problema relacionado con la eficiencia. Imaginemos que nuestra clase Individuo es una abstracción de un organismo diploide, es decir con dos series de cromosomas y que además de los cromosomas le queremos asignar un valor constante para el carácter que sea (por ejemplo color) y otros valores variables como altura, edad etc. Respecto a los cromosomas, vamos a considerar 10 pares de cromosomas cada uno con un millón de posiciones que pueden tener valor 0 o 1. Es decir, el Individuo padre consta de 20 vectores de longitud un millón de posiciones cada uno además de una constante y algunas variables. Si recordamos lo anteriormente dicho para el paso por valor, cuando llamemos a la función reproduce el primer argumento de la llamada se copiará en el primer parámetro (padre). Esta copia exigirá copiar cada una de las millones de posiciones de los vectores así como las constantes, variables y cualquier procedimiento que contengan los objetos de la clase Individuo. Total para usar, sin cambiarlos, unos valores que ya existían fuera de la función. Incurrimos por tanto en un gasto innecesario de memoria y tiempo de computación. En modelos de simulación con miles de individuos y procesos que se repiten durante muchos miles de generaciones el coste puede ser prohibitivo. ¿Cuál es la solución? Pues obviamente es el paso por referencia. La definición adecuada de la función sería
void reproduce(const Individuo& padre, const Individuo& madre, Individuo& hijo){..código de la función..}
Donde todos los argumentos se pasan por referencia de modo que no se realiza la copia de los millones de posiciones cromosómicas sino que simplemente se pasa una referencia al comienzo de la dirección del objeto Individuo correspondiente. El modificador const delante de padre y madre indica que estos parámetros no se modifican dentro de la función e impiden que se haga por error puesto que el compilador se quejará si intentamos cambiarlo en el código.
Otra alternativa que quizás genera un código más claro sería no pasar ningún argumento que se modifique y hacer que la función devuelva un nuevo individuo que será el hijo:
Individuo reproduce(const Individuo& padre, const Individuo& madre){ Indiv hijo; /* codigo....*/ return hijo; }
Sin embargo esta última opción no será posible si es necesario cambiar más de un objeto (si queremos por ejemplo generar dos hijos de cada pareja de padres) y tampoco es lo más razonable cuando queremos modificar vectores. En estos casos el paso por referencia no constante es más conveniente.
Punteros (paso por dirección)
Una última opción para el paso de argumentos en C++ es que el parámetro de la función sea un puntero. En este caso estamos pasando la dirección de la variable y cualquier cambio en el contenido de esa dirección afectará al valor de la variable. En el paso por dirección obtenemos por tanto un comportamiento igual al del paso por referencia. Pero ahora el nombre del parámetro ya no es un alias de la variable sino que corresponde a la dirección de la variable y por tanto hay que desreferenciarla para acceder al contenido de esa dirección (usando el operador *). Dejamos aquí está discusión pues queda más allá de la intención del post. Para saber más sobre referencias y punteros y cuando sería más conveniente el uso de unos u otros se puede consultar el libro de Stroustrup citado en la bibliografía. Personalmente suelo usar siempre referencias puesto que simplifica el código de la función ya que no es necesario desreferenciar la variable (una referencia equivale a un puntero que se desreferencia automáticamente). El paso de punteros puede generar problemas adicionales si ocurren llamadas a la función usando punteros nulos.
Bibliografía
Lenguajes de programación. Principios y práctica. K.C. Louden.
Programming: Principles and Practice Using C++. B. Stroustrup
Object-Oriented Programming in C++. R. Lafore.
Effective C++. S. Meyers.
Hola, bienvenido al blog de BiosDev, un biólogo contradictorio; apasionado de la naturaleza, la montaña y especialmente el mar. Pero que vive entre ordenadores, fascinado por los modelos matemáticos y la programación. En fin, un lío.
Actualmente, me dedico a la docencia e investigación en Genética y Biología Computacional. Concretamente, la Biología Computacional es un campo de estudio que considero especial. Por un lado es un ámbito multidisciplinar donde entran en juego conocimientos de biología, estadística-matemática y programación. Por otro lado es un ámbito con una gran demanda actual, y previsiblemente futura, de profesionales. Sin embargo, y aunque parezca extraño, existe en España muy poca gente formada o formándose en biología computacional, bioinformática o disciplinas afines. Las causas de esta situación no las voy a discutir aquí aunque podrían ser el tema de una futura entrada en este blog. La idea del mismo nació ya a finales de 2008 pero diversas vicisitudes personales y profesionales la fueron diluyendo. Mi intención es retomar ahora la idea original e ir compartiendo inquietudes, conocimientos e ideas que han ido surgiendo a lo largo de 15 años de programación de modelos y simulaciones en biología, especialmente en genética evolutiva.
Es también una suerte, y un privilegio, el haberme formado en dos campos a priori bien diferentes como son la Biología y la Informática. En realidad, tras un poco de estudio, se perciben muchas conexiones entre ambos mundos. Quizás por eso, siempre trato de enfocar mi trabajo desde una doble perspectiva. Por un lado, el problema biológico que subyace en un programa o modelo determinado y por otro, el lenguaje, la eficiencia y la técnica de programación utilizada. Pero ya lo dice el refrán, "quién mucho abarca, poco aprieta" y seguramente muchos programadores podrán plantear soluciones más eficientes a cuestiones de código que aquí se expongan. Y lo mismo ocurrirá desde el ámbito de la biología. Por eso toda crítica constructiva será siempre muy bien recibida.
Actualmente, me dedico a la docencia e investigación en Genética y Biología Computacional. Concretamente, la Biología Computacional es un campo de estudio que considero especial. Por un lado es un ámbito multidisciplinar donde entran en juego conocimientos de biología, estadística-matemática y programación. Por otro lado es un ámbito con una gran demanda actual, y previsiblemente futura, de profesionales. Sin embargo, y aunque parezca extraño, existe en España muy poca gente formada o formándose en biología computacional, bioinformática o disciplinas afines. Las causas de esta situación no las voy a discutir aquí aunque podrían ser el tema de una futura entrada en este blog. La idea del mismo nació ya a finales de 2008 pero diversas vicisitudes personales y profesionales la fueron diluyendo. Mi intención es retomar ahora la idea original e ir compartiendo inquietudes, conocimientos e ideas que han ido surgiendo a lo largo de 15 años de programación de modelos y simulaciones en biología, especialmente en genética evolutiva.
Es también una suerte, y un privilegio, el haberme formado en dos campos a priori bien diferentes como son la Biología y la Informática. En realidad, tras un poco de estudio, se perciben muchas conexiones entre ambos mundos. Quizás por eso, siempre trato de enfocar mi trabajo desde una doble perspectiva. Por un lado, el problema biológico que subyace en un programa o modelo determinado y por otro, el lenguaje, la eficiencia y la técnica de programación utilizada. Pero ya lo dice el refrán, "quién mucho abarca, poco aprieta" y seguramente muchos programadores podrán plantear soluciones más eficientes a cuestiones de código que aquí se expongan. Y lo mismo ocurrirá desde el ámbito de la biología. Por eso toda crítica constructiva será siempre muy bien recibida.
jueves, 25 de octubre de 2012
Biología (real e imaginaria) del número e
Siempre me fascinaron los números que se ocultan tras un seudónimo, números cuyo largo o extraño "nombre numérico" apenas podemos pronunciar y nos vemos obligados a renombrarlos con identificadores más próximos a nuestro lenguaje habitual. O bien, números que son tan importantes que los distinguimos con una denominación alternativa. Entre los más famosos tenemos, por supuesto, a π, acompañado del número e y de la unidad imaginaria i.
Aún recuerdo mi asombro, y mi deleite, al descubrir que estos tres podían relacionarse por medio del uno y el cero. El uno, el único número sólo divisible por sí mismo (y por el uno, je) que no es primo. Su primicidad sacrificada en aras del teorema fundamental de la aritmética. El cero, que separa lo positivo de lo negativo. Esta relación debida a Euler es la bella identidad:
Entre estos, el número e cobra una especial relevancia para el alumno de biología interesado en los procesos y modelos de la ecología y la genética de poblaciones. Eso pensé cuando e comenzó a cruzarse en mi camino con más frecuencia que el resto. Con exagerada, incluso cansina, frecuencia, desde el punto de vista del estudiante que sólo veía una letra sin sentido, incomprensible y reiteradamente asociada a logaritmos, derivadas e integrales. Para empeorar las cosas, a nadie parecía perturbar y tampoco nadie se molestaba en explicar la omnipresencia de la dichosa letra-número. Ni siquiera los genéticos, la familia más matemática de la tribu de biólogos, parecían tener mucho que decir sobre el tema. Soy de los que necesitan explicaciones para poder seguir avanzando. La siguiente entrada narra en unas breves pinceladas la historia de mi búsqueda para entenderme con e.
En los comienzos de mi, inacabado viaje hacia el número e, fue particularmente clarificadora lo que llamaré la analogía del bicho virtual logarítmico. Me explico. Imaginemos que simulamos en nuestro ordenador una entidad biológica cualquiera (hongo, bicho o planta virtual) muy simple y con sólo una característica cambiante en el tiempo, a saber, que aumenta de tamaño en una sola dimensión. El tamaño inicial de nuestro e-bicho será 1 unidad de tamaño virtual y programamos la regla de que su crecimiento en unidades por año es, en todo momento, equivalente a su tamaño actual. Así, inicialmente en tiempo simulación 0, crece a razón de 1 unidad anual; cuando mida 1,235 unidades, crecerá 1,235 unidades anuales, etc, etc.
Para saber cuánto va a medir nuestro e-bicho tras un breve lapso de tiempo, dividimos el año en muchas partes muy pequeñas, cada una de las cuales denominamos instante. Por tanto, para N suficientemente grande, 1/N corresponderá a un instante. Puesto que, por la regla antes definida, la ganancia anual en el instante primero es su tamaño inicial que vale 1, la ganancia instantánea (no anual) será una enésima parte de ese 1, esto es, 1/N.
Por tanto el nuevo tamaño tras el primer instante será 1 + 1/N.
En el segundo instante ganará la enésima parte de su tamaño actual, por tanto tendremos
(1 + 1/N ) + (1 + 1/N) / N = 1 + 2/N + 1/N2 = (1 + 1/N)2
y así sucesivamente para los siguientes instantes. Como hemos dividido el año en N instantes tendremos finalmente que el tamaño de nuestro e-bicho tras un año (= N instantes) será
(1 + 1/N)N
No está mal, pero ahora surge un problema práctico. ¿Cuál es el tamaño de nuestro e-bicho? ¡El público ruge pidiendo un número! Ahora bien, ¿quién conoce el número N de instantes en un año? esta no es una pregunta fácil. A lo sumo, uno puede indicar que son muchos y por tanto N será grande, pero cuánto es grande. En este punto y en el contexto de crisis actual seguro que a alguno el presente año 2012 se le está haciendo especialmente largo, como si viniera con un extra de instantes... Perdón, estoy divagando...
Desde un punto de vista formal, el número infinitamente grande de instantes infinitamente pequeños utiliza el recurso del límite, que nos dice que N tiende a tan grande que no termina (equivalentemente pero en términos del análisis no estándar, N es un número ilimitado; o el instante, un infinitesimal). En cualquier caso, y desde un punto de vista más práctico, podemos dar un valor concreto a N y asumir como aproximación bastante tosca el valor de 31,5 millones de instantes para un año. Es decir, nuestro instante vale aproximadamente un segundo. Rápidamente podemos calcular el tamaño de nuestro e-bicho que será:
(1 + 0,032 x 10-6) 31500000 = 2.7182817853117158805093639805235.
Si nuestro instante valiera un día tendríamos
(1 + 1/365)365 = 2.7145674820218743031938863066851
y si valiera un mes (1 + 1/12)12 = 2.6130352902246781602995330443549.
Que no son sino distintas aproximaciones al valor constante del número e que es precisamente el límite de (1 + 1/N)N para N cuando N tiende a infinito.
Jacob Bernoulli fue el primero que calculó este límite. Inicialmente se llamaba b y fue Euler quien, al desarrollar la función exponencial mediante series infinitas, introdujo la notación con la letra e que ahora conocemos. Nuestro e-bicho tras un año creciendo, mide e unidades. Se deja como ejercicio para el sufrido lector, el comprobar, aplicando la misma regla antes indicada, que nuestro e-bicho en el segundo año medirá e2 unidades y así sucesivamente. Es decir, crece exponencialmente con el número de años y la base es e, de ahí que digamos que el e-bicho sea logarítmico puesto que podemos medir su crecimiento anual usando los exponentes a los que está elevado e, es decir, los logaritmos.
He aquí nuestra primera pista respecto a la pregunta ¿porqué aparece e en tantos modelos de la biología? Vemos que, si medimos un crecimiento o cambio por unidad de variable, y este crecimiento vale exactamente el valor de la entidad que cambia (p.ej. crece lo mismo que pesa o mide, etc), nos acercaremos al valor del número e, tanto más, cuanto más pequeño sea el lapso de tiempo en que medimos el cambio.
¿Porqué digo que lo anterior es una pista? pues porque en Biología es muy común hallar situaciones donde los procesos de cambio dependen o son regulados en función del propio valor de la entidad que cambia, bien en sentido positivo o negativo. Por ejemplo, el crecimiento de una población será mayor cuanto mayor sea la población, al menos en tanto haya recursos disponibles. La competencia entre especies, el control de las reacciones en el metabolismo, la expansión de una epidemia, etc, etc.
El cambio de ex según varía x es ex
La derivada de una función f es su incremento por unidad de variable x. Desde un punto de vista analítico, la derivada de una función f en un punto a de su dominio es el límite en el que x tiende a a de (f(x) – f(a)) / (x –a) que también podemos expresar como el límite cuando h tiende a 0 de (f(a+h) – f(a))/h donde hemos sustituido x = a+h, Quedémonos con esta última definición.
Vamos a considerar la función f(x) = ex. Por tanto f(a) = ea. El incremento infinitesimal h utilizando nuestra notación previa será 1/N donde N tiende a infinito. Aplicando la definición dada arriba, la derivada de la función ex en el punto a será el límite cuando N tiende a infinito de
(ea+1/N – ea) / (1/N) = (ea e1/N – ea ) / (1/N) = (ea(1 + 1/N) – ea ) / (1/N) = ea
(nótese que hemos sustituido e1/N = 1 + 1/N por definición de e).
Sustituyendo a por x tenemos finalmente que
d(ex) / d(x) = ex,
es decir, la derivada de ex es ella misma, algo que ya sabíamos, pues la función ex, es "famosa" por eso mismo.
En el caso de la Biología es más común que el cambio o la evolución de la función no sea igual sino proporcional a su propio valor, que entonces sería ekx y aplicando de nuevo la definición de derivada llegamos a:
d(eka) = Keka.
Una vez, hace ya demasiados años, estaba yo en quinto de carrera y además en clase y para colmo, atendiendo a las explicaciones del profesor. Se hablaba de genética y había un silencio atento, pues el profe era bueno. Ocurrió de súbito, así ocurren siempre estas cosas en la vida, cuando no te lo esperas. Lo que aquel hombre había escrito en la pizarra se parecía a una ecuación. Pero había algo extraño, algo que no terminaba de encajar. Una ecuación, pensaba yo por aquellos tiempos, tiene que llevar su parte de y y su parte de x…, esas partes pueden ser más o menos complicadas pero ahí se acaba la cosa, o debería. A lo sumo, y ya en territorio comanche, si es derivada pues la y prima y las x. O si, peor aún, es integral pues el simbolito "ese" y dentro la parte de x con la dx. Pero no, aquel alien con forma de ecuación no encajaba en ninguna categoría conocida. En aquello coexistían la y con la x y con la derivada de y, todo revuelto.
Mi sorpresa fue inmensa cuando, sin apenas atisbo de nervios por lo que había osado plantarnos en la pizarra, nos comentó, creo recordar que con cierta sorna, que aquello no era sino era una ecuación diferencial. No había peligro, dijo, pues no era más que una ecuación que llevaba una derivada dentro. Recuerdo perfectamente aquella sensación agridulce de pensar por un lado, que efectivamente no parecía para tanto la cosa, pero por otro, todas las alarmas saltando al pensar que, si ya las derivadas se bastaban solas para complicarle la vida a uno, nada bueno podía esperarse tras el hecho de introducirlas dentro de otra ecuación. No recuerdo qué ecuación era aquella pero sí que, ese simple comentario, desnudando de parte de su oscuro misterio a las temidas ecuaciones diferenciales, me sirvió para posteriormente profundizar un poco en su estudio, aunque siempre con una “mirada” biológica.
Supe así, que había ecuaciones diferenciales de primer orden, de orden n, ecuaciones no lineales, estabilidad, ecuaciones en derivadas parciales, etc. La enfermedad contraída fue tan grave como para llegar a impartir un curso de postgrado de ecuaciones diferenciales para biólogos (dos alumnas, a día de hoy, ambas retiradas de la investigación, quiero pensar que por motivos totalmente ajenos a aquel curso). Pasó el tiempo, y hondo fue mi pesar el día que comprendí que nada de todo esto era suficiente para describir los modelos más interesantes, y todavía básicamente sencillos, en genética de poblaciones donde lo necesario era, y es, otro tipo de ecuaciones, llamadas ecuaciones diferenciales estocásticas. Pero esa es otra historia...
Si has llegado hasta aquí, querido lector, cuentas con todo mi aprecio y también con algo de mi asombro. Por tu parte, acaso estás preguntándote con cierta melancolía, sobre los caminos que te han llevado hasta una entrada de blog como esta. ¿Qué puede haber fallado? pensarás... Por otro lado, aparentemente hemos abandonado a e en los últimos párrafos. Pero, como se suele decir, no hay fallo, está aquí mismo e inmediatamente retorna al escenario en, lo prometo, el último acto.
Veamos. Un modelo matemático en que aparece una expresión relacionando una función f(x), con su primera o sucesivas derivadas, es una ecuación diferencial ordinaria. Es decir, conocida la relación entre una función f, su variable x y sus derivadas, lo que deseamos es encontrar la función incognita que satisfaga esa relación. Resolver la ecuación diferencial consiste en integrarla. Recordemos aquí, que la función que relaciona magnitudes y crecimientos (aquella que es proporcional a su derivada) es ekx y por tanto es de esperar que el número e juegue un importante papel en la interpretación de las ecuaciones diferenciales, como así es.
Un tipo de modelos sencillos y muy comunes en biología de poblaciones, ecología, epidemiología, etc. toma como variable independiente el tiempo t de modo que la función y(t) representa a los individuos de la población, por ejemplo su número, siendo la derivada el cambio del tamaño de población en el tiempo. He aquí la primera respuesta muy genérica a la pregunta de porqué aparecía el número e por todas partes en biología en general y de poblaciones en particular. El tipo de modelo que relaciona cantidad de individuos o moléculas o virus, con su evolución (cambio) en el tiempo son las ecuaciones diferenciales cuya solución involucra casi por definición al número e.
Parece que fuel el economista Tomas Malthus (1766-1834) el que primero observó que muchas poblaciones biológicas crecen de modo proporcional a la población. En genética de poblaciones es común usar la letra N para la variable del número de individuos de una población (por tanto aquí ya no corresponde al número de instantes). Por ello, la función y(t) la notaremos entonces como N(t), o más bien como N para simplificar la notación, y la razón dN/dt de cambio de la población (esto es, su derivada) es proporcional al valor actual de N. Por tanto
dN/dt = rN.
La constante r es la razón de crecimiento (si es positiva) o disminución (si negativa). Si la población está creciendo asumimos r mayor que 0 e imponemos la condición inicial de que en tiempo 0 el tamaño poblacional es N0 es decir N(0) = N0 .
La solución de esta ecuación la obtenemos si, agrupamos las N a un lado de la ecuación y la t en el otro, dN/N = rdt, e integramos ambos lados de modo que
log N = rt + C donde log es el logaritmo en base e y C es una constante. Por tanto
N = ert+C = C’ert. Y como en t=0, N = N0 tenemos que N0 = C’ y finalmente
N = N0 ert.
Que es el modelo Malthusiano de crecimiento exponencial antecesor del modelo conocido como crecimiento logístico, que surge al sustituir la constante r de razón de crecimiento por una función del tamaño poblacional N, r = f(N). El sentido de hacer esto es que, si bien en muchas poblaciones el modelo exponencial es correcto durante un tiempo definido, es obvio que los recursos son limitados y mantener el crecimiento exponencial durante tiempo ilimitado no es factible.
El modelo logístico surge a partir del Malthusiano, para tener en cuenta que la razón de crecimiento depende de la población. Para incorporar un límite al tamaño poblacional según los recursos disponibles, basta elegir un f(N) que valga aproximadamente r para N pequeño, pero de modo que f(N) decrezca según aumenta N y que a partir de un N suficientemente grande, f(N) se haga negativo. La función más sencilla que cumple estos requisitos es
f(N) = r –aN donde a es positivo.
Si sustituimos esta función por r en la ecuación de Malthus obtenemos
dN/dt = Nr(1 – N/K) donde se ha hecho el cambio a = r /K.
Esta es la ecuación logística, introducida por el matemático belga P.F. Verlhulst como modelo para el crecimiento de la población humana, con la poca fortuna de disponer de datos de censo inadecuados con lo que su modelo cayó en el olvido. Lotka volvió a obtener la ecuación en 1925 y la llamó ley del crecimiento poblacional. Al parámetro r se le conoce como razón de crecimiento intrínseco y K es la cota superior a la que tiende el tamaño poblacional N tras un tiempo infinito. Por tanto a K se le conoce como nivel de saturación o, en ecología y biología de poblaciones, como la capacidad de carga de la población (carrying capacity). La solución de la ecuación incorpora, como no, nuestra ya conocida ert.
La importancia de esta ecuación ha sido tal, que en ecología se distinguen dos estrategias de supervivencia conocidas como estrategia-r y estrategia-K en relación al conflicto entre número de descendencia versus cuidados parentales invertidos en la descendencia. La estrategia r fomenta muchos descendientes en los que se invierte escaso cuidado parental (corresponde a valores elevados de la tasa de crecimiento r y bajos de K) mientras que la estrategia K implica poca descendencia en la que se invierte mucho cuidado y por tanto tienen mayor probabilidad de supervivencia (baja r, elevada K). Si bien es necesario aclarar que a día de hoy, este modelo K-r ha sido superado por otros más completos, que incluyen distintos componentes de la historia vital de las especies teniendo en cuanta, por ejemplo, la estructura de edades.
Obviamente, queda mucho por decir sobre el papel que el número e tiene en la formalización y modelado biológicos, algo que, con cierto desenfado, he llamado biología del número e. El rastro de e lo podemos seguir a través de distribuciones probabilisticas de gran importancia en Biología como la Normal, la Poisson, etc, así como también en procesos estocásticos. Es clave en el estudio de la evolución de las secuencias de ADN, donde aparece la ecuación P = eQt donde P y Q son matrices, que permite calcular las probabilidades de transición de unas secuencias de información genética a otras. También hallaremos a e en la dinámica de sistemas ecológicos, etc etc... En cualquier caso, paremos aquí. Aunque yo no esté satisfecho, y tú ¿lector? estés harto... A veces me pregunto ¿porqué e? ¿porqué precisamente entre el 2 y el 3?.. Pero dado que llevo dos días escribiendo esta entrada sin ton ni son, será mejor para todos, terminar aquí ahora.
PD: sí, ya sé que podía haber puesto alguna figura!
Bessiere, G. "El cálculo integral fácil y atractivo" trad por A.G.T. Ed. Naval, 1946.
Boyce, WE, DiPrima, RC. "Ecuaciones diferenciales y problemas con valores en la frontera". Ed. Limusa, 1998.
Dunham, W. "Euler. El maestro de todos los matemáticos". Ed. NIVOLA, 2000.
Roughgarden, J. "Theory of population genetics and evolutionary ecology. An introduction". Ed.Prentice-Hall, 1996.
“Esta entrada participa en la edición 3.1415926 del Carnaval de Matemáticas, alojado en el blog Series divergentes.”
Aún recuerdo mi asombro, y mi deleite, al descubrir que estos tres podían relacionarse por medio del uno y el cero. El uno, el único número sólo divisible por sí mismo (y por el uno, je) que no es primo. Su primicidad sacrificada en aras del teorema fundamental de la aritmética. El cero, que separa lo positivo de lo negativo. Esta relación debida a Euler es la bella identidad:
eiπ + 1 = 0
El e-bicho logaritmico (modificado de G Bessiere, trad por A.G.T Ed. Naval, 1946)
En los comienzos de mi, inacabado viaje hacia el número e, fue particularmente clarificadora lo que llamaré la analogía del bicho virtual logarítmico. Me explico. Imaginemos que simulamos en nuestro ordenador una entidad biológica cualquiera (hongo, bicho o planta virtual) muy simple y con sólo una característica cambiante en el tiempo, a saber, que aumenta de tamaño en una sola dimensión. El tamaño inicial de nuestro e-bicho será 1 unidad de tamaño virtual y programamos la regla de que su crecimiento en unidades por año es, en todo momento, equivalente a su tamaño actual. Así, inicialmente en tiempo simulación 0, crece a razón de 1 unidad anual; cuando mida 1,235 unidades, crecerá 1,235 unidades anuales, etc, etc.
Para saber cuánto va a medir nuestro e-bicho tras un breve lapso de tiempo, dividimos el año en muchas partes muy pequeñas, cada una de las cuales denominamos instante. Por tanto, para N suficientemente grande, 1/N corresponderá a un instante. Puesto que, por la regla antes definida, la ganancia anual en el instante primero es su tamaño inicial que vale 1, la ganancia instantánea (no anual) será una enésima parte de ese 1, esto es, 1/N.
Por tanto el nuevo tamaño tras el primer instante será 1 + 1/N.
En el segundo instante ganará la enésima parte de su tamaño actual, por tanto tendremos
(1 + 1/N ) + (1 + 1/N) / N = 1 + 2/N + 1/N2 = (1 + 1/N)2
y así sucesivamente para los siguientes instantes. Como hemos dividido el año en N instantes tendremos finalmente que el tamaño de nuestro e-bicho tras un año (= N instantes) será
(1 + 1/N)N
No está mal, pero ahora surge un problema práctico. ¿Cuál es el tamaño de nuestro e-bicho? ¡El público ruge pidiendo un número! Ahora bien, ¿quién conoce el número N de instantes en un año? esta no es una pregunta fácil. A lo sumo, uno puede indicar que son muchos y por tanto N será grande, pero cuánto es grande. En este punto y en el contexto de crisis actual seguro que a alguno el presente año 2012 se le está haciendo especialmente largo, como si viniera con un extra de instantes... Perdón, estoy divagando...
Desde un punto de vista formal, el número infinitamente grande de instantes infinitamente pequeños utiliza el recurso del límite, que nos dice que N tiende a tan grande que no termina (equivalentemente pero en términos del análisis no estándar, N es un número ilimitado; o el instante, un infinitesimal). En cualquier caso, y desde un punto de vista más práctico, podemos dar un valor concreto a N y asumir como aproximación bastante tosca el valor de 31,5 millones de instantes para un año. Es decir, nuestro instante vale aproximadamente un segundo. Rápidamente podemos calcular el tamaño de nuestro e-bicho que será:
(1 + 0,032 x 10-6) 31500000 = 2.7182817853117158805093639805235.
Si nuestro instante valiera un día tendríamos
(1 + 1/365)365 = 2.7145674820218743031938863066851
y si valiera un mes (1 + 1/12)12 = 2.6130352902246781602995330443549.
Que no son sino distintas aproximaciones al valor constante del número e que es precisamente el límite de (1 + 1/N)N para N cuando N tiende a infinito.
Jacob Bernoulli fue el primero que calculó este límite. Inicialmente se llamaba b y fue Euler quien, al desarrollar la función exponencial mediante series infinitas, introdujo la notación con la letra e que ahora conocemos. Nuestro e-bicho tras un año creciendo, mide e unidades. Se deja como ejercicio para el sufrido lector, el comprobar, aplicando la misma regla antes indicada, que nuestro e-bicho en el segundo año medirá e2 unidades y así sucesivamente. Es decir, crece exponencialmente con el número de años y la base es e, de ahí que digamos que el e-bicho sea logarítmico puesto que podemos medir su crecimiento anual usando los exponentes a los que está elevado e, es decir, los logaritmos.
Primeros atisbos de respuesta a la pregunta de un biólogo enfrentado a e
He aquí nuestra primera pista respecto a la pregunta ¿porqué aparece e en tantos modelos de la biología? Vemos que, si medimos un crecimiento o cambio por unidad de variable, y este crecimiento vale exactamente el valor de la entidad que cambia (p.ej. crece lo mismo que pesa o mide, etc), nos acercaremos al valor del número e, tanto más, cuanto más pequeño sea el lapso de tiempo en que medimos el cambio.
¿Porqué digo que lo anterior es una pista? pues porque en Biología es muy común hallar situaciones donde los procesos de cambio dependen o son regulados en función del propio valor de la entidad que cambia, bien en sentido positivo o negativo. Por ejemplo, el crecimiento de una población será mayor cuanto mayor sea la población, al menos en tanto haya recursos disponibles. La competencia entre especies, el control de las reacciones en el metabolismo, la expansión de una epidemia, etc, etc.
El cambio de ex según varía x es ex
La derivada de una función f es su incremento por unidad de variable x. Desde un punto de vista analítico, la derivada de una función f en un punto a de su dominio es el límite en el que x tiende a a de (f(x) – f(a)) / (x –a) que también podemos expresar como el límite cuando h tiende a 0 de (f(a+h) – f(a))/h donde hemos sustituido x = a+h, Quedémonos con esta última definición.
Vamos a considerar la función f(x) = ex. Por tanto f(a) = ea. El incremento infinitesimal h utilizando nuestra notación previa será 1/N donde N tiende a infinito. Aplicando la definición dada arriba, la derivada de la función ex en el punto a será el límite cuando N tiende a infinito de
(ea+1/N – ea) / (1/N) = (ea e1/N – ea ) / (1/N) = (ea(1 + 1/N) – ea ) / (1/N) = ea
(nótese que hemos sustituido e1/N = 1 + 1/N por definición de e).
Sustituyendo a por x tenemos finalmente que
d(ex) / d(x) = ex,
es decir, la derivada de ex es ella misma, algo que ya sabíamos, pues la función ex, es "famosa" por eso mismo.
En el caso de la Biología es más común que el cambio o la evolución de la función no sea igual sino proporcional a su propio valor, que entonces sería ekx y aplicando de nuevo la definición de derivada llegamos a:
d(eka) = Keka.
Mi primera ecuación diferencial
Una vez, hace ya demasiados años, estaba yo en quinto de carrera y además en clase y para colmo, atendiendo a las explicaciones del profesor. Se hablaba de genética y había un silencio atento, pues el profe era bueno. Ocurrió de súbito, así ocurren siempre estas cosas en la vida, cuando no te lo esperas. Lo que aquel hombre había escrito en la pizarra se parecía a una ecuación. Pero había algo extraño, algo que no terminaba de encajar. Una ecuación, pensaba yo por aquellos tiempos, tiene que llevar su parte de y y su parte de x…, esas partes pueden ser más o menos complicadas pero ahí se acaba la cosa, o debería. A lo sumo, y ya en territorio comanche, si es derivada pues la y prima y las x. O si, peor aún, es integral pues el simbolito "ese" y dentro la parte de x con la dx. Pero no, aquel alien con forma de ecuación no encajaba en ninguna categoría conocida. En aquello coexistían la y con la x y con la derivada de y, todo revuelto.
Mi sorpresa fue inmensa cuando, sin apenas atisbo de nervios por lo que había osado plantarnos en la pizarra, nos comentó, creo recordar que con cierta sorna, que aquello no era sino era una ecuación diferencial. No había peligro, dijo, pues no era más que una ecuación que llevaba una derivada dentro. Recuerdo perfectamente aquella sensación agridulce de pensar por un lado, que efectivamente no parecía para tanto la cosa, pero por otro, todas las alarmas saltando al pensar que, si ya las derivadas se bastaban solas para complicarle la vida a uno, nada bueno podía esperarse tras el hecho de introducirlas dentro de otra ecuación. No recuerdo qué ecuación era aquella pero sí que, ese simple comentario, desnudando de parte de su oscuro misterio a las temidas ecuaciones diferenciales, me sirvió para posteriormente profundizar un poco en su estudio, aunque siempre con una “mirada” biológica.
Supe así, que había ecuaciones diferenciales de primer orden, de orden n, ecuaciones no lineales, estabilidad, ecuaciones en derivadas parciales, etc. La enfermedad contraída fue tan grave como para llegar a impartir un curso de postgrado de ecuaciones diferenciales para biólogos (dos alumnas, a día de hoy, ambas retiradas de la investigación, quiero pensar que por motivos totalmente ajenos a aquel curso). Pasó el tiempo, y hondo fue mi pesar el día que comprendí que nada de todo esto era suficiente para describir los modelos más interesantes, y todavía básicamente sencillos, en genética de poblaciones donde lo necesario era, y es, otro tipo de ecuaciones, llamadas ecuaciones diferenciales estocásticas. Pero esa es otra historia...
Si has llegado hasta aquí, querido lector, cuentas con todo mi aprecio y también con algo de mi asombro. Por tu parte, acaso estás preguntándote con cierta melancolía, sobre los caminos que te han llevado hasta una entrada de blog como esta. ¿Qué puede haber fallado? pensarás... Por otro lado, aparentemente hemos abandonado a e en los últimos párrafos. Pero, como se suele decir, no hay fallo, está aquí mismo e inmediatamente retorna al escenario en, lo prometo, el último acto.
El número e, las ecuaciones diferenciales y la evolución de poblaciones varias
Veamos. Un modelo matemático en que aparece una expresión relacionando una función f(x), con su primera o sucesivas derivadas, es una ecuación diferencial ordinaria. Es decir, conocida la relación entre una función f, su variable x y sus derivadas, lo que deseamos es encontrar la función incognita que satisfaga esa relación. Resolver la ecuación diferencial consiste en integrarla. Recordemos aquí, que la función que relaciona magnitudes y crecimientos (aquella que es proporcional a su derivada) es ekx y por tanto es de esperar que el número e juegue un importante papel en la interpretación de las ecuaciones diferenciales, como así es.
Un tipo de modelos sencillos y muy comunes en biología de poblaciones, ecología, epidemiología, etc. toma como variable independiente el tiempo t de modo que la función y(t) representa a los individuos de la población, por ejemplo su número, siendo la derivada el cambio del tamaño de población en el tiempo. He aquí la primera respuesta muy genérica a la pregunta de porqué aparecía el número e por todas partes en biología en general y de poblaciones en particular. El tipo de modelo que relaciona cantidad de individuos o moléculas o virus, con su evolución (cambio) en el tiempo son las ecuaciones diferenciales cuya solución involucra casi por definición al número e.
Parece que fuel el economista Tomas Malthus (1766-1834) el que primero observó que muchas poblaciones biológicas crecen de modo proporcional a la población. En genética de poblaciones es común usar la letra N para la variable del número de individuos de una población (por tanto aquí ya no corresponde al número de instantes). Por ello, la función y(t) la notaremos entonces como N(t), o más bien como N para simplificar la notación, y la razón dN/dt de cambio de la población (esto es, su derivada) es proporcional al valor actual de N. Por tanto
dN/dt = rN.
La constante r es la razón de crecimiento (si es positiva) o disminución (si negativa). Si la población está creciendo asumimos r mayor que 0 e imponemos la condición inicial de que en tiempo 0 el tamaño poblacional es N0 es decir N(0) = N0 .
La solución de esta ecuación la obtenemos si, agrupamos las N a un lado de la ecuación y la t en el otro, dN/N = rdt, e integramos ambos lados de modo que
log N = rt + C donde log es el logaritmo en base e y C es una constante. Por tanto
N = ert+C = C’ert. Y como en t=0, N = N0 tenemos que N0 = C’ y finalmente
N = N0 ert.
Que es el modelo Malthusiano de crecimiento exponencial antecesor del modelo conocido como crecimiento logístico, que surge al sustituir la constante r de razón de crecimiento por una función del tamaño poblacional N, r = f(N). El sentido de hacer esto es que, si bien en muchas poblaciones el modelo exponencial es correcto durante un tiempo definido, es obvio que los recursos son limitados y mantener el crecimiento exponencial durante tiempo ilimitado no es factible.
El modelo logístico surge a partir del Malthusiano, para tener en cuenta que la razón de crecimiento depende de la población. Para incorporar un límite al tamaño poblacional según los recursos disponibles, basta elegir un f(N) que valga aproximadamente r para N pequeño, pero de modo que f(N) decrezca según aumenta N y que a partir de un N suficientemente grande, f(N) se haga negativo. La función más sencilla que cumple estos requisitos es
f(N) = r –aN donde a es positivo.
Si sustituimos esta función por r en la ecuación de Malthus obtenemos
dN/dt = Nr(1 – N/K) donde se ha hecho el cambio a = r /K.
Esta es la ecuación logística, introducida por el matemático belga P.F. Verlhulst como modelo para el crecimiento de la población humana, con la poca fortuna de disponer de datos de censo inadecuados con lo que su modelo cayó en el olvido. Lotka volvió a obtener la ecuación en 1925 y la llamó ley del crecimiento poblacional. Al parámetro r se le conoce como razón de crecimiento intrínseco y K es la cota superior a la que tiende el tamaño poblacional N tras un tiempo infinito. Por tanto a K se le conoce como nivel de saturación o, en ecología y biología de poblaciones, como la capacidad de carga de la población (carrying capacity). La solución de la ecuación incorpora, como no, nuestra ya conocida ert.
La importancia de esta ecuación ha sido tal, que en ecología se distinguen dos estrategias de supervivencia conocidas como estrategia-r y estrategia-K en relación al conflicto entre número de descendencia versus cuidados parentales invertidos en la descendencia. La estrategia r fomenta muchos descendientes en los que se invierte escaso cuidado parental (corresponde a valores elevados de la tasa de crecimiento r y bajos de K) mientras que la estrategia K implica poca descendencia en la que se invierte mucho cuidado y por tanto tienen mayor probabilidad de supervivencia (baja r, elevada K). Si bien es necesario aclarar que a día de hoy, este modelo K-r ha sido superado por otros más completos, que incluyen distintos componentes de la historia vital de las especies teniendo en cuanta, por ejemplo, la estructura de edades.
Obviamente, queda mucho por decir sobre el papel que el número e tiene en la formalización y modelado biológicos, algo que, con cierto desenfado, he llamado biología del número e. El rastro de e lo podemos seguir a través de distribuciones probabilisticas de gran importancia en Biología como la Normal, la Poisson, etc, así como también en procesos estocásticos. Es clave en el estudio de la evolución de las secuencias de ADN, donde aparece la ecuación P = eQt donde P y Q son matrices, que permite calcular las probabilidades de transición de unas secuencias de información genética a otras. También hallaremos a e en la dinámica de sistemas ecológicos, etc etc... En cualquier caso, paremos aquí. Aunque yo no esté satisfecho, y tú ¿lector? estés harto... A veces me pregunto ¿porqué e? ¿porqué precisamente entre el 2 y el 3?.. Pero dado que llevo dos días escribiendo esta entrada sin ton ni son, será mejor para todos, terminar aquí ahora.
PD: sí, ya sé que podía haber puesto alguna figura!
Bibliografía
Bessiere, G. "El cálculo integral fácil y atractivo" trad por A.G.T. Ed. Naval, 1946.
Boyce, WE, DiPrima, RC. "Ecuaciones diferenciales y problemas con valores en la frontera". Ed. Limusa, 1998.
Dunham, W. "Euler. El maestro de todos los matemáticos". Ed. NIVOLA, 2000.
Roughgarden, J. "Theory of population genetics and evolutionary ecology. An introduction". Ed.Prentice-Hall, 1996.
“Esta entrada participa en la edición 3.1415926 del Carnaval de Matemáticas, alojado en el blog Series divergentes.”
Suscribirse a:
Entradas (Atom)