/***********************************************************\ | A didactic script - Generalized Optimization | | by P.A. Goloboff (2021) | | | | Options: | | | | -v skip mapping for the characters | | | | The script assumes you have read the data set and have | | trees in memory; it operates always on the first tree | | in memory. It can handle only additive and nonadditive | | characters (Sankoff characters trigger an error). | | This script uses one of the worst possible methods for | | finding optimal character state reconstructions; it is | | intended only to illustrate some general principles. | | | | Feel free to modify and use. | | | \***********************************************************/ int charnum , total_length ; unsigned int state_sets [ nnodes( 0 ) + 1 ] ; int tmpcost [ 32 ] ; /******* Non-Additive optimization ***************\ | (trying states one-by-one) | \**************************************************/ void nonadditive_down ( int node ) { unsigned int x , y , z , newstate ; int i , thiscost , best = 10e6 ; y = state_sets [ firstdes( 0 , node ) ] ; z = state_sets [ sister ( 0 , firstdes( 0 , node ) ) ] ; x = 1 ; for ( i = 0 ; ( x <= y || x <= z ) && i < 31 ; x <<= 1 ) { ++ i ; thiscost = 0 ; if ( !( x & y ) ) thiscost = 1 ; if ( !( x & z ) ) ++ thiscost ; if ( thiscost < best ) { newstate = x ; best = thiscost ; } else if ( thiscost == best ) newstate |= x ; } state_sets [ node ] = newstate ; total_length += best * weight ( 0 ) ; } void nonadditive_up ( int node ) { unsigned int w , x , xa , y , z , newstate ; int i , j , thiscost , best ; if ( node == root ) { tagset ( node , name ( "bitset" , state_sets [ node ] ) ) ; return 1 ; } y = state_sets [ firstdes( 0 , node ) ] ; z = state_sets [ sister ( 0 , firstdes( 0 , node ) ) ] ; w = state_sets [ anc ( 0 , node ) ] ; xa = 1 ; state_sets [ node ] = 0 ; for ( j = 0 ; xa <= w && j < 31 ; xa <<= 1 ) { // Try every optimal state in ancestor ++ j ; if ( ( xa & w ) == 0 ) continue ; best = 10e6 ; x = 1 ; for ( i = 0 ; ( x <= y || x <= z || x <= xa ) && i < 31 ; x <<= 1 ) { ++ i ; thiscost = 0 ; if ( !( x & y ) ) thiscost = 1 ; if ( !( x & z ) ) ++ thiscost ; if ( !( x & xa ) ) ++ thiscost ; if ( thiscost < best ) { newstate = x ; best = thiscost ; } else if ( thiscost == best ) newstate |= x ; } state_sets [ node ] |= newstate ; // and add all optimal states to the final state sets } tagset ( node , name ( "bitset" , state_sets [ node ] ) ) ; } /********* Additive optimization ****************\ | (trying states one-by-one) | \**************************************************/ int distance_between ( unsigned int x , unsigned int y ) // Use the naive method (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 additive_down ( int node ) { unsigned int x , y , z , n , newstate ; int i , d , thiscost , best = 10e6 , at ; y = state_sets [ firstdes( 0 , node ) ] ; z = state_sets [ sister ( 0 , firstdes( 0 , node ) ) ] ; x = 1 ; for ( i = 0 ; ( x <= y || x <= z ) && i < 31 ; x <<= 1 ) { ++ i ; at = firstdes( 0 , node ) ; thiscost = 0 ; if ( !( x & y ) ) thiscost = distance_between ( x , y ) ; if ( !( x & z ) ) thiscost += distance_between ( x , z ) ; if ( thiscost < best ) { best = thiscost ; newstate = x ; } else if ( thiscost == best ) newstate |= x ; } total_length += best * weight ( 0 ) ; state_sets [ node ] = newstate ; } void additive_up ( int node ) { unsigned int w , x , xa , y , z , n , newstate ; int i , j , thiscost , best ; if ( node == root ) { tagset ( node , name ( "bitset" , state_sets [ node ] ) ) ; return 1 ; } y = state_sets [ firstdes( 0 , node ) ] ; z = state_sets [ sister ( 0 , firstdes( 0 , node ) ) ] ; w = state_sets [ anc ( 0 , node ) ] ; xa = 1 ; state_sets [ node ] = 0 ; for ( j = 0 ; xa <= w && j < 31 ; xa <<= 1 ) { // Try every optimal state in ancestor ++ j ; if ( ( xa & w ) == 0 ) continue ; best = 10e6 ; x = 1 ; for ( i = 0 ; ( x <= y || x <= z || x <= xa ) && i < 31 ; x <<= 1 ) { ++ i ; thiscost = 0 ; if ( !( x & y ) ) thiscost += distance_between ( x , y ) ; if ( !( x & z ) ) thiscost += distance_between ( x , z ) ; if ( !( x & xa ) ) thiscost += distance_between ( x , xa ) ; if ( thiscost < best ) { newstate = x ; best = thiscost ; } else if ( thiscost == best ) newstate |= x ; } state_sets [ node ] |= newstate ; // and add all optimal states to the final state sets } tagset ( node , name ( "bitset" , state_sets [ node ] ) ) ; } void copy_term_states ( void ) // Label terminal branches with their corresponding states { int i ; unsigned int x ; for ( i = 0 ; i < Ntax ; ++ i ) { states ( &x , -1 , charnum , i ) ; if ( x == missing ) tagset ( i , "?" ) ; else tagset ( i , name ( "bitset" , x ) ) ; }} void main ( int argc , char ** argv ) { int domaps = 1 ; int start , finish ; time ( &start ) ; if ( argc > 1 ) if ( !strcmp ( argv [ 1 ] , "-v" ) ) domaps = 0 ; total_length = 0 ; if ( domaps ) tagtree ( 0 ) ; for ( charnum = 0 ; charnum < Nchar ; ++ charnum ) { if ( domaps ) copy_term_states () ; states ( state_sets , -1 , charnum ) ; if ( issank ( charnum ) ) error ( "Cannot optimize Sankoff characters (%i)\n" , charnum ) ; if ( isadd ( charnum ) ) { travtree ( "down" , "additive_down" , 0 ) ; if ( domaps ) travtree ( "up" , "additive_up" , 0 ) ; } else { travtree ( "down" , "nonadditive_down" , 0 ) ; if ( domaps ) travtree ( "up" , "nonadditive_up" , 0 ) ; } if ( domaps ) { tnt ( "ttag>MAPPING FOR CHAR. %i;\n" , charnum ) ; tnt ( "ttag;" ) ; }} time ( &finish ) ; tntprintf ( "Length %i\nTime: %.1f sec\n" , total_length , difftime ( finish , start ) ) ; }