/************************************************************\ | A didactic script - Most Parsimonious Reconstructions | | by P.A. Goloboff (2021) | | | | Options: | | | | -v make run non-verbose (which saves a looot | | of output) | | | | The code is designed for legibility and clarity more | | than efficiency. | | | | Feel free to use and modify. | | | \************************************************************/ int state_sets[ Nnodes(0) ] , chosen_state[ Nnodes(0) ] ; int list_of_nodes [ Nnodes(0) ] , verbose = 1 ; int transf_cost [ 32 ] [ 32 ] , num_reconsts = 0 , total_length = 0 , target_length ; void fill_transf_costs ( int i ) { int j , k ; if ( iscont ( i ) || islmark ( i ) ) error ( "This code is for discrete characters only!!!" ) ; if ( issank ( i ) ) for ( j = 0 ; j <= maxstate( i ) ; ++ j ) for ( k = 0 ; k <= maxstate( i ) ; ++ k ) transf_cost [ j ] [ k ] = cost ( i , j , k ) ; else if ( isadd ( i ) ) for ( j = 0 ; j <= maxstate( i ) ; ++ j ) for ( k = j ; k <= maxstate( i ) ; ++ k ) transf_cost [ j ] [ k ] = transf_cost [ k ] [ j ] = k - j ; else for ( j = 0 ; j <= maxstate( i ) ; ++ j ) for ( k = j ; k <= maxstate( i ) ; ++ k ) transf_cost [ j ] [ k ] = transf_cost [ k ] [ j ] = 1 ; } int cost_from_to ( int state , int left , int right ) { int w , x , y , z , i , this , lbest , rbest ; i = -1 ; x = chosen_state [ left ] ; z = chosen_state [ right ] ; w = 1 << state ; if ( ( w & x ) && ( w & z ) ) this = 0 ; else { lbest = rbest = 10e6 ; if ( ( w & x ) ) lbest = 0 ; if ( ( w & z ) ) rbest = 0 ; for ( y = 1 ; ( y <= x || y <= z ) && ++ i < 32 ; y <<= 1 ) { if ( ( y & w ) ) continue ; if ( ( y & x ) ) if ( transf_cost [ state ] [ i ] < lbest ) lbest = transf_cost [ state ] [ i ] ; if ( ( y & z ) ) if ( transf_cost [ state ] [ i ] < rbest ) rbest = transf_cost [ state ] [ i ] ; } this = rbest + lbest ; } return this ; } void recurse_reconstructions ( int at ) { int curnod , total_lengthwas = total_length ; int i , x , y , left , right ; if ( at == listsize ) { if ( verbose ) { tnt ( "ttag>Rec. %i; ttag; " , num_reconsts ) ; chkbreak () ; } num_reconsts ++ ; return ; } curnod = list_of_nodes [ at ] ; left = firstdes( 0 , curnod ) ; right = sister ( 0 , firstdes ( 0 , curnod ) ) ; x = state_sets [ curnod ] ; if ( numbits ( x ) == 1 ) { for ( i = 0 ; i < 32 ; ++ i ) if ( ( x & ( 1 << i ) ) ) break ; total_length += cost_from_to ( i , left , right ) ; if ( total_length > target_length ) { total_length = total_lengthwas ; continue ; } chosen_state [ curnod ] = state_sets [ curnod ] ; tagset ( curnod , name ( "bitset" , x ) ) ; recurse_reconstructions ( at + 1 ) ; total_length = total_lengthwas ; return ; } i = -1 ; for ( y = 1 ; y <= x && ++ i < 32 ; y <<= 1 ) { if ( ! ( x & y ) ) continue ; // state does not occur in state set total_length += cost_from_to ( i , left , right ) ; if ( total_length > target_length ) { total_length = total_lengthwas ; continue ; } chosen_state [ curnod ] = y ; tagset ( curnod , name ( "bitset" , y ) ) ; recurse_reconstructions ( at + 1 ) ; total_length = total_lengthwas ; } } void main ( int argc , char ** argv ) { int i , begin , finish ; if ( argc < 2 ) error ( "Must indicate character to reconstruct (always on tree 0)" ) ; i = atoi ( argv [ 1 ] ) ; target_length = length ( 0 , i ) ; // get target length and full state sets from TNT... easiest. states ( state_sets , 0 , i ) ; fill_transf_costs ( i ) ; if ( argc > 2 ) if ( !strcmp ( argv [ 2 ] , "-v" ) ) verbose = 0 ; downlist ( list_of_nodes , 0 ) ; tagtree ( 0 ) ; for ( i = 0 ; i < Ntax ; ++ i ) { chosen_state [ i ] = state_sets [ i ] ; // terminals never modified! if ( state_sets [ i ] == missing ) tagset ( i , "?" ) ; else tagset ( i , name ( "bitset" , state_sets [ i ] ) ) ; } time ( &begin ) ; recurse_reconstructions ( 0 ) ; time ( &finish ) ; tntprintf ( "TOTAL: %i reconstructions at %i steps (%.1f sec).\n" , num_reconsts , target_length , difftime ( finish, begin ) ) ; }