viernes, octubre 23, 2009

TNT macros: Variables (2)


In a previous entry I wrote about variable declaration under TNT macro language, in this post I show the most important aspect of variables: their access.

Variables can be accessed in two ways, for reading, in which the value stored on the variable is retrieved. The second access way is for writing, in which a value is assigned to a variable.

General variables

In TNT general variables (variables declared with var keyword) are accessed by its name inside single quotations (').
quote The number stored on variable is 'value' ;
hold 'number' ;
In this example, the number stored on value is printed and the number stored in number is used as the maximum tree hold. This is an important feature of the language: you can replace any command parameter by the value of variables. This gives an extraordinary flexibility to TNT's macros.

Sometimes, specially at linux version, it is necessary to use parenthesis to read successfully the variable:
hold ('numero') ;
To access arrays, the same format is used, with an index inside squared brackets ([]).
quote the value of the fourth element is  'vector [4]' ;
It can be a more complex (and useful) way for this example
quote the value of the 'i' element is 'vector ['i']' ;
Note that the variable i is also inside quotations. If you are somewhat familiar with some computer programming language, you note that quotations are a way to "dereference" of a variable.

This is an example for a multidimensional array
quote the value of the element 'i' , 'j' is 'mat ['i' , 'j']' ;
Using parenthesis (()) you can use math operations as index values. Here is an example with a two dimensional array, but that is exactly equal with any kind of array:
quote the value of the element 4, 'j' is 'mat [ (2 + 2 ) , 'j' ]' ;
Another option is the use of series
quote the value of the elements 3 to 8 is 'list [ 3 - 8 ]' ;
In this case, note that there are no parenthesis limiting the scope! If you put parenthesis, TNT interpret it as a mathematical operation.

To write values on general variables, the keword set is used with the name of the variable to assing, without quotations, and followed with the assigned value
set val 4 ;           /* Assigns 4 to val */
set res (3 + 4) ; /* Assigns a math operation */
set num 'val' ; /* Assigns the content of val to num */
Note that as is the value of the variable that is assigned to num, then val must be inside quotations. To assign arrays, it is possible to use the same format, just using one element at time.
set vec [4] 5 ;                 /* Assigns 5 to the 4th element of the array*/
set arreglo [ (2 + 3) ] 8 ; /* Assigns 8 to the 5th element of the array */
set arreglo ['i'] 10 ; /* Assigns 10 to the i-element of the array */
set arreglo [ ('j' + 'k') ] 7 ; /* Assigns 7 to the element j + k of the array */
Note that indexes must be dereferenced!

Many times, an operation to a variable is just a modification of its value, for example increase its value by one
set val 'val' + 1 ;
It is more clear and intuitive, modifying directly the variable without dereference it, this can be done using a C-like sintaxys
set val ++ ;              /* increase by 1 */
set arr [3] -- ; /* decrease the content of element 3 in 1 */
set mat [2, 4] += 'val' ; /* adds the value of val to the element 2, 4 of mat */
set num *= (3 + 'j') ; /* multiplies the content of num by (3 plus j)-times */
set arr [2] /= 3 ; /* divides the content of element 2 of arr by 3 */
set val -= 4 ; /* decreases the content of val in 4 */
In some cases, you want to store all array elements at the same time. This can be done with the keyword setarray. It is important that array dimensions coincide with the dimensions used in the declaration. The format of setarray is the name of the array followed by the name of the array an the elements.
var:
lista [5, 4]
;
setarray 5, 4 lista 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 2 2 2 2 ;
In this example, the following matrix is assigned to a 5x4 array
0 1 1 1
1 0 1 1
1 1 1 0
2 2 2 2
The order of the elements follows the dimensions from left to right, then the first four elements assigned to the array are the array elements [0, 0], [0, 1], [0, 2] and [0, 3].

Variables in loops

Usually, loop variables are control variables, so the best option is not modifying them. So in principle, write the code as if loop variables are read-only. If a loop variable needs to be modified TNT can do it, but that usually shows a design problem, and it is not a good practice to make that modifications.

Loop variables can be accessed using using names or number. Is a good practice to use names to identify loop variables.

To read a loop variable the number or name of the cicle must be preceded by number character (#). For example to save a simple search of several k values using implied weights [1]
loop =kval 1 10
keep 0 ;
piwe = #kval ;
mult ;
tsave* resu#kval.tre ; save; tsave /;
stop
As #kval is the name of the first cicle, then is possible to use #1, then to save the tree it would be
tsave* resu#1..tre ;
Note the double point after #1, that is because a point is interpreted by TNT as a decimal fraction, and assumes that the point is part of the name. In any case, it is better to use names ;). The numbering of loops is assigned in relation with its nestedness, starting with #1.

To modify a loop variable the keyword setloop is used, that change the most nested loop (i.e. the loop in which the instruction is used). As changing the value of loops distorts their sequence, you must be aware of infinite-loops. This is the reason that makes the change of loop variables unsafe and unrecommended. But, some times it is necessary a non-stop loop that finish only when a specified condition is meet. In that case, using endloop coupled with setloop can be very useful
var:
i
;

set i 0 ;
loop =non 0 1
/* several instructions */
if ('i' == 1)
endloop ;
else
setloop 0 ;
end ;
stop
It is assumed that somewhere inside the loop, the value of i is changed to 1.

Command line variables

Variables from command line are read-only. Then, we only can access their value. To access a command line variable it is used the percent sign (%).
set j %1 ;
Assigns j to the first parameter from the command line (by the way, %0 is the name of the macro). Trying to read more parameters than parameters actually provided, is an error and stops the execution. When you write a macro, keep in mind that command line parameters is the only way to communicate with the user, so it is necessary to do a good choice of reading order and default values.

As always, do not forgive to keep your TNT copy updated. And check the TNT wiki or join the TNT google group to deal with any question!

References
[1] Goloboff, P.A. 1993. Estimating character weights during tree search. Cladistics 9: 83-91. DOI: 10.1111/j.1096-0031.1993.tb00209.x

Previous post on TNT's macros

Macros de TNT: Variables (2)


En la pasada entrega (hace como mil años :P) escribí sobre como declarar las variables en los macros de TNT, en el post de hoy, voy a mostrar lo verdaderamente importante de las variables: su acceso.

Las variables se pueden acceder de dos maneras, una es la lectura, en donde el valor guardado en la variable es leído de alguna manera. La segunda forma de acceso es la escritura, en donde se asigna un valor a una variable.

Variables generales

En TNT para leer a las variables generales (las declaradas con la palabra clave var) se utiliza el nombre de la variable entre comillas simples (').
quote El valor de la variable es 'valor' ;
hold 'numero' ;
En este ejemplo, se imprime el número contenido en la variable valor. Y se retienen solo la cantidad de árboles que indica numero. Esta cualidad es muy importante, porque podemos reemplazar los parámetros de cualquier comando de TNT por variables, lo que le da una flexibilidad enorme a los macros.

A veces, en especial en la versión de linux, es necesario usar un paréntesis para que se lea correctamente la variable:
hold ('numero') ;
Para acceder a los arreglos, se utiliza el mismo formato, con el indice entre paréntesis cuadrados ([])
quote valor del elemento 4 del arreglo es 'arreglo [4]' ;
Es posible ser un poco más complejo es este caso
quote valor del elemento 'i' del arreglo es 'arreglo ['i']' ;
En este caso, se muestra el valor del elemento i dentro del arreglo. Véase que i va dentro de comillas. Para los que estén familiarizados con los lenguajes de programación pueden ver que al acceder a el valor de una variables este se "dereferencia" con las comillas.

Este es el ejemplo con arreglos multidimencionales
quote valor del elemento 'i', 'j' es 'dosd ['i', 'j' ]' ;
El indice del arreglo puede ser una operación entre paréntesis típicos (()), el ejemplo va con un array de dos dimensiones, pero es igualmente valido para arreglos de cualquier dimensión.
quote valor del elemento 4, 'j' es 'dosd [ (2 + 2), 'j' ] ;
Otra posibilidad muy útil en arreglos es el uso se series
quote el valor de los elementos 3 al 8 son 'lista [ 3 - 8 ];
En este caso es importante ¡no colocar los valores entre paréntesis! De lo contrario TNT interpreta una operación matemática.

Para escribir variables en TNT se utiliza la palabra clave set, usando el nombre de la variable a asignar, sin comillas, y seguida del valor asignado.
set valor 4 ;                   /* Asigna 4 a valor */
set resultado (3 + 4) ; /* Asigna una operación */
set numero 'valor' ; /* Asigna el contenido de valor a numero */
Véase que cuando lo que se asigna es el valor dentro de una variable, esta debe ir entre comillas. Para asignar arreglos se puede utilizar el mismo formato, asignando un elemento a la vez
set arreglo [4] 5 ;             /* Asigna 5 al elemento 4 del arreglo */
set arreglo [ (2 + 3) ] 8 ; /* Asigna 8 al elemento 5 del arreglo */
set arreglo ['i'] 10 ; /* Asigna 10 al i-elemento de arreglo */
set arreglo [ ('j' + 'k') ] 7 ; /* Asigna 7 al elemento j + k del arreglo */
Es importante anotar que los indices del arreglo deben dereferenciarse!

Muchas veces, las operaciones sobre las variables consisten en modificar su valor, por ejemplo, incrementar su valor en 1.
set valor 'valor' + 1 ;
Es mucho más claro e intuitivo, modificar directamente la variable sin dereferenciarla, esto se logra mediante operadores que usan una sintaxis como la de C en esos casos.
set valor ++ ;                 /* incrementa valor en 1 */
set arreglo [3] -- ; /* disminuye el contenido del elemento 3 de arreglo en 1 */
set matriz [2, 4] += 'valor' ; /* adiciona el contenido de valor al elemento 2, 4 de matriz */
set numero *= (3 + 'j') ; /* multiplica el contenido de numero por 3 más j */
set arreglo [2] /= 3 ; /* divide el contenido de arreglo en 3 */
set valor -= 4 ; /* disminuye el contenido de valor en 4 */
Para los arreglos puede ser más útil asignar a la vez todos los elementos, para ello se utiliza setarray. Es muy importante que las dimensiones del arreglo coincidan con las de su declaración. El formato de set array es las dimensiones del arreglo, seguido por el nombre del arreglo, seguido por los elementos.
var:
lista [5, 4]
;

setarray 5, 4 lista 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 2 2 2 2 ;
En este ejemplo se creo un array de 5 x 4 que tiene la siguiente matriz
0 1 1 1
1 0 1 1
1 1 0 1
1 1 1 0
2 2 2 2
El orden de los elementos es siguiendo el orden de las dimensiones de izquierda a derecha, así en el ejemplo los primeros 4 elementos son los elementos [0, 0], [0, 1], [0, 2], y [0, 3] del arreglo.

Variables en ciclos

En general, las variables de los ciclos son de control, y por ello, lo mejor es no modificarlas. Así en principio, las variables de los ciclos son de solo lectura. Si se necesita modificar el valor de la variable de un ciclo, TNT también puede hacerlo, aunque es importante recalcar que las modificaciones de las variables en un ciclo, suelen reflejar algún problema de diseño, y no es una buena practica recurrir a esas modificaciones.

Las variables de ciclos se pueden leer usando nombres o números. Es una buena practica el utilizar nombres para identificar las variables de los ciclos.

Para leer una variable de un ciclo se utiliza el número del ciclo antecedido por el carácter numeral (#). Por ejemplo, para recorrer varios valores de k usando pesos implicados [1]
loop =kval 1 10
keep 0 ;
piwe = #kval ;
mult ;
tsave* resu#kval.tre ; save; tsave /;
stop
El #kval se refiere al primer ciclo, también podría utlizarse #1, en ese caso, para guardar el arbol debería usarse
tsave* resu#1..tre
Véase que al guardar el nombre del árbol se usa un doble punto despues del número. Eso es porque con un solo punto, TNT interpreta que forma parte del número (es decir, que se trata de un real con decimales). De todas formas, mejor utilizar el nombre ;). Los números de los ciclos se asignan de acuerdo a su anidamiento (empezando desde #1).

Para modificar el valor de una variable de un ciclo se utiliza la palabra clave setloop, que modifica el ciclo más anidado en el que se encuentra la instrucción. Dado que el cambiar el valor de los ciclos distorsiona la secuencia del ciclo, hay que tener cuidado con los ciclos infinitos. Es por eso que su uso no es recomendable. Sin embargo, a veces necesitamos de un ciclo que corra indefinidamente hasta que se cumpla alguna condición. En ese caso, utilizar setloop junto a endloop puede ser muy útil. Por ejemplo
var:
i
;

set i 0 ;
loop =non 0 1
/* varias instrucciones */
if ('i' == 1)
endloop ;
else
setloop 0 ;
end ;
stop
Aquí se asume que en alguna condición dentro del ciclo el valor de i se cambia a 1, lo que permite que el ciclo finalice.

Variables de linea de comando

Las variables que vienen desde la linea de comando son de solo lectura. Así solo podemos accesar su valor. Para acceder al valor de una variable de linea de comando se usa el signo (%)
set i %1 ;
asigna a i el valor del primer parámetro pasado por la linea de comando (por interés %0 es el nombre del macro). Tratar de leer más parámetros de los que se pasaron por linea de comando, es un error, y detiene la ejecución. A la hora de escribir un macro, hay que tener en cuenta que son estos parámetros de la línea de comando la única manera de comunicarse con el usuario, así que hay que hacer una buena elección de su orden, y de los valores por omisión.

Como siempre, no se olviden de chequear que su versión de TNT esta actualizada, así mismo no duden en consultar la wiki de TNT o suscribirse al grupo google de TNT!

Referencias
[1] Goloboff, P.A. 1993. Estimating character weights during tree search. Cladistics 9: 83-91. DOI: 10.1111/j.1096-0031.1993.tb00209.x

Post anteriores en macros de TNT

lunes, octubre 19, 2009

The behemot! for free! / El behemot gratuito!

Santí noted that the "behemot"'s paper [1] is now available for free from Wiley/Blackwell web-page, this is the link:


Or you can jump directly to the paper abstract and download the pdf:


If you want to look the results, Pablo post a quick guide to manage the trees and extract info from them. If you want to check you favorite group, or just play for a moment, there are no excuses ;)!

And the best of it, as they are macros, you can use for your own trees! and don't forget to check out the TNT wiki!


Santí me hizo notar que el paper donde el "behemot" vio la luz [1] se encuentra disponible de manera gratuita en el sitio de blackwell. El enlace es este:


O pueden ir directamente a la pagina del artículo y descargar el pdf:

Para quienes quieran jugar con los resultados, Pablo publicó una guía rápida para manejar los árboles y poder extraer info de ellos. Así que si tienen curiosidad por ver como quedo el grupo que ustedes trabajan, o simplemente quieren pasar un rato, ya no hay excusas ;)

Y lo mejor de todo, como son macros, se pueden usar en sus propios resultados :D! No se olviden de consultar la wiki de TNT ;)

[1] Goloboff, P.A. et al. 2009. Phylogenetic analysis of 73 060 taxa corroborates major eukaryotic groups. Cladistics 25: 211-230. DOI: 10.1111/j.1096-0031.2009.00255.x

sábado, agosto 29, 2009

TNT macros: variables (1)


I forget an important feature. The commentaries.

TNT's macro language support the same multiline comment style of C (but not the single line comment style), enclosing comments within /* and */. Comments has no effect on the script, but are very useful to document the code!

/* This is a comment */

Comments are really important, no only to made code legible to other people, but also to our own understanding (specially, if it is code written long time away!).

Comments are valid after the macro mode is open (using macro =).

Following a well knew tradition I put this simple macro ;)

