macro=; 
/*
   This macro performs the support weighting of Farris (2001)

*/ 


/* SOURCE: http://www.lillo.org.ar/phylogeny/tnt/scripts/supwt.run */


if ( argnumber < 2 )
  errmsg Must give two numbers:
      1   N  N = num. replications
      2   N  N = max. weighting rounds
     [3]  +  weight against homoplasy 
             (instead of support weighting)
     ; end
report- ;
macfloat 5 ;
var :
  avgwt [ ( nchar + 1 ) ] 
  oldwts [ ( nchar + 1 ) ] 
  mkup [ nchar ]
  mkdn [ nchar ]
  baswts [ ( nchar + 1 ) ] 
  numup
  numdn
  i j k
  search [ 50 ]
  finalsearch [ 50 ] 
  thischgs
  numchgs
  bufsta [ ( 2 * ( ntax + 1 ) ) ]
  ints
  firsround
  stabilized
  nunos
  nutrees
  nuwt
  numreplics
  maxwtrounds
  dosucwt
  maxwtis
; 
/****
     Here, set search and final search algorithms
*****/
set search                $ mult 5 = hold 3 .,;
set finalsearch           $ ho200.,mu20=hol10.,bb., ; 


set dosucwt 0 ;
if ( argnumber > 2 )
   if ( eqstring [ %3 + ] )
      set dosucwt 1 ;
   else 
      errmsg Wrong argument: %3 ; 
      end
   end 

set numreplics              %1 ;
set maxwtrounds             %2 ;
set maxwtis                100 ;   /*  maximum possible weight */
                                   /*  in the range 1 - maxwtis */ 

if ( 'numreplics' > 1 ) 
   tsave * wtsup.tre ; 
   end 
collapse rule 3 ;
lquote [ ;
ckeep ;

loop 0 nchar 
   set avgwt[ #1 ] 0 ; 
   set baswts [ #1 ] weight [ #1 ] ; 
stop 

loop 1 'numreplics'
     rseed * ; 
     macfloat 0 ;
     /*   First, resample... */
     set numup 0 ;
     set numdn 0 ;
     loop 0 nchar
        if ( ! isact [ #2 ] ) continue ; end 
        if ( weight [ #2 ] == 0 ) continue ; end
        if ( ! isinfo [ #2 ] ) continue ; end 
        set i getrandom [ 1 100 ] ;
        if ( 'i' > 66 ) continue ; end
        set i getrandom [ 0 1 ] ;
        if ( 'i' )
             set mkup [ 'numup' ] #2 ; set numup ++ ; 
        else set mkdn [ 'numdn' ] #2 ; set numdn ++ ; end
        stop
     ccode ! ;   /*  ...start from original ccode  */
     if ( 'numup' )
         set numup -- ; 
         ccode /2 'mkup[0-'numup']' ;
         end
     if ( 'numdn' )
         set numdn -- ; 
         ccode ] 'mkdn[0-'numdn']' ;
         end

     /*  Then search until stabilization...    */
     set firsround 1 ; 
     loop 1 'maxwtrounds' 

       set stabilized 1 ;
       if ( 'firsround' )
          set stabilized 0 ;
          end 
       $search ;  
       if ( 'firsround' )
           ccode [ 'mkdn[0-'numdn']' ;
           end

       if ( !'dosucwt' )
          if ( ( ntrees + 1 ) == maxtrees )
             set i ntrees ;
             keep 'i' ;
             end 
          nelsen * ;
          set nunos tnodes[ ntrees ] ;
          set nutrees ntrees ;
          keep 'nutrees' ;
          end

       if ( ntrees > 9 )
          tgroup =0 *10 ;
          tchoose { 0 }  ;
          end
                /*   Nesting:
                   1 - outer iteration (resampling)
                   2 - inner iteration (rounds of succ wt)
                   3 - char 
                   4 - tree
                   5 - node    */ 

       loop 0 nchar set oldwts[ #3 ] weight[ #3 ] ; stop 
       ccode ! ; 

       loop 0 nchar
          if ( ! isact [ #3 ] ) continue ; end
          if ( 'oldwts [ #3 ]' == 0 ) continue ; end
          if ( ! isinfo [ #3 ] ) continue ; end 
         
          macfloat 5 ;

          if ( !'dosucwt' ) 
                            /*  weight with support weighting */
             set numchgs 0 ;
             loop 0 ntrees
                set thischgs 0 ;
                set bufsta states[ #3 . #4 ] ; 
                loop ( ntax + 2 ) nnodes[#4]
                   if ( ( 'bufsta[ #5 ]' & 'bufsta[ anc [ #4 #5 ] ]' ) == 0 )
                      set thischgs ++ ;
                      end
                   stop
                if ( 'thischgs' > 'numchgs' )
                   set numchgs 'thischgs' ;
                   end
                stop
             set nuwt ( 'maxwtis' * 'numchgs' ) / 'nunos' ;  

          else 
                             /*  re-weight with inverse of homoplasy,
                                 with concavity 5  */
             ccode /1 #3 ; 
             loop 0 ntrees
                set numchgs 10000000 ;
                set thischgs length[ #4 #3 ] ;
                if ( 'thischgs' < 'numchgs' )
                    set numchgs 'thischgs' ;
                    end
                stop 
             set numchgs -= minsteps [ #3 ] ; 
             set nuwt ( 'maxwtis' * 2.5 ) / ( 2.5 + 'numchgs' ) ;

          end

          macfloat 0 ;
          if ( 'nuwt' == 0 ) set nuwt 1 ; end
          if ( 'nuwt' > 'maxwtis' ) set nuwt 'maxwtis' ; end 
          set nuwt *= 'baswts [ #3 ]' ;       /*  take into account   */
          if ( 'nuwt' != 'oldwts[ #3 ]' )     /*  user-defined weights!  */ 
              set stabilized 0 ;
              end 
          ccode /'nuwt' #3 ;   /*  re-set weight!  */ 
          stop

       if ( 'stabilized' ) endloop ; end 
       set firsround 0 ;
       stop

     /*  Now, search under final weights  */ 
     $finalsearch ; 
     if ( ( ntrees + 1 ) == maxtrees )
          set i ntrees ;
          keep 'i' ;
          end 
     nelsen * ;
     if ( 'numreplics' > 1 ) 
          save / ;
     else tchoose / ; 
          end 

     loop 0 nchar
       set avgwt [ #2 ] += weight [ #2 ] ;
       stop 

     stop
progress/;
ccode ! ; 
if ( 'numreplics' == 1 ) proc/; end 
keep 0 ;
tsave/; 
proc wtsup.tre ;
macfloat 2 ; 
collapse notemp ; 
loop 0 nchar
   set avgwt[ #1 ] /= 'numreplics' ;
   stop
if ( !'dosucwt' ) quote RESULTS OF SUPPORT WEIGHTING : ;
else quote RESULTS OF WEIGHTING AGAINST HOMOPLASY ; end 
major = 60 ;
major = 60 * ;
maketable avgwt Average character weights:; 
tchoose/; 
set j 0 ; 
set k 0 ; 
loop 0 nchar
   set j += 'avgwt [#1 ]' ; 
stop 
quote Average weight of all chars. = 'j' ; 
proc/;