macro=;
if ( argnumber < 2 )
  errmsg
  This scripts merges two trees (mixing groups of as
  even a size as possible).  Must specify two trees
  to use as source for new tree. ;
end
var:
   i j 
   thenod
   nodlist[2] 
   reftree
   daroot
   nodXinA
   nodYinB
   nodYinA
   nodXinB
   tmp
   mindiff
   treeA
   treeB
   partX
   partY
;

/****   First, find a shared node that
        partitions both trees evenly ****/

set thenod 0 ;
set mindiff ntax ; 
collapse notemp ;
tgroup =0 (nothing) ;  /*  just make sure consensus goes to group 1 */
nelsen * %1 %2 ;
loop (ntax+2) nnodes[ntrees]
    set tmp numdes[ntrees #1] - ( ( ntax + 1 ) / 2 ) ;
    if ( 'tmp' < 0 ) set tmp ( - 'tmp' ) ; end
    if ( 'tmp' >= 'mindiff' ) continue ; end
    set mindiff 'tmp' ; 
    set thenod #1 ;
    set partX numdes[ntrees #1] ; 
stop
macfloat 0 ; 
set partY ( ntax+1 ) - 'partX' ; 
quote --------------------------------------------------------- ; 
quote DIVIDING TREES AT NODE 'thenod' OF CONSENSUS ('partX'-'partY') ;
quote --------------------------------------------------------- ; 


/**** Second, reroot source trees (as
      temporary copies) in such a way
      that both partitions are monophyletic
      and identify equivalent nodes   *******/

set reftree ntrees ;
set daroot root ;
set nodXinA eqgroup[ 'reftree' 'thenod' %1 ] ; 
set nodXinB eqgroup[ 'reftree' 'thenod' %2 ] ; 
edit ] 'reftree' 'daroot' 'daroot' 'thenod' ;
copytree %1 ; 
set treeA ntrees ;
copytree %2 ; 
set treeB ntrees ; 
edit ] 'treeA' 'daroot' 'daroot' 'nodXinA' ;
edit ] 'treeB' 'daroot' 'daroot' 'nodXinB' ;
set nodlist deslist[ 'reftree' 'daroot' ] ;
set nodXinA eqgroup[ 'reftree' 'nodlist[0]' 'treeA' ] ;
set nodYinA eqgroup[ 'reftree' 'nodlist[1]' 'treeA' ] ; 
set nodXinB eqgroup[ 'reftree' 'nodlist[0]' 'treeB' ] ; 
set nodYinB eqgroup[ 'reftree' 'nodlist[1]' 'treeB' ] ; 


/****  Now, mix the source trees,
       placing the new trees in a
       group of its own, and source
       trees in a separate group.  ********/
tgroup noauto ; 
tread
   ( +'treeA' 'nodXinA' +'treeB' 'nodYinB' ) * 
   ( +'treeA' 'nodYinA' +'treeB' 'nodXinB' ) ;
set i ntrees ;
set j ntrees - 1 ; 
tgroup =0 (trees_%1_%2_mixed ) 'i' 'j' ;
tgroup =1 (source_trees) %1 %2 ; 
tgroup auto ;


/****  Finally, get rid of the
       temporary trees and reroot
       all trees uniformly **********/

tchoose +. - 'reftree' 'treeA' 'treeB' ;
reroot ; 

tsave pis ;
save { source_trees } ;
tsave/;
tsave caca ;
save { trees_%1_%2_mixed } ;
tsave/;

tsave pis + ; 
keep 0 ;
sho caca ;
tc 0 ; bb = nom ;  sav ; 
keep 0 ;
sho caca ;
tc 1 ; bb = nom ; sav ;

tsave/;
keep 0 ; 
sho pis ;
le ; 

if ( linux ) 
   system rm pis caca^;
else
   system del pis caca^;
end 

proc/;