macro =;
/* My first TNT macro */

quote Hello world! ;
p/ ;

Now, we can enter on the subject of variables...

Variables are objects that the program manipulates by an specific objective given by the user. Variables store values, that can be modified or read. Usually we are interested in that value. In other circumstances the variable value is useful for controlling the program flow. Also they are used to store values given by the user.

When macro mode is open, TNT sets a default number of variables (1000). They can be accessed by a number, starting with 0. It is possible to change the amount of variables using macro* N K ; where N is the number of loops, and K the number of variables.

General variables

Using numbers to access variables can be useful for small scripts, but with more complex macros, the management of code can be very difficult. Then, is a good practice to use names for variables. Variables can be named in any part of the code. It is a good practice naming it just after starting the macro's file (just after macro = ). Variables names are declared using the keyword var.

There are two ways to name variables, the first is explicit, that is useful to backward compatibility with scripts written in old versions of TNT.

Var =
0 variable1
1 variable2
5 variable 3
;

In this format, using var =, each variable is declared using the internal number of TNT. The declaration finish with a semi colon. The principal problem with this style, is that it is an invitation to use the numbers instead of the names. This can be dangerous when parts of the code change (for example, adding more variables).

A more elegant way, and more secure (protects against usage of unnamed variables), is declaring directly the names, using var : instead of var =

