sec : rec 1 ; macro[15000 ; macro * 10 250000 ; macro=; lquote [ ; if ( !argnumber ) errmsg ----------------------------------------------- This script runs (serially) a sectorial approximmation of supports, for large data sets. It automatically breaks the tree in chunks of 200-166 nodes. The groups are those in a tree found by a quick search (see below). You have to give the name of the data set as first argument (name only, extension *must* be &34.tnt&34). Results are written to screen (or log file, iif opened before running this routine), and saved to a file with the same name as the data set but extension &34.sup&34. ---------------------------------------------- ; end /** Read the data set (with memory-saving options) *******/ rseed ! ; tshrink ! ; cost < ; p %1.tnt ; /** Variables and some settings... ****/ rseed ! ; ra 1 ; var : wrktodo wrkdone damin damax cycles num_fixed_nodes timesect bigvals[(nnodes[0]+1)] redvals[(root+1)] GCfullvals[(nnodes[0]+1)] GCbigvals[(nnodes[0]+1)] GCredvals[(root+1)] nodlist[(root+1)] i j k ; if ( ntax > 5000 ) mult : wc 1000 ; bb : clu 40 ; end txts 10000 ; goto =%0 ; keep 0 ; hold 105 ; set num_fixed_nodes 0 ; set damin ntax / 200 ; /* minimum number of pieces in which to break tree */ set damax ntax / 166 ; /* maximum number of pieces in which to break tree */ set cycles 3 ; /* num of cycles, set to either 2 or 3 */ if ( ( 'damin' < 4 ) || ( 'damax' < 6 ) ) set damin 5 ; set damax 8 ; end if ( 'damax' > bittype ) set damin (bittype - 8) ; set damax bittype ; end if ( 'cycles' == 3 ) set wrktodo 2 * 'damin' + ( ( 'damax' - 'damin' ) / 2 ) + 'damax' ; else set wrktodo 'damin' + 'damax' ; end set wrkdone 0 ; if ( !windows ) if ( !bground ) report= ; end ; else report= ; end /****** Search a starting tree ... CHANGE THESE ROUTINES IF YOU WANT TO SEARCH DIFFERENTLY, OR MAKE IT READ A TREE FROM A FILE IF YOU WANT TO USE YOUR OWN TREE... *****/ rseed 0 ; rseed * ; mu1 =ho1 ; sec = rss godr 100000 maxs 100 mins 40 ; sec = xss10-15+3-1 gocomb 10 combst 5 fuse 4 drift 6 ; drift = iter 6 ; keep 1 ; /*** Begin Sectorial Resampling ****/ loop root nnodes[0] set bigvals[ #1 ] 200 ; /* i.e. an impossible value! */ set GCbigvals[ #1 ] 200 ; stop report- ; resettime ; sil = con ; bb: clu 15 ; sec = track tree 0 chkroot xss '/.0damax'-'/.0damin'+'cycles' [ set wrkdone ++ ; set i time ; set j 'i' * 'wrktodo' / 'wrkdone' ; progress 'wrkdone' 'wrktodo' Predicted time '/.0j's ; coll tbr ; resample rep 100 [ mu1=ho1 ] from 0 savetree ; tc 0 2. ; coll none ; set redvals freqlist [0] 1. ; set GCredvals freqdlist [0] 1. ; set nodlist biginsect[0] ; loop 0 (listsize-1) set i 'nodlist[#1]' ; set j nodtosect[ 'i' 0 ] ; if ( 'j' <= ntax ) continue ; end if ( 'bigvals[ 'i' ]' > 'redvals[ 'j' ]' ) set bigvals[ 'i' ] 'redvals[ 'j' ]' ; end if ( 'GCbigvals[ 'i' ]' > 'GCredvals[ 'j' ]' ) set GCbigvals[ 'i' ] 'GCredvals[ 'j' ]' ; end stop ] ; progress/; /*** Make sure no nodes are left-over ******/ loop (root+1) nnodes[0] if ( anc [ 0 #1 ] == root ) continue ; end if ( 'bigvals[#1]' < 200 ) continue ; end sil - buff ; quote Fixing values for node #1 ... ; sil = buff ; set num_fixed_nodes ++ ; sec = track tree 0 dss #1 8 [ coll tbr ; resample rep 100 [ mu1=ho1 ] from 0 savetree ; tc 0 2. ; coll none ; set redvals freqlist [0] 1. ; set GCredvals freqdlist [0] 1. ; set nodlist biginsect[0] ; loop 0 (listsize-1) set i 'nodlist[#1]' ; set j nodtosect[ 'i' 0 ] ; if ( 'j' <= ntax ) continue ; end if ( 'bigvals[ 'i' ]' > 'redvals[ 'j' ]' ) set bigvals[ 'i' ] 'redvals[ 'j' ]' ; end if ( 'GCbigvals[ 'i' ]' > 'GCredvals[ 'j' ]' ) set GCbigvals[ 'i' ] 'GCredvals[ 'j' ]' ; end stop ] ; stop /*** Time for complete estimation ***/ set timesect time ; sil - buffer file ; sil = con ; report=; /*** Write labels to tree ****/ set bigvals[ 0 ] 100 ; set GCbigvals[ 0 ] 100 ; ttag - ; ttag = ; nak = ; tp 0 ; loop (root+1) nnodes[0] if ( anc [ 0 #1 ] == root ) ttag +#1 100/100; continue ; end ttag +#1 '/.0bigvals[ #1 ]'/'/.0GCbigvals[ #1 ]'; stop /*** Output results ****/ set wrkdone root ; quote RESULTS FOR %1, 'wrkdone' taxa (Freq/GC); quote Time: 'timesect' sec., additional pass over 'num_fixed_nodes' nodes. ; ttag ; /** Save values *****/ tsav * %1.sup; sav * ; tsav / ; proc/;