sil=all ; 
screen 39x80 ;

if ( ntrees < 0 )   
   errmsg Must read data (and trees) before running ;
   end 
macro *10 ( ( ( 2*ntax+2 ) * 180 ) + 1000 ) ; 
macro [ ( ( ( ( ( 2*ntax+2 ) * 400 ) * 4 + 1000 ) / 1000 ) + 5000 ) ; 
macro=;


/***

   This script is a refinement of the original script of
   Giannini & Goloboff (2010), delcor.run.

   It is intended for use in large data sets where there
   are too many reconstructions to be handled by the
   original script.

   This version also adds the option of saving the reconstructions
   for later use (so that they do not have to be generated again),
   and the option of using fewer significant digits and/or
   averages for the terminal taxa (which may speed up enumeration
   of reconstructions).

   The original script is retained for documentation
   purposes.

***/


macfloat 6 ;

/************************************************************\
/************************************************************\
|                                                            |
|        M A I N    B O D Y    O F    S C R I P T            |
|                                                            |
|      This contains the main set of instructions, and       |
|        (with "goto") calls the rest of the routines        |
|                that are used in the script.                |
|                                                            |
\************************************************************/

var:
   numsigdigs averagit
   rdsvtype savefilename[ 100 ]  
   signif matches logged dograph userlengths 
   reccycles numstordx numstordy notest 
   reclo rechi testdown maxstorrays 
   twotailed incx [ (2*ntax+2) ] incy [ (2*ntax+2) ]
   visitd [ (2*ntax+2) ] distosrc [ (2*ntax+2) ]
   radius randwts lor hir 
   reconum sample numrepls 
   outnod innod tmpradius
   usertree ispos
   userx  worsr 
   usery
   minincx minincy 
   myx myy
   drag [ (2*ntax+2) ]
   xpairs [ (2*ntax+2) ]
   ypairs [ (2*ntax+2) ]
   blength [ (2*ntax+2) ]
   numpairs 
   wts [ (2*ntax+2) ]
   delfac didpair stasis dax day i ;

/********************

   Explanation of main variables used in script: 

	dograph : when set to N, plots pairs for set of reconstructions number N
	maxstorrays : maximum number of reconstructions to store for X and Y
	stordx[] & stordy[] : stored arrays for increments (to facilitate randomization,
		and to make sure randomization uses a random sample of reconstructions)  
	notest : skip the randomization
	testdown : test as tree-down.  That is, check whether change Y is always
		preceded by a change in X.  Alternative (default) is checking 
		whether a change in X always implies a change in Y. 
   	twotailed : make test two-tailed
   	radius : how many nodes to go up (or down) to search for a match
	randwts : whether to use a randomized radius value 
   	reccycles : if using randwts, how many cycles to do for observed correlation
	lor & hir : value of lowest/highest correlation 
   	reclo & rechi : number of reconstructions which produces lowest/highest correlation
   	usertree : number of tree to use for testing
   	userx & usery : number of characters, X and Y, to test  
   	minincx & minincy : threshold of change in X and Y to be considered 
   	blength [] : array to store branch lengths provided by user
	userlengths :  use branch lengths provided by user (otherwise, all unity) 
   	delfac : delay factor, in formula to set weight as 1 / ( 1 + ( Dist x delfac ) ) 
	stasis : 

**********/

goto = %0 ;
goto READPARMS %1 %2 %3 %4 %5 %6 %7 %8 %9 %(10) %(11)
   %(12) %(13) %(14) %(15) %(16) %(17) %(18) %(19) %(20) ; 

goto CHECKPARMS ; 

set maxstorrays 'sample' ; 
var: 
   stordx [ 'maxstorrays' (2*ntax+2) ]
   stordy [ 'maxstorrays' (2*ntax+2) ]
;

report- ;
lquote=;
lquote[ ;
naked=;

if ( 'userlengths' )
    quote &10SETTING USER BRANCH LENGTHS ...&10 ; 
                 /********************************************************** 
                     Read user-supplied branch lengths, when so requested  */
    loop 0 nnodes['usertree']
       if ( #1 == root ) continue ; end 
       if ( eqstring[ x$ttag #1 x ] )
          errmsg Sorry, branch nr. #1 has no length defined... &10;
          end 
       set blength[#1]  ( $ttag #1 ) ;
       stop
    quote &32    ...all branch lengths set successfully.&10; 
else
           /*************************************************************
              When no user-supplied branch lengths, all branches unity  */
    loop 0 nnodes['usertree'] set blength[#1] 1 ; stop
    end 

if ( !exstatus ) proc/; end   

if ( !'randwts' ) set reccycles 1 ; end 

if ( !'testdown' ) 
     goto CORRELATE_UP ;
else goto CORRELATE_DOWN ; end

if ( 'dograph' || 'notest' ) proc/; end
if ( 'numrepls' == 0 ) proc/; end 

if ( ( 'lor' <= 0 ) && ( 'hir' >= 0 ) )
   sil - all ; 
   quote &10NO NEED TO RANDOMIZE ANYTHING!!&10 ;
   proc/;
end 

if ( 'lor' > 0 ) set ispos 1 ;
else set ispos 0 ; end 

if ( !'testdown' )
     goto RANDOMIZE_UP ;
else goto RANDOMIZE_DOWN ; end 

proc/;


/************************************************************\ 
/************************************************************\ 
|                                                            |
|             I N I T I A L I Z E   A N D                    |
|            R E A D   P A R A M E T E R S                   |
|                                                            |
\************************************************************/

label READPARMS ;
set usertree 0 ;       /*  the tree for the test   */ 
set userx 0 ;          /*  primary (X) variable    */ 
set usery 1 ;          /*  secondary (Y) variable  */ 
set minincx 0.0001 ;   /*  minimum change in X     */
set minincy 0.0001 ;   /*  minimum change in Y     */
set radius 4 ;         /*  max. nr. of nodes to visit seeking
                           a change in y            */
set numpairs (-1) ;    /*  initial number of XY pairs */
set delfac 1 ;         /*  factor to downweight pairs,
                           wt = 1 / ( 1 + ( B*factor) ),
                           where B=branches in between  */ 
set stasis 0 ;         /*  are cases of (0,0) or (0,1) considered?  */ 
set randwts 0 ;        /*  randomize weights for cases of stasis?  */
set sample 50 ;        /*  number of combinations of reconsts. to look at   */
set numrepls 100 ;     /*  num replications */
set dograph 0 ;        /*  graphically display results and quit */
set testdown 0 ;       /*  does Y change only if X changed? */ 
set reccycles 5 ;      /*  number of cycles for each reconst
                           (randomized radius)               */
set userlengths 0 ;    /*  Dont use supplied branch lengths  */ 
set twotailed 0 ;      /*  Make test two-tailed  */
set notest 0 ;         /*  notest: skip randomization  */ 
set numsigdigs 3 ;     /*  number of significant digits to use */
set averagit  0 ;      /*  use average of min/max for terminals */
set rdsvtype 0 ;       /*  dont use files; 1 save to a file; 2 read from file */


var: i j k ;

if ( argnumber == 1 )
  if ( eqstring [ %1 help ] )
    goto DISPLAY_HELP ; 
    return 0 ;
    end
end 

if ( argnumber )

 loop 1 argnumber

   set i #1 ;
   set j #1 + 1 ;
   set k #1 + 2 ;

   if ( eqstring [ %('i') tree ] )
      set usertree %('j') ; var usertree ; 
      setloop 'k' ;
   else
   if ( eqstring [ %('i') chars ] )
      set userx %('j') ; var userx ; 
      set usery %('k') ; var usery ; 
      setloop ( #1 + 3 ) ; 
   else
   if ( eqstring [ %('i') minincx ] )
      set minincx %('j') ; var minincx ; 
      setloop ( #1 + 2 ) ;
   else
   if ( eqstring [ %('i') minincy ] )
      set minincy %('j') ; var minincy ; 
      setloop ( #1 + 2 ) ;
   else
   if ( eqstring [ %('i') radius ] )
      set radius %('j') ; var radius ; 
      setloop ( #1 + 2 ) ;
   else 
   if ( eqstring [ %('i') stasis ] )
      set stasis 1 ; var stasis ; 
   else
   if ( eqstring [ %('i') delfac ] )
      set delfac %('j') ; var delfac ; 
      setloop ( #1 + 2 ) ;
   else
   if ( eqstring [ %('i') sample ] )
      set sample %('j') ; var sample ; 
      setloop ( #1 + 2 ) ;
   else
   if ( eqstring [ %('i') cycles ] )
      set reccycles %('j') ; var reccycles ; 
      setloop ( #1 + 2 ) ;
   else
   if ( eqstring [ %('i') repls ] )
      set numrepls %('j') ; var numrepls ; 
      setloop ( #1 + 2 ) ;
   else
   if ( eqstring [ %('i') dograph ] )
      set dograph %('j') ; var dograph ;
      setloop ( #1 + 2 ) ; 
   else
   if ( eqstring [ %('i') twotailed ] )
      set twotailed 1; var twotailed ; 
   else
   if ( eqstring [ %('i') userlengths ] )
      set userlengths 1; var userlengths ; 
   else
   if ( eqstring [ %('i') testdown ] )
      set testdown 1; var testdown ; 
   else
   if ( eqstring [ %('i') notest ] )
      set notest 1; var notest ;
   else
   if ( eqstring [ %('i') randwts ] )
      set randwts 1; var randwts ; 
   else
   if ( eqstring [ %('i') saveto ] ) 
      set rdsvtype 1 ; 
      set savefilename $ %('j'); 
      setloop ( #1 + 2 ) ; 
   else 
   if ( eqstring [ %('i') rdfrom ] ) 
      set rdsvtype 2 ; 
      set savefilename $ %('j'); 
      setloop ( #1 + 2 ) ; 
   else 
   if ( eqstring [ %('i') average ] )
      set averagit 1 ; var averagit ; 
   else 
   if (eqstring [ %('i') sigdigs ] )
      set numsigdigs %('j') ; var numsigdigs ; 
      setloop ( #1 + 2 ) ;
   else
      errmsg Wrong argument: %('i'); 
      end end end end end end end end end end
      end end end end end end end end end end
 stop

else

 if ( !windows ) return 1 ; end
 var : 
   stgetminx[4]
   stgetminy[4] 
   stgetradius[4]
   stgetdelfac[4]
   stgetnumrepls[4]
   havtreetags ; 
 set stgetminx $'/.3minincx' ; 
 set stgetminy $'/.3minincy' ; 
 set stgetradius $'/.3radius' ; 
 set stgetdelfac $'/.3delfac' ;
 set stgetnumrepls $'/.0numrepls' ;
 set havtreetags ( eqstring[ x$ttag 0 x ] == 0 ) ; 

 opendlg 15 15 570 360 Measure and test correlation...

     frame 10 10 160 120 Choose chars. and tree ;
           spin 0 ntrees usertree 15  30 75 20 tree nr.; 
           spin 0 nchar userx     15  55 75 20 char. X; 
           spin 0 nchar usery     15  80 75 20 char. Y;
           subdlg                 50 105 110 20  H  e  l  p ;
              10 10 220 140 Choose chars. and tree ;
               showtxt 10 10 190 60
                    @@ CHOOSE_CHARS ; 
              closedlg

     frame 10 140 160 145 General settings ;
           gettxt  stgetminx   15 160 50 20 ;  showtxt 75 160 100 18 Min.incr.X;
           gettxt  stgetminy   15 185 50 20 ;  showtxt 75 185 100 18 Min.incr.Y;
           gettxt  stgetradius 15 210 50 20 ;  showtxt 75 210 100 18 (Max)radius;
           gettxt  stgetdelfac 15 235 50 20 ;  showtxt 75 235 100 18 Delay factor;
           subdlg              50 260 110 20  H  e  l  p ;
              10 140 500 100 General settings...;
               showtxt 10 10 475 300
                     @@ GENERAL_SETTINGS ; 
              closedlg

     frame 190 10 355 120 Observed R ;
           showtxt                   205 35 85 20 Sample up to ;
           spin    10 10000 sample   290 35 180 20 reconsts. for each character ;
           subdlg 460 55 80 20 Files... ;
               60 90 340 220 Save/read reconstructions... ;
               choose rdsvtype 
                   20 20 260 20 Do not use files 
                   20 40 260 20 Save reconstructions to a file for later use 
                   20 60 280 20 Use reconstructions previously saved to a file ;
               showtxt 20 85 280 20 NOTE!
               showtxt 20 105 300 20 Please use file names without spaces in their path!!;
               closedlg ;
           check randwts             205 70 120 20 randomize radius ;
              +
                spin 1 100 reccycles 325 70 40 20 times ;
              = 
           subdlg                    420 100 110 20  H  e  l  p ;
              190 10 500 100 Observed r;
               showtxt 10 10 475 300
                    @@ OBSERVED_R ; 
              closedlg

     frame 190 140 120 145 Randomized R ;
           gettxt  stgetnumrepls  220 160 50 20 ;  showtxt 220 180 70 30 Number of replications;
           check twotailed 220 225 70 20 two-tailed ; 
           subdlg              205 260 95 20  H  e  l  p ;
              190 140 500 100 Randomization ...;
               showtxt 10 10 475 300
                     @@ RANDOMIZATION ; 
              closedlg


     frame 330 140 215 145 Type of test ;
           choose testdown
               340 155 170 18 Does X make Y change?
                 + 
                    check stasis 380 175 150 18 Consider stasis ;
                 = 
               340 200 180 18 Does Y change only with X? ;
           check [('havtreetags')] userlengths 340 230 180 18 Use supplied branch lengths ; 
           subdlg 425 260 110 20  H  e  l  p ;
               300 140 300 120 Type of test;
               showtxt 10 10 275 300
                     @@ TYPE_OF_TEST ; 
              closedlg

      subdlg              40 295 190 20  Significant digits and ranges .... ; 
           60 90 280 140 Change to no more than... ; 
              spin   0 3 numsigdigs 20 20 180 20 significant digits;
              check averagit 20 45 190 20 Average min/max for terminals; 
      closedlg

 closedlg

  if ( 'rdsvtype' )
    if ( 'rdsvtype' == 1 )
         getfname savefilename write Save reconstructions to ... ;
    else getfname savefilename read Read reconstructions from ... ; end
    end 

  if ( ( 'numsigdigs' < 3 ) || 'averagit' ) 
    var : q r ;
    macfloat 'numsigdigs' ; 
    report- ; 
    if ( 'averagit' ) 
      loop 0 ntax
        if ( contmins['userx' #1] != 65.535 ) 
          set q contmins['userx' #1] + contmaxs['userx' #1] / 2 ;
          xread ='userx' #1 'q';
          end 
        if ( contmins['usery' #1] != 65.535 ) 
          set q contmins['usery' #1] + contmaxs['usery' #1] / 2 ;
          xread ='usery' #1 'q'; 
          end 
        stop 
    else 
    loop 0 ntax
      if ( contmins['userx' #1] != 65.535 ) 
        set q contmins['userx' #1] ;
        set r contmaxs['userx' #1] ;
        xread ='userx' #1 'q'-'r' ; 
        end 
      if ( contmins['usery' #1] != 65.535 ) 
        set q contmins['usery' #1] ;
        set r contmaxs['usery' #1] ;
        xread ='usery' #1 'q'-'r' ; 
        end 
     stop 
    end 
    macfloat 6 ; 
    end

 if ( exstatus )
     set minincx ( $stgetminx ) ;
     set minincy ( $stgetminy ) ; 
     set radius ( $stgetradius ) ; 
     set delfac ( $stgetdelfac ) ; 
     set numrepls ( $stgetnumrepls ) ;
     if ( 'minincy' == 0 ) set minincy 0.0001 ; end
     if ( 'minincx' == 0 ) set minincx 0.0001 ; end
     return 1 ;
 else return 0 ; end 
end

if ( ( 'numsigdigs' < 3 ) || 'averagit' )
  var : q r ;
  macfloat 'numsigdigs' ; 
  report- ; 
  if ( 'averagit' )
    loop 0 ntax
      if ( contmins['userx' #1] != 65.535 ) 
        set q contmins['userx' #1] + contmaxs['userx' #1] / 2 ;
        xread ='userx' #1 'q'; 
        end 
      if ( contmins['usery' #1] != 65.535 ) 
         set q contmins['usery' #1] + contmaxs['usery' #1] / 2 ;
         xread ='usery' #1 'q'; 
         end 
      stop 
  else 
  loop 0 ntax
    if ( contmins['userx' #1] != 65.535 ) 
      set q contmins['userx' #1] ;
      set r contmaxs['userx' #1] ;
      xread ='userx' #1 'q'-'r' ; 
      end 
    if ( contmins['usery' #1] != 65.535 ) 
      set q contmins['usery' #1] ;
      set r contmaxs['usery' #1] ;
      xread ='usery' #1 'q'-'r' ; 
      end 
   stop 
  end 
  macfloat 6 ; 
  end

return 1 ;

/********************************************************************\
/********************************************************************\
|                                                                    |
|                  C O R R E L A T E   U P                           |
|                                                                    |
|   This uses the command "iterrecs" to generate reconstructions,    |
|               and travtree to move along tree                      |
|                                                                    |
\********************************************************************/

label CORRELATE_UP ;
set reconum 0 ; 
silent = all ;
var : numrecsx numrecsy numrecs daprob ;
if ( 'dograph' ) 
  var: comesfrom[ (nnodes['usertree']+1) 2 ] ;
  end 

if ( 'rdsvtype' == 2 )
    proc $savefilename ;
else 
  /*********  Randomly select some reconstrucions for X
              (for subsequent testing)   ***************/
  set numstordx (-1) ;
  iterrecs 'usertree' 'userx' +/'maxstorrays' incx
      set numstordx ++ ;
      loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop 
      endrecs
  /*************** Likewise, for Y  ***************/
  set numstordy (-1) ;
  iterrecs 'usertree' 'usery' +/'maxstorrays' incy
      set numstordy ++ ;
      loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop 
      endrecs
  end

/*******************************************************
  Save to a file for later use, if so requested...
*/
if ( 'rdsvtype' == 1 )
    var : i ; 
    macfloat 3 ; 
    log $savefilename ;
    sil - file ;
    quote sil - all ., 
       if ( &39maxstorrays&39 != 'maxstorrays ) 
          errmsg Sorry, must be using the same number of reconstructions! (you saved 'maxstorrays', not &39maxstorrays&39) &38 10., 
          end ., ;
    set i ntax + 1 ; 
    quote if ( ( ntax + 1 ) != 'i' ) 
             errmsg Sorry, must be using the same number of taxa! (you saved reconstructions for 'i')  &38 10., 
             end ., ;
    quote &13set numstordx 'numstordx' ., ;
    quote &13set numstordy 'numstordy' ., ;
    var stordx * ;
    var stordy * ;
    quote sil = all., proc/., ;
    log / ;
    macfloat 6 ; 
    end 

/**************************************************************
    Note: if not using a random radius, a constant value of
          'tmpradius' will be used in all cases
*/
set tmpradius 'radius' ;


/**************************************************************  
    Now, begin enumerating combinations of reconsts. 
*/

set numrecsx ('numstordx'+1) ;
set numrecsy ('numstordy'+1) ; 
set numrecs 'numrecsx' * 'numrecsy' ; 

sil-all ;
quote TOTAL RECONSTRUCTIONS TO LOOK AT: 'numrecs' ('numrecsx'x'numrecsy')&10;
sil = all ;


set lor 100 ;
set hir (-100) ; 


loop =recx 0 'numstordx'

  loop =noxnum 0 nnodes['usertree']
    set incx[ #noxnum ] 'stordx[ #recx #noxnum ]' ;  
    stop 

  loop =recy 0 'numstordy' 

    loop =noynum 0 nnodes['usertree']
      set incy[ #noynum ] 'stordy[ #recy #noynum ]' ;  
      stop 

    set reconum ++ ;

    if ( 'dograph' ) 
       if ( 'reconum' != 'dograph' ) proc/; end
       end 

    progress 'reconum' 'numrecs' Checking... ;

    loop =cycles 1 'reccycles' 

      set numpairs (-1) ;
      loop =node 0 nnodes['usertree'] set visitd[#node] 0 ; stop 
  
      travtree up terms 'usertree' root - outnod
  
         if ( 'visitd['outnod']' ) proc/; end
         set myx 'incx['outnod']';
         if ( 'randwts' )
            set tmpradius getrandom[ 0 'radius' ] ;
            end
         if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx') ) )
  
             set didpair 0 ;
             set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ;
             set drag[anc['usertree' 'outnod']] 0 ; 

             travtree up terms 'usertree' 'outnod' innod

                set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ;
                set day 'incy['innod']' ;
                set drag['innod'] 'incx['innod']' + 'drag[anc['usertree' 'innod']]' ;

                if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) || ( 'innod' < ntax ) )
                   set didpair 1 ; 
                   set numpairs ++ ;
                   set xpairs['numpairs'] 'drag['innod']' ; 
                   set ypairs['numpairs'] 'day' ;
                   if ( 'dograph' )
                      set comesfrom['numpairs' 0] 'outnod' ; 
                      set comesfrom['numpairs' 1] 'innod' ;
                      end 
                   set wts['numpairs'] 1 / ( 1 + ( 'distosrc['innod']' * 'delfac' ) ) ;
                   travtree path 'usertree' 'innod' 'outnod' i  set visitd['i'] 1 ; endtrav
                   skipdes ;  
                   end 
             endtrav 
             if ( !'didpair' )
                 set numpairs ++ ;
                 set xpairs['numpairs'] 'myx' ;
                 set ypairs['numpairs'] 'incy['outnod']' ;
                 set wts['numpairs'] 1 ; 
                 if ( 'dograph' ) 
                    set comesfrom['numpairs' 0] 'outnod' ;
                    set comesfrom['numpairs' 1] 'innod' ;
                    end 
                 proc/; end
  
         else
                 
             if ( !'stasis' ) proc/; end      
             set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ;
             travtree up terms 'usertree' 'outnod' innod
                set dax 'incx['innod']' ; 
                if ( ( 'dax' >= 'minincx' ) || ( 'dax' <= (-'minincx' ) ) )
                   skipdes ;
                   proc/;
                   end 
                set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ;
                set day 'incy['innod']' ;
                if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) )
                   set numpairs ++ ;
                   set xpairs['numpairs'] 'myx' ; 
                   set ypairs['numpairs'] 'day' ;
                   if ( 'dograph' )
                      set comesfrom['numpairs' 0] 'outnod' ;
                      set comesfrom['numpairs' 1] 'innod' ;
                      end 
                   set wts['numpairs'] 1 ; 
                   travtree path 'usertree' 'innod' 'outnod' i  set visitd['i'] 1 ; endtrav
                   skipdes ;  
                   end 
             endtrav 
  
           end
  
      endtrav
 
    set numpairs ++ ;
    if ( 'numpairs' > 3 )
       if ( 'dograph' )    /******  Graphing ...  *****/ 
          sil - all ;
          quote &10HAVE 'numpairs' PAIRS&10;
          nak-; tp; 
          nak=; ttag=; tp; 
          var: daan dades ; 
          loop =node 0 nnodes['usertree'] set visitd[#node] 0 ; stop 
          loop =pair 0 ('numpairs'-1)
             set daan 'comesfrom[#pair 0]' ;
             set dades 'comesfrom[#pair 1]' ;
             if ( 'dades' != 'daan' ) 
                ttag +'dades' <#pair:'/+.3incx['dades']','/+.3incy['dades']'('/.3wts[#pair]');
                travtree path 'usertree' 'dades' 'daan' - i
                    if ( ( !'visitd['i']' ) && ( 'incx['i']' != 0 ) ) 
                         ttag +'i' '/+.3incx['i']'; 
                         set visitd['i'] 1 ; 
                         end 
                    if ( 'i' == 'daan' ) 
                         ttag +'i' |;
                    else ttag +'i' <; end
                    endtrav
             else
                ttag +'dades' X#pair:'/+.3incx[#pair]'/'/+.3incy[#pair]'('/.3wts[#pair]');
                end 
             stop
          ttag;
          ttag - ; 
          loop =pair 0 ('numpairs'-1)
             quote Pair #pair ('/3.0comesfrom[#pair 0]','/3.0comesfrom[#pair 1]'): '/+.3xpairs[#pair]','/+.3ypairs[#pair]'&10;
             stop
          quote &10;
          var: tmp['numpairs'];
          loop =pair 0 ('numpairs'-1) set tmp[#pair] 'xpairs[#pair]' ; stop 
          var + 60 25 tmp ypairs ;
          end   
       var & xpairs ypairs *wts 'numpairs' ;
       if ( 'dograph' ) killrecs ; end 
       if ( regr < 'lor' ) set lor regr ; set reclo 'reconum' ; end
       if ( regr > 'hir' ) set hir regr ; set rechi 'reconum' ; end
       end
 
     stop

  stop

stop

progress/;

if ( !'dograph' ) 
 if ( 'lor' > 10 )  /*  in this case, never did any pairs!!  */ 
    sil - all ;
    set lor 0 ;
    set hir 0 ; 
    quote CANNOT ESTABLISH CORRELATION.&10(TOO FEW PAIRS)&10;
    proc/; 
    end
 if ( ( 'lor' <= 0 ) && ( 'hir' >= 0 ) )
    sil - all ;
    quote RANGE OF OBSERVED R VALUES INCLUDES ZERO! ;
    proc/;
    end
 if ( 'lor' < 0 )
      set worsr 'hir' ;
 else set worsr 'lor' ; end 
 silent - all ;
 quote OBSERVED RANGE OF R: '/.4lor' - '/.4hir' (recs. 'reclo' - 'rechi')&10 ;
 sil = all ; 
end 

proc/; 


/**********************************************************\   
/**********************************************************\   
|                                                          |
|               R A N D O M I Z E    U P                   |
|                                                          |
\**********************************************************/

label RANDOMIZE_UP ;
set reconum 0 ; 
silent = all ;
var : lapstime nocalcs i j numrecs daprob xpick ypick 
      rali[ (nnodes['usertree']+1)] ;
set nocalcs 0 ; 
sil-all ;
quote RANDOMIZING RECONSTRUCTIONS: &10;
sil = all ; 
set matches 1 ; 
set logged 1 ; 
set tmpradius 'radius' ;

loop 1 'numrepls' 

    set numpairs (-1) ;
    loop 0 nnodes['usertree'] set visitd[#2] 0 ; stop

               /*****************************
                     Pick a reconstruction of X, at random, out of the 'numstordx' 
                     reconstructions stored, and copy it to the 'incx'. 
                     Also, pick a reconstruction of Y, at random, out of the 'numstordy' 
                     reconstructions stored, and randomly redistribute it over the tree 
                     (using rali, a random list of the branches)                        ***/

    set rali randomlist[ (nnodes['usertree']+1) ] ;
    set xpick getrandom[ 0 'numstordx'] ; 
    set ypick getrandom[ 0 'numstordy'] ; 
    loop 0 nnodes['usertree']
        set incx[#2] 'stordx['xpick' #2]' ;
        set incy['rali[#2]'] 'stordy['ypick' #2]' ;
        stop
                   /************************************************************
                        Now, establish correlation for those reconstructions, 
                        just like done for the observed correlation, but using 
                        the randomized Y-increments                             */
    travtree up terms 'usertree' root - outnod
       if ( 'visitd['outnod']' ) proc/; end
       set myx 'incx['outnod']'; 

       if ( 'randwts' )
          set tmpradius getrandom[ 0 'radius'] ;
          end

       if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx') ) )

           set didpair 0 ;
           set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ;
           set drag[anc['usertree' 'outnod']] 0 ; 
           travtree up terms 'usertree' 'outnod' innod
              set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ;
              set day 'incy['innod']' ;
              set drag['innod'] 'incx['innod']' + 'drag[anc['usertree' 'innod']]' ;
              if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) || ( 'innod' < ntax ) )
                 set didpair 1 ; 
                 set numpairs ++ ;
                 set xpairs['numpairs'] 'drag['innod']' ; 
                 set ypairs['numpairs'] 'day' ; 
                 set wts['numpairs'] 1 / ( 1 + ( 'distosrc['innod']' * 'delfac' ) ) ; 
                 travtree path 'usertree' 'innod' 'outnod' i  set visitd['i'] 1 ; endtrav
                 skipdes ;  
                 end 
           endtrav 
           if ( !'didpair' )
               set numpairs ++ ;
               set xpairs['numpairs'] 'myx' ;
               set ypairs['numpairs'] 'incy['outnod']' ;
               set wts['numpairs'] 1 ; 
               proc/; end

       else

           if ( !'stasis' ) proc/; end 
           set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ;
           travtree up terms 'usertree' 'outnod' innod
              set dax 'incx['innod']' ; 
              if ( ( 'dax' >= 'minincx' ) || ( 'dax' <= (-'minincx' ) ) )
                 skipdes ;
                 proc/;
                 end 
              set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ;
              set day 'incy['innod']' ;
              if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) )
                 set numpairs ++ ;
                 set xpairs['numpairs'] 'myx' ; 
                 set ypairs['numpairs'] 'day' ;
                 set wts['numpairs'] 1 ; 
                 travtree path 'usertree' 'innod' 'outnod' i  set visitd['i'] 1 ; endtrav
                 skipdes ;  
                 end 
           endtrav 

         end

    endtrav 

  set numpairs ++ ; 
  if ( 'numpairs' > 3 )
     var & xpairs ypairs *wts 'numpairs' ;
               /*******************************************************
                  Establish regression for this randomization, and 
                  compare with observed value of r                     */
     if ( 'twotailed' ) 
        if ( regr < 0 )
            if ( ( 'ispos' && ( (-regr) >= 'lor' ) ) ||
                 ( (!'ispos') && ( regr <= 'hir' ) ) )
                 set matches ++ ; end
        else
            if ( ( (!'ispos') && ( (-regr) <= 'hir' ) ) ||
                 ( 'ispos' && ( regr >= 'lor' ) ) )
                 set matches ++ ; end
        end
     else 
       if ( 'ispos' && ( regr >= 'lor' ) ) set matches ++ ;
       else if ( (!'ispos') && ( regr <= 'hir' ) )
                set matches ++ ;
                end
            end
       end 
     set signif ( 'matches' / ( #1 + 1 ) ) ; 
     set i #1 ; 
     set j regr ;
     set lapstime time ; 
     if ( 'logged' == 10 ) sil - all ; 
     else if ( 'lapstime' >= 1 ) resettime ; sil - con ; end end 
     if ( !windows ) 
        quote &0 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j';
        if ( 'logged' == 10 ) quote &10 ; set logged 0 ; end 
     else 
        progress 'i' 'numrepls' 'i' REPLS - P(inst) '/.3signif' ('matches' matches) ; 
        if ( 'logged' == 10 ) 
            quote &32 &32 &32 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j'&10;
            set logged 0 ; end 
         end
     set logged ++ ; 
     sil = all ; 
  else
     set nocalcs ++ ; 
     setloop #1 ;
     end

stop 

if ( windows ) progress/; end
lquote - ; 
sil - all ;
if ( !windows ) 
    quote &0 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10;
else 
    quote &32 &32 &32 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10;
end 

proc/; 


/****************************************************************\ 
/****************************************************************\ 
|                                                                |
|               C O R R E L A T E    D O W N                     |
|                                                                |
\****************************************************************/

label CORRELATE_DOWN ;

set reconum 0 ;
silent = all ;
var : distoy numrecsx numrecsy foundanx numrecs daprob ;
if ( 'dograph' ) 
  var: comesfrom[ (nnodes['usertree']+1) 2 ] ;
  end 

if ( 'rdsvtype == 2 )
    proc $savefilename ;
else 
  /*********  Randomly select some reconstrucions for X
              (for subsequent testing)   ***************/
  set numstordx (-1) ;
  iterrecs 'usertree' 'userx' +/'maxstorrays' incx
      set numstordx ++ ;
      loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop 
      endrecs
  /*************** Likewise, for Y  ***************/
  set numstordy (-1) ;
  iterrecs 'usertree' 'usery' +/'maxstorrays' incy
      set numstordy ++ ;
      loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop 
      endrecs
  end

/*******************************************************
  Save to a file for later use, if so requested...
*/
if ( 'rdsvtype' == 1 )
    var : i ; 
    macfloat 3 ; 
    log $savefilename ;
    sil - file ;
    quote sil - all ., 
       if ( &39maxstorrays&39 != 'maxstorrays ) 
          errmsg Sorry, must be using the same number of reconstructions! (you saved 'maxstorrays', not &39maxstorrays&39)  &38 10., 
          end ., ;
    set i ntax + 1 ; 
    quote if ( ( ntax + 1 ) != 'i' ) 
             errmsg Sorry, must be using the same number of taxa! (you saved reconstructions for 'i')  &38 10., 
             end ., ;
    quote &13set numstordx 'numstordx' ., ;
    quote &13set numstordy 'numstordy' ., ;
    var stordx * ;
    var stordy * ;
    quote sil = all., proc/., ;
    log / ;
    macfloat 6 ; 
    end 

/**************************************************************
    Note: if not using a random radius, a constant value of
          'tmpradius' will be used in all cases
*/
set tmpradius 'radius' ;


/**************************************************************  
    Now, begin enumerating combinations of reconsts. 
*/

set numrecsx ('numstordx'+1) ;
set numrecsy ('numstordy'+1) ; 
set numrecs 'numrecsx' * 'numrecsy' ; 

sil-all ;
quote TOTAL RECONSTRUCTIONS TO LOOK AT: 'numrecs' ('numrecsx'x'numrecsy')&10;
sil = all ;

set lor 100 ;
set hir (-100) ; 

loop =recx 0 'numstordx'

  loop =nonum 0 nnodes['usertree']
    set incx[ #nonum ] 'stordx[ #recx #nonum ]' ;  
    stop 

  loop =recy 0 'numstordy' 

    loop =nonum 0 nnodes['usertree']
      set incy[ #nonum ] 'stordy[ #recy #nonum ]' ;  
      stop 

    set reconum ++ ;

    if ( 'dograph' )
       if ( 'reconum' != 'dograph' ) proc/; end
       end 

    progress 'reconum' 'numrecs' Checking... ;

    loop =cycle 1 'reccycles' 
 
      set numpairs (-1) ;
      set tmpradius 'radius' ;
  
      loop =node 0 nnodes['usertree']
  
        if ( #node == root ) continue ; end
  
        if ( 'randwts' ) set tmpradius getrandom [ 0 'radius' ] ; end 
        set day 'incy[ #node ]' ;
        set outnod #node ;
        if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) ) 
            set i 0 ;
            set dax 0 ;
            set foundanx 0 ;
            set distoy 0 ; 
            travtree path 'usertree' 'outnod' root innod
                set myx 'incx['innod']';
                set dax += 'myx' ;
                if ( !'foundanx' ) 
                  if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx' ) ) )
                     set distoy 'i' ;
                     set foundanx 1 ;
                     end
                  end 
                set i ++ ;
                if ( 'i' > 'tmpradius' ) killtrav ; end 
               endtrav 
            set numpairs ++ ;
            set xpairs [ 'numpairs' ] 'dax' ;
            set ypairs [ 'numpairs' ] 'day' ;
            if ( 'dograph' )
                set comesfrom[ 'numpairs' 0 ] 'innod' ;
                set comesfrom[ 'numpairs' 1 ] 'outnod' ; 
                end
            if ( 'foundanx' ) 
                 set wts['numpairs'] 1 / ( 1 + ( 'distoy' * 'delfac' ) ) ;
            else set wts['numpairs'] 1 ; end 
            end
  
        stop
  
      set numpairs ++ ;
      if ( 'numpairs' > 3 )
        if ( 'dograph' )     /***   Graphing ...   **/ 
          sil - all ;
          quote &10HAVE 'numpairs' PAIRS&10;
          loop =pair 0 ('numpairs'-1)
             quote Pair #pair ('/3.0comesfrom[#pair 0]','/3.0comesfrom[#pair 1]'): '/+.3xpairs[#pair]','/+.3ypairs[#pair]'&10;
             stop
          quote &10;
          var: tmp['numpairs'];
          loop =pair 0 ('numpairs'-1) set tmp[#pair] 'xpairs[#pair]' ; stop 
          var + 60 25 tmp ypairs ;
          end   
        var & xpairs ypairs *wts 'numpairs' ;
        if ( 'dograph' ) killrecs ; end 
        if ( regr < 'lor' ) set lor regr ; set reclo 'reconum' ; end
        if ( regr > 'hir' ) set hir regr ; set rechi 'reconum' ; end
        end

      stop

  stop

stop

progress/;

if ( !'dograph' ) 
 if ( 'lor' > 10 )  /*  in this case, never did any pairs!!  */ 
    sil - all ;
    set lor 0 ;
    set hir 0 ; 
    quote CANNOT ESTABLISH CORRELATION.&10(TOO FEW PAIRS)&10;
    proc/; 
    end
 if ( ( 'lor' <= 0 ) && ( 'hir' >= 0 ) )
    sil - all ;
    quote RANGE OF OBSERVED R VALUES INCLUDES ZERO! ;
    proc/;
    end
 if ( 'lor' < 0 )
      set worsr 'hir' ;
 else set worsr 'lor' ; end 
 silent - all ;
 quote OBSERVED RANGE OF R: '/.4lor' - '/.4hir' (recs. 'reclo' - 'rechi')&10 ;
 sil = all ; 
end 

proc/;

/**************************************************************\ 
/**************************************************************\ 
|                                                              |
|               R A N D O M I Z E    D O W N                   |
|                                                              |
\**************************************************************/

label RANDOMIZE_DOWN ;
set reconum 0 ; 
silent = all ;
var : lapstime nocalcs i j numrecs daprob xpick ypick 
      rali[ (nnodes['usertree']+1)] foundanx distoy ;
set nocalcs 0 ; 
sil-all ;
quote RANDOMIZING RECONSTRUCTIONS: &10;
sil = all ; 
set matches 1 ; 
set logged 1 ; 
set tmpradius 'radius' ;

loop 1 'numrepls' 

    set numpairs (-1) ;
    set tmpradius 'radius' ;

    set rali randomlist[ (nnodes['usertree']+1) ] ;
    set xpick getrandom[ 0 'numstordx'] ; 
    set ypick getrandom[ 0 'numstordy'] ; 
    loop 0 nnodes['usertree']
        set incx[#2] 'stordx['xpick' #2]' ;
        set incy['rali[#2]'] 'stordy['ypick' #2]' ;
        stop

    loop =node 0 nnodes['usertree']

      if ( #node == root ) continue ; end

      if ( 'randwts' ) set tmpradius getrandom [ 0 'radius' ] ; end 
      set day 'incy[ #node ]' ;
      set outnod #node ;
      if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) ) 
          set i 0 ;
          set dax 0 ;
          set foundanx 0 ;
          set distoy 0 ; 
          travtree path 'usertree' 'outnod' root innod
              set myx 'incx['innod']';
              set dax += 'myx' ;
              if ( !'foundanx' ) 
                if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx' ) ) )
                   set distoy 'i' ;
                   set foundanx 1 ;
                   end
                end 
              set i ++ ;
              if ( 'i' > 'tmpradius' ) killtrav ; end 
             endtrav 
          set numpairs ++ ;
          set xpairs [ 'numpairs' ] 'dax' ;
          set ypairs [ 'numpairs' ] 'day' ;
          if ( 'dograph' )
              set comesfrom[ 'numpairs' 0 ] 'innod' ;
              set comesfrom[ 'numpairs' 1 ] 'outnod' ; 
              end
          if ( 'foundanx' ) 
               set wts['numpairs'] 1 / ( 1 + ( 'distoy' * 'delfac' ) ) ;
          else set wts['numpairs'] 1 ; end 
          end

      stop
  
   set numpairs ++ ;
   if ( 'numpairs' > 3 )
     var & xpairs ypairs *wts 'numpairs' ;
     if ( 'twotailed' ) 
        if ( regr < 0 )
            if ( ( 'ispos' && ( (-regr) >= 'lor' ) ) ||
                 ( (!'ispos') && ( regr <= 'hir' ) ) )
                 set matches ++ ; end
        else
            if ( ( (!'ispos') && ( (-regr) <= 'hir' ) ) ||
                 ( 'ispos' && ( regr >= 'lor' ) ) )
                 set matches ++ ; end
        end
     else 
        if ( 'ispos' && ( regr >= 'lor' ) ) set matches ++ ;
        else if ( (!'ispos') && ( regr <= 'hir' ) )
                set matches ++ ;
                end
             end
        end 
     set signif ( 'matches' / ( #1 + 1 ) ) ; 
     set i #1 ; 
     set j regr ;
     set lapstime time ; 
     if ( 'logged' == 10 ) sil - all ; 
     else if ( 'lapstime' >= 1 ) resettime ; sil - con ; end end 
     if ( !windows ) 
        quote &0 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j';
        if ( 'logged' == 10 ) quote &10 ; set logged 0 ; end 
     else 
        progress 'i' 'numrepls' 'i' REPLS - P(inst) '/.3signif' ('matches' matches) ; 
        if ( 'logged' == 10 ) 
            quote &32 &32 &32 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j'&10;
            set logged 0 ; end 
         end
     set logged ++ ; 
     sil = all ; 
  else
     set nocalcs ++ ; 
     setloop #1 ;
     end

  stop 

if ( windows ) progress/; end 
lquote - ; 
sil - all ; 
if ( !windows ) 
    quote &0 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10;
else 
    quote &32 &32 &32 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10;
end 

proc/;

/************************************************************\
/************************************************************\
|                                                            |
|                  C H E C K P A R M S                       |
|           (it makes sure parameters are ok)                |
|                                                            |
\************************************************************/

label CHECKPARMS ;

if ( (!iscont['userx']) || (!iscont['usery']) ) 
   errmsg Sorry, both characters X and Y must be continuous ;
   end

if ( tsize['usertree'] != root )
   errmsg Sorry, all taxa in the matrix must be included in the tree!!;
   end 

proc/; 


/***************************************************************\
/***************************************************************\
|                                                               |
|                   D I S P L A Y   H E L P                     |
|                                                               |
\***************************************************************/

label DISPLAY_HELP 
lquote [ ; 
sil - all ; 
quote * 

This script tests character correlation
between continuous characters.

Argument:           default:

   tree N                0    
   chars X Y            0,1  
   minincx N             0    
   minincy N             0    
   radius N              4    
   stasis                no   
   delfac N              1    
   sample N             50
   randwts               no 
   cycles N          1, if randwts: 5
   repls N              100
   userlengths           no
   twotailed             no
   testdown              up
   dograph N             no
   notest               test 
   average              no
   sigdigs N            3 
   saveto filename
   readfrom filename
&10; 

proc/; 


/***************************************************************\
/***************************************************************\
|                                                               |
|             T E X T   F O R   D I A L O G S                   |
|                                                               |
\***************************************************************/

label CHOOSE_CHARS
Note that the tree must be complete (i.e. have all taxa).  Both X and Y must be continuous characters ;
@@

label GENERAL_SETTINGS
Changes below Min/Max increases are ignored. Max radius is the number of nodes beyond which the dependent variable is considered to not have changed. As the delay factor increases, changes in the dependent variable occuring more nodes away are considered to be less compelling evidence of correlation.
@@

label OBSERVED_R
The observed R is calculated by using a sample of reconstructions for the dependent and independent variables (the lowest R is used).  The radius can be randomized for each combination of reconstructions, so that the conclusion does not depend on a specific radius value.
@@ 

label RANDOMIZATION
The randomization re-distributes observed changes (for different alternative reconstructions for X and Y) on the tree, and recalculates R.
@@

label TYPE_OF_TEST
Two types of test can be performed.  Branch-lengths must be provided as tree-tags.
@@