var :
variable1
variable2
variable3
;

Here, I always use this style, that produce a more legible code.

TNT variables can be simple, that is, just with one value, or can be arrays, a vector of several consecutive values. TNT arrays are static, that is, they cannot be resized in running time. Then declare them only after you know the size of the array. For example, just after reading the data, and we known the number of nodes in a tree.

An array declaration is simple, just as a normal variable, but dimensions inside brackets ([]), if it is a multidimensional matrix, use a comma (,) to separate the values, in a similar way that arrays are declared in pascal.

Var:
anArray [25]
aMatrix [3, 4]
;

Variables in loops

In TNT loop variables are independent from general variables. This variables are managed by TNT and not directly by the programmer. As general variables they can be named, o accessed using a number that start from 1, and increase in each nesting level. I generally use numbers. Maybe a better practice is to use names, as that freed the loop from its nesting level, an that's the only way to access some variables (for example, if you access it from another file).

Loop variables are declared using an equal before the name:

loop =cicle 1 10

Variables in the command line

When a macro is called, it is possible to assign some initial values from the command line associated with the script. This variables are sometimes known as arguments or parameters. For example, a macro called dojac can receive as parameters the number of replicates and the number of iterations per replicate

dojac 1000 20 ;

Then, dojac receives 1000 replicates and each replicate with 20 iterations.

As TNT's scripts are not interactive, this is the only way to the final user to modify the behavior of the scripts. Then, if the idea is to distribute the macro, calling the macro without parameters shows a help screen explaining the objective, the parameters, and the conditions necessary for using the macro in a correct way.

The function argnumber returns the number of arguments used when calling the macro.

Not forget to check out for TNT updates, and of course, several scripts and useful documentation at the TNT wiki :D

lunes, julio 13, 2009

Macros de TNT: Variables (1)


En el post anterior olvide algo muy importante. Se trata de los comentarios.

El lenguaje de macros de TNT permite el uso de comentarios en la misma forma que el tradicional C, es decir entre un /* y un */. Los comentarios no tienen ningún efecto sobre el script, pero son muy útiles para documentar el código.

/* Este es un comentario */

Es muy importante colocar comentarios, no solo para que otras personas puedan entender el código, sino para que uno mismo pueda entenderlo más rápidamente cuando vuelve a el (en especial, si fue escrito hace tiempo!).

Los comentarios son validos, una vez se a iniciado el modo de macros (usando macro=)

Por tradición, coloco este macro simple

macro =;
/* Mi primer macro de TNT */

quote Hola mundo!;
p/;

Ahora entrando en materia con el tema de este post...


Las variables son los objetos que el programa manipula con un propósito especificado por el usuario. Las variables guardan valores, y estos valores pueden ser modificados o leídos. Muchas veces, es el valor contenido en una variable lo que nos interesa. En otras ocasiones, el contenido de una variable sirve para controlar el flujo del programa. También pueden usarse variables para guardar valores dados por el usuario.

Al arrancar un macro, TNT define por omisión un número de variables (1000). Estas se pueden acceder usando un número, que inicia desde 0. Es posible cambiar el numero de variables usando macro* N K; donde N es el número de loops, y K el número de variables.

Variables generales

Usar los números de las variables puede ser cómodo en ,macros pequeños, pero al incrementar la complejidad de los macros, se hace más difícil el mantenimiento del código, por lo que es una buena práctica usar nombres para las variables. Las variables se pueden nombrar en cualquier parte del programa. Pero es una buena práctica nombrarlas justo después de iniciar el archivo con el macro (después del macro=). Para su declaración se utiliza la palabra clave var.

Hay dos formas de nombrar las variables, una es una forma explicita, que sirve, más que todo, para mantener la compatibilidad con macros anteriores.

var =
0 variable1
1 variable2
5 variable3
;

En este formato, usando var= se tiene que nombrar cada variable con el número interno (de TNT) que la identifica. La declaración termina con un punto y coma. El principal problema con esta forma, es que invita al uso intercambiable de números y nombres en las variables. Eso puede tener un efecto nosivo a la hora de hacer mantenimiento al código (por ejemplo, querer agregar más variables).

Una forma mucho más cómoda, y más recomendable (evita el uso de variables sin nombre), es declarar directamente los nombres, usando dos punto (:) en vez de igual

var:
variable 1
variable 2
variable 3
;

De aquí en adelante, utilizo esta forma, que me parece, mejora la legibilidad del código escrito.

Las variables en TNT pueden ser simples, es decir, que contienen un solo valor, o pueden ser arreglos, series de valores consecutivos. Los arreglos de TNT son estáticos, es decir que no pueden redimensionarse en tiempo de ejecución. Por eso es bueno declararlos solo cuando ya sabemos el tamaño que puede tener. Por ejemplo, luego de leer los datos, ya sabemos cuantos nodos puede tener un árbol.

Declarar un arreglo es muy simple, es como una variable normal, pero con las dimensiones entre corchetes ([]), si es una matriz multidimensional, los valores se separan con comas (,) de forma similar a como se declaran los arreglos en pascal.

var:
arreglo [25]
matriz [3, 4]
;


Variables en ciclos

En TNT las variables de los ciclos son independientes de las variables generales. Estas variables son manejadas por TNT y no directamente por el usuario. Al igual que las variables tradicionales, pueden nombrarse, o usar los números que vienen por omisión que comienzan desde 1, y se incrementan al irse anidando los ciclos.

Yo en general utilizo los números. Quizá sea una mejor práctica utilizar nombres puesto que de esa manera se independiza de el anidamiento del ciclo, y además, es la única forma de acceder variables de ciclos desde otro archivo.

La forma de declarar una variable en un ciclo es anteponiendo un igual (=) antes del nombre:

loop =ciclo 1 10


Variables de linea de comando

Al llamar a un macro, es posible asignarle valores iniciales usando la línea de comandos asociada al llamado del macro. Estas variables se conocen como "argumentos" o "parámetros". Por ejemplo, uno podría tener un macro que se llame dojac, y que reciba como parámetros el número de replicas, y un número de iteraciones en cada replica, así

dojac 1000 20 ;

En este caso, se reciben 1000 replicas y cada replica tiene 20 iteraciones.

En general, cuando uno distribuye un macro, si no se da ningún parámetro, se ofrece una especie de salida de ayuda, donde explica el objetivo y las condiciones necesarias para utilizar el macro de forma correcta.

La función argnumber retorna el número de argumentos usados al invocar el macro.

No se olviden de revisar las novedades de TNT, y por supuesto, consultar la wiki de TNT ;), hay varios scripts y otra documentación muy útil :D

sábado, mayo 16, 2009

Algorithms for phylogenetics 0b: Trees


