/************************************************************\ | A didactic script - Optimization by (almost) Brute | | Force Enumeration of Character-State Reconstructions | | by P.A. Goloboff (2021) | | | | Options: call it without arguments | | | | This shows how optimization had to be done before Farris | | and Fitch invented the first algorithms for optimization. | | | | It produces correct results, but the time to run it | | increases very rapidly with taxa and number of states | | | | Feel free to use and modify. | | | \************************************************************/ int state_sets[ nnodes(0) + 1 ] ; int final_state_sets [ nnodes(0) + 1 ] ; int node_list [ Ntax ] ; int curnode , curlength , best_length , chartodo , number_of_states ; int distance_between ( unsigned int x , unsigned int y ) // Use naive method to do distances (this is supossed to be naive code) { int dis = 0 ; if ( x < y ) while ( !( x & y ) ) { x <<= 1 ; ++ dis ; } else while ( !( x & y ) ) { y <<= 1 ; ++ dis ; } return dis ; } void recurse_nodes ( int at ) // the meat of the algorithm { int curnod = node_list [ at ] , initial_length = curlength ; unsigned int x , l , r ; int i , j ; l = state_sets [ firstdes ( 0 , curnod ) ] ; r = state_sets [ sister ( 0 , firstdes ( 0 , curnod ) ) ] ; x = 1 ; for ( i = 0 ; i < number_of_states ; ++ i ) { curlength = initial_length ; if ( isadd ( chartodo ) ) { if ( !( x & l ) ) curlength += distance_between ( x , l ) ; if ( !( x & r ) ) curlength += distance_between ( x , r ) ; } else { if ( !( x & l ) ) ++ curlength ; if ( !( x & r ) ) ++ curlength ; } state_sets [ curnod ] = x ; if ( curlength <= best_length ) // continue with rest of reconstruction only if still within best if ( curnod == root ) { // reconstruction is complete; store it! if ( curlength < best_length ) { best_length = curlength ; for ( j = nnodes(0) + 1 ; j -- > root ; ) final_state_sets [ j ] = 0 ; } for ( j = nnodes(0) + 1 ; j -- > root ; ) final_state_sets [ j ] |= state_sets [ j ] ; } else recurse_nodes ( at + 1 ) ; // reconstruction still incomplete; go to next node x <<= 1 ; }} void report_and_map ( int start ) { int i , finish ; tntprintf ( "Best length: %i\n" , best_length ) ; tagtree ( 0 ) ; for ( i = 0 ; i < Ntax ; ++ i ) if ( state_sets[ i ] == missing ) tagset ( i , "?" ) ; else tagset ( i , name ( "bitset" , state_sets [ i ] ) ) ; for ( i = root ; i <= nnodes(0) ; ++ i ) tagset ( i , name ( "bitset" , final_state_sets [ i ] ) ) ; time ( &finish) ; tnt("ttag>MAPPING FOR CHAR. %i; ttag;" , chartodo ); tntprintf ( "Time: %.1f sec.\n" , difftime ( finish , start ) ) ; } void main ( int argc , char ** argv ) { int chartodo ; int start ; time ( &start ) ; if ( argc < 2 ) error ( "\n ------------------------------------------------------------------------------- Script to optimize trees by (almost) brute force enumeration\n Time to complete optimization increases with ca. S^T (where S = number of states, and T = number of taxa). Branch-and-bound used to cut-down reconstructions, which saves some time, but not much. Must indicate character to map (always mapped on tree 0); second argument [optional] can be \"-v\" (skip mapping) -------------------------------------------------------------------------------\n" ) ; chartodo = atoi ( argv [ 1 ] ) ; if ( issank ( chartodo ) ) error ( "Sorry, only Fitch/Farris characters can be optimized with this script" ) ; downlist ( node_list , 0 ) ; states ( state_sets , -1 , chartodo ) ; number_of_states = maxstate ( chartodo ) + 1 ; curlength = 0 ; best_length = 10e6 ; recurse_nodes ( 0 ) ; report_and_map ( start ) ; }