Another basic component of phylogenetic analysis' programs are the trees. Trees can be the result of an analysis (as in cladistic analysis' programs) or part of the data used (as in biogeography analysis' programs).

Trees are one of the most well known data structures of computer science. In computers science a tree is a collection of elements (nodes), nodes are related by a "kindship" relation that imposes a hierarchy on the nodes. Kinship is a pairwise relationship: one node is the ancestor and the other the descendant. A node can have an indefinite number of descendants, but, at maximum, only one ancestor.

In general, in phylogenetics, the term node is used only for elements with descendants, and terminal of leafs to elements without descendants (in computer science the used terms are internal node and terminal node, or leaf node. Here I will use the terminology from computer science, and when using node I refer to any element of the tree).

A trees has at maximum a single root node. A root node is an internal node without ancestor. Different from most abstract trees of computer science, the only elements with have the user data are the terminals (from the data matrix). Data stored by internal nodes is inferred through algoritmic means.

There are many ways to represent a tree. For me, the most natural one is using structures and pointers. Here is the basic node
typedef struct tagPHYLONODE PHYLONODE;

struct tagPHYLONODE {
PHYLONODE* anc;
PHYLONODE* left;
PHYLONODE* right;
STATES* recons;
};
This structure is a node from a binary tree, that is, each internal node have two descendants, the pointers left and right, note that they are pointers to PHYLONODE elements. The pointer to ancestor is anc. In this structure, the root node has anc as NULL (i.e. points to no elements), and terminal nodes has left and right as NULL. The array recons store the data (as bitfields).

I usually store the information from a terminal (extracted from data matrix) into an independent structure, for example
typedef struct tagTERMINAL TERMINAL;
struct tagTERMINAL {
char name [32];
STATES* chars;
};
And include a pointer to TERMINAL as part of PHYLONODE
struct tagPHYLONODE {
... //Otros campos aquí
TERMINAL* term;
};
Internal nodes has term as NULL.

A nice property of trees, is that subtrees have the same properties of a tree, then recursive algorithms a very natural way to work with trees. Nevertheless, I always prefer to have an independent structure to store the whole tree
typedef struct tagPHYLOTREE PHYLOTREE;
struct tagPHYLOTREE {
unsigned numNodes;
PHYLONODE* root;
PHYLONODE* nodeArray;
DATAMATRIX* matrix;
};
This structure store a pointer to the root node (root), an array to the nodes of the tree (nodeArray), a pointer to a data matrix (matix) and the number of nodes.

This approach has the advantage of having several independent sets of trees, each tree with his on characterisitcs (id, name, data matrix), that are no properties of any node, but of the whole set of nodes (A typical case: in biogeography, several cladograms form different organism are used).

Different from structures it is possible to store trees as coordinate arrays. This option is used in [1] for the example code, and it is also the way in which TNT macro language uses trees (the language does not support structures). The declarations of arrays is simple
int* anc;
int* left;
int* right;
STATES** chars;
In this case, all arrays are dynamic (they are declared as pointers, but as soon as they are assigned, they can be treated as arrays), and instead of have a pointer, they store the number that is the index of the array. As it is a coordinate array, every element with the same index is the same node. For example anc [17] = 20, assigns the node 20 as the ancestor of node 17. In this case, left [20] or right [20] must be equal to 17. chars [20] has the sates assigned to node 20. If you look at the tree entry in wikipedia (specially the entry on heaps) you will find that is the way used to keep trees.

For me, the bad side of that approximation is that requires an strict bookkeeping (take note that in TNT macro language, TNT automatically updates the information, so that is not a problem!). This is my own experience, I worked with tree reconciliation, and work with someone that does not know working with pointers. So I try a coordinate array approach. But we not know how many nodes would be on the tree after reconciliation, then the code turns messy and extremely difficult to update.

That experience does not show if some alternative is best than the other, just that sometimes the nature of the problem might constraint the way to resolve it, and that some styles of programming are better than others for some programmers ;). In this series I always will use structures.

Just to have some code, and to mark the recursivity of trees, here are a function to store a tree in a parenthetical notation (as in Hennig86, NONA and TNT)
#include < stdio.h >

void PrintTree (FILE* out, PHYLOTREE* tree) {
fprintf (out, "tread \'Simple tree output\'\n");
PrintNode (out, tree - > root)
fprintf (out, ";\np/;");
}

void PrintNode (FILE* out, PHYLONODE* node) {

if (node - > term != NULL) {
fprinf ("%s ", node - > term - > name);
return;
}
fprintf (out, "(");
PrintNode (out, node - > left);
PrintNode (out, node- > right);
fprintf (out, ") ");
}
This functions used standard IO operations. The function PrintTree puts the tree header for Hennig86/NONA/TNT trees. And call PrintNode on the root node. PrintNode checks if the node is a terminal, if this is true, it writes the terminal name, if not, the node is an internal one, so it prints the parenthesis and call recursively PrintNode on each one of its descendants.

I think that some time ago Mike Keesey write a post about trees under a object oriented paradigm. I'm unable to find the post, this one is the most closely alternative :P.

References
[1] Goloboff, P. A. 1997. Self weighted optimization: tree searches and character state reconstruction under implied transformation costs. Cladistics 13: 225-245. doi: 10.1111/j.1096-0031.1997.tb00317.x

Previous post on Algorithms for phylogenetics
0a. Bitfields

viernes, mayo 15, 2009

Algoritmos en filogenia 0b: Árboles


Otro de los componentes básicos de un programa de análisis filogenético son los árboles. Los árboles pueden ser el resultado del análisis (como en los programas de análisis cladístico), o pueden ser parte de los datos usados (como en los programas de análisis biogeográfico).

Los árboles, son además, una de las estructuras de datos más conocidas y usadas en programación. En programación un árbol es una colección de elementos (nodos), los nodos están relacionados entre sí por un "parentesco" que les impone una jerarquía. El parentesco es una relación de pares: un nodo es el ancestro y el otro es el descendiente. Un nodo puede tener un número indefinido de descendientes, más como máximo solo puede tener un ancestro.

En general, en filogenía se llama nodo, únicamente a los elementos que poseen descendientes, y hojas o terminales, a los elementos sin descendientes (en ciencias de computadores se utiliza más nodo interno y nodo terminal o nodo hoja, aquí yo sigo la notación de ciencias de computadores, y utilizo nodos cuando me refiero a todos los elementos del árbol).

Un árbol tiene como máximo un nodo raíz, que es un nodo interno sin ancestro. A diferencia de los árboles más abstractos de ciencias de computadores, los únicos elementos con datos incluidos por el usuario son los terminales (que entran vía la matriz de datos). Los datos guardados por los nodos internos son inferidos mediante algún algoritmo.

Existen varias formas de representar un árbol. Para mi, la más natural es usando estructuras y apuntadores. He aquí un nodo básico
typedef struct tagPHYLONODE PHYLONODE;

struct tagPHYLONODE {
PHYLONODE* anc;
PHYLONODE* left;
PHYLONODE* right;
STATES* recons;
};
Esta estructura es el nodo de un árbol binario, es decir, que cada nodo interno tiene dos descendientes, los apuntadores left, y right, que como pueden notar apuntan a elementos del tipo PHYLONODE. anc es el puntero a el ancestro. En esta estructura el nodo raíz tiene a anc como NULL (es decir que no apunta a ningún elemento), y los nodos terminales tienen left y right como NULL. recons es un arreglo con los caracteres (en forma de campo de bits).

Yo suelo guardar la información del terminal (extraída de la matriz de datos) en una estructura aparte, por ejemplo
typedef struct tagTERMINAL TERMINAL;
struct tagTERMINAL {
char name [32];
STATES* chars;
};
E incluyo un apuntador a TERMINAL como parte de PHYLONODE
struct tagPHYLONODE {
... //Otros campos aquí
TERMINAL* term;
};
Si es un nodo interno, term es NULL.

Una de las ventajas de los árboles es que un subárbol tiene las mismas propiedades de un árbol, por lo que los algoritmos recursivos son muy naturales al trabajarlos con árboles. Aún así, yo prefiero tener siempre una estructura que guarda el árbol
typedef struct tagPHYLOTREE PHYLOTREE;
struct tagPHYLOTREE {
unsigned numNodes;
PHYLONODE* root;
PHYLONODE* nodeArray;
DATAMATRIX* matrix;
};
En esta estructura se guarda un puntero al nodo raíz (root), un array con los nodos del árbol (nodeArray), un apuntador a la matriz de datos (matrix) y el número de nodos.

Esta aproximación tiene la ventaja de que uno bien puede tener varios conjuntos de árboles de forma independiente, por ejemplo, cada árbol con su propio id, o nombre, que no es una característica del nodo raíz, sino del conjunto completo de los nodos.

Aparte de las estructuras es posible guardar los árboles como arrays coordinados. Esta opción, por ejemplo, es usada en [1] para el código que aparece en los ejemplos, y también es la manera de manejar los árboles en el lenguaje de macros de TNT (que no maneja estructuras). La declaración de los arreglos coordinados es más sencilla
int* anc;
int* left;
int* right;
STATES** chars;
En este caso anc, left y right son arreglos dinámicos (por eso están declarados como apuntadores, pero se pueden tratar como arrays), y en vez de contener un apuntador, contienen un numero que es el índice del array. Como es un arreglo coordinado, todos los elementos con el mismo índice son parte del mismo nodo. Así, por ejemplo anc [17] = 20; dice que el ancestro del nodo 17, es el nodo 20. En este caso left [20] o right [20] debe ser igual a 17. chars [20] contiene los estados asignados al nodo 20. Si revisan las entradas sobre árboles en wikipedia (en especial las entradas sobre heaps) notaran que esa es la manera con la que tratan los árboles.

Para mi, lo malo de esa aproximación, es que requiere un control más fuerte sobre el mantenimiento de la información (cabe anotar que en el lenguaje de macros TNT, TNT mantiene todo por nosotros, así que eso no es ningún problema!). En mi propia experiencia, yo trabaje en reconciliación de árboles, y trabajaba con alguien que no sabia manejar apuntadores. Así que trate de hacer la implementación con arreglos coordinados. pero en árboles reconciliados no siempre sabemos cuantos nodos nuevos va a tener el árbol, por lo que mantener el código con los arreglos coordinados se hizo bien problemático.

Lo que muestra esa experiencia, no es que una alternativa es mejor que la otra, solo que a veces, la naturaleza del problema puede complicar la forma de resolverlo, y que además algunos estilos de programación se dan mejor que otros, según quien este haciendo el código ;). De aquí en adelante en otros posts, yo voy a usar estructuras.

Para tener algo de código y resaltar la recursividad de los árboles he aquí una función que guarda un árbol en formato parentética (como Hennig86, NONA y TNT)
#include < stdio.h >

void PrintTree (FILE* out, PHYLOTREE* tree) {
fprintf (out, "tread \'Simple tree output\'\n");
PrintNode (out, tree - > root)
fprintf (out, ";\np/;");
}

void PrintNode (FILE* out, PHYLONODE* node) {

if (node - > term != NULL) {
fprinf ("%s ", node - > term - > name);
return;
}
fprintf (out, "(");
PrintNode (out, node - > left);
PrintNode (out, node- > right);
fprintf (out, ") ");
}
Estas funciones usas las operaciones IO estándar. La función PrintTree coloca la cabecera de un archivo de árboles para Hennig86/NONA/TNT. Y llama a PrintNode. PrintNode chequea si el nodo actual es un terminal, si es así escribe el nombre del terminal, de lo contrario, como en un nodo interno, imprime los paréntesis y llama recirsivamente a PrintNode en cada uno de sus descendientes.

Me parece que Mike Keesey escribió un post sobre árboles bajo programación orientada a objetos, pero no pude encontrarlo, esta es la alternativa más parecida :P

Referencias
[1] Goloboff, P. A. 1997. Self weighted optimization: tree searches and character state reconstruction under implied transformation costs. Cladistics 13: 225-245. doi: 10.1111/j.1096-0031.1997.tb00317.x


Post anteriores en Algoritmos en filogenia
0a. Campos de bits

miércoles, mayo 06, 2009

TNT macros: Intro


One of the most powerful utilities of TNT is its macro's language. This language allows access to TNT internal variables, and tree and data manipulation. This coupled with the computer power of TNT makes a terrific combination.

Unfortunatelly, maybe because of a lack of an extensive manual, or because there are few papers that explicitly use TNT macros (but they are growing!), then this capability is ignored for several users. In this series of posts, I want to give an introduction of the many possibilities allowed by the usage of macros.

Part of this idea born after I see the book "phylogenetics with R" [1], and the constant exchange with Santi Catalano about the TNT's wiki. I hope to post part of this series on the TNT's wiki, but for the moment, just to get the writing feeling, I post it first on my blog ;).

The style is somewhat orientated to programming--I'm powerfully influenced by Kenighan and Ritchie :)--because I think it is the best way to learn the language.

Notation

The TNT macro language is an interpreted language, this is, each instruction is executed as it is readed, then it is possible to “write” the programs in real time (just like the old BASIC, lisp, or the fashionable python).

This can be very nice, for example, to make some simple mathematical operations directly on TNT.

An example:
> macro = ;
Macro language is ON
> var: dest ;
> set dest 4 + 5 ;
> quote 'dest' ;
9
This simple calculator can introduce us in some particularities of the language.

To activate macros, the command macro = ; must be invoked, macro - ; deactivates the macros.

Each instruction finish with a semi colon (;). Although not necessary, for clarity, an space before the semicolon is a good practice.

The keyword var is used to declare variables.

The keyword set assigns a value to a variable. In the example it is a sum. Set accepts the four basic operations, as well as more complex combinations using parenthesis.

The command quote prints on the screen. In TNT to access the value stored into a variable, the name of the variable must be between simple quotes ('). Note that set assigns a value, so the name is written without quotations, but if the value assigned is stored into a variable, then that values must be in quotations:
set dest 'first' * 'second' ;
This operation assigns to dest the value of first times second.

But more important, is that TNT macros can be written on separate file, and executing like any commando of the program.

If the macro is saved with the ".run" extention, and is on the same working directory of the program (for example c:\bin\tnt) it is possible to call the macro just typing the name of the file. For example, TNT is distributed with the macro stats.run, that calculates consistency and retention index of the trees in memory. To execute it, only is necessary to type:
> stats ;
in the program command line, and the macro runs automatically.

As any file of TNT it can be accessed using proc command (but you lost the parameter list) of the macro. It is better to use the command run (or just write the macro's name). As any command macros can be accessed trough the OS command line.

Exercises

As exercises try to use the different ways to call a macro, as well as accessing values with the commands, think on scripts like this:
macro = ;
var:
numIts
numTrees
itTrees;

set numIts 100 ;
set itTrees 10 ;
set numTrees 'numIts' * 'itTrees' ;

rseed 0;
hold 'numTrees' ;
mult = replic 'numIts' hold 'itTrees' ;
nel * ;
p/ ;
That does noting that you can do more easily by hand, but allows you to get familiarity with the notation.

References
[1] Paradis, E. 2006. Analysis of phylogenetics and evolution with R. Springer

Macros de TNT: Intro


Una de las utilidades más poderosas de TNT es su lenguaje de macros. Este lenguaje permite acceder variables internas de TNT, manipular árboles y datos. Esto en conjunto con el poder de computo de TNT es una excelente combinación.

Infortunadamente, bien sea porque no hay un manual extenso, o porque aún son pocos los que han usado explicitamente macros en sus artículos, esta utilidad suele pasar más o menos desapercibida para los usuarios. En esta serie de posts quiero dar un poco de remedio a esta situación, y dar una especie de introducción a las muchas posibilidades que permite el uso de macros.

Parte de la idea nació a que una vez vi el libro de "R para filogenia" [1], y pues a que he estado en contacto con Santi Catalano, que ha estado trabajando en la wiki de TNT. Espero poder postear después parte de esto a la wiki de TNT, pero por ahora, como para ver como me siento escribiendolo, lo coloco primero en mi blog ;).

El estilo, es bien orientado a la programación--mi influencia más poderosa es Kernighan y Ritchie :)--pues la idea es que se le saque el máximo provecho al lenguaje!

Notación

El lenguaje de macros de TNT es un lenguaje interpretado, esto quiere decir que cada instrucción se ejecuta a medida que se va leyendo, por lo que es posible ir "escribiendo" el programa en tiempo real desde la linea de comando.

Esto puede ser muy útil, por ejemplo, para realizar cálculos matemáticos directamente en TNT.

Por ejemplo:
> macro = ;
Macro language is ON
> var: dest ;
> set dest 4 + 5 ;
> quote 'dest' ;
9
Al escribir esa secuencia uno puede tener una calculadora simple. Esto sirve para introducirnos a las cosas particulares del lenguaje.

Para activar los macros, se usa el comando macro = ; con macro - ; se desactivan los macros.

Al terminar cada instrucción se coloca un punto y coma (;). Aunque no es necesario, por claridad es aconsejable dejar un espacio antes del punto y coma.

La palabra clave var se utiliza para declarar las variables.

Con la palabra clave set se asigna un valor a una variable. En este ejemplo el valor de una suma. Set acepta las cuatro operaciones básicas, así como operaciones más complejas usando paréntesis.

El comando quote imprime en pantalla. En TNT, para acceder al valor guardado en las variables, se coloca el nombre de la variable entre comillas simples. Así en este ejemplo, se imprime lo que esta dentro de dest.

Es importante diferenciar entre el valor guardado en la variable y su nombre. Set asigna a una variable, por ello la variable no va dentro de comillas, pero si lo que es asignado es otra variable, entonces como es un valor, debe estar entre comillas:
set dest 'primero' * 'segundo' ;
Esta instrucción asigna a dest la multiplicación de los valores de primero y segundo.

Sin embargo, la aplicación más importante de los macros de TNT es que pueden guardarse en un archivo y ejecutarse como si se tratara de un comando del programa.

Si los macros se guardan con la extensión ".run", y son colocados en el mismo directorio de trabajo del programa (por ejemplo c:\bin\tnt) es posible llamarlo simplemente usando el nombre del archivo. Por ejemplo, en la instalación de TNT esta incluido stats.run, que calcula el indice de consistencia y retención de los árboles en memoria. Para ejecutarlo, simplemente se escribe
> stats ;
en la linea de comando, y el macro es ejecutado automáticamente.

Como todos los archivos de TNT se pueden acceder usando el comando proc, pero se pierde la lista de parámetros (si los hay) del macro. Por eso es mejor usar el commando run (o simplemente escribir el nombre). Como cualquier comando se puede acceder directamente en la lista de parámetros del programa.

Como ejercicios de esta introducción se puede probar las diferentes maneras de invocar los macros, así como acceder valores desde comandos. Por ejemplo scripts como este:
macro = ;
var:
numIts
numTrees
itTrees;

set numIts 100 ;
set itTrees 10 ;
set numTrees 'numIts' * 'itTrees' ;

rseed 0;
hold 'numTrees' ;
mult = replic 'numIts' hold 'itTrees' ;
nel * ;
p/ ;
Que no hacen nada que no se puede hacer directamente (y segura más fácil!), pero sirven para familiarizarse con el lenguaje ;).

Referencias
[1] Paradis, E. 2006. Analysis of phylogenetics and evolution with R. Springer

viernes, mayo 01, 2009

The TNT wiki

Just few days ago, the TNT wiki was released. The project is led by my friend Santi Catalano, Matt Yoder and Ximo Mengual. They meet at the Hennig Meeting, in Tucumán, with the idea of provide a friendly environment for TNT users.

This is the link: http://tnt.insectmuseum.org/

For the most part, the wiki is an on-line version of the program help. But as it is a wiki, it is possible to provide several examples and tricks.

I happy to help with the pages relevant to implied weights: IW and IW quick tutorial pages. And I hope to write somethings for the scripting pages :).

Pablo is not directly involved with the project. Nonetheless, from time to time he reviewed some content, and of course, he always encourage its development. :).

Also, there is a TNT google group: http://groups.google.com/group/tnt-tree-analysis-using-new-technology, open to questions from any user!

A behemoth analysis, and a wiki, this is a great week for TNT :).

Just in case, this is the TNT is a free software for phylogenetic analysis at any scale. It is sponsored by the Willi Hennig society, and can be downloaded here: http://www.zmuc.dk/public/phylogeny/TNT/

La wiki de TNT

Hace unos días, la wiki de TNT fue lanzada al publico. El proyecto es liderado por mi amigo Santi Catalano, Matt Yoder y Ximo Mengual. Ellos se encontraron en el Hennig meeting, en Tucumán, con la idea de proporcionar una ayuda amigable a los usuarios de TNT.

Este es el enlace: http://tnt.insectmuseum.org/

En su mayor parte, la wiki es una versión en linea de la ayuda del programa. Pero como es una wiki es posible agrandarla con muchos ejemplos y trucos.

Yo colabore con las páginas de pesado implicado: la pagina básica y un tutorial rápido. Y espero escribir algunas cosas para la sección de scripts :).

Pablo no esta directamente asociado con el proyecto. Pero, cada cierto tiempo, el revisa algo del contenido, y por supuesto, siempre a apoyado el desarrollo del sitio :).

Hay además un grupo de google de TNT: http://groups.google.com/group/tnt-tree-analysis-using-new-technology abierto a las preguntas de cualquier usuario!

El análisis de un behemot, y una wiki, esta es una gran semana para TNT :).

Por si acaso, TNT es un software gratuito para análisis filogenético a cualquier escala. Esta patrocinado por la Willi Hennig society, y puede descargarse aquí: http://www.zmuc.dk/public/phylogeny/TNT/

martes, abril 28, 2009

Filogenia de 73060 eucariotas


Finalmente, el Behemot vio la luz. Nuestro paper de un análisis de parsimonia de 73060 especies de eucariotas (y 7800 caracteres mol+morf) acaba de ser puesto en linea en Cladistics [doi: 10.1111/j.1096-0031.2009.00255.x].



Pablo hizo un trabajo brutal optimizando cada cosa en las búsquedas de arboles con TNT. Y todos aquí trabajaron intensamente para poder manejar toda esa cantidad de datos!

Al principio me sorprendió la enorme exactitud de los árboles encontrados, porque el set de datos esta lleno de entradas faltantes. Además, me puse muy contento porque al incluir los datos morfológicos, aún a esta escala tan grande, los resultados fueron mejores que usando únicamente datos moleculares!

Hace uno meses, esto fue posteado en dechronization:
[Cassy Dunn] hizo un llamado convincente acerca de la necesidad de técnicas análiticas revolucionarias ahora que entramos en una era en que los alcances de la computación serán más limitantes que la disponibilidad de datos.
Yo creo que nuestro trabajo muestra exactamente lo contrario: nuestras herramientas de búsqueda son lo suficientemente buenas, pero no hay datos suficientes (el conjunto de datos más grande es SSU con 20000 especies, y solo un puñado de genes tiene más de 10000 especies).

Otra lección?.... No necesitamos los super-árboles!

lunes, abril 27, 2009

A phylogeny of 73060 eukaryotes


Finally, the behemoth has seen the light :). Our paper with a parsimony analysis of 73060 eukariotic species (and 7800 mol+morf characters) was just published (as “online early”) in Cladistics [doi:10.1111/j.1096-0031.2009.00255.x].


Pablo does a wonderful work optimizing every aspect of the tree-searches in TNT. And all of the guys worked really hard to manage that amount of data!

At first I was surprised with the high accuracy of the trees founded, because the data set is full of missing entries. Also, I fill happy because the inclusion of morphological data, even at this huge scale, produce better results than molecules alone!

Just few months ago, this was posted in dechronization:
He [Cassey Dunn] makes a convincing case for the idea that a revolution in analytical techniques will be needed as we enter an era during which computational capabilities will be more limiting than data availability.
I think our study shows exactly the inverse: that our actual search capabilities are good enough, but we do not have sufficient data (the largest gene set is SSU with 20000 species, and a handful of genes has more than 10000 species).

The second lesson?... We do not need super-trees!

viernes, abril 24, 2009

Algorithms for phylogenetics 0a: Bitfields


Intro

I want to start a series of post about the algorithms used in phylogenetic analyses. I feel a bit disappointed with the entry on computational phylogenetics (and related subjects) in wikipedia, that are somewhat biased towards model-based methods and "bioinformatics", with a poor representation of algorithms for phylogenetics.

There are two outstanding examples of algorithms for phylogenetics on the web. The book of Wheeler et al. [1], and, in a similar vein, a manuscript by DeLaet [2]. But their presentation of the algorithms is somewhat general, that is be nice for educational purpouse, but in some cases, far away from the actual implementation in computer software.

As the presentation of the algorithms in [1] and [2] is excellent, I want to focus on the implementation. The idea is to go from simpler to more complex algorithms (as far as I can comprehend them xD).

I want to acknowledge Pablo Goloboff. He always have the time to discuss with me about algorithms (and any topic on methodological phylogenetics) :D. Of course, any error in my implementations are my own errors!

As this series progress, I will post the source code on google source or sourceforge, if someone wants to check it, requires an knowledge of C programming language. Thanks to Rubén who showme some tips with the HTML :).

Bitfields

Characters are a basic component of phylogenetic analysis. For simplicity assume that we only have a single character in each terminal. Each character can be stored into a variable, assigning to them 0, 1, 2... as indicated by the character state. This idea is implicit in the original description of Wagner's "groundpland divergence" [3], and in Farris' formalization [4, see 5]. Then operations between characters would be operations of ranges. But it is more simpler to think about characters as a set of states [4, 6]. Then operations between characters would be set operations. This set has a peculiar property: they are perfectly boolean, that is, a taxon has or has not a particular state. This can be translated in two positive consequences: (i) character states can be stored as a bitfield: (ii) character operations are bifield operations.

Computers store numbers as binary numbers. A bitfields use this characteristic to represent a set of positioned bits. This example shows it better:
State    binary representation       "Number" (human)
0 0001 1
1 0010 2
2 0100 4
3 1000 8
0, 1 0011 3
1, 3 0101 5
As can be seen, each state has his on bit position, and the states combination is just the union of its states. Many programming languages include operations to work directly on bits (not, and, or and xor), then the translation from set operations to bit operations is direct.

In this series I will use char to store character states. In C, char has 8 bits, then only 8 states can be stored. I define the bitfield with a typedef, so changing the number of bits on the bitfield (and then the number of states) is easy.
#define BITS_ON_STATE 8
#define ALL_BITS_ON 255
typedef char STATES;
This small example reads the caracters of a taxon, from a file (in), and store it into an array of characters (chars). The number of characters is nc, and the functions getc (in stdio.h) reads a single character from the file, SkipSpaces ignores the blank characters, and isdigit (ctype.h) detects if the character is a number.
for (j = 0; j < nc; ++ j) {
error = SkipSpaces (in);
if (error != NO_ERROR) return error;
c = getc (in);
if (isdigit (c))
chars [j] = 1 < < (c - '0');
else if ((c == '?') || (c == '-'))
chars [j] = ALL_BITS_ON;
else return BAD_FORMAT;
}
The operation chars [j] 1 < < (c – '0') substract the ASCII value of '0' (30) from the ASCII value stored in c (a number, so it is form 30 to 39). The resulting value is used to shift the bits of 1. For example if a '0' is readed the operations are:
chars [j] = 1 < < (30 – 30)
chars [j] = 1 < < 0
[0000 0001 < < 0 = 0000 0001]
chars [j] = 1
If a 5 is readed
chars [j] = 1 < < (35 – 30)
chars [j] = 1 < < 5
[0000 0001 < < 5 = 0010 0000]
chars [j] = 32
The 1 is shifted five bits to the left.

Note that inapplicables and unknowns are coded with all bits on.

This work nicely with non additive characters. For additive characters, it is possible to “recode it” as several binary characters [7][8], then it is possible to use the character as independent binary characters.

There are more sophisticated usages of bitfields, packing several characters in a single variable (for example, eight 4-bit characters can be stored into a singel 32-bits variable) [8][9]. This packing allows an speed up during searches, because several characters can be examined simultaneously.

References
[1] Wheeler, W. et al. 2006. Dynamic homology and phylogenetic systematics: a unified approach using POY. American Mus. of Natural History, published in cooperation with NASA. online: http://research.amnh.org/scicomp/pdfs/wheeler/Wheeler_etal2006b.pdf
[2] DeLaet, J. 2005. Pseudocode for some tree search algorithms in phylogenetics. Manuscript online: http://www.plantsystematics.org/publications/jdelaet./algora.pdf
[3] Wagner, W.H. 1961. Problems in the classification of ferns. In: Recent advances in botany. Toronto Univ. Press. pp. 841-844.
[4] Farris, J.S. 1970. Methods for computing Wagner trees. Syst. Zool. 19: 83-92.
[5] Goloboff, P.A. et al. 2006. Continuos characters analyzed as such. Cladistics 22: 589-601. doi: 10.1111/j.1096-0031.2006.00122.x
[6] Fitch, W.M. 1971. Toward defining the course of evolution: minimum change for a specific tree topology. Syst. Zool. 20: 406-416.
[7] Farris, J.S. et al. 1970. A numerical approach to phylogenetic systematics. Syst. Zool. 19: 172-189.
[8] Moilanen, A. 1999. Searching for most parsimonious trees with simulated evolutionary optimization. Cladistics 15: 39-50. doi: 10.1111/j.1096-0031.1999.tb00393.x
[9] Goloboff, P.A. 2002. Optimization of polytomies: state and parallel set operations. Mol Phyl. Evol. 22: 269-275. doi: 10.1006/mpev.2001.1049

Algoritmos en filogenia 0a: Campos de bits


Intro

Quiero iniciar una serie de posts sobre los algoritmos usados en análisis filogenético. El motivo principal fue ver el árticulo sobre “filogenética computacional” (y los tópicos de sistemática) en wikipedia, que pues me parecen algo sesgados hacía los métodos basados en modelos, y en general los algoritmos de filogenia no están muy representados.

Hay un par de ejemplos muy buenos en i-net de algoritmos de filogenia, uno es el libro de Wheeler et al. [1], y en una vena similar, un manuscrito de DeLaet [2]. Sin embargo, los algoritmos que presentan son más bien generales, lo cual es muy bueno desde el punto de vista didáctico, pero en varios casos, alejado de lo que hacen los programas en la actualidad!

Como esas fuentes son muy buenas, espero enfocarme más en el aspecto de la implementación de los algoritmos. La idea es pasar de los más sencillos a los más complejos (en la medida en la que yo pueda comprenderlos xD).

Antes de empezar, pues quiero agradecer a Pablo Goloboff, que siempre ha estado dispuesto a discutir de algoritmos conmigo :). Por supuesto, los errores que tengan mis implementaciones son mis errores!

A medida que avance, iré posteando el código fuente en google code o sourceforge, por si desean examinarlo, requiere que sepan algo de C. El HTML tiene problemas con los símbolos mayor-que y menor-que, por los que los mostrare en estos posts como texto, en vez de usar el símbolo. Gracias a Rubén por el tip del HTML :)!

Campos de bits

Un componente básico del análisis filogenéticos son los caracteres. Asumamos que solo tenemos un carácter para cada terminal. Uno bien podría guardar cada estado en una variable, por ejemplo asignandole 0, 1, 2... según corresponda. Una idea así esta implícita en la descripción inicial del “groundpland divergence” de Wagner [3], y más o menos, en la formalización de Farris [4, ver 5]. Las operaciones entre caracteres serían operaciones de intervalos. Pero es mucho más sencillo ver a los caracteres como un conjunto de estados [4, 6]. Las operaciones entre caracteres serian operaciones de conjuntos. Esos conjuntos tienen además una característica particular: son perfectamente booleanos, es decir, un taxon tiene o no un determinado estado. Esto tiene dos consecuencias positivas: (i) podemos guardar el conjunto de estados como un campo de bits; (ii) las operaciones entre caracteres son operaciones de los campos de bits.

Las compus guardan los datos como números binarios. En un campo de bits, usamos esa particularidad para representar conjuntos de bits. Un ejemplo, aclara mejor la situación:
Estado    Representación binaria    "Número" (humano)
0 0001 1
1 0010 2
2 0100 4
3 1000 8
0, 1 0011 3
1, 3 0101 5
Como se ve, cada estado individual tiene su propio bit, y la combinación de estados es simplemente la unión los estados. Muchos lenguajes de programación incorporan operaciones para trabajar directamente con bits (not, and, or y xor) con lo cual la traducción de operaciones de conjuntos a operaciones de bits es directa.

En esta serie voy a utilizar un char para guardar los estados de carácter. En C, char es de 8 bits. Yo lo defino con un typedef, de manera que después es posible cambiar el número de bits en el bitfield.
#define BITS_ON_STATE 8
#define ALL_BITS_ON 255
typedef char STATES;
En este pequeño ejemplo, se leen los caractes de un taxon, desde un archivo (in), y se guardan en un array de caracteres (chars). El número de caracteres es nc, y se usan las funciones getc para leer caracteres de texto, SkipSpaces para ignorar los espacios de texto, y isdigit para reconocer si el carácter de texto es un número.
for (j = 0; j < nc; ++ j) {
error = SkipSpaces (in);
if (error != NO_ERROR) return error;
c = getc (in);
if (isdigit (c))
chars [j] = 1 < < (c - '0');
else if ((c == '?') || (c == '-'))
chars [j] = ALL_BITS_ON;
else return BAD_FORMAT;
}
En la operación chars [j] = 1 < < (c – '0') se resta el valor ASCII de '0' (30) al valor ASCII guardado en c (es un número, va de 30 a 39). El valor de esa resta lo usamos para correr los bits de 1. Así por ejemplo, si se lee un '0', la operación seria:
chars [j] = 1 < < (30 – 30)
chars [j] = 1 < < 0
[0000 0001 < < 0 = 0000 0001]
chars [j] = 1
Si leemos un 5
chars [j] = 1 < < (35 – 30)
chars [j] = 1 < < 5
[0000 0001 < < 5 = 0010 0000]
chars [j] = 32
Corrimos el 1 cinco bits a la izquierda.

Observen que los no aplicables o desconocidos, se codifican con todos los bits encendidos.

Esto funciona bastante bien para caracteres no aditivos. Para caracteres aditivos, es posible "recodificarlos" como varios caracteres binarios [7][8], una vez recodificados es posible utilizarlos como si fueran caracteres no aditivos independientes.

Hay formas más sofisticas de usar los campos de bits, uniendo varios caracteres en una sola variable (por ejemplo, en una variable de 32 bits, se pueden usar 8 caracteres de 4 bits, como nucleotidos) [8][9]. Ese empaquetado permite una evaluación muchisimo más rápida de los caracteres durante las búsquedas, puesto que en cada operación se pueden revisar simultáneamente varios caracteres.

Referencias
[1] Wheeler, W. et al. 2006. Dynamic homology and phylogenetic systematics: a unified approach using POY. American Mus. of Natural History, published in cooperation with NASA. online: http://research.amnh.org/scicomp/pdfs/wheeler/Wheeler_etal2006b.pdf
[2] DeLaet, J. 2005. Pseudocode for some tree search algorithms in phylogenetics. Manuscrito online: http://www.plantsystematics.org/publications/jdelaet./algora.pdf
[3] Wagner, W.H. 1961. Problems in the classification of ferns. En: Recent advances in botany. Toronto Univ. Press. pp. 841-844.
[4] Farris, J.S. 1970. Methods for computing Wagner trees. Syst. Zool. 19: 83-92.
[5] Goloboff, P.A. et al. 2006. Continuos characters analyzed as such. Cladistics 22: 589-601. doi: 10.1111/j.1096-0031.2006.00122.x
[6] Fitch, W.M. 1971. Toward defining the course of evolution: minimum change for a specific tree topology. Syst. Zool. 20: 406-416.
[7] Farris, J.S. et al. 1970. A numerical approach to phylogenetic systematics. Syst. Zool. 19: 172-189.
[8] Moilanen, A. 1999. Searching for most parsimonious trees with simulated evolutionary optimization. Cladistics 15: 39-50. doi: 10.1111/j.1096-0031.1999.tb00393.x
[9] Goloboff, P.A. 2002. Optimization of polytomies: state and parallel set operations. Mol Phyl. Evol. 22: 269-275. doi: 10.1006/mpev.2001.1049

lunes, marzo 30, 2009

doptores

Esta semana vi la defensa de las tesis doctorales de dos chicos aquí de Tucumán. Guillermo Suarez que trabaja con musgos, y pues mostró un enorme trabajo sobre musgos. Es tremendo lo que hay que hacer para poder sacar caracteres morfológicos en plantas! Además, pues un excelente trabajo de taxonomía y filogenética combinado :). También se presento Sebastian Barrionuevo, el trabaja con unas ranas especialmente acuáticas, muchas de ellas viven en grandes alturas, por lo que es muy complicado muestrearlas. Una excelente labor de morfología, y de unión de la filogenia con otros campos, como la ontogenía y de alguna manera la ecología.

También se presentaba Santiago Catalano en BsAs, pero pues no me pude comunicar con el, además de que estaba muy cansado, pues fue al otro día del toque de Radiohead (me hubiera quedado dormido! xP). Y también supe de un chico, Ignacio Escapa, de Trelew (muuuy lejos xD) que presento uno de los posters que más me gusto en el meeting de la Willi Hennig.

Tanto Santí como Sebastian ganaron premios en el meeting ;).

Esperemos que cada día sea mayor la cantidad de gente que trabaja en cladística :D, y que el interes por la morfología siga siendo un área fuerte!

sábado, marzo 21, 2009

Just some links

There are some nice things on line lately :).

First, Rod Page talk at NHM of London, he puts his slide-show on line, there is also a video and a pod-cast of the talk (I do not see them :P). The slide-show can give a nice idea about the talk. There are two important things on the talk: (a) the importance of the availability of data! I just agree with him: data must be freely available, and previously published data must be easy accessible; and (b) Scientific data must be readily usable, for example the “Encyclopedia of life” is just a fun site just like wikipedia or the web tree of life, they are full of nice pics and info, but cientifically irrelevant: no hard data are attached to it (i.e. morfological info, character matrices), this contrast with the most simpler ispaces of Page!

Although I do not like some of the comments, Malte Ebach made a good point against the revival of a “pragmatic” classification (for horticulture, for example). I think that the only rule of classification is the phylogeny. If an immutable classification is the objective, the only way I can thing is an “alphabetic” or “numeric” system, but the preservation of “traditional” names is not the solution. Maybe, it is better a change of the rules...

And speaking of change of the rules, I recommend the continuos checking of several of Mike Keesey posts. I do not think that phylocode is the solution (in fact I think that it is even worse than traditional Linnean classification), but he writes interesting things about databasing!

Also, it is seems to change some rules about the publication of names for taxonomy (via Evolving thoughts). I thing that only peer reviewed publication must count, at least from the last fifty years or so. Also I thing that new names, at least for species must be published on journals instead of books...